API Reference

BlackbodyRadiation(; name=:BlackbodyRadiation)

Implements Planck's blackbody radiation law (Eq. 4.2) from Seinfeld & Pandis.

F_B(λ) = 2πc²hλ⁻⁵ / (exp(ch/kλT) - 1)

This equation gives the monochromatic emissive power of a blackbody at temperature T for a given wavelength λ.

Outputs (Variables)

• F_B_λ: Monochromatic emissive power (W m⁻³)

Inputs (Parameters)

• T: Temperature (K)
• λ: Wavelength (m)

Constants

• h: Planck's constant = 6.626 × 10⁻³⁴ J·s
• c: Speed of light = 2.9979 × 10⁸ m/s
• k: Boltzmann constant = 1.381 × 10⁻²³ J/K

BrOxCycle(; name=:BrOxCycle)

Create a ModelingToolkit System for the BrOx catalytic ozone destruction cycle.

Bromine is approximately 50 times more effective than chlorine in destroying ozone on an atom-for-atom basis (Page 169).

Key Reactions (Section 5.5.2, Page 166):

Br + O3 → BrO + O2 (k = 1.7 × 10⁻¹¹ exp(-800/T) cm³/molec/s) BrO + O → Br + O2 BrO + ClO → Br + Cl + O2 (or BrCl + O2) BrO + HO2 → HOBr + O2

Note on Reservoir Species (Page 169):

Unlike chlorine, bromine does not form a stable HBr reservoir because Br + CH4 is endothermic and extremely slow.

COOxidation(; name)

ModelingToolkit System for CO oxidation chemistry diagnostics.

Implements diagnostic calculations for the CO oxidation mechanism (Eqs. 6.9-6.17) from Seinfeld & Pandis Chapter 6, including HOx cycling and steady-state relationships.

Input Variables

• CO, OH, HO2, NO, NO2, O3: Species concentrations [m⁻³]

Output Variables

• P_O3: Net ozone production rate [m⁻³ s⁻¹]
• L_HOx: Total HOx loss rate [m⁻³ s⁻¹]
• chain_length: HOx chain length (cycles before termination)
• HO2_ss: HO₂ steady-state concentration (high NOx limit) [m⁻³]

Rate Constants at 298 K

• kCOOH = 2.4 × 10⁻¹³ cm³ molecule⁻¹ s⁻¹ = 2.4 × 10⁻¹⁹ m³ s⁻¹
• kHO2NO = 8.1 × 10⁻¹² cm³ molecule⁻¹ s⁻¹ = 8.1 × 10⁻¹⁸ m³ s⁻¹
• kHO2HO2 = 2.9 × 10⁻¹² cm³ molecule⁻¹ s⁻¹ = 2.9 × 10⁻¹⁸ m³ s⁻¹
• kOHNO2 = 1.0 × 10⁻¹¹ cm³ molecule⁻¹ s⁻¹ = 1.0 × 10⁻¹⁷ m³ s⁻¹
• kHO2O3 = 2.0 × 10⁻¹⁵ cm³ molecule⁻¹ s⁻¹ = 2.0 × 10⁻²¹ m³ s⁻¹
• kOHO3 = 7.3 × 10⁻¹⁴ cm³ molecule⁻¹ s⁻¹ = 7.3 × 10⁻²⁰ m³ s⁻¹

ChapmanMechanism(; name=:ChapmanMechanism)

Create a ModelingToolkit System for the Chapman mechanism.

The Chapman mechanism describes the basic production and destruction of ozone in the stratosphere through photolysis of O2 and subsequent reactions.

Reactions (Section 5.2, Seinfeld & Pandis 2006):

1. O2 + hν → O + O (j_O2)
2. O + O2 + M → O3 + M (k2)
3. O3 + hν → O + O2 (j_O3)
4. O + O3 → O2 + O2 (k4)

Rate Equations (Equations 5.1-5.2):

d[O]/dt = 2jO2[O2] - k2[O][O2][M] + jO3[O3] - k4[O][O3] d[O3]/dt = k2[O][O2][M] - j_O3[O3] - k4[O][O3]

Steady-State Ozone (Equation 5.13):

[O3]ss = 0.21 × (k2 × jO2 / (k4 × j_O3))^(1/2) × [M]^(3/2)

ClOxCycle(; name=:ClOxCycle)

Create a ModelingToolkit System for the ClOx catalytic ozone destruction cycle.

ClOx Cycle 1 (Section 5.5.1, Page 162):

Cl + O3 → ClO + O2 (k1 = 2.3 × 10⁻¹¹ exp(-200/T) cm³/molec/s) ClO + O → Cl + O2 (k2 = 3.0 × 10⁻¹¹ exp(70/T) cm³/molec/s) Net: O3 + O → O2 + O2

Rate of Odd Oxygen Destruction (Equation 5.29):

d[Ox]/dt = -2 k2[ClO][O]

Steady-State [Cl]/[ClO] Ratio (Equation 5.30):

[Cl]/[ClO] = (kClO+O[O] + kClO+NO[NO]) / (k_Cl+O3[O3])

Reservoir Species (Section 5.6):

• HCl formed by Cl + CH4 → HCl + CH3
• ClONO2 formed by ClO + NO2 + M → ClONO2 + M

ClimateFeedback(; name) -> ModelingToolkitBase.System

Climate feedback model relating radiative forcing to equilibrium temperature change with feedbacks.

This model implements the climate sensitivity equations with feedback from Seinfeld & Pandis (2006), Chapter 23 "Climate and Chemical Composition of the Atmosphere", Section 23.4.

Note: This extends the basic no-feedback climate sensitivity from Chapter 4 (see ClimateSensitivity) by incorporating climate feedback processes.

The model includes:

