Using EarthSciML

As we saw in Using ModelingToolkit, keeping track of variables and equations when working with complex models can be a challenge. EarthSciML is designed to streamline this process in the case of geoscientific model by providing utilities for coupling model components together and running them, and also by providing reference implementations of standard model components that can be used as building blocks to create larger-scale geoscientific models.

Model Components

EarthSciML model components are ModelingToolkit ODESystems which contain metadata which specifies how each component should be coupled to other components.

So, for example, we can create an instance of the SuperFast gas-phase chemistry model model from the GasChem package like this:

using GasChem

chem = SuperFast()

\[ \begin{align} \frac{\mathrm{d} \mathtt{O3}\left( t \right)}{\mathrm{d}t} &= \mathtt{jNO2} \mathtt{NO2}\left( t \right) - \mathtt{jO32OH} \mathtt{O3}\left( t \right) - \mathtt{ISOP}\left( t \right) \mathtt{O3}\left( t \right) \mathtt{arrhenius21.k}\left( t \right) - \mathtt{NO}\left( t \right) \mathtt{O3}\left( t \right) \mathtt{arrhenius6.k}\left( t \right) - \mathtt{O3}\left( t \right) \mathtt{arrhenius2.k}\left( t \right) \mathtt{HO2}\left( t \right) - \mathtt{O3}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{arrhenius1.k}\left( t \right) \\ \frac{\mathrm{d} \mathtt{OH}\left( t \right)}{\mathrm{d}t} &= \mathtt{jCH3OOH} \mathtt{CH3OOH}\left( t \right) + 2 \mathtt{jH2O2} \mathtt{H2O2}\left( t \right) + 2 \mathtt{jO32OH} \mathtt{O3}\left( t \right) - \mathtt{CH4} \mathtt{arrhenius9.k}\left( t \right) \mathtt{OH}\left( t \right) - \mathtt{ISOP}\left( t \right) \mathtt{arrhenius18.k}\left( t \right) \mathtt{OH}\left( t \right) - 0.5 \mathtt{ISOP}\left( t \right) \mathtt{arrhenius20.k}\left( t \right) \mathtt{OH}\left( t \right) - \mathtt{ISOP}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{arrhenius19.k}\left( t \right) + \mathtt{arrhenius7.k}\left( t \right) \mathtt{NO}\left( t \right) \mathtt{HO2}\left( t \right) - \mathtt{rate\_OHCO10.k}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{CO}\left( t \right) - \mathtt{H2O2}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{arrhenius5.k}\left( t \right) - \mathtt{NO2}\left( t \right) \mathtt{arr\_3rdbody8.k}\left( t \right) \mathtt{OH}\left( t \right) - \mathtt{arrhenius11.k}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{CH2O}\left( t \right) - \mathtt{CH3OOH}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{arrhenius13.k}\left( t \right) + \mathtt{O3}\left( t \right) \mathtt{arrhenius2.k}\left( t \right) \mathtt{HO2}\left( t \right) - \mathtt{O3}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{arrhenius1.k}\left( t \right) - \mathtt{arrhenius3.k}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{HO2}\left( t \right) \\ \frac{\mathrm{d} \mathtt{HO2}\left( t \right)}{\mathrm{d}t} &= 2 \mathtt{jCH2Oa} \mathtt{CH2O}\left( t \right) + \mathtt{jCH3OOH} \mathtt{CH3OOH}\left( t \right) + 0.06 \mathtt{ISOP}\left( t \right) \mathtt{O3}\left( t \right) \mathtt{arrhenius21.k}\left( t \right) - \mathtt{arrhenius7.k}\left( t \right) \mathtt{NO}\left( t \right) \mathtt{HO2}\left( t \right) + \mathtt{rate\_OHCO10.k}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{CO}\left( t \right) + 0.8 \left( \mathtt{CH3O2}\left( t \right) \right)^{2} \mathtt{arrhenius16.k}\left( t \right) + \mathtt{CH3O2}\left( t \right) \mathtt{NO}\left( t \right) \mathtt{arrhenius15.k}\left( t \right) - \mathtt{CH3O2}\left( t \right) \mathtt{HO2}\left( t \right) \mathtt{arrhenius12.k}\left( t \right) + \mathtt{H2O2}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{arrhenius5.k}\left( t \right) + \mathtt{arrhenius11.k}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{CH2O}\left( t \right) - \mathtt{O3}\left( t \right) \mathtt{arrhenius2.k}\left( t \right) \mathtt{HO2}\left( t \right) + \mathtt{O3}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{arrhenius1.k}\left( t \right) - \mathtt{arrhenius3.k}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{HO2}\left( t \right) - 2 \left( \mathtt{HO2}\left( t \right) \right)^{2} \mathtt{rate\_HO2HO24.k}\left( t \right) \\ \frac{\mathrm{d} \mathtt{NO}\left( t \right)}{\mathrm{d}t} &= \mathtt{jNO2} \mathtt{NO2}\left( t \right) - \mathtt{arrhenius7.k}\left( t \right) \mathtt{NO}\left( t \right) \mathtt{HO2}\left( t \right) - \mathtt{CH3O2}\left( t \right) \mathtt{NO}\left( t \right) \mathtt{arrhenius15.k}\left( t \right) - \mathtt{NO}\left( t \right) \mathtt{O3}\left( t \right) \mathtt{arrhenius6.k}\left( t \right) \\ \frac{\mathrm{d} \mathtt{NO2}\left( t \right)}{\mathrm{d}t} &= - \mathtt{jNO2} \mathtt{NO2}\left( t \right) - \mathtt{H2O} \mathtt{NO2}\left( t \right) \mathtt{rateconvert17.k}\left( t \right) + \mathtt{arrhenius7.k}\left( t \right) \mathtt{NO}\left( t \right) \mathtt{HO2}\left( t \right) + \mathtt{CH3O2}\left( t \right) \mathtt{NO}\left( t \right) \mathtt{arrhenius15.k}\left( t \right) + \mathtt{NO}\left( t \right) \mathtt{O3}\left( t \right) \mathtt{arrhenius6.k}\left( t \right) - \mathtt{NO2}\left( t \right) \mathtt{arr\_3rdbody8.k}\left( t \right) \mathtt{OH}\left( t \right) \\ \frac{\mathrm{d} \mathtt{CH3O2}\left( t \right)}{\mathrm{d}t} &= \mathtt{CH4} \mathtt{arrhenius9.k}\left( t \right) \mathtt{OH}\left( t \right) + 2 \mathtt{ISOP}\left( t \right) \mathtt{arrhenius18.k}\left( t \right) \mathtt{OH}\left( t \right) + 1.86 \mathtt{ISOP}\left( t \right) \mathtt{O3}\left( t \right) \mathtt{arrhenius21.k}\left( t \right) - 2 \left( \mathtt{CH3O2}\left( t \right) \right)^{2} \mathtt{arrhenius16.k}\left( t \right) - \mathtt{CH3O2}\left( t \right) \mathtt{NO}\left( t \right) \mathtt{arrhenius15.k}\left( t \right) - \mathtt{CH3O2}\left( t \right) \mathtt{HO2}\left( t \right) \mathtt{arrhenius12.k}\left( t \right) + \mathtt{CH3OOH}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{arrhenius13.k}\left( t \right) \\ \frac{\mathrm{d} \mathtt{CH2O}\left( t \right)}{\mathrm{d}t} &= - \mathtt{jCH2Oa} \mathtt{CH2O}\left( t \right) - \mathtt{jCH2Ob} \mathtt{CH2O}\left( t \right) + \mathtt{jCH3OOH} \mathtt{CH3OOH}\left( t \right) + 0.87 \mathtt{ISOP}\left( t \right) \mathtt{O3}\left( t \right) \mathtt{arrhenius21.k}\left( t \right) + 2 \left( \mathtt{CH3O2}\left( t \right) \right)^{2} \mathtt{arrhenius16.k}\left( t \right) + \mathtt{CH3O2}\left( t \right) \mathtt{NO}\left( t \right) \mathtt{arrhenius15.k}\left( t \right) - \mathtt{arrhenius11.k}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{CH2O}\left( t \right) + \mathtt{arrhenius14.k}\left( t \right) \mathtt{CH3OOH}\left( t \right) \mathtt{OH}\left( t \right) \\ \frac{\mathrm{d} \mathtt{CO}\left( t \right)}{\mathrm{d}t} &= \mathtt{jCH2Oa} \mathtt{CH2O}\left( t \right) + \mathtt{jCH2Ob} \mathtt{CH2O}\left( t \right) + 0.05 \mathtt{ISOP}\left( t \right) \mathtt{O3}\left( t \right) \mathtt{arrhenius21.k}\left( t \right) - \mathtt{rate\_OHCO10.k}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{CO}\left( t \right) + \mathtt{arrhenius11.k}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{CH2O}\left( t \right) \\ \frac{\mathrm{d} \mathtt{CH3OOH}\left( t \right)}{\mathrm{d}t} &= - \mathtt{jCH3OOH} \mathtt{CH3OOH}\left( t \right) + \mathtt{CH3O2}\left( t \right) \mathtt{HO2}\left( t \right) \mathtt{arrhenius12.k}\left( t \right) - \mathtt{arrhenius14.k}\left( t \right) \mathtt{CH3OOH}\left( t \right) \mathtt{OH}\left( t \right) - \mathtt{CH3OOH}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{arrhenius13.k}\left( t \right) \\ \frac{\mathrm{d} \mathtt{ISOP}\left( t \right)}{\mathrm{d}t} &= - \mathtt{ISOP}\left( t \right) \mathtt{arrhenius18.k}\left( t \right) \mathtt{OH}\left( t \right) - \mathtt{ISOP}\left( t \right) \mathtt{O3}\left( t \right) \mathtt{arrhenius21.k}\left( t \right) \\ \frac{\mathrm{d} \mathtt{H2O2}\left( t \right)}{\mathrm{d}t} &= - \mathtt{jH2O2} \mathtt{H2O2}\left( t \right) - \mathtt{H2O2}\left( t \right) \mathtt{OH}\left( t \right) \mathtt{arrhenius5.k}\left( t \right) + \left( \mathtt{HO2}\left( t \right) \right)^{2} \mathtt{rate\_HO2HO24.k}\left( t \right) \\ \frac{\mathrm{d} \mathtt{HNO3}\left( t \right)}{\mathrm{d}t} &= 0.5 \mathtt{H2O} \mathtt{NO2}\left( t \right) \mathtt{rateconvert17.k}\left( t \right) + \mathtt{NO2}\left( t \right) \mathtt{arr\_3rdbody8.k}\left( t \right) \mathtt{OH}\left( t \right) \\ \mathtt{arrhenius1.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius1.K\_300}}{T} \right)^{\mathtt{arrhenius1.b0}} P \mathtt{arrhenius1.A} \mathtt{arrhenius1.a0} \mathtt{arrhenius1.ppb\_unit} e^{\frac{\mathtt{arrhenius1.c0}}{T}}}{T \mathtt{arrhenius1.R}} \\ \mathtt{arrhenius2.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius2.K\_300}}{T} \right)^{\mathtt{arrhenius2.b0}} P \mathtt{arrhenius2.A} \mathtt{arrhenius2.a0} \mathtt{arrhenius2.ppb\_unit} e^{\frac{\mathtt{arrhenius2.c0}}{T}}}{T \mathtt{arrhenius2.R}} \\ \mathtt{arrhenius3.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius3.K\_300}}{T} \right)^{\mathtt{arrhenius3.b0}} P \mathtt{arrhenius3.A} \mathtt{arrhenius3.a0} \mathtt{arrhenius3.ppb\_unit} e^{\frac{\mathtt{arrhenius3.c0}}{T}}}{T \mathtt{arrhenius3.R}} \\ \mathtt{rate\_HO2HO24.k}\left( t \right) &= \left( 1 + \frac{1.4 \cdot 10^{-30} \mathtt{H2O}^{2} P \mathtt{rate\_HO2HO24.A} \mathtt{rate\_HO2HO24.num\_density\_unit\_inv} \mathtt{rate\_HO2HO24.ppb\_inv} e^{\frac{\mathtt{rate\_HO2HO24.T\_0}}{T}}}{T \mathtt{rate\_HO2HO24.R}} \right) \left( \frac{P \mathtt{rate\_HO2HO24.A} \mathtt{rate\_HO2HO24.num\_density\_unit\_inv} \mathtt{rate\_HO2HO24.k1.k}\left( t \right)}{T \mathtt{rate\_HO2HO24.R}} + \mathtt{rate\_HO2HO24.k0.k}\left( t \right) \right) \\ \mathtt{rate\_HO2HO24.k0.k}\left( t \right) &= \frac{\left( \frac{\mathtt{rate\_HO2HO24.k0.K\_300}}{T} \right)^{\mathtt{rate\_HO2HO24.k0.b0}} P \mathtt{rate\_HO2HO24.k0.A} \mathtt{rate\_HO2HO24.k0.a0} \mathtt{rate\_HO2HO24.k0.ppb\_unit} e^{\frac{\mathtt{rate\_HO2HO24.k0.c0}}{T}}}{T \mathtt{rate\_HO2HO24.k0.R}} \\ \mathtt{rate\_HO2HO24.k1.k}\left( t \right) &= \frac{\left( \frac{\mathtt{rate\_HO2HO24.k1.K\_300}}{T} \right)^{\mathtt{rate\_HO2HO24.k1.b0}} P \mathtt{rate\_HO2HO24.k1.A} \mathtt{rate\_HO2HO24.k1.a0} \mathtt{rate\_HO2HO24.k1.ppb\_unit} e^{\frac{\mathtt{rate\_HO2HO24.k1.c0}}{T}}}{T \mathtt{rate\_HO2HO24.k1.R}} \\ \mathtt{arrhenius5.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius5.K\_300}}{T} \right)^{\mathtt{arrhenius5.b0}} P \mathtt{arrhenius5.A} \mathtt{arrhenius5.a0} \mathtt{arrhenius5.ppb\_unit} e^{\frac{\mathtt{arrhenius5.c0}}{T}}}{T \mathtt{arrhenius5.R}} \\ \mathtt{arrhenius6.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius6.K\_300}}{T} \right)^{\mathtt{arrhenius6.b0}} P \mathtt{arrhenius6.A} \mathtt{arrhenius6.a0} \mathtt{arrhenius6.ppb\_unit} e^{\frac{\mathtt{arrhenius6.c0}}{T}}}{T \mathtt{arrhenius6.R}} \\ \mathtt{arrhenius7.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius7.K\_300}}{T} \right)^{\mathtt{arrhenius7.b0}} P \mathtt{arrhenius7.A} \mathtt{arrhenius7.a0} \mathtt{arrhenius7.ppb\_unit} e^{\frac{\mathtt{arrhenius7.c0}}{T}}}{T \mathtt{arrhenius7.R}} \\ \mathtt{arr\_3rdbody8.k}\left( t \right) &= \frac{0.6^{\frac{1}{1 + \left( \log_{10}\left( \frac{P \mathtt{arr\_3rdbody8.A} \mathtt{arr\_3rdbody8.num\_density\_unit\_inv} \mathtt{arr\_3rdbody8.alow.k}\left( t \right)}{T \mathtt{arr\_3rdbody8.R} \mathtt{arr\_3rdbody8.ahigh.k}\left( t \right)} \right) \right)^{2}}} P \mathtt{arr\_3rdbody8.A} \mathtt{arr\_3rdbody8.num\_density\_unit\_inv} \mathtt{arr\_3rdbody8.alow.k}\left( t \right)}{T \mathtt{arr\_3rdbody8.R} \left( 1 + \frac{P \mathtt{arr\_3rdbody8.A} \mathtt{arr\_3rdbody8.num\_density\_unit\_inv} \mathtt{arr\_3rdbody8.alow.k}\left( t \right)}{T \mathtt{arr\_3rdbody8.R} \mathtt{arr\_3rdbody8.ahigh.k}\left( t \right)} \right)} \\ \mathtt{arr\_3rdbody8.alow.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arr\_3rdbody8.alow.K\_300}}{T} \right)^{\mathtt{arr\_3rdbody8.alow.b0}} P \mathtt{arr\_3rdbody8.alow.A} \mathtt{arr\_3rdbody8.alow.a0} \mathtt{arr\_3rdbody8.alow.ppb\_unit} e^{\frac{\mathtt{arr\_3rdbody8.alow.c0}}{T}}}{T \mathtt{arr\_3rdbody8.alow.R}} \\ \mathtt{arr\_3rdbody8.ahigh.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arr\_3rdbody8.ahigh.K\_300}}{T} \right)^{\mathtt{arr\_3rdbody8.ahigh.b0}} P \mathtt{arr\_3rdbody8.ahigh.A} \mathtt{arr\_3rdbody8.ahigh.a0} \mathtt{arr\_3rdbody8.ahigh.ppb\_unit} e^{\frac{\mathtt{arr\_3rdbody8.ahigh.c0}}{T}}}{T \mathtt{arr\_3rdbody8.ahigh.R}} \\ \mathtt{arrhenius9.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius9.K\_300}}{T} \right)^{\mathtt{arrhenius9.b0}} P \mathtt{arrhenius9.A} \mathtt{arrhenius9.a0} \mathtt{arrhenius9.ppb\_unit} e^{\frac{\mathtt{arrhenius9.c0}}{T}}}{T \mathtt{arrhenius9.R}} \\ \mathtt{rate\_OHCO10.k}\left( t \right) &= \frac{0.6^{\frac{1}{1 + \left( \log_{10}\left( \frac{P \mathtt{rate\_OHCO10.A} \mathtt{rate\_OHCO10.num\_density\_unit\_inv} \mathtt{rate\_OHCO10.klo2.k}\left( t \right)}{T \mathtt{rate\_OHCO10.R} \mathtt{rate\_OHCO10.khi2.k}\left( t \right)} \right) \right)^{2}}} \mathtt{rate\_OHCO10.klo2.k}\left( t \right)}{1 + \frac{P \mathtt{rate\_OHCO10.A} \mathtt{rate\_OHCO10.num\_density\_unit\_inv} \mathtt{rate\_OHCO10.klo2.k}\left( t \right)}{T \mathtt{rate\_OHCO10.R} \mathtt{rate\_OHCO10.khi2.k}\left( t \right)}} + \frac{0.6^{\frac{1}{1 + \left( \log_{10}\left( \frac{P \mathtt{rate\_OHCO10.A} \mathtt{rate\_OHCO10.num\_density\_unit\_inv} \mathtt{rate\_OHCO10.klo1.k}\left( t \right)}{T \mathtt{rate\_OHCO10.R} \mathtt{rate\_OHCO10.khi1.k}\left( t \right)} \right) \right)^{2}}} P \mathtt{rate\_OHCO10.A} \mathtt{rate\_OHCO10.num\_density\_unit\_inv} \mathtt{rate\_OHCO10.klo1.k}\left( t \right)}{T \mathtt{rate\_OHCO10.R} \left( 1 + \frac{P \mathtt{rate\_OHCO10.A} \mathtt{rate\_OHCO10.num\_density\_unit\_inv} \mathtt{rate\_OHCO10.klo1.k}\left( t \right)}{T \mathtt{rate\_OHCO10.R} \mathtt{rate\_OHCO10.khi1.k}\left( t \right)} \right)} \\ \mathtt{rate\_OHCO10.klo1.k}\left( t \right) &= \frac{\left( \frac{\mathtt{rate\_OHCO10.klo1.K\_300}}{T} \right)^{\mathtt{rate\_OHCO10.klo1.b0}} P \mathtt{rate\_OHCO10.klo1.A} \mathtt{rate\_OHCO10.klo1.a0} \mathtt{rate\_OHCO10.klo1.ppb\_unit} e^{\frac{\mathtt{rate\_OHCO10.klo1.c0}}{T}}}{T \mathtt{rate\_OHCO10.klo1.R}} \\ \mathtt{rate\_OHCO10.khi1.k}\left( t \right) &= \frac{\left( \frac{\mathtt{rate\_OHCO10.khi1.K\_300}}{T} \right)^{\mathtt{rate\_OHCO10.khi1.b0}} P \mathtt{rate\_OHCO10.khi1.A} \mathtt{rate\_OHCO10.khi1.a0} \mathtt{rate\_OHCO10.khi1.ppb\_unit} e^{\frac{\mathtt{rate\_OHCO10.khi1.c0}}{T}}}{T \mathtt{rate\_OHCO10.khi1.R}} \\ \mathtt{rate\_OHCO10.klo2.k}\left( t \right) &= \frac{\left( \frac{\mathtt{rate\_OHCO10.klo2.K\_300}}{T} \right)^{\mathtt{rate\_OHCO10.klo2.b0}} P \mathtt{rate\_OHCO10.klo2.A} \mathtt{rate\_OHCO10.klo2.a0} \mathtt{rate\_OHCO10.klo2.ppb\_unit} e^{\frac{\mathtt{rate\_OHCO10.klo2.c0}}{T}}}{T \mathtt{rate\_OHCO10.klo2.R}} \\ \mathtt{rate\_OHCO10.khi2.k}\left( t \right) &= \frac{\left( \frac{\mathtt{rate\_OHCO10.khi2.K\_300}}{T} \right)^{\mathtt{rate\_OHCO10.khi2.b0}} P \mathtt{rate\_OHCO10.khi2.A} \mathtt{rate\_OHCO10.khi2.a0} \mathtt{rate\_OHCO10.khi2.ppb\_unit} e^{\frac{\mathtt{rate\_OHCO10.khi2.c0}}{T}}}{T \mathtt{rate\_OHCO10.khi2.R}} \\ \mathtt{arrhenius11.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius11.K\_300}}{T} \right)^{\mathtt{arrhenius11.b0}} P \mathtt{arrhenius11.A} \mathtt{arrhenius11.a0} \mathtt{arrhenius11.ppb\_unit} e^{\frac{\mathtt{arrhenius11.c0}}{T}}}{T \mathtt{arrhenius11.R}} \\ \mathtt{arrhenius12.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius12.K\_300}}{T} \right)^{\mathtt{arrhenius12.b0}} P \mathtt{arrhenius12.A} \mathtt{arrhenius12.a0} \mathtt{arrhenius12.ppb\_unit} e^{\frac{\mathtt{arrhenius12.c0}}{T}}}{T \mathtt{arrhenius12.R}} \\ \mathtt{arrhenius13.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius13.K\_300}}{T} \right)^{\mathtt{arrhenius13.b0}} P \mathtt{arrhenius13.A} \mathtt{arrhenius13.a0} \mathtt{arrhenius13.ppb\_unit} e^{\frac{\mathtt{arrhenius13.c0}}{T}}}{T \mathtt{arrhenius13.R}} \\ \mathtt{arrhenius14.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius14.K\_300}}{T} \right)^{\mathtt{arrhenius14.b0}} P \mathtt{arrhenius14.A} \mathtt{arrhenius14.a0} \mathtt{arrhenius14.ppb\_unit} e^{\frac{\mathtt{arrhenius14.c0}}{T}}}{T \mathtt{arrhenius14.R}} \\ \mathtt{arrhenius15.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius15.K\_300}}{T} \right)^{\mathtt{arrhenius15.b0}} P \mathtt{arrhenius15.A} \mathtt{arrhenius15.a0} \mathtt{arrhenius15.ppb\_unit} e^{\frac{\mathtt{arrhenius15.c0}}{T}}}{T \mathtt{arrhenius15.R}} \\ \mathtt{arrhenius16.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius16.K\_300}}{T} \right)^{\mathtt{arrhenius16.b0}} P \mathtt{arrhenius16.A} \mathtt{arrhenius16.a0} \mathtt{arrhenius16.ppb\_unit} e^{\frac{\mathtt{arrhenius16.c0}}{T}}}{T \mathtt{arrhenius16.R}} \\ \mathtt{rateconvert17.k}\left( t \right) &= \frac{P \mathtt{rateconvert17.A} \mathtt{rateconvert17.a0} \mathtt{rateconvert17.ppb\_unit}}{T \mathtt{rateconvert17.R}} \\ \mathtt{arrhenius18.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius18.K\_300}}{T} \right)^{\mathtt{arrhenius18.b0}} P \mathtt{arrhenius18.A} \mathtt{arrhenius18.a0} \mathtt{arrhenius18.ppb\_unit} e^{\frac{\mathtt{arrhenius18.c0}}{T}}}{T \mathtt{arrhenius18.R}} \\ \mathtt{arrhenius19.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius19.K\_300}}{T} \right)^{\mathtt{arrhenius19.b0}} P \mathtt{arrhenius19.A} \mathtt{arrhenius19.a0} \mathtt{arrhenius19.ppb\_unit} e^{\frac{\mathtt{arrhenius19.c0}}{T}}}{T \mathtt{arrhenius19.R}} \\ \mathtt{arrhenius20.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius20.K\_300}}{T} \right)^{\mathtt{arrhenius20.b0}} P \mathtt{arrhenius20.A} \mathtt{arrhenius20.a0} \mathtt{arrhenius20.ppb\_unit} e^{\frac{\mathtt{arrhenius20.c0}}{T}}}{T \mathtt{arrhenius20.R}} \\ \mathtt{arrhenius21.k}\left( t \right) &= \frac{\left( \frac{\mathtt{arrhenius21.K\_300}}{T} \right)^{\mathtt{arrhenius21.b0}} P \mathtt{arrhenius21.A} \mathtt{arrhenius21.a0} \mathtt{arrhenius21.ppb\_unit} e^{\frac{\mathtt{arrhenius21.c0}}{T}}}{T \mathtt{arrhenius21.R}} \end{align} \]

