1. Benchmark

• Grab the all-in-one Jupyter Notebook (YAMCHA_tutorial_ch8) 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 #5, let's code up the analytical solution of the condensational growth of aerosols in Chapter 13 of the Seinfeld & Pandis book, which we'll use later for evaluation. IMPORTANT: this simplified solution (13.2.2) has two assumptions: (i) continuum regime, i.e. Kn is very small or f(Kn, alpha)=1; and (ii) unity mass accommodation coefficient, i.e. alpha =1.

• Firstly, let's try to reproduce Figure 13.3 in the Seinfeld & Pandis book. Two key functions here:

# --- this is not dN/dlogDp, but rather dN/Dp, so the unit is a bit counter intuitive

#     but integrate this over a Dp range should get N

def lognormal_dist_dN_dDp(in_Dp, in_Dp_mean, in_sigma, in_N0):

    return in_N0/np.log(in_sigma)/in_Dp/np.sqrt(2.*np.pi) * np.exp(-0.5*(np.log(in_Dp/in_Dp_mean)**2.)/np.log(in_sigma)/np.log(in_sigma))

# --- growth of a lognormal aerosol distribution: analytical

#     note this is dN/dDp, rather than dN/dlogDp

def lognormal_aerosol_growth_func(Dp_um, t_s, p_Pa, peq_Pa, Dg_cm_s,Mw_g_mol,N0, in_Dp_mean, in_sigma):

    A_cm2_s = 4.*Dg_cm_s*Mw_g_mol*(p_Pa-peq_Pa)/8.314/temp/(rho_g_cm3*1000000.)

    Dp2_2At = ((Dp_um*1e-4)**2.-2.*A_cm2_s*t_s)/1e-4/1e-4            # unit of this term: um2 (converted from cm2)

    fk = np.log(np.sqrt(Dp2_2At)/in_Dp_mean)**2.                     # this term is dimensionless

    out = Dp_um/Dp2_2At * N0/np.sqrt(2.*np.pi)/np.log(in_sigma) * np.exp(-0.5*fk/np.log(in_sigma)/np.log(in_sigma))

    return out

The first function, lognormal_dist_dN_dDp, is for dN/dDp (rather than dN/dlogDp). If you integrate the dN/dDp over aerosol diameter (Dp), you get the aerosol number concentration. The second function, lognormal_aerosol_growth_func, is just the analytical solution at time t, which is given in the book.

• Once other aerosol parameters are set based on the book, run this cell and you should get the figure on the left, while the figure on the right is a screenshot from the book. Looks good to me!

5. Mass Transfer Coefficient

• Now proceed to the main model cell. In the very beginning, a function is defined to calculate the mass transfer coefficient. Note that to calculate Kn (Knudsen number), we need mean free path, which has several formulations. The formulation used here is simple but probably is not the best suited for "larger" molecules in the air. More sophisticated formulations are available but some will require specific information about the molecule itself. Overall there's probably not a big differenc and mean free path is on the order of 10^-8 to 10^-7 m.

def kmt_gas_aerosol(in_radius_cm, in_Dg_cm2_s, in_N_cm3, in_alpha,

                    in_temp, in_press, in_MW_gas_kg_mol):

    MW_air_kg_mol = 0.029

    DynamicViscosity_kg_m_s = 1.715747771E-5 + 4.722402075E-8 * (in_temp-273.15) \

                            - 3.663027156E-10 * ((in_temp-273.15)**2.0)    \

                            + 1.873236686E-12 * ((in_temp-273.15)**3.0)    \

                            - 8.050218737E-14 * ((in_temp-273.15)**4.0)

    thermalspeed_m_s = np.sqrt(8.0*8.314*in_temp/np.pi/in_MW_gas_kg_mol)

    rho_air_kg_m3 = in_press*MW_air_kg_mol/8.314/in_temp

    MeanFreePath_m = 2.*DynamicViscosity_kg_m_s/rho_air_kg_m3/thermalspeed_m_s

    Kn = (MeanFreePath_m*100.) / in_radius_cm

    f = 0.75*in_alpha*(1.+Kn)/( Kn*(1.+Kn) + 0.283*in_alpha*Kn + 0.75*in_alpha )

    # --- NOTE: if evaluate the model using the Seinfeld & Pandis book, Chapter 13.2.2

    #           there are two conditions: (i) continuum regime, i.e. Kn is small and f=1

    #           (ii) unity mass accommodation coefficient, i.e. alpha=1

    f = 1.

    return 4.*np.pi*in_radius_cm*in_Dg_cm2_s*in_N_cm3*f