• Eq. 23.1: ΔTs = λ ΔF (climate sensitivity with feedbacks)
• Eq. 23.2: ΔT₀ = λ₀ ΔF (no-feedback temperature response)
• Eq. 23.3: Eᵢ = λᵢ / λ_CO₂ (efficacy of forcing agents)
• Eq. 23.4: ΔFₑ = ΔFᵢ Eᵢ (effective forcing)
• Eq. 23.7: ΔTunrealized = (ΔF - ΔFr) λ (unrealized warming, p. 1045)

Reference: Seinfeld, J.H. and Pandis, S.N. (2006) Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 2nd Edition, Chapter 23, John Wiley & Sons.

Example

using GasChem, ModelingToolkit
sys = ClimateFeedback()

ClimateSensitivity(; name=:ClimateSensitivity)

Implements the climate sensitivity equations (Eqs. 4.8-4.10) from Seinfeld & Pandis.

ΔF_net = ΔF_S - ΔF_L       (Eq. 4.8 - Net radiative energy change)
ΔT_e = λ₀ × ΔF_net         (Eq. 4.9 - Temperature response)
λ₀ = 1/(4σT_e³) = T_e/(4F_L) (Eq. 4.10 - Climate sensitivity factor)

Outputs (Variables)

• ΔF_net: Net radiative forcing (W/m²)
• ΔT_e: Change in equilibrium temperature (K)
• λ_0: Climate sensitivity factor (K m²/W)
• F_L: Reference emitted longwave flux (W/m²)

Inputs (Parameters)

• T_e: Reference equilibrium temperature (K)
• ΔF_S: Change in absorbed solar flux (W/m²)
• ΔF_L: Change in emitted longwave flux (W/m²)

Constants

• σ: Stefan-Boltzmann constant = 5.671 × 10⁻⁸ W m⁻² K⁻⁴

FastJX(
    t_ref::AbstractFloat;
    name
) -> ModelingToolkitBase.System

Description: This is a box model used to calculate the photolysis reaction rate constant using the Fast-JX scheme (Neu, J. L., Prather, M. J., and Penner, J. E. (2007), Global atmospheric chemistry: Integrating over fractional cloud cover, J. Geophys. Res., 112, D11306, doi:10.1029/2006JD008007.)

Argument:

• t_ref: Reference time for the model, can be a DateTime or a Unix timestamp (in seconds).

Example

Build Fast-JX model:

fj = FastJX(DateTime(2000, 1, 1))

Description: This is a box model used to calculate the photolysis reaction rate constant using the Fast-JX scheme (Neu, J. L., Prather, M. J., and Penner, J. E. (2007), Global atmospheric chemistry: Integrating over fractional cloud cover, J. Geophys. Res., 112, D11306, doi:10.1029/2006JD008007.)

Argument:

• t_ref: Reference time for the model, can be a DateTime or a Unix timestamp (in seconds).

Example

Build Fast-JX model:

fj = FastJX(DateTime(2000, 1, 1))

GEOS-Chem full-chem mechanism adapted from GEOS-Chem version 14.1.1

• Adapted from file https://github.com/geoschem/geos-chem/blob/4722f288e90291ba904222f4bbe4fc216d17c34a/KPP/fullchem/fullchem.eqn
• The GEOS-Chem license applies: https://github.com/geoschem/geos-chem/blob/main/LICENSE.txt

If the keyword argument rxn_sys is set to true, the function will return a reaction system instead of an ODE system.

===============================================================================

REFERENCES (alphabetical order)