It's just an ODESystem, so we can run it like any other ModelingToolkit model:

using ModelingToolkit, OrdinaryDiffEq, Plots

sol = solve(ODEProblem(structural_simplify(chem)), tspan=(0,2*3600))

plot(sol, legend=:topright, xlabel="Time (s)", ylabel="Concentration (ppb)",
    title="SuperFast Chemistry")

Coupling Model Components

EarthSciML model components differ from ordinary ModelingToolkit models in that they contain metadata that specifies how they should be coupled to other components. (In particular, this is done by specifying :coupletype metadata and the EarthSciMLBase.couple2 function, which is covered in more detail here.) This allows us to use the EarthSciMLBase.couple function to couple model components together into a single model that can be run as a single unit, for example as shown below:

using EarthSciMLBase

chem_phot = couple(
    SuperFast(),
    FastJX()
)

CoupledSystem containing 2 system(s), 0 operator(s), and 0 callback(s).

The above example couples the SuperFast with the Fast-JX photolysis model, and when we run them together as a coupled model, we can see a characteristic diurnal cycle in the O₃ concentration, were concentrations are lower at night and higher during the day owing to photochemistry:

chem_phot_sys = convert(ODESystem, chem_phot)
sol = solve(ODEProblem(chem_phot_sys), tspan=(0,3*24*3600))