• Much of this cell should look pretty familiar to you! Several important things:

  1. In this chapter we're not covering nucleation or coagulation, or the deposition of aerosols, therefore the number concentration of aerosols remain unchanged.

  2. Kelvin effect accounts for the impact of the curvature of aerosols on the vapor pressure. It matters more for small aerosols. Normally this should be turned on (will require surface tension), but again for the sake of evaluating using the Seinfeld and Pandis book, let's turn if off for now.

  3. Mass accommodation coefficient: 0.1 should work for most organic compounds (e.g. Zaveri et al. 2008). But the analytical solution in the Seinfeld and Pandis book assumes unity accommodation coefficient, so for now let's set it to 1.

  4. Another important aspect of the analytical solution in the Seinfeld & Pandis book is that, the "driving force", i.e. the difference between the partial pressure of the condensing gas and the equilibrium vapor pressure (p - peq) is set to a constant value. This is equivalent to the gas-phase concentration of the condensing gas is held constant. This is done by adding one line in the derivative function (dC[i]/dt = 0, see the line highlighted in red below).

def dcodt(t,conc):

    rate = np.zeros(n_reaction)

    dcodt_out = np.zeros(n_spc)

    for n in range(n_reaction): rate[n] = eval(rates_2eval_compiled[n])

    for n in range(n_spc): dcodt_out[n] = eval(derivatives_2eval_compiled[n])

    dcodt_out[species_list.index('C_gas')] = 0.0

    return dcodt_out

Under the hood, if the derivatives are generated by the model, then the Jacobian matrix is automatically taken care of. Now that you *touched* derivatives, in theory you should edit the Jacobian matrix accordingly. Quiz: in this case, do you need to manually edit the Jacobian matrix?

You can try running the model with and without this line. Without this line, the model runs normally and you should see the gas-phase concentratinon (C_gas) decays away as the aerosols grow. But with this line, the gas-phase concentration will remain constant, so its time series will remain flat!

• 5. Aerosol growth is accounted for by adding the "new mass" (due to condensation) to the existing aerosols. This is done in the following code snippet (towards the end of the big time loop):

for ibin in range(len(init_bin_center_um)):

        # --- get "new" mass per particle

        if i==0: dC_aer[ibin] = conc[i,species_list.index('C_aer'+str(ibin))]

        else: dC_aer[ibin] = (conc[i,species_list.index('C_aer'+str(ibin))]-conc[i-1,species_list.index('C_aer'+str(ibin))])

        m_per_particle_g[ibin] = ((4./3.)*np.pi*(0.5*new_Dp_um[ibin]*1e-4)**3.)*rho_g_cm3 + dC_aer[ibin]/init_bin_N_cm3[ibin]/6.0232e+23*MW_g_mol

        # --- get new aerosol size

        new_Dp_um[ibin] = 2.*1e+4*((m_per_particle_g[ibin]/rho_g_cm3)*(3./4.)/np.pi)**(1./3.)

        # --- correct for Kelvin effect

        if if_Kelvin: KelvinFactor = np.exp(4.*MW_g_mol*surface_tension_water_N_m/8.314/temp/(rho_g_cm3*1000000.)/(new_Dp_um[ibin]*1e-6))

        else: KelvinFactor = 1.

        SatVapConc_molec_cm3[ibin] = KelvinFactor * SatVapPres_Pa/8.314/temp *6.0232e+23/1000000.

        # --- update kinetics

        kmt[ibin] = kmt_gas_aerosol(0.5*new_Dp_um[ibin]*1e-4, 0.1, init_bin_N_cm3[ibin], mass_accommodation_coeff, temp, press, MW_g_mol/1000.)

        kmt_Cstar_Cgas[ibin] = kmt[ibin] * SatVapConc_molec_cm3[ibin]

        kmt_Cstar_Caer[ibin] = -1.*kmt[ibin] * SatVapConc_molec_cm3[ibin]

You may wonder why this is done in the end of the time loop. This is because the "new mass" is calculated as the difference of the concentrations of a given aerosol bin between two time steps, i.e. conc[i,n]-conc[i-1,n], where n is the species index and i is the time index. If you do this before calling the ODE solver, the "new concentration" (conc[i,n]) has not been computed yet.

NOTE: One important caveat of this approach is that, we're essentially assuming the aerosol sizes remain unchanged during the entire time step. Technically it shouldn't be like this, i.e. kmt is a function of size and hence the concentrations too. You can move these to dcodt, i.e. calculating all the rate coefficients (e.g. kmt) using the "local" concentrations at each solver internal step. But this simplification doesn't seem to affect the result much, as long as the time step is smaller than tens of seconds. Let's stick to this simple solution for now.

• Let's run the model! In the end a few plots are generated. Also the analytical solution is obtained by using the funtion that we coded before (lognormal_aerosol_growth_func), and the inputs are the same as the model. You can see that this model can reproduce the condensational growth in the Seinfeld and Pandis book very well! The dashed lines are the "ground truth", i.e. the analytical solutions calculated using Equation 13.27 in the book, and the solid lines are the results from our kinetic model!