• Atkinson2003: Atkinson and Arey, Chem. Rev., doi:10.1021/cr0206420, 2003.
• Bates2014: Bates et al., J. Phys. Chem A, 118, doi:10.1021/jp4107958, 2014.
• Bates2019: Bates and Jacob, Atmos. Chem. Phys., doi:10.5194/acp-19-9613-2019, 2019.
• Bates2021a: Bates et al, JGR, https://doi.org/10.1029/2020JD033439, 2021.
• Bates2021b: Bates et al, ACP, https://doi.org/10.5194/acp-2021-605, 2021.
• Browne2011: Browne et al., Atmos. Chem. Phys., doi:10.5194/acp-11-4209-2011, 2011.
• Browne2014: Browne et al., Atmos. Chem. Phys., doi:10.5194/acp-14-1225-2014, 2014.
• Chen2017: Chen et al., Geophys. Res. Lett., doi:10.1002/2017GL073812, 2017.
• Crounse2012: Crounse et al., J. Phys. Chem. A, doi:10.1021/jp211560u, 2012.
• Eastham2014: Eastham et al., Atmos. Env., doi:10.1016/j.atmosenv.2014.02.001, 2014.
• Fischer2014: Fischer et al., Atmos. Chem. Phys., doi:10.5194/acp-14-2679-2014, 2014.
• Fisher2016: Fisher et al., Atmos. Chem. Phys., doi:10.5194/acp-16-5969-2016, 2016.
• Fisher2018: Fisher et al., J. Geophys. Res., doi:10.1029/2018JD029046, 2018.
• Fry2014: Fry et al. Environ. Sci. Technol., doi:10.1021/es502204x, 2014.
• Gill2002: Gill and Hites, J. Phys. Chem. A, doi:10.1021/jp013532, 2002.
• Goliff2013: Goliff et al., Atmos. Environ., doi:10.1016/j.atmosenv.2012.11.038, 2013.
• Jacobs2014: Jacobs et al., Atmos. Chem. Phys., doi:10.5194/acp-14-8933-2014, 2014.
• Jenkin2015: Jenkin et al., Atmos. Chem. Phys., doi:10.5194/acp-15-11433-2015, 2015.
• Kasibhatla2018: Kasibhatla et al., Atmos. Chem. Phys., doi:10.5194/acp-18-11185-2018, 2018
• IUPAC ROO19: https://iupac-aeris.ipsl.fr/htdocs/datasheets/pdf/ROO19CH3O2NO3.pdf
• JPL 10-6: JPL Publication 10-6, https://jpldataeval.jpl.nasa.gov/previous_evaluations.html, 2011.
• JPL 15-10: JPL Publication 15-10, https://jpldataeval.jpl.nasa.gov, 2015.
• Kwon2020: Kwon et al, Elementa, https://doi.org/10.1525/elementa.2021.00109, 2020.
• Lee2014: Lee et al., J. Phys. Chem. A, doi:10.1021/jp4107603, 2014.
• Marais2016: Marais et al., Atmos. Chem. Phys, doi:10.5194/acp-16-1603-2016, 2016.
• Miller2017: Miller et al., Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-1042, 2017.
• Millet2015: Millet et al., Atmos. Chem. Phys., doi:10.5194/acp-15-6283-2015, 2015. Moch et al, JGR, https, * Moch2020 # //doi.org/10.1029/2020JD032706, 2020.
• Muller2014: Muller et al., Atmos. Chem. Phys., doi:10.5194/acp-14-2497-2014, 2014.
• Parrella2012: Parrella et al. Atmos. Chem. Phys, doi:10.5194/acp-12-6723-2012, 2012.
• Paulot2009: Paulot et al., Atmos. Chem. Phys., doi:10.5194/acp-9-1479-2009, 2009a and Paulot et al., Science, doi:10.1126/science.1172910, 2009b.
• Peeters2010: Peeters and Muller, Phys. Chem. Chem. Phys., doi:10.1039/C0CP00811G, 2010.
• Peeters2014: Peeters et al., J. Phys. Chem. A, doi:10.1021/jp5033146, 2014.
• Pye2010: Pye et al., Atmos. Chem. Phys., doi:10.5194/acp-10-11261-2010, 2010.
• Roberts1992: Roberts and Bertman, Int. J. Chem. Kinet., doi:10.1002/kin.550240307, 1992.
• Sherwen2016b: Sherwen et al., Atmos. Chem. Phys., doi:10.5194/acp-16-12239-2016, 2016b.
• Sherwen2017: Sherwen et al., Faraday Discuss., doi:10.1039/C7FD00026J, 2017.
• StClair2016: St. Clair et al., J. Phys. Chem. A, doi:10.1021/acs.jpca.5b065322016, 2016.
• Travis2016: Travis et al., Atmos. Chem. Phys., doi:10.5194/acp-16-13561-2016, 2016.
• Wolfe2012: Wolfe et al., Phys. Chem. Chem. Phys., doi: 10.1039/C2CP40388A, 2012.
• Xie2013: Xie et al., Atmos. Chem. Phys., doi:10.5194/acp-13-8439-2013, 2013.

GHGForcing(; name) -> ModelingToolkitBase.System

Greenhouse gas radiative forcing model.

Computes the radiative forcing from well-mixed greenhouse gases using values from IPCC (2001) as presented in Seinfeld & Pandis (2006), Chapter 23, Section 23.3.

Reference: Seinfeld, J.H. and Pandis, S.N. (2006) Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 2nd Edition, Chapter 23, John Wiley & Sons.

Example

using GasChem, ModelingToolkit
sys = GHGForcing()

GWP_exponential(τ_A, a_A, t_f; τ_CO2=150.0)

Calculate the Global Warming Potential using the closed-form solution for exponential decay.

This implements the formula from Problem 23.3 (p. 1050) of Seinfeld & Pandis (2006):

\[GWP = \frac{a_A \tau_A (1 - e^{-t_f/\tau_A})}{a_{CO_2} \tau_{CO_2} (1 - e^{-t_f/\tau_{CO_2}})}\]

Arguments

• τ_A: Atmospheric lifetime of species A [years]
• a_A: Radiative efficiency of species A relative to CO₂ (dimensionless)
• t_f: Time horizon for GWP calculation [years]
• τ_CO2: Effective CO₂ lifetime (default: 150 years)

Returns

• GWP: Global Warming Potential relative to CO₂ (dimensionless)

Example

# Calculate 100-year GWP for CH₄ (τ = 12 yr, fitted relative efficiency ≈ 140)
gwp_ch4 = GWP_exponential(12.0, 140.0, 100.0)

Reference: Seinfeld, J.H. and Pandis, S.N. (2006) Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 2nd Edition, Chapter 23, Problem 23.3.

GlobalWarmingPotential(; name) -> ModelingToolkitBase.System

Global Warming Potential (GWP) calculation for exponentially decaying species.

Implements Eq. 23.5 and 23.6 from Seinfeld & Pandis (2006), Chapter 23, for the special case of exponential decay of atmospheric perturbations. Also implements the closed-form GWP formula from Problem 23.3 (p. 1050).

The GWP is defined as the ratio of the time-integrated radiative forcing from a pulse emission of 1 kg of compound A to that of 1 kg of CO₂ over a specified time horizon.

Equations implemented:

• Eq. 23.5: GWP = ∫₀^{tf} aₐ[A(t)]dt / ∫₀^{tf} a_R[R(t)]dt
• Eq. 23.6: AGWP_A = ∫₀^{tf} aₐ[A(t)]dt
• Problem 23.3: GWP = [aₐ τₐ (1 - e^{-tf/τₐ})] / [aCO₂ τCO₂ (1 - e^{-tf/τ_CO₂})]

Reference: Seinfeld, J.H. and Pandis, S.N. (2006) Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 2nd Edition, Chapter 23, John Wiley & Sons.

Example

using GasChem, ModelingToolkit
sys = GlobalWarmingPotential()

HOxCycle(; name=:HOxCycle)

Create a ModelingToolkit System for the HOx catalytic ozone destruction cycle.

HOx Cycle 1 (Page 159):

OH + O3 → HO2 + O2 HO2 + O → OH + O2 Net: O3 + O → O2 + O2

HOx Cycle 2 (Page 159):

OH + O3 → HO2 + O2 HO2 + O3 → OH + O2 + O2 Net: O3 + O3 → O2 + O2 + O2