plot(sol.t, sol[chem_phot_sys.SuperFast₊O3], legend=:topright, xlabel="Time (s)",
    ylabel="O₃ Concentration (ppb)", label=:none,
    title="SuperFast Chemistry + Fast-JX Photolysis")

External Data

We can also bring in external data to provide forcing to our model using EarthSciData. For example, below we load GEOS-FP meteorological data and pollutant emissions from the U.S. National Emissions Inventory.

using EarthSciData
using Dates

domain = DomainInfo(
    DateTime(2016, 5, 1),
    DateTime(2016, 5, 2);
    lonrange = deg2rad(-115):deg2rad(2.5):deg2rad(-68.75),
    latrange = deg2rad(25):deg2rad(2):deg2rad(53.7),
    levrange = 1:15,
    dtype = Float64)

emis = NEI2016MonthlyEmis("mrggrid_withbeis_withrwc", domain; stream=false)
emis = EarthSciMLBase.copy_with_change(emis, discrete_events=[]) # Workaround for bug.

equations(emis)[1:5]

\[ \begin{align} \mathtt{ACET}\left( t \right) &= ifelse\left( \mathtt{lev} < 2, \frac{\mathtt{ACET\_itp}\left( t, \mathtt{lon}, \mathtt{lat} \right)}{\mathtt{{\Delta}z}}, \mathtt{zero\_emis} \right) \\ \mathtt{ACROLEIN}\left( t \right) &= ifelse\left( \mathtt{lev} < 2, \frac{\mathtt{ACROLEIN\_itp}\left( t, \mathtt{lon}, \mathtt{lat} \right)}{\mathtt{{\Delta}z}}, \mathtt{zero\_emis} \right) \\ \mathtt{ALD2}\left( t \right) &= ifelse\left( \mathtt{lev} < 2, \frac{\mathtt{ALD2\_itp}\left( t, \mathtt{lon}, \mathtt{lat} \right)}{\mathtt{{\Delta}z}}, \mathtt{zero\_emis} \right) \\ \mathtt{ALD2\_PRIMARY}\left( t \right) &= ifelse\left( \mathtt{lev} < 2, \frac{\mathtt{ALD2\_PRIMARY\_itp}\left( t, \mathtt{lon}, \mathtt{lat} \right)}{\mathtt{{\Delta}z}}, \mathtt{zero\_emis} \right) \\ \mathtt{ALDX}\left( t \right) &= ifelse\left( \mathtt{lev} < 2, \frac{\mathtt{ALDX\_itp}\left( t, \mathtt{lon}, \mathtt{lat} \right)}{\mathtt{{\Delta}z}}, \mathtt{zero\_emis} \right) \end{align} \]

