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 highly 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)
Grab the all-in-one Jupyter Notebook (YAMCHA_tutorial_ch7) from GitHub.
Run the first 3 Cells. Comment out or skip Cell #4. In this Chapter, we'll manually setup the chemical mechanism.
Proceed to Cell #4, and manually setup the mechanism:
ChemMech_MultiPhase = pd.DataFrame({'R0': ['SO2 = SIV_aq', 'kt_g2a_SO2'], # gas-to-aqueous transfer of SO2
'R1': ['SIV_aq = SO2', 'kt_a2g_SO2'], # aqueous-to-gas transfer of SO2
'R2': ['H2O2 = H2O2_aq', 'kt_g2a_H2O2'], # gas-to-aqueous transfer of H2O2
'R3': ['H2O2_aq = H2O2', 'kt_a2g_H2O2'], # aqueous-to-gas transfer of H2O2
'R4': ['SIV_aq + H2O2_aq = SVI_aq', 'k_SIV_H2O2'], # aqueous-phase rate coefficient
})
ChemMech_MultiPhase.index = ['reaction', 'rate_coefficient']
ChemMech_MultiPhase.T
here "_aq" stands for stuff in the aqueous-phase. SIV_aq represent all forms of S(IV) in the aqueous-phase, including SO2-H2O, HSO3-, and SO3^2-. This simple mechanism has 5 reactions and 5 species. Again if loaded correctly, it'll print a summary of the mechanism under this cell.
Also we need to provide a few user-defined functions. The first one calculates the phase-transfer coefficients:
def PhaseTransfer_kt_s(in_Heff_M_atm, in_lwc_vol_vol, in_T_K, in_radius_cm, in_dg_cm2, in_MW_gas_kg_mol, in_alpha):
R = 0.08205
thermalspeed_cm_s = 100.0 * np.sqrt(8.0 * 8.314 * in_T_K / np.pi / in_MW_gas_kg_mol)
kt = 1.0/(in_radius_cm*in_radius_cm/3./in_dg_cm2 + (4./3.)*in_radius_cm/thermalspeed_cm_s/in_alpha)
k_g2a, k_a2g = kt * in_lwc_vol_vol, kt / in_Heff_M_atm / R / in_T_K
return k_g2a,k_a2g
This function returns the gas-to-aqueous (k_g2a) and aqueous-to-gas (k_a2g) transfer coefficients. Units are both in s^-1.
We also need a function to calculate the pH-dependent effective Henry's law constant:
def get_Heff_M_atm(in_T_K, in_pH, in_kh298, in_dhr, in_k1_298, in_dh1r, in_k2_298, in_dh2r, in_if_acid):
aH = 10.0**(-1.*in_pH)
kh = in_kh298 * np.exp(in_dhr*(1./in_T_K-1./298.))
k1 = in_k1_298 * np.exp(in_dh1r*(1./in_T_K-1./298.))
k2 = in_k2_298 * np.exp(in_dh2r*(1./in_T_K-1./298.))
if in_if_acid: heff = kh * (1. + k1/aH * (1. + k2/aH)) # acid formulation
else: heff = kh * (1. + k1*aH/k2) # base formulation
return heff
Lastly, the reaction rate coefficient of R4: S(IV, aq) + H2O2 (aq) is a strong function of pH. This is given in the Seinfeld and Pandis (2016) book. 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.
def get_k_SIV_H2O2(in_pH, in_T_K):
Hp = 10**(-1*in_pH)
K1_SO2 = 1.3E-2 * np.exp(1960*(1./in_T_K - 1./298.))
K2_SO2 = 6.6E-8 * np.exp(1500*(1./in_T_K - 1./298.))
f_HSO3m = K1_SO2*Hp / (Hp*Hp + K1_SO2*Hp + K1_SO2*K2_SO2)
return 7.45E+7 * Hp / (1 + 13*Hp) * f_HSO3m
Before going too far, let's do a quick sanity check: plot the pH dependency of the S(IV,aq) + H2O2 (aq) reaction. The plot is attached to the right. You can compare that to the plot on Page 311 in Seinfeld and Pandis (2016).
At this point you're already familiar with the model.
# --- time stamps
total_run_time_days = 3./24. # how long you want to run the model?
output_frequency_seconds = 60. # how often do you save the output?
...
# ================================
# Which chemical mechanism to use?
# ================================
in_ChemMech = ChemMech_MultiPhase
...
# --- this is the solar zenith angle for the j-value parameterization (MCM)
SZA = 64. * np.pi / 180.
# --- environmental conditions
temp = 298 # Kelvin
press = 101325. # Pa
RH = 60.
...
# --- initial condition. need to convert to molec/cm3
conc_init['SO2'] = (5.) *press*(6.0232E+8)/(8.314*temp)
conc_init['H2O2'] = (1.) *press*(6.0232E+8)/(8.314*temp)
Note that SZA and RH are not really used in this work..
Now set a few parameters that are needed for the phase-transfer and aqueous reactions:
# --- set aqueous conditions
a_pH = 2. # pH of the aqueous-phase
a_LWC = 1e-7 # liquid water content
a_radius_cm = 0.005 # radius of the droplets
dg_cm_s = 0.1 # gas diffusivity
aq_rate_conv = 1.e3/(6.0232e+23*a_LWC)
# --- set Henry's law
Heff_H2O2_M_atm = get_Heff_M_atm(temp, a_pH, 8.70e+04, 7320.,2.2e-12,-3730., 0., 0., True)
Heff_SO2_M_atm = get_Heff_M_atm(temp, a_pH, 1.23, 0., 1.3E-2, 1960, 6.6E-8, 1500, True)
# --- set phase-transfer coefficients
kt_g2a_H2O2, kt_a2g_H2O2 = PhaseTransfer_kt_s(Heff_H2O2_M_atm, a_LWC, temp, a_radius_cm, dg_cm_s, 34.0147e-3, 0.18)
kt_g2a_SO2, kt_a2g_SO2 = PhaseTransfer_kt_s(Heff_SO2_M_atm, a_LWC, temp, a_radius_cm, dg_cm_s, 64.066e-3, 0.12)
# --- set aqueous-phase rate coefficients
k_SIV_H2O2 = get_k_SIV_H2O2(a_pH, temp) *aq_rate_conv
A few notes: a_LWC (liquid water content) and a_radius_cm (droplet radius) are typical values for cloud droplets.
pH affects a number of things here. 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.
Unit conversion (VERY IMPORTANT): Note the term aq_rate_conv here (highlighted in red). 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?
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 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.
Firstly run the model with a_pH = 2 (very acidic). Then change a_pH to 8 (basic) and run the model again. The results are shown on the right.
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.
You may have noticed that at pH = 2, gaseous SO2 drops "slower" than H2O2 does. But when pH = 8, 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. In Cell #8, let's shut off the H2O2+S(IV) reaction:
k_SIV_H2O2 = 0
Set the output time interval to 10 seconds, and the total run time to 0.02 day:
total_run_time_days = 0.02 # how long you want to run the model?
output_frequency_seconds = 10. # how often do you save the output?
Then change initial concentrations of SO2 and H2O2 both to 2 ppb
conc_init['SO2'] = (2.) *press*(6.0232E+8)/(8.314*temp)
conc_init['H2O2'] = (2.) *press*(6.0232E+8)/(8.314*temp)
Now in this case we don't have aqueous-phase reaction, so phase-transfer only. The reults are shown on the right. 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).
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.
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.