Ch 7. Multiphase SO2-H2O2 Chemistry

OVERVIEW

    • In previous chapters we learnt how to model simple gas-phase and heterogeneous chemistry using a kinetic box model. Water is a vital component of the Earth's atmosphere, and liquid water exists in a variety of forms in the atmosphere, such as cloud, fog, liquid aerosols, etc. Chemical reactions may occur in the atmospheric waters and affect the physical and chemical properties of the atmosphere.
    • In this chapter we'll learn how to model simple multiphase chemistry. We'll take SO2 + H2O2 chemistry as an example. SO2 is emitted from a variety of natural and anthropogenic sources, it is also produced in the atmosphere from chemical reactions (e.g., DMS oxidation). SO2 is somewhat soluble, and its solubility depends on pH. H2O2 is an important oxidant in the atmosphere, which is very soluble. SO2 oxidation in the atmosphere produces sulfate, a major component of the atmospheric aerosol particles. SO2 oxidation in the gas-phase is slow, but SO2 may undergo cloud-processing & oxidized by a number of oxidants in the aqueous-phase, including H2O2.
    • In this chapter we'll model the explicit phase-transfer of SO2 and H2O2, as well as the reaction between S(IV) - the dissolved product of SO2 and H2O2 in the aqueous-phase.
    • The phase-transfer is largely based on Stephen E. Schwartz's legendary 1986 paper (1), which is one of my all time favorate. Many key concepts are elaborated in THE BIG BOOK by John H. Seinfeld and Spyros N. Pandis (2). Rate coefficients and parameters are mostly from the latest JPL compilation (3)

References 1. Schwartz, 1986. Mass-Transport Considerations Pertinent to Aqueous Phase Reactions of Gases in Liquid-Water Clouds 2. Seinfeld and Pandis, 2006. Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. 3. Burkholder et al., 2015. Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies

1. Multiphase chemical mechanism

  • The chemical mechanism in this demo is super simple, with 5 species and 3 reactions. See the screenshot of the BIG TABLE. The subscript "_aq" stands for "aqueous-phase"
  • For simplicity, we lump all aqueous-phase S(IV) products into SIV_aq, which is the sum of [SO2-H2O], [HSO3-], and [SO3^2-].
  • You may manually type in these...

2. Generate ODE

  • Click the ODE scratchboard button. You know the trick - if there is already an ODE_scratchboard pops out, kill it.
  • Then click the Make New ODE button, and generate your own ODEs based on the chemical mechanism you just entered.
  • It should look something like this screenshot.

3. Set up phase-transfer

  • Now the fun part: scroll to the bottom of the ODE_scratchboard.ipf, copy-paste the following snippet:
Variable/G LWC, TEMP Variable/G kt_SO2, kt_H2O2 Variable/G Heff_SO2, Heff_H2O2 // ---------------------------- // Phase transfer of SO2 // ---------------------------- YDOT[0] = YDOT[0] - kt_SO2 * (Conc[0]*LWC - Conc[2]/Heff_SO2/0.082/TEMP) YDOT[2] = YDOT[2] + kt_SO2 * (Conc[0]*LWC - Conc[2]/Heff_SO2/0.082/TEMP) // ------------------------------ // Phase transfer of H2O2 // ------------------------------ YDOT[1] = YDOT[1] - kt_H2O2 * (Conc[1]*LWC - Conc[4]/Heff_H2O2/0.082/TEMP) YDOT[4] = YDOT[4] + kt_H2O2 * (Conc[1]*LWC - Conc[4]/Heff_H2O2/0.082/TEMP)
  • Let's go thru this: firstly we declare a number of global variables:
LWC: liquid water content, in vol-water per vol-air TEMP: temperature, in K kt_SO2 and kt_H2O2: phase-transfer coefficients for SO2 and H2O2, respectively, in s^-1 Heff_SO2 and Heff_H2O2: (effective) Henry's law constants for SO2 and H2O2, respectively, in M/atm.
  • The values of these are defined in the main procedure, and we'll go back to this later. I'm lazy so I pass these parameters as global variables. You may pack them into one or more waves and pass them into this function too (see this example in Chapter 4). Be careful when using global variables! You may unintentially change the value!
  • For each species, we need two YDOT terms, describing its gas-to-aqueous and aqueous-to-gas transfer. In this example, species 0 is SO2, and species 2 is its aqueous counterpart S(IV). See Schwartz (1986) and Seinfeld and Pandis (2016) for details.

4. Main procedure: set up conditions

  • Now hit ctrl+M and jump to the main procedure that you're already familiar with. Copy-paste the following snippet into your Main Procedure and overwrite the corresponding sections:
// =========== // House keeping // =========== Variable/G TEMP = 298, Pressure = 101325, LWC = 1e-7 Variable Radius_nm = 50000 Variable MassAccommodationCoeff_SO2 = 0.12 Variable MassAccommodationCoeff_H2O2 = 0.18 Variable Dg_SO2_cm2_s = 0.15 Variable Dg_H2O2_cm2_s = 0.12 Variable Hp = 10^(-1*2) Variable H_SO2, K1_SO2, K2_SO2 Variable/G kt_SO2, kt_H2O2 Variable/G Heff_SO2, Heff_H2O2 // ================================================= // Initial concentrations // If set here, will overwrite the intial concentrations in the BIG TABLE // ================================================= Concentration_Matrix[0][0] = (5)*Pressure*(6.0232E+8)/(8.314*TEMP) // SO4 Concentration_Matrix[0][1] = (1)*Pressure*(6.0232E+8)/(8.314*TEMP) // H2O2
  • Everything should be straightforward. We set to 298 K, 101325 Pa, LWC = 1e+7 isn't a bad estimate for cloud droplet. We assume all cloud droplets are 50 microns in radius. Then we set mass accommodation coefficients, gas-phase diffusion coefficients. Then we declare a number of variables, Henry's law constants, phase-transfer coefficients, etc.
  • pH affects the H2O2-S(IV) rate coefficient. Usually in cloud/fog droplets, aerosol liquid waters, pH is controlled/buffered by a number of things. To make our life easier (for now), we set constant pH. In this demo, pH is set to 2. You may changes to your favorate pH later.
  • Then we set initial concentrations: 5 ppb for SO2 and 2 ppb for H2O2.

5. Main procedure: set up equilibrium, rates, etc

  • Now scroll a little lower, go into the for-loop (where we loop thru each time step), and copy-paste the following:
H_SO2 = 1.23 K1_SO2 = 1.3E-2 * exp(1960*(1/TEMP - 1/298)) K2_SO2 = 6.6E-8 * exp(1500*(1/TEMP - 1/298))
 kt_SO2 = (((Radius_nm*1E-7)^2)/3/(Dg_SO2_cm2_s) + 4*(Radius_nm*1E-7)/3/(100*sqrt(8*8.314*TEMP*1000/3.14159/64.066))/MassAccommodationCoeff_SO2)^-1 kt_H2O2 = (((Radius_nm*1E-7)^2)/3/(Dg_H2O2_cm2_s) + 4*(Radius_nm*1E-7)/3/(100*sqrt(8*8.314*TEMP*1000/3.14159/64.066))/MassAccommodationCoeff_SO2)^-1 Heff_SO2 = H_SO2 * (1 + K1_SO2/Hp + K1_SO2*K2_SO2/Hp/Hp) Heff_H2O2 = 1E+5 // ================================ // Define time-dependent rate-coefficients here // ================================ kr = 0 kr[2] = k_H2O2_SO2_M_s(Hp, TEMP) / (6.0232e+23 / 1000 * LWC)
  • Everything should be straightforward. We firstly set the Henry's law constant of SO2 (not effective), and it's dissociation constants (K1_SO2 and K2_SO2). We then setup the phase-transfer coefficients and the effective Henry's law constants (these will be passed into the ODE_scratchboard.ipf). The formulation of kt can be found in See Schwartz (1986) and Seinfeld and Pandis (2016) for details. Note the effective Henry's law constant of SO2 depends on pH.
  • Then we setup the rate coefficients. We have three reactions: R0 and R1 and phase-transfer reactions, and their "rates" are setup in the ODE_Scratchboard.ipf, so we don't need to define them here. kr[2] is the aqueous-phase reactions of H2O2 + S(IV). Its 2nd order rate coefficient is given in a function k_H2O2_SO2_M_s(Hp, TEMP), in M^-1 s^-1. We'll get back to this later.
  • Unit conversion (VERY IMPORTANT): Note the term 1/(6.0232e+23 / 1000 * LWC) here. This converts the unit of the rate coefficient from M^-1 s^-1 to (molec/cm3)^-1 s^-1. That is, although both the reactants are in the aqueous-phase, they're still solved in the model as if they're in the gas-phase, so the units of all aqueous-phase reactions have to be consistent.
  • The unit conversion from aqueous-phase to gas-phase depends on the order of the reaction. If it's a first-order reaction (in s^-1), then no unit conversion is needed (yay). If 2nd order, you already know how to do it. How about 3rd order?

Bottle neck effect: mass transfer considerations