geosfp = GEOSFP("0.5x0.625_NA", domain; stream=false)
geosfp = EarthSciMLBase.copy_with_change(geosfp, discrete_events=[]) # Workaround for bug.

equations(geosfp)[1:5]

\[ \begin{align} \mathtt{A3dyn.DTRAIN}\left( t \right) &= \mathtt{A3dyn.DTRAIN\_itp}\left( t, \mathtt{lon}, \mathtt{lat}, \mathtt{lev} \right) \\ \mathtt{A3dyn.OMEGA}\left( t \right) &= \mathtt{A3dyn.OMEGA\_itp}\left( t, \mathtt{lon}, \mathtt{lat}, \mathtt{lev} \right) \\ \mathtt{A3dyn.RH}\left( t \right) &= \mathtt{A3dyn.RH\_itp}\left( t, \mathtt{lon}, \mathtt{lat}, \mathtt{lev} \right) \\ \mathtt{A3dyn.U}\left( t \right) &= \mathtt{A3dyn.U\_itp}\left( t, \mathtt{lon}, \mathtt{lat}, \mathtt{lev} \right) \\ \mathtt{A3dyn.V}\left( t \right) &= \mathtt{A3dyn.V\_itp}\left( t, \mathtt{lon}, \mathtt{lat}, \mathtt{lev} \right) \end{align} \]