Steady-State Ratio (Equation 5.28):

[HO2]/[OH] = kOH+O3[O3] / (kHO2+NO[NO])

Rate of Odd Oxygen Destruction (Equation 5.27):

d[Ox]/dt = -2 kHO2+O3[HO2][O3] - 2 kHO2+O[HO2][O]

MethaneOxidation(; name)

ModelingToolkit System implementing reaction rate diagnostics for the methane oxidation mechanism from Table 6.1 of Seinfeld & Pandis Chapter 6.

This system computes individual reaction rates and diagnostic production/loss terms.

Species (Input Variables)

• CH4, CH3, CH3O2, CH3O, CH3OOH: Methane chain species [m⁻³]
• HCHO, HCO, CO: Formaldehyde and products [m⁻³]
• OH, HO2, H: HOx species [m⁻³]
• NO, NO2, O, O2, O3: NOx and oxygen species [m⁻³]
• M: Total air density [m⁻³]

Diagnostics (Output Variables)

• R1-R17: Individual reaction rates [m⁻³ s⁻¹]
• PO3gross: Gross O₃ production from peroxy+NO reactions [m⁻³ s⁻¹]
• P_HCHO: HCHO production rate [m⁻³ s⁻¹]
• L_CH4: CH₄ loss rate [m⁻³ s⁻¹]

Rate Constants

All rate constants from Table 6.1 at 298 K are implemented as parameters. Bimolecular rate constants converted from cm³/molec/s to m³/s (×10⁻⁶). Termolecular rate constants converted from cm⁶/molec²/s to m⁶/s (×10⁻¹²). # Parameters - Rate constants at 298 K from Table 6.1 (converted to SI)

MethaneOxidationODE(; name)

Full ODE system for methane oxidation with species time derivatives.

This uses Catalyst.jl's reaction network DSL to define the 17 reactions from Table 6.1 plus auxiliary reactions (CO+OH, OH+NO₂, HO₂+HO₂, NO+O₃) and an external OH source. The reaction network is converted to an ODE system.

Note: This is a stiff system due to the wide range of timescales (radicals have lifetimes of seconds, while CH₄ has a lifetime of years).

O₂ and M (total air density) are treated as parameters with default values for termolecular reactions (reactions 2, 6, 13, 14, 17).

NOxCycle(; name=:NOxCycle)

Create a ModelingToolkit System for the NOx catalytic ozone destruction cycle.

NOx Cycle 1 (Section 5.3.2, Page 154):

NO + O3 → NO2 + O2 (k1 = 3.0 × 10⁻¹² exp(-1500/T) cm³/molec/s) NO2 + O → NO + O2 (k2 = 5.6 × 10⁻¹² exp(180/T) cm³/molec/s) Net: O3 + O → O2 + O2

Steady-State (Equation 5.20):

0 = k1[NO][O3] - j_NO2[NO2] - k2[NO2][O]

Rate of Odd Oxygen Destruction (Equation 5.22):

d[Ox]/dt = -2 k2[NO2][O]

NOxPhotochemistry(; name)

ModelingToolkit System for the NOx photochemical cycle.

Implements the photostationary state relationships (Eqs. 6.5-6.8) from Seinfeld & Pandis Chapter 6.

Input Variables

• NO: Nitric oxide concentration [m⁻³]
• NO2: Nitrogen dioxide concentration [m⁻³]
• O3: Ozone concentration [m⁻³]
• O2: Molecular oxygen concentration [m⁻³]
• M: Total air number density [m⁻³]

Output Variables

• O: Ground-state oxygen atom concentration [m⁻³]
• O3_pss: Photostationary state ozone concentration [m⁻³]
• Φ: Photostationary state parameter [dimensionless]
• P_O3: Net ozone production rate [m⁻³ s⁻¹]

Parameters

• j_NO2: NO₂ photolysis rate [s⁻¹]
• k_O_O2_M: Rate constant for O + O₂ + M → O₃ [m⁶ s⁻¹]
• k_NO_O3: Rate constant for NO + O₃ → NO₂ + O₂ [m³ s⁻¹]

Rate Constants at 298 K (from Table B.1)

• j_NO2 ≈ 8 × 10⁻³ s⁻¹ (typical midday value)
• kOO2_M = 6.0 × 10⁻³⁴ cm⁶ molecule⁻² s⁻¹ = 6.0 × 10⁻⁴⁶ m⁶ s⁻¹ # Parameters (rate constants converted to SI)
• kNOO3 = 1.9 × 10⁻¹⁴ cm³ molecule⁻¹ s⁻¹ = 1.9 × 10⁻²⁰ m³ s⁻¹ (p. 211)

OHProduction(; name)

ModelingToolkit System for OH radical production from O₃ photolysis.

Implements the steady-state O(¹D) concentration (Eq. 6.1), OH production rate (Eq. 6.3), and OH yield (Eq. 6.4) from Seinfeld & Pandis Chapter 6.

This is an algebraic system that computes diagnostic quantities from input concentrations.

Input Variables (must be provided)

• O3: Ozone concentration [m⁻³]
• H2O: Water vapor concentration [m⁻³]
• M: Total air number density [m⁻³]

Output Variables (computed)

• O1D: Excited oxygen O(¹D) steady-state concentration [m⁻³]
• P_OH: OH production rate [m⁻³ s⁻¹]
• ε_OH: OH yield (fraction of O(¹D) producing OH) [dimensionless]

Parameters

• j_O3: O₃ photolysis rate producing O(¹D) [s⁻¹] (default: 6×10⁻⁵ s⁻¹ at surface, solar zenith 0°)
• k3_N2: Rate constant for O(¹D) + N₂ quenching [m³ s⁻¹]
• k3_O2: Rate constant for O(¹D) + O₂ quenching [m³ s⁻¹]
• k4: Rate constant for O(¹D) + H₂O → 2OH [m³ s⁻¹]
• f_N2: Fraction of M that is N₂ [dimensionless]
• f_O2: Fraction of M that is O₂ [dimensionless]

