Computer exercise 2

Problem statement

Consider a simplified $\ce{O3}$-$\ce{NO_x}$-VOC system in which the VOC is represented by a single molecule, propene ($\ce{C3H6}$). Writing this as a generalized mechanism (propene written as $\ce{RH}$), reactions and rate constants are shown in the following table (Seinfeld and Pandis, 2006):

\begin{align} \newcommand{\conc}{\text{cm$^{3}$ molec$^{-1}$ s$^{-1}$}} & \text{Reaction} & \text{Rate constant (298 K)} \\ \hline 1\quad & \ce{RH}+\ce{OH} \stackrel{\ce{O2}}{\to} \ce{RO2}+\ce{H2O} & 26.3 \times 10^{-12} \conc\\ 2\quad & \ce{RO2}+\ce{NO} \stackrel{\ce{O2}}{\to} \ce{NO2}+\ce{R^'CHO}+\ce{HO2} & 7.7 \times 10^{-12} \conc\\ 3\quad & \ce{HO2}+\ce{NO} \to \ce{NO2}+\ce{OH} & 8.1 \times 10^{-12} \conc\\ 4\quad & \ce{OH}+\ce{NO2} \stackrel{\ce{M}}{\to} \ce{HNO3} & 1.1 \times 10^{-11} \conc \text{(at 1 atm)}\\ 5\quad & \ce{HO2}+\ce{HO2}\to \ce{H2O2}+\ce{O2} & 2.9 \times 10^{-12} \conc\\ 6\quad & \ce{RO2}+\ce{HO2}\to \ce{ROOH}+\ce{O2} & 5.2 \times 10^{-12} \conc\\ 7\quad & \ce{NO2}+h\nu \stackrel{\ce{O2}}{\to} \ce{NO}+\ce{O3} & \text{typical} \approx 0.015 \text{s$^{-1}$} \\ 8\quad & \ce{O3}+\ce{NO} \to \ce{NO2}+\ce{O2} & 1.9 \times 10^{-14} \conc \end{align}

We will define maximum ozone concentrations from a 10-hour simulation specified from intial values of $\ce{RH}$, $\ce{NO}$, and $\ce{NO2}$ (constant photolysis rate).


  1. State which reactions are responsible for:
    • a. Formation of $\ce{O3}$
    • b. Destruction of $\ce{O3}$
    • c. Production of the precursor in (a)
  2. Derive analytical expressions for $[\ce{OH}]$, $[\ce{HO2}]$, $[\ce{RO2}]$ as a function of $[\ce{RH}]$, $[\ce{NO}]$, and $[\ce{NO2}]$ concentrations based on pseudo steady-state assumptions (see below).
  3. Plot evolving concentrations of $\ce{RH}$, $\ce{NO}$, $\ce{NO2}$, and $\ce{O3}$ for initial values of $[\ce{RH}]_0$ = 500 ppb and $[\ce{NO_x}]_0$ = 50 ppb for your selected time step.
  4. Evaluate maximum ozone concentrations for $[\ce{RH}]_0$ = {50, 60, ..., 100, 200, ..., 500} ppb and $[\ce{NO_x}]_0$ = {1, 2, ..., 10, 20, ..., 50} ppb.
    • a. Plot an ozone isopleth for these concentrations.
    • b. Plot the maximum ozone concentration (ppb) as a function of $\ce{NO_x}$ (ppb) for the different initial RH concentrations (ppb) as separate lines in a single plot.
    • c. Plot the percent change in maximum ozone concentrations at each initial concentration of $\ce{RH}$ when initial $\ce{NO_x}$ concentrations are reduced from 50 ppb to 30 ppb (a 40% reduction in $\ce{NO_x}$).
    • d. What is the minimum and maximum percent change in maximum ozone concentrations over the range simulated in Problem 4c?
  5. Select the cases where RH = 50 ppb (minimum). Plot the values below and discuss whether these trends are consistent with O3 concentrations in low and high $\ce{NO_x}$ regimes.
    • a. the ratio $k_5[\ce{HO2}][\ce{HO2}]$/$(k_4[\ce{OH}][\ce{NO2}])$ at $t=0$ as a function of $[\ce{NO_x}]_0$
    • b. the ratio $k_3[\ce{HO2}][\ce{NO}]$/$(k_4[\ce{OH}][\ce{NO2}])$ at $t=0$ as a function of $[\ce{NO_x}]_0$


  • Constant production rate of [$\ce{OH}$] is given by $P_{\ce{HO_x}}$ = 0.1 ppt s$^{-1}$.
  • Initial $[\ce{NO}]_0/[\ce{NO2}]_0 = 2$.
  • $T = 298$ K and $P = 1$ atm.


  • Initial concentration of $\ce{O3}$ can be calculated from photostationary state for initial values of $\ce{NO}$ and $\ce{NO2}$. You will need to use rate constants $k_7$ and $k_8$.
  • Since we define $[\ce{HO_x}] = [\ce{OH}] + [\ce{HO2}] + [\ce{RO2}]$ for this problem, $$ \frac{d[\ce{HO_x}]}{dt} = \frac{d[\ce{OH}]}{dt} + \frac{d[\ce{HO2}]}{dt} + \frac{d[\ce{RO2}]}{dt} $$ Note that $d[\ce{OH}]/dt$ includes $P_{\ce{HO_x}}$ in addition to the reactions listed above.
  • Pseudo steady-state relations for the following:
    • a. $d[\ce{HO_x}]/dt$ derived from balance of reactions 4-6 and $P_{\ce{HO_x}}$. (Combining all reactions 1-8 for each of $d[\ce{OH}]/dt$, $d[\ce{HO2}]/dt$, and $d[\ce{RO2}]/dt$ will result in an expression that leaves reactions 4-6 after cancellations of rates in production and degradation.) $$ 0 \approx P_{\ce{HO_x}} - r_4 - 2 r_5 - 2 r_6 $$
    • b. $d[\ce{RO2}]/dt$ derived from balance of reactions 1 and 2 (reaction 6 is an order of magnitude slower). I.e., production and destruction rates of $\ce{RO2}$ from reactions 1 and 2 are equal to each other. $$ 0 \approx r_1 - r_2 $$
    • c. $d[\ce{HO2}]/dt$ derived from balance of reactions 2 and 3 (reactions 5 and 6 are orders of magnitude slower). I.e., production and destruction rates of $\ce{HO2}$ from reactions 2 and 3 are equal to each other. $$ 0 \approx r_2 - r_3 $$
  • Further reaction of products $\ce{R^'CHO}$, $\ce{H2O2}$, and $\ce{ROOH}$ can be neglected.

Solution strategy (hints)

  • Write helper functions to convert between mixing ratio (ppb) and molecular concentrations.
    • Atmospheric concentrations are given in and should be reported in mixing ratio
    • Compute rates in molecular concentrations (the units of the kinetic rate constants) - this is a typical convention
  • Use PSSA relations to express radical concentrations as functions of $\ce{RH}$, $\ce{NO}$, and $\ce{NO2}$. There will be a quadratic solution involved - analyze which root makes sense based on the expectation that physical quantities are non-negative.
  • Write function to pass to ODE solver. Function should return vector of rate changes of concentrations: $$ \frac{dx_i}{dt} = \sum_{j\ \in\ \text{reactions}}\!\!\!\!\! \nu_{ji} r_j, \ \text{where}\ r_j = k_j\!\!\!\!\! \prod_{i\ \in\ \text{reactants}}\!\!\!\!\! x_i^{\nu_{ji}} $$ $x$ denotes concentration, $k$ denotes reaction rate constant, $\nu$ denotes stoichiometric number (negative for reactants and positive for products), $i$ denotes species, and $j$ denotes reaction number.
    • For the $\ce{HO2}$ self-reaction (reaction 5), two molecules of $\ce{HO2}$ are consumed in the reaction ($\nu_{5,\ce{HO2}}=2$). So the contribution to $d[\ce{HO2}]/dt$ is $-2k_5[\ce{HO_2}]^2$.
    • In contrast, reaction 3 involves only one molecule of $\ce{HO2}$ ($\nu_{3,\ce{HO2}}=1$). So its contribution to $d[\ce{HO2}]/dt$ is $-k_3[\ce{HO_2}][\ce{NO}]$.
  • Solve the ODEs for a single set of initial conditions (Problem 3). Do the numbers make sense? Evaluate if order of magnitude is sensible based on tropospheric concentrations or regulatory limits of ozone.
  • Write a loop to evaluate max ozone concentrations for a large number of initial conditions (Problem 4).

Regarding PSSA

We can write the rate equation above with initial conditions in vector form:

$$ \frac{d \mathbf{x}}{dt} = \boldsymbol{\nu} \mathbf{r} = \boldsymbol{f}\left(\mathbf{x}, \mathbf{k}, \boldsymbol{\nu}\right), \quad \mathbf{x}(t=0) = \mathbf{x}_0 \\ $$

Essentially, the rates of change in concentrations are a functions $\boldsymbol{f}$ of concentrations, rate coefficients $\mathbf{k}$, and stoichiometric numbers $\boldsymbol{\nu}$ relevant for each reaction.

The PSSA (pseudo steady-state assumption) (also referred to as the quasi steady-state assumption, or QSSA) assumes the existence of a set of slow and fast species $\mathbf{x} = \{\mathbf{x}^{(1)}, \mathbf{x}^{(2)}\}$; the dynamics for the fast species drive them toward their equilibrium distribution while the concentration of the slow species vary little over the same timescale. Mathematically, PSSA leads to the assumption that the dynamics of these species can be decoupled and solved by different methods to obtain the same chemical trajectories obtained by solving the full set of equations:

\begin{align} \frac{d \mathbf{x}^{(1)}}{dt} &= \boldsymbol{f}^{(1)}\left(\mathbf{x}, \mathbf{k}, \boldsymbol{\nu}\right), \quad \mathbf{x}^{(1)}(t=0) = \mathbf{x}^{(1)}_0 \\ \frac{d \mathbf{x}^{(2)}}{dt} &= \boldsymbol{f}^{(2)}\left(\mathbf{x}, \mathbf{k}, \boldsymbol{\nu}\right) \approx \boldsymbol{0} \\ \end{align}

Clarifications about PSSA:

  • PSSA does not mean that the concentrations of the fast species $\mathbf{x}^{(2)}$ is constant over time.
  • Therefore, $d\mathbf{x}^{(2)}(t)/dt \neq \boldsymbol{0}$ for all $t$, and $\mathbf{x}^{(2)}(t) \neq \mathbf{x}^{(2)}_0$.

Turányi and Tomlin (2014) can be consulted for more information (with historical perspective), if desired.

For this chemical system, $\mathbf{x}^{(1)} = \{\ce{[RH]}, \ce{[NO]}, \ce{[NO2]}, \ce{[O3]}, \ldots\}$, $\mathbf{x}^{(2)} = \{\ce{[OH]}, \ce{[HO2]}, \ce{[RO2]}\}$. The problem statement determines which specific reactions to consider for the steady-state expressions for each of the species in $\mathbf{x}^{(2)}$. $\{\ldots\}$ for $\mathbf{x}^{(1)}$ includes other products which are formed in the reactions that can be tracked but can also be neglected as they do not participate in further reactions leading to ozone production (for the purpose of this exercise). To obtain solutions:

  • $\mathbf{x}^{(1)}(t)$: solve differential equations
  • $\mathbf{x}^{(2)}(t)$: solve analytical expression given by $\boldsymbol{f}^{(2)}\left(\{\mathbf{x}^{(1)}(t), \mathbf{x}^{(2)}(t)\}, \mathbf{k}, \boldsymbol{\nu}\right) = \boldsymbol{0}$

Note that the function of rate equations you will want to integrate with the numerical solver should accept a vector of concentrations, $\mathbf{x}^{(1)}$ and return a vector of rates $d\mathbf{x}^{(1)}/dt$. Concentrations for $\mathbf{x}^{(2)}$ should also be calculated (through the PSSA) within the same function as they are necessary for calculating $d\mathbf{x}^{(1)}/dt$ at every time step.

The numerical solver recommended for GNU Octave is lsode and for MATLAB is ode45 or ode15s (for further information, see Beers, 2007).

Coding hints (implementation)

General strategy

Good practice:

  1. Define constants only once in your code.
  2. Write small functions that perform specific tasks. Assemble your program from the small functions rather than writing a single, monolithic script.
  3. Strive for readability.

For instance, functions you will want to have:

  1. Unit conversions. ppb $\to$ molec/cm$^{3}$, and molec/cm$^{3}$ $\to$ ppb.
  2. A function which returns the PSSA solution for $\mathbf{x}^{(2)}$: $\ce{OH}$, $\ce{HO2}$, and $\ce{RO2}$.
  3. A function for $d\mathbf{x}^{(1)}/dt$, which calls function defined in (2) to pass to the solver. (See example below.)


It is often the case that you have a vector of concentrations. To access individual elements of the vector, it may be convenient to pass around a vector of indices so the species referred to is clear. For example, you can define a structure like the following in your main script to hold values of vector indices:

ind = struct('RH', 1, 'NO', 2, 'NO2', 3, 'O3', 4);

If x is vector of concentrations at a given time step, you can access the $\ce{NO_2}$ concentrations using the syntax, x(ind.NO2).


Rate of change in concentrations. Here is an example of the function to be integrated - function #3 described above. You fill in the {...}:

function xdot = rateeqns(x, t, k, P_HOx, ind)

  %% x is concentration
  %% t is time
  %% k are rate constants
  %% P_HOx is HOx production rate
  %% ind is the structure of indices

  %% PSSA

  [OH, HO2, RO2] = pssasoln(x, k, P_HOx, ind);

  %% rates

  r = zeros(size(k));
  r(1) = k(1)*x(ind.RH)*OH;
  r(2) = {...}
  r(3) = {...}
  r(4) = {...}
  r(5) = {...}
  r(6) = {...}
  r(7) = {...}
  r(8) = k(8)*x(ind.O3)*x(ind.NO);

  %% output

  xdot = zeros(size(x)); % {RH, NO, NO2, O3}
  xdot(ind.RH) = -r(1);
  xdot(ind.NO) = {...}
  xdot(ind.NO2) = {...}
  xdot(ind.O3) = r(7) - r(8);


Plotting contours. For Problem 4a, you should obtain a matrix which contains maximum ozone concentrations for each initial value of $[\ce{RH}]_0$ and $[\ce{NO_x}]_0$. To plot a contour plot, you can find some help on these pages: 1 2 rather than in the help documentation. You will likely have a vector of initial RH and a vector of initial NOx concentrations. If the vectors are called init_RH and init_NOx with ozone matrix max_O3, you will call meshgrid prior to contour as shown in the following example:

[grid_RH, grid_NOx] = meshgrid(init_RH, init_NOx);
contour(grid_RH, grid_NOx, max_O3)
scale(gca, 'xscale', 'log', 'yscale', 'log');

Check your solutions

Problem 3:


Problem 4a (generated with the contour function):



Beers, K. (2007): Numerical Methods for Chemical Engineering: Applications in MATLAB. Cambridge University Press, New York, NY.

Seinfeld, J. H., and S. N. Pandis (2006): Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. 2nd edition. John Wiley & Sons, New York, NY.

Turányi, T. and A. S. Tomlin (2014): Analysis of Kinetic Reaction Mechanisms. Springer, New York, NY.