In order to load data products, we will usually need to specify a spatiotemporal domain as demonstrated above. Simulations that don't include a spatial component will use data from the centerpoint of the provided spatial domain.

We can add our data providers to our model from above. We'll also add in other processes important to a simplified model of air quality, such as deposition. The coupling process will automatically make the appropriate connections between the model components to build a unified model.

using AtmosphericDeposition

model = couple(
    SuperFast(),
    FastJX(),
    #DrydepositionG(), Not currently working
    Wetdeposition(),
    emis,
    geosfp,
)

CoupledSystem containing 5 system(s), 0 operator(s), and 0 callback(s).

Box Model Simulation

Once we specify our model, we can try running it in a zero-dimensional "box model" format. This can allow us study the dynamics or diagnose any programs without using a lot of computing power.

model_sys = convert(ODESystem, model)
sol = solve(ODEProblem(model_sys), Rosenbrock23(), tspan=get_tspan(domain))

plot(unix2datetime.(sol.t), sol[model_sys.SuperFast₊O3],
    ylabel="O₃ Concentration (ppb)",
    xlabel="Time", label=:none,
    title="Air Quality Box Model")

3D Model Simulation

Once we're satisfied with our box model, we can move on to a full 3D simulation. To do this, we'll add in any spatial operators that we'll need, for example advection from the EnvironmentalTransport package. We can also add in a model component that will write the model results to the disk, which can be useful for large simulations where the results may not all fit in the memory of your computer. Finally, we also need to add our spatiotemporal domain to the model.