Rate Constants at 298 K (from Table B.1, Seinfeld & Pandis)

• j_O3 ≈ 6 × 10⁻⁵ s⁻¹ (at surface, solar zenith angle 0°, Figure 4.15)
• k3_N2 = 2.6 × 10⁻¹¹ cm³ molecule⁻¹ s⁻¹ = 2.6 × 10⁻¹⁷ m³ s⁻¹
• k3_O2 = 4.0 × 10⁻¹¹ cm³ molecule⁻¹ s⁻¹ = 4.0 × 10⁻¹⁷ m³ s⁻¹
• k4 = 2.2 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹ = 2.2 × 10⁻¹⁶ m³ s⁻¹

OzoneProductionEfficiency(; name)

ModelingToolkit System for calculating Ozone Production Efficiency (OPE).

OPE = Δ[O₃] / Δ[NOx] represents the number of ozone molecules produced per NOx molecule consumed. This is a key metric for understanding ozone-NOx-VOC chemistry in different regimes.

From Equations 6.21-6.24:

• At high NOx (VOC-limited): OPE is low
• At low NOx (NOx-limited): OPE is high
• The transition occurs near the "ridge line" in O₃-NOx-VOC space

Input Variables

• OH, HO2, RO2, NO, NO2: Species concentrations [m⁻³]

Output Variables

• P_O3: Gross ozone production rate [m⁻³ s⁻¹]
• L_NOx: NOx loss rate [m⁻³ s⁻¹]
• OPE: Ozone production efficiency [dimensionless]

PhotonEnergy(; name=:PhotonEnergy)

Implements the photon energy equation (Eq. 4.1) from Seinfeld & Pandis.

Δε = hν = hc/λ

The energy of a photon is related to its frequency (ν) or wavelength (λ) through Planck's constant (h) and the speed of light (c).

Outputs (Variables)

• Δε: Photon energy (J)
• ν: Frequency (Hz = s⁻¹)

Inputs (Parameters)

• λ: Wavelength (m) - input parameter

Constants

• h: Planck's constant = 6.626 × 10⁻³⁴ J·s
• c: Speed of light = 2.9979 × 10⁸ m/s

PhotostationaryState(; name)

Simplified model for analyzing photostationary state deviations.

This system calculates the deviation from photostationary state (Φ - 1), which indicates the net effect of peroxy radical chemistry on the NO-NO₂-O₃ cycle.

When Φ > 1: Additional oxidants (HO₂, RO₂) are converting NO to NO₂ When Φ < 1: Additional reductants are present When Φ = 1: Pure photostationary state (no net ozone production)

PlanetaryEnergyBalance(; name=:PlanetaryEnergyBalance)

Implements the planetary energy balance equations (Eqs. 4.5-4.7) from Seinfeld & Pandis.

F_S = S₀/4 × (1 - R_p)     (Eq. 4.5 - Absorbed solar flux)
F_L = σT_e⁴                 (Eq. 4.6 - Emitted longwave flux)
T_e = ((1 - R_p)S₀/4σ)^(1/4) (Eq. 4.7 - Equilibrium temperature)

At equilibrium: FS = FL

Outputs (Variables)

• F_S: Absorbed solar flux (W/m²)
• F_L: Emitted longwave flux (W/m²)
• T_e: Planetary equilibrium temperature (K)

Inputs (Parameters)

• S_0: Solar constant = 1370 W/m² (can be varied)
• R_p: Planetary albedo ~ 0.3 (can be varied)

Constants

• σ: Stefan-Boltzmann constant = 5.671 × 10⁻⁸ W m⁻² K⁻⁴

Pollu()

This atmospheric chemical system model is built based on this paper:

"GAUSS-SEIDEL ITERATION FOR STIFF ODES FROM CHEMICAL KINETICS" (1994), J. G. VERWER.

The input of the function is Temperature, concentrations of all chemicals, and reaction rates of photolysis reactions

If the keyword argument rxn_sys is set to true, the function will return a reaction system instead of an ODE system.

Example

using GasChem, EarthSciMLBase, DifferentialEquations, Plots
rs = Pollu()
sol = solve(ODEProblem(mtkcompile(rs), [], (0,360), [], combinatoric_ratelaws=false), AutoTsit5(Rosenbrock23()), saveat=10.0)
plot(sol)

RadiationFundamentals(; name=:RadiationFundamentals)

Combined system implementing all radiation fundamentals equations from Seinfeld & Pandis Chapter 4 (Eqs. 4.1-4.11).

This is a composed system that includes:

• PhotonEnergy: Energy-frequency-wavelength relations (Eq. 4.1)
• BlackbodyRadiation: Planck's law (Eq. 4.2)
• WienDisplacement: Peak emission wavelength (Eq. 4.3)
• StefanBoltzmann: Total emissive power (Eq. 4.4)
• PlanetaryEnergyBalance: Earth's energy balance (Eqs. 4.5-4.7)
• ClimateSensitivity: Temperature response to forcing (Eqs. 4.8-4.10)
• TOARadiativeForcing: Net flux at top of atmosphere (Eq. 4.11)

RemoteConditions(; name)

ModelingToolkit System representing typical remote/background conditions with low NOx.

All concentrations are parameters in SI units (m⁻³).

StefanBoltzmann(; name=:StefanBoltzmann)

Implements the Stefan-Boltzmann law (Eq. 4.4) from Seinfeld & Pandis.

F_B = σT⁴

The total emissive power of a blackbody integrated over all wavelengths.

Outputs (Variables)

• F_B: Total blackbody emissive power (W/m²)

Inputs (Parameters)

• T: Temperature (K)

Constants

