1. Prep

• You already got the KPP file (mcm_export_isoprene.eqn) from the previous Chapter (in case not, here's the GitHub link). We can use the same file in this Chapter.

• Grab the all-in-one Jupyter Notebook (YAMCHA_tutorial_ch4) from GitHub. Dump the KPP file in Colab.

• Run the first 3 Cells to load the libraries & utility functions. Then setup Cell #4 to load the KPP file. If it's loaded successfully, a summary of the mechanism will be displayed under Cell #4.

2. Mixing Layer, emissions & depositions

• From the box model perspective: we put all chemical species in a box and, we assume everything is well mixed inside this box. In real world if you randomly draw an imaginary box, the chance that things are well mixed inside the box is slim. But sometimes it's not the worst assumption. For instance, during daytime when the ground is heated, convection and turbulence drive the mixing of tracers, extending to a certain height depending on stability. Within this layer, tracers (water vapor, trace gases, etc) are somewhat well mixed, and this layer is often referred to as "mixing layer", and the depth is called "mixing layer depth", which often has clear diurnal variabilities. Well this is my "layman's explanation". Consult a proper textbook if you want to learn more.

• Now come back to the box model: The *size* of the box didn't matter in our previous examples, because the model was operating on the volume basis. But if we want to mimic the emission and deposition processes, the size of the box matters. Emissions and depositions are often expressed as fluxes (in molec/cm2/s), i.e. the number of molecules crossing the interface (ground) and enter / leaving the box. If we (have to) assume everything is well mixed inside the box, then it makes sense to set the height of the box to the "mixing layer depth", so the concentration changes caused by emissions/depositions are evenly distributed vertically. In other words: fluxes (in molec/cm2/s) need to be converted to volume rates (molec/cm3/s).

• Therefore, for emissions, we can do:

EmissionRate (molec/cm3/s) = EmissionFlux (molec/cm2/s) / MixingLayerHeight (cm)

From time to time, you may find emission fluxes measurements for different species, under different conditions.

• Deposition is treated slightly different, as often times deposition velocity (in cm/s) is reported. Therefore:

DepositionRate (molec/cm3/s) = (DepositionVelocity (cm/s) / MixingLayerHeight (cm)) * Concentration (molec/cm3)

• The height of a well mixed mixing layer varies: usually lower at night and higher during daytime. Let's assume it's 500 m at night, and reaches 1500 m at local noon. Again this is vaguely representative for a good sunny day in the mid/low-latitude.

4. Dry Deposition: Direct Approach (Slightly More Complicated)

• Dry deposition velocities (in cm/s) are often reported to describe the dry deposition of gas-phase species. Divide the dry deposition velocity by the mixing layer depth, you got a first-order loss dry deposition rate (in s^-1) which often is a good estimation of dry deposition.

• There are several ways to implement dry deposition in a box model. For instance, you can directly modify the derivatives and add the deposition terms. Since this box model also provides Jacobian matrix, if you decide to manually modify the derivatives, you also need to modify the Jacobian matrix accordingly. I'll very briefly describe the method here. Let's say you want to add dry deposition for compound X, and it's species index is n.

dcodt[n] = dcodt[n] - DryDepVelocity_cm_s/(MLH_m * 100) * conc[n]

jac[n,n] = jac[n,n] - DryDepVelocity_cm_s/(MLH_m * 100)

where dcodt is its derivative and jac[n,n] is its Jacobian matrix term, DryDepVelocity_cm_s is its dry deposition velocity, MLH_M is the mixing layer height in meters, and conc[n] is the local concentration of this species. The dcodt term should be fairly straightforward. The jac terms are essentially partial derivatives but in this case only jac[n,n] is non-zero.

• As you can see this is a little complicated, so I'm going to introduce a simpler approach. Since this dry deposition implementation operates as a first-order loss, we can add a *fake* first-order loss reaction in the reaction list to represent dry deposition, and the mechanism parser will automatically process it and get the dcodt and jac terms behind the scene!

5. Dry Deposition: Add A Fake Reaction (Easier)

• As you can see this is a little complicated, so I'm going to introduce a simpler approach. Since this dry deposition implementation operates as a first-order loss, we can add a *fake* first-order loss reaction in the reaction list to represent dry deposition, and the mechanism parser will automatically process it and get the dcodt and jac terms behind the scene!

• Proceed to Cell #7. I'm adding three fake chemical reactions to represent the dry deposition of O3, NO2, and HNO3, and their first-order dry deposition rates are provided as rate coefficients. We will get the values for these rate coefficients later.

ChemMech_dep = pd.DataFrame({'D0': ['O3 = ',  'dep_rate_O3_sec', None],

                             'D1': ['NO2 = ', 'dep_rate_NO2_sec', None],

                             'D2': ['HNO3 = ', 'dep_rate_HNO3_sec', None],

                             'E0': [' = C5H8', 'emission_rate_isoprene_molec_cm3_s', None],

                             })

ChemMech_dep.index = ['reaction', 'rate_coefficient', 'RO2']

ChemMech_dep.T

Ignore the "None" and "RO2" stuff here which are useless but will make our life easier later. Same as any user-defined rate coefficients, you have to define "dep_rate_O3_sec" etc later.

• Run this cell. If all goes right, you should see a small table printed under this cell, with these three *fake* chemical reactions.

• Now proceed to Cell #8, which merges these three *fake* reactions with the main chemical mechanism. 

ChemMech_merge = pd.concat([ChemMech,ChemMech_dep],axis=1)

ChemMech_merge.T

If all okay, you should see the summary of the merged mechanism printed below, which is the original mechanism with the three "dry deposition reactions" listed in the end. 

• Don't forget: when running the box model, you must load this merged mechanism into the model! The relevant code snippets in Cell #9 look like these:

# ===================

# Model configuration

# ===================

# --- time stamps

total_run_time_days = 1.               # how long you want to run the model?

output_frequency_seconds = 1800        # how often do you save the output?

tout_s = np.arange(0,total_run_time_days*86400,output_frequency_seconds)  # output time stamps, unit: second

# ================================

# Which chemical mechanism to use?

# ================================

# in_ChemMech = ChemMech

in_ChemMech = ChemMech_merge          # make sure you load the merged mechanism

species_list, reactions, rates_2eval, derivatives_2eval, rate_coeff_express_2eval, jac_2eval, RO2_sum_2eval = pre_process(in_ChemMech)

...

As you noticed, we also set the total run time to 1 day, and output results every 1800 seconds (half an hour).

6. Emissions

• Emission rates are different from dry deposition rates. Emission rates are independent of concentrations, therefore are similar to zero-order reactions. Zero-order reactions are not very common in atmospheric chemistry. But the parser can certainly handle zero-order reactions!

ChemMech_dep = pd.DataFrame({'D0': ['O3 = ',  'dep_rate_O3_sec', None],

                             'D1': ['NO2 = ', 'dep_rate_NO2_sec', None],

                             'D2': ['HNO3 = ', 'dep_rate_HNO3_sec', None],

                             'E0': [' = C5H8', 'emission_rate_isoprene_molec_cm3_s', None],

                             })

ChemMech_dep.index = ['reaction', 'rate_coefficient', 'RO2']

ChemMech_dep.T

E0 here is the fake zero-order reaction that is added to represent the emissions. In this reaction, there is no reactant, and the "product" is just isoprene, and the rate coefficient has the unit of molec/cm3/s.

• Alternatively, you can also implementing emissions directly, by adding an emission term in the derivatives.

      dcodt[n] = dcodt[n] + Emission_rate_molec_cm3_s

where Emission_rate_molec_cm3_s is just the emission rate (in molec/cm3/s) for a compound with species index of n. You may wonder why do we not need a Jacobian term as the deposition does? Again Emission_rate_molec_cm3_s is independent of concentration, so all its Jacobian terms will be zero. If you want to implement emissions in this eay, you can touch dcodt(...):

# ===================================================================

# this calculates the derivatives: dC/dt

# I'm using eval() to get the derivatives from the string expressions

# ===================================================================

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('C5H8')] += emission_rate_isoprene_molec_cm3_s

    return dcodt_out