using EnvironmentalTransport

dt = 300.0 # Operator splitting timestep
advection = AdvectionOperator(dt, upwind1_stencil, ZeroGradBC())

outfile = ("RUNNER_TEMP" ∈ keys(ENV) ? ENV["RUNNER_TEMP"] : tempname()) * "out.nc" # This is just a location to save the output.
output = NetCDFOutputter(outfile, 3600.0) # Save output every 3600 sec. of simulation time,

model_3d = couple(model, advection, output, domain)

CoupledSystem containing 5 system(s), 1 operator(s), and 1 callback(s).

Graph Visualization

Once we've created our model, we can also check a graph representation of its components. This shows each of the components in the model as nodes and couplings between the components as edges.

using MetaGraphsNext
using CairoMakie, GraphMakie

g = graph(model_3d)
f = Figure(backgroundcolor = :transparent)
ax = Axis(f[1, 1], backgroundcolor = :transparent)
p = graphplot!(ax, g; ilabels=labels(g), edge_color=:gray)
CairoMakie.xlims!(ax, -2.4, 1.5)
CairoMakie.ylims!(ax, -2.3, 2.5)
hidedecorations!(ax); hidespines!(ax)
f

So for example in the figure above we can see that WetDeposition is coupled to SuperFast (because it is depositing the SuperFast chemical species) and to GEOS-FP (because it gets information such as Temperature and cloud fraction from GEOS-FP).

