HSSModel.jl Documentation

HSSModel implements a basic NOx-HOx steady state model that can be used to compute NOx lifetime, ozone production efficency, and other quantities for theoretical conditions.

Installation

Currently, this package is only available to install from GitHub. From the Julia REPL:

julia> ]add https://github.com/joshua-laughner/HSSModel.git

Since this will track the master branch, you may also wish to pin this package so that it does not update until you want it to:

julia> ]pin HSSModel

Basic Walkthrough

Importing HSSModel exposes the main driver function, hox_ss_solver from the SteadyStateSolvers module. This function requires 5 inputs:

NO concentration (in molec. cm$^{-3}$)
NO2 concentration (in molec. cm$^{-3}$)
HOx production rate (in molec. cm$^{-3}$ s$^{-1}$)
Total VOC reactivity (in s$^{-1}$)
The NO + RO2 branching ratio, $\alpha$ (unitless).

For this example, we'll be using values from Table 4 of Murphy et al. 2006 for P(HOx), VOC reactivity, and $\alpha$, plus we will assume 5 ppb of NOx. Since our inputs need to be in number density, rather than mixing ratio, we'll need the number density of air to convert.

# gas constant in Pa cm^3 K^-1 molec^-1
R = 8.314 * 1e6 / 6.022e23

# number density of air at STP
nair = 101325 / (R * 298)

# Convert from NOx in ppb to NO and NO2 in molec. cm^-3
# assuming a 4:1 NO2:NO ratio
nox_ppb = 5;
no_nd = nox_ppb*1e-9 * 0.2 * nair;
no2_nd = nox_ppb*1e-9 * 0.8 * nair;

# Other default values from Murphy et al. 2006
phox = 2.5e-9 * nair / 3600;  # convert ppb hr^-1 -> molec. cm^-3 s^-1
vocr = 5.8;
alpha = 0.04;

With these conditions, we call the steady state solver, which will compute the steady state concentrations of OH, HO2, and RO2 given these fixed values.

using HSSModel;
result = hox_ss_solver(no_nd, no2_nd, phox, vocr, alpha);

println("OH = $(result.oh), HO2 = $(result.ho2), RO2 = $(result.ro2)");

OH = 1.3288062945776744e7, HO2 = 3.555827622868582e8, RO2 = 3.8550423326615673e8

This result structure contains information about the steady state concentrations, as well as the model configuration. Generally, most information you would need to plot the model results is contained in this structure.

These result structure can also be use to calculate secondary quantities. The DerivedQuantities module exports two functions to do so: nox_lifetime and ozone_production_efficiency. Calling these with your results produces:

nox_lifetime(result)

Dict{String,Float64} with 3 entries:
  "total" => 1.96583
  "hno3"  => 2.37549
  "ans"   => 11.3992

This returns a dictionary with three different lifetimes, in hours:

total is the overall NOx lifetime.
hno3 is the lifetime with respect to loss to HNO3
ans is the lifetime with respect to loss to alkyl nitrates

ozone_prod_efficiency(result)

10.279835377647869

Unlike nox_lifetime, this returns a single value. This represents the ratio of ozone production to NOx loss.

Running an ensemble

The steady state model function works well with Julia's built in broadcasting behavior. For a simple example, let's look at running a range of NOx concentrations:

# Compute NO and NO2 number densities for logarithmically-spaced
# NOx mixing ratios between 0.01 and 10 ppb, still assuming a
# 4:1 NO2:NO ratio.
nox_nd = exp10.(range(-2., stop=1., length=20)) .* 1e-9 .* nair;
no2_nd = 0.8 .* nox_nd;
no_nd = 0.2 .* nox_nd;

# Note the broadcasting operator (the period following `hox_ss_solver`)
results = hox_ss_solver.(no_nd, no2_nd, phox, vocr, alpha);
size(results)

(20,)

If you want to test an ensemble with multiple variables, you can use broadcasting in multiple dimensions. In this example, we'll use the same NOx vector (a 1D, 75 element vector) along with a vector of VOC reactivity values. By making the VOCR values a row (1-by-$N$) vector, the broadcasting will result in a 2D array of results.

# Create a 1-by-10 vector of VOC reactivities
vocr_vec = collect(transpose(range(1., stop=10., length=10)));
results = hox_ss_solver.(no_nd, no2_nd, phox, vocr_vec, alpha);
size(results)

(20, 10)

A useful trick for quickly extracting values from these results arrays (for plotting or downstread analysis) is to broadcast the getfield function:

getfield.(results, :oh)

20×10 Array{Float64,2}:
 1.70933e7  8.55371e6  5.70404e6  …  2.13975e6  1.90204e6  1.71187e6
 1.71007e7  8.56047e6  5.70923e6     2.14202e6  1.90408e6  1.71372e6
 1.7123e7   8.5761e6   5.72065e6     2.14677e6  1.90833e6  1.71756e6
 1.71781e7  8.6101e6   5.74476e6     2.15649e6  1.91701e6  1.7254e6 
 1.72999e7  8.6804e6   5.79374e6     2.17585e6  1.93428e6  1.74098e6
 1.7547e7   8.81798e6  5.88861e6  …  2.21291e6  1.96731e6  1.77078e6
 1.80091e7  9.07057e6  6.06186e6     2.28018e6  2.02724e6  1.82481e6
 1.88068e7  9.50386e6  6.35855e6     2.39515e6  2.12966e6  1.91715e6
 2.00859e7  1.02016e7  6.83704e6     2.58095e6  2.29519e6  2.06641e6
 2.20083e7  1.12664e7  7.57094e6     2.86772e6  2.5508e6   2.29696e6
 2.47215e7  1.28168e7  8.65044e6  …  3.29486e6  2.93183e6  2.64086e6
 2.82704e7  1.49677e7  1.01764e7     3.91265e6  3.48376e6  3.13961e6
 3.23832e7  1.77659e7  1.22325e7     4.78088e6  4.26158e6  3.84403e6
 3.60763e7  2.10258e7  1.48023e7     5.95756e6  5.32119e6  4.80761e6
 3.73347e7  2.40196e7  1.75806e7     7.46232e6  6.69014e6  6.06262e6
 3.41351e7  2.52282e7  1.96835e7  …  9.18573e6  8.29248e6  7.55685e6
 2.73379e7  2.3181e7   1.96983e7     1.07277e7  9.80613e6  9.02775e6
 2.01573e7  1.86061e7  1.70493e7     1.13205e7  1.05569e7  9.88305e6
 1.43693e7  1.3774e7   1.31629e7     1.03748e7  9.90906e6  9.47509e6
 1.01194e7  9.86258e6  9.60545e6     8.37773e6  8.15128e6  7.93244e6