• σ: Stefan-Boltzmann constant = 5.671 × 10⁻⁸ W m⁻² K⁻⁴

StratosphericOzoneSystem(; name=:StratosphericOzoneSystem)

Create a comprehensive ModelingToolkit System combining all stratospheric ozone chemistry cycles.

This system includes:

• Chapman mechanism (O, O3 production and loss)
• O(¹D) photochemistry including O(¹D) + H2O and O(¹D) + CH4
• N2O + O(¹D) → 2NO as the NOx source (Section 5.3.1)
• NOx cycle (catalytic O3 destruction)
• HOx cycle (catalytic O3 destruction)
• ClOx cycle (catalytic O3 destruction)
• BrOx cycle (catalytic O3 destruction)

All rate coefficients are temperature-dependent, computed from Arrhenius parameters defined as @constants.

Key Equations from Seinfeld & Pandis Chapter 5:

Odd Oxygen Balance (Equation 5.9):

d[Ox]/dt = 2j_O2[O2] - 2k4[O][O3]

Steady-State O3 (Equation 5.13):

[O3]ss = 0.21 × (k2 × jO2 / (k4 × j_O3))^(1/2) × [M]^(3/2)

[O]/[O3] Ratio (Equation 5.7):

[O]/[O3] = j_O3 / (k2[O2][M])

Time to Steady State (Equation 5.17):

τO3^ss = (1/4) × (k2[M] / (k4 × jO2 × j_O3))^(1/2)

NO Yield from N2O (Equation 5.18, Section 5.3.1):

N2O + O(¹D) → NO + NO (k = 6.7 × 10⁻¹¹ cm³/molec/s) N2O + O(¹D) → N2 + O2 (k = 4.9 × 10⁻¹¹ cm³/molec/s)

SuperFast()

This atmospheric chemical system model is built based on the Super Fast Chemical Mechanism, which is one of the simplest representations of atmospheric chemistry. It can efficiently simulate background tropheric ozone chemistry and perform well for those species included in the mechanism. The chemical equations we used is included in the supporting table S2 of the paper:

"Evaluating simplified chemical mechanisms within present-day simulations of the Community Earth System Model version 1.2 with CAM4 (CESM1.2 CAM-chem): MOZART-4 vs. Reduced Hydrocarbon vs. Super-Fast chemistry" (2018), Benjamin Brown-Steiner, Noelle E. Selin, Ronald G. Prinn, Simone Tilmes, Louisa Emmons, Jean-François Lamarque, and Philip Cameron-Smith.

The input of the function is Temperature, concentrations of all chemicals, and reaction rates of photolysis reactions

If the keyword argument rxn_sys is set to true, the function will return a reaction system instead of an ODE system.

Example

using GasChem, EarthSciMLBase, DifferentialEquations, Plots
rs = SuperFast()
sol = solve(ODEProblem(mtkcompile(rs), [], (0,360), [], combinatoric_ratelaws=false), AutoTsit5(Rosenbrock23()), saveat=10.0)
plot(sol)

TOARadiativeForcing(; name=:TOARadiativeForcing)

Implements the top of atmosphere radiative forcing based on Eq. 4.11 from Seinfeld & Pandis.

F_net = S₀/4 × (1 - R_p) - F_L

Note: Seinfeld & Pandis define Eq. 4.11 as -Fnet = S₀/4(1-Rp) - FL, using the convention that -Fnet represents net downward flux. Here we define Fnet as net incoming flux directly, so positive Fnet indicates the planet is gaining energy (warming).

Outputs (Variables)

• F_net: Net radiative flux at TOA (W/m²) - positive = energy gain
• F_S: Absorbed solar flux (W/m²)

Inputs (Parameters)

• S_0: Solar constant = 1370 W/m² (can be varied)
• R_p: Planetary albedo ~ 0.3 (can be varied)
• F_L: Emitted longwave flux (W/m²) - input parameter

TroposphericChemistrySystem(; name)

Combined ModelingToolkit System for tropospheric chemistry diagnostics.

This system couples the OH production, NOx cycling, and CO/CH₄ oxidation mechanisms to create a comprehensive diagnostic model of tropospheric photochemistry.

Key coupling points:

• OH is produced from O₃ photolysis and HO₂ + NO reaction
• OH is consumed by CO, CH₄, and other species
• HO₂ is produced from CO/VOC oxidation
• HO₂ is lost to NO (producing OH) and self-reaction
• NOx cycles between NO and NO₂ through O₃ and peroxy radicals
• O₃ is produced when peroxy radicals oxidize NO to NO₂

Subsystem Composition:

• oh: OHProduction subsystem (Section 6.1)
• nox: NOxPhotochemistry subsystem (Section 6.2)
• co: COOxidation subsystem (Section 6.3)

Input Variables

All species concentrations must be provided as inputs [m⁻³].

Diagnostic outputs:

• PO3net: Net ozone production rate [m⁻³ s⁻¹]
• OPE: Ozone production efficiency [dimensionless]
• HOx: Total HOx (OH + HO₂) [m⁻³] # =========================================================================
• chain_length: HOx chain length [dimensionless] # Subsystems (composed from individual component functions)

TypicalConditions(; name)

ModelingToolkit System representing typical lower troposphere conditions.

All concentrations are parameters in SI units (m⁻³).

Based on values from Seinfeld & Pandis Chapter 6.

UrbanConditions(; name)

ModelingToolkit System representing typical urban conditions with elevated NOx.

All concentrations are parameters in SI units (m⁻³).

WienDisplacement(; name=:WienDisplacement)

Implements Wien's displacement law (Eq. 4.3) from Seinfeld & Pandis.

λ_max = 2.897 × 10⁻³ / T  (SI: λ_max in m, T in K)

Note: S&P express this as λmax = 2.897 × 10⁶ / T when λmax is in nm.

This gives the wavelength at which the blackbody emission spectrum peaks.