The line highlighted in red is what you need. Later we'll asign values to this variable (emission_rate_isoprene_molec_cm3_s).

• IMPORTANT: if you add the emissions as "fake" zero-order reactions, you don't need to touch dcodt! Otherwise you're double-counting emissions! In this Chapter, isoprene emission will be added using the simpler way, i.e. by adding a "fake" reaction.

7. Emission Fluxes and Deposition Velocities

• Emission fluxes (molec/cm2/s) are often reported which is also commonly used in models. Depositions, on the other hand, are often represented by deposition velocities (cm/s). We talked before how these are implemented in the box model. The following chunk of codes sets a number of user-defined inputs, with emission fluxes and dry deposition velocities defined in the end (highlighted in red).

# --- environmental conditions

temp = 298      # Kelvin

press = 101325. # Pa

RH = 60.

M = press*6.0232E+17/8.314/temp  # air density. molec/cm3

N2,O2 = 0.78*M, 0.21*M

H2O = ((RH/100.)*6.1078*np.exp(17.269*(temp-273.3)/temp)*100./8.314/temp)*6.0232e+23/1000000.

RO2 = 0         # total RO2 concentration (molec/cm3). MCM needs this. not really saved

RO2_ts = []     # this is a MCM thing

# --- initial condition. need to convert to molec/cm3

conc_init['O3'] = (20.)     *press*(6.0232E+8)/(8.314*temp)

conc_init['NO'] = (5.)      *press*(6.0232E+8)/(8.314*temp)

conc_init['NO2'] = (10.)    *press*(6.0232E+8)/(8.314*temp)

conc_init['CO'] = (100.)    *press*(6.0232E+8)/(8.314*temp)

conc_init['CH4'] = (1800.)  *press*(6.0232E+8)/(8.314*temp)

# --- set emissions and depositions

