4. Now let's scroll all the way to the top of the Notebook, and run the very first cell:

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import time

from scipy.integrate import solve_ivp,cumulative_trapezoid

This cell imports Python packages/libraries that will be used throughout this Notebook. Remember how to run a cell? You either hit CTRL + ENTER, or click that button at the upper left corner of the cell. Importing these will not generate any outputs, so you won't see anything except a green checkmark at the upper left corder indicating the cell was run successfully.

5. Now move down to the second cell (Cell #2), which has several utility functions to load/process the KPP files generated from the MCM website. We don't need these at this moment but let's run this cell anyway.

6. Move further down to Cell #3, which has a bunch of utility functions that can parse the chemical mechanism and translate it into a format that is needed by the ODE solver. This is the heart of this box model:

• pre_process(...): this processes the mechanism, then generate the derivatives (dC/dt). This is required by the ODE solver.

• get_jac(...): this calculates the Jacobian matrix. Many ODE solvers can estimate Jacobian matrix but this will slow down the solver dramatically. So this function here will make the solver must faster!

• mcm_jval(...): photolysis frequencies, which are based on MCM (parameterized as a function of solar zenith angle)

• pd_mcm_kinetics(...): generic rate coefficients, also from MCM.

Run this cell!

Provide the model with a chemical mechanism:

7. Move down to Cell #4. Note that this entire cell is commented out (by adding # at the beginning a line). So even if you run this cell, nothing will happen. But this cell demonstrates how to load a KPP file, which we will use this in the next Chapter! So for now let's skip this cell. Or you can run it.

8. Now on Cell #5. A Pandas Dataframe is defined here, with a list reactions and the corresponding rate coefficients. In this Chapter we will play around with a simple nighttime NOx chemistry, consisting of five species (NO, NO2, NO3, N2O5, and O3). 

ChemMech_simple = pd.DataFrame({'R0': ['NO + O3 = NO2', '1.4e-12*EXP(-1310/temp)'],

                                'R1': ['NO2 + O3 = NO3', '3.5e-17'],

                                'R2': ['NO + NO3 = 2*NO2', '2.6e-11'],

                                'R3': ['NO2 + NO3 = N2O5', '1.2e-12'],

                                'R4': ['N2O5 = NO2 + NO3', '4.46e-2'],

                                })

ChemMech_simple.index = ['reaction', 'rate_coefficient']

ChemMech_simple.T

R0-R5 are the reaction indices. The rate coefficients are number or generic math expressions in string format (the model can recognize most of them). You can also provide user-defined funtions to calculate these rate expressions, which we will cover in later Chapters.

Now run this cell as well! You should see the formatted DataFrame being printed as a table right below the cell.

10. Now Cell #6 is the actual model, which is a bit long but I'll walk you thru.

Firstly, there are two functions defined here, and you don't need to touch these at this moment.

def dcodt(t,conc):

    ...

def jac(t,conc):

    ...

Then you'll configure a few important things about the model.

total_run_time_days = 12./24.          # how long you want to run the model?

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

The comments tell you what they mean. The units are also given. So as written, the model runs for 12 hours (0.5 day), and we ask the model to output results every 600 seconds.

Then a bunch of other stuff you don't need to worry about. Then you'll see a bunch of user inputs.

# --- this is the solar zenith angle for the j-value parameterization (MCM)

SZA = 90. * np.pi / 180.

# --- environmental conditions

temp = 298      # Kelvin

press = 101325. # Pa

RH = 6.

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['NO'] = (10.)  *press*(6.0232E+8)/(8.314*temp)

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

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

SZA is solar zenith angle, and here it's set to 90 degrees (and converted to radians). Note that the chemical mechanism you defined in Cell #5 (Step 8 above) does not involve any photolysis reactions, therefore SZA there will have zero impact on the model.

conc_init is an array that carries initial concentrations of all species (initial condition), and the unit must be molec/cm3. Here I'm setting the initial concentrations of NO, NO2, and O3 to 10 ppb, 5 ppb, and 30 ppb, respectively, which are then converted to the required unit which is molec/cm3.

As you may tell, this experiment is configured for 298 K and 101325 Pa, with 6% relative humidity (RH).

11. Scroll down and find a big for loop, which loops over all time steps.

for i,t in enumerate(tout_s):

    # --- update photolysis frequencies

    jval = mcm_jval(SZA)

    # --- update generic rate coefficients

    if ('HO2' in species_list) & ('NO' in species_list):

        mcm_kinetics = pd_mcm_kinetics(temp,M,N2,O2,H2O,

                                       conc[max(i-1,0),species_list.index('NO')], # only needed for K16ISOM

                                       conc[max(i-1,0),species_list.index('HO2')], RO2)

    else: mcm_kinetics = pd_mcm_kinetics(temp,M,N2,O2,H2O)

    for k in k_label: exec("%s = np.array(mcm_kinetics['%s'])[0]" % (k,k))

    # --- need to update all environmental variables, mcm kinetics, jval before eval()

    for n in range(n_reaction): rate_coeff[n] = eval(rate_coeff_express_2eval[n])

    # --- now solve the ODE. Use BDF since it's mostly stiff

    #     accuracy controlled by relative/absolute tolerance (rtol,atol)

    #     integrate step by step, since some env variables may change over the course of time

    if i==0: conc[i,:] = conc_init

    else: conc[i,:] = solve_ivp(dcodt,tout_s[0:2],conc[i-1,:],t_eval=[output_frequency_seconds],

                                jac=jac, method='BDF', rtol=1e-5,atol=1e+2).y[:,0]

As you can see, at each time step, the model firstly updates the photolysis and all other rate coefficients, then update solve the ODE using solve_ivp (ivp stands for initial value problem). Note that in this simple example, the environmental conditions remain unchanged. If you want to have different environmental conditions, e.g., temperature, pressure, RH, and water vapor (H2O), you need to update these at the very beginning of each time step, because these may affect the rate coefficients.

13. Now Cell #7 provides some high-level overview for post analysis. A few things you need to know:

• species_list: exactly what it looks like.. All species are indexed, and if you want to find the species index of a species, e.g. NO, this is what you need: species_list.index('NO')

• reactions: a list of all reactions. This is equivalent to ChemMech_simple.T['reaction']. 

• conc: this is a Pandas Dataframe carrying the time series of all species, and the unit is molec/cm3. You may access the time series of a species by specifying its species name, e.g. conc['O3'] is the time series of O3 concentration.

• reaction_rates: this is a Pandas DataFrame for the time series of all reaction rates, and the unit is molec/cm3/s. Reaction rate is defined as the product of rate coefficient and the concentrations of all reactants.

And if you're particularly nerdy like I am, there are a few other things that might be of interest to you:

• rates_2eval: this is a list of strings, with each string represent the math expression of the rate of a reaction.

• derivatives_2eval: this is a list of strings, with each string represent the math expression of the derivative of a species, e.g. dC/dt.

• jac_2eval: this is a list of strings, with each string represent the math expression of the jacobian matrix:

This is it! You now have a 0D box model with simple nighttime NOx chemistry!

High five!