This multiphase chemical system consists of several processes:

  • Diffusion of SO2 and H2O2 from "far far away" to near the droplet
  • SO2 and H2O2 attach onto the surface of the droplet
  • SO2 and H2O2 undergo processes in the droplet (e.g., dissociation, hydrolysis, chemical reactions)

The phase-transfer velocity term, kt, consists of two components:

  • r^2/3/Dg: diffusion in the gas-phase, where r is the radius of the droplet and Dg is the gas-phase diffusion coefficient
  • 4r/3/c/a: surface accommodation, where c is the thermal speed, and a is the mass accommodation coefficient

Overall the slowest process limits the overall speed. We discussed this in more detail in Chapter 5. For highly viscous solutions, aqueous-phase reactions may be limited by diffusion in the bulk aqueous-phase. We will not cover this in this chapter. Usually for dilute conditions, many chemical reactions probably won't be affected by bulk diffusion.

6. pH-dependent rate coefficient of H2O2+S(IV)

  • Now insert the following function into your procedure window:
Function k_H2O2_SO2_M_s(Hp, TEMP) Variable Hp, TEMP Variable K1_SO2 = 1.3E-2 * exp(1960*(1/TEMP - 1/298)) Variable K2_SO2 = 6.6E-8 * exp(1500*(1/TEMP - 1/298)) Variable f_HSO3m = K1_SO2*Hp / (Hp*Hp + K1_SO2*Hp + K1_SO2*K2_SO2) Return 7.45E+7 * Hp / (1 + 13*Hp) * f_HSO3m End
  • The original pH-dependent rate expression, kII (M^-1 s^-1) = 7.45E+7 * Hp / (1 + 13*Hp) can be found in Seinfeld and Pandis (2016). Note that this rate coefficient is for H2O2 + HSO3-, but here we write our reaction as H2O2 + S(IV), therefore we have to multiply it by the fraction of HSO3-, i.e., the f_HSO3m term.
  • Not sure if we did this right? Make a plot of kII vs pH (see the plot) and compare that to the plot on Page 311 in Seinfeld and Pandis (2016).

7. Now let's roll!

  • Set the total run time = 0.1 day and output time interval to 60 seconds, and hit the run button! See SO2 and H2O2 both decay and aqueous S(VI)!
  • Now change pH from 2 to 8, and run the model again. See the difference?
  • pH dependency: pH affects two things here: firstly, the H2O2+S(IV) rate, and secondly, the SO2 solubility. As seen in the plot, when change pH from 2 to 8, the formation of S(VI) decreased dramatically because H2O2+S(IV) is much slower. But SO2 drops faster, since the effective Henry's law constant of SO2 increased with increasing pH.

8. Equilibrium timescale

  • You may have noticed that at pH = 2, gaseous SO2 drops "faster" than H2O2 does. But when pH = 7, it's the other way around. This is because, it takes longer for soluble species to reach equilibrium.
  • To further test this, let's do a few more things: shut off the H2O2+S(IV) reaction, by setting kr[2] = 0; set the output time interval to 20 seconds, and the total run time to 0.02 day; change initial concentrations of SO2 and H2O2 both to 2 ppb. Now we don't have any aqueous-phase reactions, we only have phase-transfer. You'll see the decrease of gas-phase species and the increase of aqueous-phase counter species both hit a plateau, i.e., reaching the phase-transfer equilibrium, after which there is no net flux between the gas-phase and aqueous-phase. When pH = 2, H2O2 is more soluble, and SO2 reaches equilibrium almost immediately. When pH = 7, SO2 effective Henry's law constant is greatly enhanced, now it takes longer for SO2 to reach equilibrium.
  • The equilibrium timescale is elaborated in both Schwartz (1986) and Seinfeld and Pandis (2016).

9. Summary

  • In this chapter, we learnt how to model simple multiphase chemistry, with explicit phase-transfer and aqueous-phase chemistry. Modeling the aqueous-phase chemistry is essentially the same idea as the gas-phase chemistry, except that we may have to convert the unit accordingly.
  • Phase-transfer is a different story though. You need to setup the explicit gas-to-aqueous and aqueous-to-gas terms for the differential equations. I have provided template. If you need to do this for a large number of species, I recommend you write a script to "generate" scripts rather than writing them by hand. You can do this in Microsoft Excel ;)
  • Lots of things you can play around with! For example, the effect of surface accommodation, size of the droplet, etc.
  • Generally speaking, a chemical system with aqueous-phase processes are computationally more "stiff", since the phase-transfer (s^-1) is usually much faster than many gas-phase reactions. As a result, sometimes the solver may crash. If this happens, you may change the tolerance setting, or decrease the step size.