max_emis_flux_isop_molec_cm2_s = 1e+12

dep_vel_O3_cm_s = 0.2

dep_vel_NO2_cm_s = 0.2

dep_vel_HNO3_cm_s = 1.

Here we'll set emission for isoprene only, and the maximum daily emission flux of isoprene is set to 1e+12 molec/cm2/s. This is a reasonable guestimate, although real measurements may be higher or lower. We'll also set deposition velocities for O3, NO2, and HNO3. For O3 and NO2, 0.2 cm/s is a reasonable guestimate as well. For HNO3, it's a bit more "sticky" so it has a stronger dry deposition velocity, so 1 cm/s is not a bad guess. Feel free to use other values if you want!

8. Time-Dependent Parameters

• Okay now we'll take care of the time-dependent parameters: SZA, mixing layer height, also emissions and depositions which are affected by mixing layer height. These are done inside the bit time-loop of course.

for i,t in enumerate(tout_s):

    # --- set time of day

    hour_of_day = (t%86400)/3600.

    # --- set time-dependent variables

    SZA_deg = (90. + (20.-90.) * diurnal_profile_sine(hour_of_day)) 

    SZA = SZA_deg * np.pi / 180.

    MLH_m = 500. + (1500.-500.) * diurnal_profile_sine(hour_of_day)

    # --- set time-dependent emissions and depositions

    emission_rate_isoprene_molec_cm3_s = diurnal_profile_sine(hour_of_day)*max_emis_flux_isop_molec_cm2_s/(MLH_m*100.)

    dep_rate_O3_sec = dep_vel_O3_cm_s/(MLH_m*100.)

    dep_rate_NO2_sec = dep_vel_NO2_cm_s/(MLH_m*100.)

    dep_rate_HNO3_sec = dep_vel_HNO3_cm_s/(MLH_m*100.)

    ...

In Cell #9, locate the big time loop. The chunk of codes above firstly sets the "hour of day", and the sine-shaped function we defined before is used to scale SZA and MLH_m (mixing layer height, in meters). 

• SZA is set to 90 degrees at night and 20 degrees at local noontime.

• MLH_m is set to 500 m at night and a daily maximum of 1500 m at noon. 

• Isoprene has a daily max defined by max_emis_flux_isop_molec_cm2_s, which is then converted to isoprene emission rate (molec/cm3/s) by scaling to MLH_m, and the emission rate (again in molec/cm3/s) is passed into dcodt and eventually the ODE solver.

• Dry deposition velocities are converted into 1st-order dry deposition rates (in s^-1), which are then passed into the *fake* deposition reaction rates and eventually sent into ODE solver.

10. Effect of Deposition?

• One of the coolest things we can do with a model is that, we can turn a process on and off and see how it affects the results.

• Let's firstly plot O3, NO, NO2, and HNO3 with deposition on, which is shown below (on the left). We can see that in the very beginning (at night), O3 rapidly reacts with NO producing NO2. NO2 and O3 then start to gradually decrease until sunrise at 6AM. Question: where does all the nitrogen go during this period of time?

• Starting from 6AM, with the presence of sunlight, photolysis kicks in. That's why O3 starts to be produced, and NO2 is converted to NO as well. In the meantime, NOx is oxidized into HNO3 which is accumulated in the model.

• Now if you turn off the dry depositions, by setting the dry deposition velocities to zero:

dep_vel_O3_cm_s = 0      #0.2

dep_vel_NO2_cm_s = 0     #0.2

dep_vel_HNO3_cm_s = 0    #1.

and run the model again. You should get the other plot below (on the right). You can see everything is pretty similar, except O3, NO2, and HNO3 are slightly higher. That's how dry deposition affects these species!

11. A few notes: spherical chicken in the vacuum

• This is so far the longest Chapter, and thanks for sticking around. But this Chapter also demonstrates a number of important things:

1. How emissions and depositions are implemented in a box model

2. How to set time-dependent parameters

3. How to manually edit chemical mechanism: you can always directly edit the KPP file. But here we show that you can create a small mechanism and merge that with the original mechanism.

• This model (as is) is only conceptually useful, as the real ambient environment is way more complicated. Many physical and chemical processes are not covered in this chapter. 

• For example, in the evening, the mixing layer shrinks, forming a inversion layer that is essentially de-coupled with the surface, this is often called "residual layer". Chemistry in the residual layer can be quite different from the near-surface. In the early morning, the increase of mixing layer height would lead to the enhanced exchange with this residual layer, this processes is often referred to as "entrainment". To describe this process, one has to either model the boundary layer dynamics as well as the chemistry in the residual layer, or, to prescribe the chemical characteristics in the residual layer (if you have measurements). Technique-wise, it's fairly easy to add this (check the Seinfeld and Pandis textbook).

• Advection also plays an important role in the real world. Again to do this, a thorough understanding of the "background air" is needed. Equations are out there and it's quite straightforward.

• Last but not least, people usually would "spin-up" the model for a certain period of time. You can simply run the model for longer periods of time and ditch the first few days.