Outputs (Variables)

• λ_max: Peak wavelength (m)

Inputs (Parameters)

• T: Temperature (K)

Constants

• b: Wien's displacement constant = 2.897 × 10⁻³ m·K

adjust_j_O31D(T, P, H2O)

Adjust the photolysis rate of O3 -> O2 + O(1D) to represent the effective rate for O3 -> 2OH. This adjustment is based on the fraction of O(1D) that reacts with H2O to produce 2 OH.

Third body effect for pressure dependence of rate coefficients. a1, b1, c1 are the Arrhenius parameters for the lower-limit rate. a2, b2, c2 are the Arrhenius parameters for the upper-limit rate. fv is the falloff curve parameter, (see ATKINSON ET. AL (1992) J. Phys. Chem. Ref. Data 21, P. 1145). Usually fv = 0.6.

Third body effect for pressure dependence of rate coefficients. a1, b1, c1 are the Arrhenius parameters for the lower-limit rate. a2, b2, c2 are the Arrhenius parameters for the upper-limit rate. fv is the falloff curve parameter, (see ATKINSON ET. AL (1992) J. Phys. Chem. Ref. Data 21, P. 1145). Usually fv = 0.6.

Third body effect for pressure dependence of first-order rate coefficients. a1, b1, c1 are the Arrhenius parameters for the lower-limit rate. a2, b2, c2 are the Arrhenius parameters for the upper-limit rate. fv is the falloff curve parameter, (see ATKINSON ET. AL (1992) J. Phys. Chem. Ref. Data 21, P. 1145). Usually fv = 0.6.

Arrhenius equation:

\[ k = a0 * exp( c0 / T ) * (T/300)^b0\]

Compute a first-order reaction rate constant using a modified Arrhenius equation, with output units in s⁻¹.

The formula used is:

\[k(t) = a₀ * exp(c₀ / T) * (300 / T)^b₀\]

Compute a second-order reaction rate constant using a modified Arrhenius equation, with output in SI units m³·molec⁻¹·s⁻¹.

The formula used is:

\[k(t) = a₀ * exp(c₀ / T) * (300 / T)^b₀\]

Arrhenius equation:

\[ k = a0 * exp( c0 / T ) * (T/300)^b0\]

Compute a third-order reaction rate constant using a modified Arrhenius equation, with output units in ppb⁻²·s⁻¹.

This function is appropriate when reactant concentrations are expressed in ppb (mol/mol_air × 10⁹).

The formula used is:

\[k(t) = a₀ * exp(c₀ / T) * (300 / T)^b₀\]

Modified Arrhenius law with output in units s⁻¹.

Modified Arrhenius law with output in units ppb·s⁻¹.

Function to create a constant rate coefficient for second-order reactions

Function to create a constant rate coefficient for first-order reactions

cos_solar_zenith_angle(lat, t, long)

This function is to compute the cosine of the solar zenith angle, given the unixtime, latitude and longitude The input variables: lat=latitude(°), long=longitude(°), t=unixtime(s) the cosine of the solar zenith angle (SZA) is given by: . cos(SZA) = sin(LAT)sin(DEC) + cos(LAT)cos(DEC)*cos(AHR)

       where LAT = the latitude angle,
             DEC = the solar declination angle,
             AHR = the hour angle, all in radians.  All in radians

Create a vector of interpolators to interpolate the cross sections σ (TODO: What are the units?) for different wavelengths (in nm) and temperatures (in K).

We use use linear interpolation with flat extrapolation.

Calculates the equilibrium constant for second-order reactions Find the backwards reaction by K=kforward/kbackwards Calculates the rate constant of the forward reaction

Used to compute the rate for these reactions: PPN = RCO3 + NO2 PAN = MCO3 + NO2

Calculates the equilibrium constant for first-order reactions Find the backwards reaction by K=kforward/kbackwards Calculates the rate constant of the forward reaction

Used to compute the rate for these reactions: ClOO {+M} = Cl + O2 {+M} Cl2O2 {+M} = 2ClO {+M}

get_conditions_dict(sys)

Helper to extract default parameter values from a conditions System as a Dict{Symbol,Float64}. This is a convenience function for use in tests and analysis code.

Get mean photolysis rates at different times

Used to compute the rate for this reactions: HO2 + HO2 = H2O2 + O2

Used to compute the rate for these reactions: IHOO1 + NO = NO2 + ... IHOO4 + NO = NO2 + ... IHP001 + NO = NO2 + ... IHP002 + NO = NO2 + ... IHP003 + NO = NO2 + ... IEPOXAOO + NO = NO2 + ... IEPOXBOO + NO = NO2 + ... ICHOO + NO = NO2 + ... ISOPNOO1 + NO = 1.728NO2 + ... ISOPNOO2 + NO = NO2 + ... IDHNDOO1 + NO = NO2 + ... IDHNDOO2 + NO = NO2 + ... IDHNBOO + NO = NO2 + ... IDHNDOO + NO = NO2 + ... INO2B + NO = 2.000NO2 + ... INO2D + NO = NO2 + ... IHPNBOO + NO = 1.065NO2 + ... IHPNDOO + NO = NO2 + ... MVKOHOO + NO = NO2 + ... MCROHOO + NO = NO2 + ...

Reaction rate for: DMS + OH = 0.750SO2 + 0.250MSA + MO2

Used to compute the rate for these reactions: RIPA + OH = 0.67IEPOXA + 0.33IEPOXB + OH + 0.005LVOC RIPB + OH = 0.68IEPOXA + 0.321IEPOB + OH + 0.005LVOC IEPOXA + OH = 0.67IEPOXA00 + 0.33IEPOXB00 IEPOXB + OH = 0.81IEPOXA00 + 0.19IEPOXB00 IHN2 + OH = 0.67IEPOXA + 0.33IEPOXB + NO2 IHN3 + OH = 0.67IEPOXA + 0.33IEPOXB + NO2 IHN1 + OH = IEPOXD + NO2 IHN4 + OH = IEPOXD + NO2 INPB + OH = OH + ITHN INPD + OH = OH + ITHN INPD + OH = NO2 + ICHE ICN + OH = NO2 + ICHE