Running 3D simulation

Running a 3D simulation is similar to running the box model simulation, except we need to specify a SolverStrategy which will tell the solver how to combine the non-spatial and spatial components of the model. Here we will use SolverStrangSerial. We also add in a callback to ensure that the concentrations of the chemical species remain positive.

We then use our solver strategy as an argument to ODEProblem. We also specify some options in solve show a progress bar (because this simulation can take a while) and to prevent the output data from being saved in memory because we have already specified that we want to save it to disk instead.

using DiffEqCallbacks

st = SolverStrangSerial(Rosenbrock23(), dt, callback=PositiveDomain(save=false))

prob = ODEProblem(model_3d, st)

using ProgressLogging # Needed for progress bar. Use `TerminalLoggers` if in a terminal.

sol = solve(prob, SSPRK22(); dt=dt, progress=true, progress_steps=1,
    save_on=false, save_start=false, save_end=false, initialize_save=false)

retcode: Success
Interpolation: 3rd order Hermite
t: Float64[]
u: Vector{Float64}[]

Visualizing 3D model results

Because we have saved our results to disk instead of in the sol variable, we need to visualize differently than we have above. To do so, we can use the NCDatasets package.

using NCDatasets

ds = NCDataset(outfile, "r");

cross_section_lat = size(ds["SuperFast₊O3"], 2) ÷ 2
anim = @animate for i ∈ 1:size(ds["SuperFast₊O3"])[4]
    Plots.plot(
        Plots.heatmap(ds["SuperFast₊O3"][:, :, 1, i]', title="Ground-Level", clim=(0,200)),
        Plots.heatmap(ds["SuperFast₊O3"][:, cross_section_lat, :, i]', clim=(0,200),
            title="Vertical Cross-Section"),
        size=(1200, 400)
    )
end
gif(anim, fps = 5)

Next Steps

Now that you have a basic understanding of how to use EarthSciML, you can take a look at the model components and features available in the different Available Software Libraries that we provide.