Used to compute the rate for this reaction: GLYC + OH = 0.732CH2O + 0.361CO2 + 0.505CO + 0.227OH

• 0.773HO2 + 0.134GLYX + 0.134HCOOH which is the "A" branch of GLYC + OH.

For this reaction, these Arrhenius law terms evaluate to 1: (300/T)^b0 * exp(c0/T) Because b0 = c0 = 0.

Used to compute the rate for this reaction: GLYC + OH = HCOOH + OH + CO which is the "B" branch of GLYC + OH.

For this reaction, these Arrhenius law terms evaluate to 1: (300/T)^b0 * exp(c0/T) Because b0 = c0 = 0.

Reaction rate for: GLYX + NO3 = HNO3 + HO2 + 2CO i.e. the HO2 + 2*CO branch

Used to compute the rate for this reaction: HAC + OH = MGLY + HO2 which is the "A" branch of HAC + OH.

Used to compute the rate for this reaction: HAC + OH = 0.5HCOOH + OH + 0.5ACTA + 0.5CO2 + 0.5CO + 0.5MO2 which is the "B" branch of HAC + OH.

Used to compute the rate for this reactions: HO2 + HO2 = H2O2 + O2

Used to compute the rate for these reactions: ISOP + OH = LISOPOH + IHOO1 ISOP + OH = LISOPOH + IHOO4

Used to compute the rate for these reactions: ISOP + OH = 0.3MCO3 + 0.3MGLY + 0.3CH2O

• 0.15HPALD3 + 0.25HPALD1 + 0.4HO2
• 0.6CO + 1.5OH + 0.3HPETHNL + LISOPOH ISOP + OH = 0.3CH2O + 0.15HPALD4 + 0.25HPALD2
• 1.5OH + 0.9CO + 0.7HO2 + 0.3MGLY
• 0.3ATOOH + LISOPOH

Used to compute the rate for these reactions: IHOO1 + NO = IHN2 IHOO4 + NO = IHN4 IHPOO1 + NO = IHTN IHPOO2 + NO = IHTN IHPOO2 + NO = IHTN IEPOXAOO + NO = IHTN IEPOXBOO + NO = IHTN IHCOO + NO = IHTN ISOPNOO1 + NO = IDN ISOPNOO2 + NO = IDN IDHNDOO1 + NO = IDN IDHNDOO2 + NO = IDN INO2B + NO = IDN INO2D + NO = IDN IHPNBOO + NO = IDN IHPNDOO + NO = IDN MVK0HOO + NO = 0.438MVKN MCROHOO + NO = MCRHN

Reaction rate for: OH + CO = HO2 + CO2 (cf. JPL 15-10)

Used to compute the rate for these reactions: HNO3 + OH = H2O + NO3 HONIT + OH = NO3 + HAC

Reaction rate for: OH + CO = HO2 + CO2 (cf. JPL 15-10)

Used to compute the rate for reaction: BZPAN –> BZCO3 + NO2

Used to compute the rate for reaction: MACR1OO + NO2 –> MPAN MACRNO2 + NO2 –> MPAN + NO2 BZCO3 + NO2 –> BZPAN

Carbon Dependence of RO2+HO2, used in these reactions: A3O2 + HO2 = RA3P PO2 + HO2 = PP KO2 + HO2 = 0.150OH + 0.150ALD2 + 0.150MCO3 + 0.850ATOOH B3O2 + HO2 = RB3P PRN1 + HO2 = PRPN

Reaction rate for the "A" branch of these RO2 + NO reactions: MO2 + NO = MENO3 in which the "a1" parameter equals exactly 1.

For these reactions, these Arrhenius law terms evaluate to 1: (300/T)^b0 (300/T)^b1 * exp(c1/T) because b0 = b1 = c1 = 0.

Special treatment for methyl nitrate based on observations as Carter and Atkinson formulation does not apply to C1. Value based on upper limit of Flocke et al. 1998 as applied in Fisher et al. 2018

Reaction rate for the "A" branch of these RO2 + NO reactions, ETO2 + NO = ETNO3 A3O2 + NO = NPRNO3 R4O2 + NO = R4N2 B3O2 + NO = IPRNO3 in which the "a1" parameter is greater than 1.0.

Reaction rate for the "B" branch of these RO2 + NO reactions: MO2 + NO = CH2O + NO2 + HO2 in which the "a1" parameter equals exactly 1.

For these reactions, these Arrhenius law terms evaluate to 1: (300/T)^b0 (300/T)^b1 * exp(c1/T) because b0 = c0 = c1 = 0.

Reaction rate for the "B" branch of these RO2 + NO reactions: ETO2 + NO = NO2 + HO2 + ... A3O2 + NO = NO2 + HO2 + ... R4O2 + NO = NO2 + 0.27HO2 + ... B3O2 + NO = NO2 + HO2 + ... in which the "a1" parameter is greater than 1.0.

Use this function when a1 input argument is greater than 1.0.

Convert reaction rate value in unit of cm^3/molec/s to (s*ppb)^-1

Temperature Dependent Branching Ratio

Computes a temperature-dependent reaction rate constant using a modified Arrhenius expression with additional tunneling effect terms. Used to compute the rate for these reactions with output in units s⁻¹: IHOO1 = 1.5OH + ... IHOO4 = 1.5OH + ...

Computes a temperature-dependent reaction rate constant using a modified Arrhenius expression with additional tunneling effect terms. Used to compute the rate for these reactions with output in units ppb·s⁻¹: IHOO1 = 1.5OH + ... IHOO4 = 1.5OH + ...