A mass- and energy-conserving framework for using machine learning to speed computations: a photochemistry example

Large air quality models and large climate models simulate the physical and chemical properties of the ocean, land surface, and/or atmosphere to predict atmospheric composition, energy balance and the future of our planet. All of these models employ some form of operator splitting, also called the method of fractional steps, in their structure, which enables each physical or chemical process to be simulated in a separate operator or module within the overall model. In this structure, each of the modules calculates property changes for a fixed period of time; that is, property values are passed into the module, which calculates how they change for a period of time and then returns the new property values, all in round-robin between the various modules of the model. Some of these modules require the vast majority of the computer resources consumed by the entire model, so increasing their computational efficiency can either improve the model's computational performance, enable more realistic physical or chemical representations in the module, or a combination of these two. Recent efforts have attempted to replace these modules with ones that use machine learning tools to memorize the input–output relationships of the most time-consuming modules. One shortcoming of some of the original modules and their machine-learned replacements is lack of adherence to conservation principles that are essential to model performance. In this work, we derive a mathematical framework for machine-learned replacements that conserves properties – say mass, atoms, or energy – to machine precision. This framework can be used to develop machine-learned operator replacements in environmental models.

1 Introduction

Complex systems require large models that simulate the wide range of physical and chemical properties that govern their performance. In the air quality realm, models include CMAQ (Foley et al., 2010), CAMx (Yarwood et al., 2007), WRF-Chem (Grell et al., 2005), and GEOS-Chem (Eastham et al., 2014). In the climate change arena, models include HadCM3 (Jones et al., 2005), GFDL CM2 (Delworth et al., 2012), ARPEGE-Climat (Somot et al., 2008), CESM (Kay et al., 2015), and E3SM (Golaz et al., 2019). These models employ operator splitting, also called the method of fractional steps (Janenko, 1971), in their structure so that each module can be tasked with representing one or a small number of physical and/or chemical processes. This modular structure enhances model maintenance and sustainability while enabling diverse physical and chemical processes to interact. Each module is tasked with simulating its processes over a fixed period of time, each module called in turn until they have all returned their results. Usually, the computational performance of these models is governed by one or two modules that consume the vast majority of the computer resources. In air quality models, this is usually the photochemistry and/or aerosol dynamics modules. In climate models, this is usually the radiative energy transport module.

Machine learning has been used to improve the computational efficiency of modules in atmospheric models for decades (Potukuchi and Wexler, 1997). As machine learning algorithms have improved, these efforts have matured (Hsieh, 2009; Kelp et al., 2018; Rasp et al., 2018; Keller and Evans, 2019; Pal et al., 2019). But the effort to replace physical and chemical operators with machine-learned modules is challenging because small systematic errors can build. For instance, a 0.1 % error over a 1 h time step could lead to a 72 % error after a month of simulation. This problem is compounded if the replacement module does not conserve quantities that are essential to model accuracy, such as atoms in a photochemical module, molecules and mass in an aerosol dynamics module, or energy in a radiative transfer module.

Recent efforts at developing and using machine-learned replacement modules has focused on memorizing how the quantities change. Some have also explored enforcing physical constraints when memorizing these quantities via postprediction balancing approaches (Krasnopolsky et al., 2010), introducing a penalty into the cost function (Beucler et al., 2019) or incorporating hard constraints on a subset of the output in neural network architecture (Beucler et al., 2019). All of these approaches focus on memorizing how quantities change and incorporate some correction strategy after all or a portion of the quantities have been predicted. If instead we focus on how the fluxes between quantities change, we can guarantee adherence to conservation principles to machine precision without a postprediction correction. In photochemical modules, the fluxes are how atoms move between chemical species as reactions progress. In aerosol dynamics, the fluxes are the condensation and evaporation or coagulation processes that move material between the gas and particle phases or between particle sizes. In radiative transfer, the fluxes are the energy movements between spatial domains.

In this work, we derive a mathematical framework that enables the use of machine learning tools to memorize these fluxes. We focus this work on atmospheric photochemistry and provide an example for a simple photochemical reaction mechanism because the number of species and the complexity of the problem exercises many aspects of the framework.

2 Derivation of the framework for photochemistry

In general, the atmospheric chemistry operator solves

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{C}}{\partial t}}=\mathit{F}\left(\mathit{C},T,\text{RH},\mathrm{actinic}\mathrm{flux},\mathrm{stability},\mathrm{etc}.\right),\end{array}$$

where C is a vector containing the current concentration of the chemical species, T is temperature, and RH is the relative humidity. A full list of symbols can be found in Appendix A. The right-hand side can be written as

$$\begin{array}{}\text{(2)}& \mathit{F}=\mathbf{A}\mathit{R},\end{array}$$

where A is a matrix describing the stoichiometry, and R is a vector of reactions. The form of the right-hand side assures mass balance because it is composed of reactions that destroy one species while creating one or more other ones, all in balance, described by A. The R terms take forms such as kC_iC_j, where k is the rate constant for a reaction between species; JC_i, where J is the photolysis reaction rate; or k(C_i(x₁)−C_i(x₂)), where k is a diffusion or mass transport rate constant between two spatial locations or between the gas and particle phases.

In the method of fractional steps, all modules integrate their equations forward for a fixed time step, Δt, that we call the operator-splitting time step. Combining these two equations and integrating gives

$$\begin{array}{}\text{(3a)}& \mathrm{\Delta }{C}_i=\sum _j{A}_{i,j}\underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}{R}_j\left(t\right)\text{d}t=\sum _j{A}_{i,j}{S}_j,\end{array}$$

or in matrix form

$$\begin{array}{}\text{(3b)}& \mathbf{\Delta }\mathit{C}=\mathbf{A}\mathit{S},\end{array}$$

where S_j equals ${\int }_t^{t+\mathrm{\Delta }t}{R}_j\left(t\right)\text{d}t$. S_j is the flux integral. For atmospheric photochemistry, it is the flux of atoms between molecules. For aerosol dynamics, it is the flux of molecules condensing on or evaporating from particles or the flux of small particles coagulating on large particles. For radiative transfer, S_j is the energy between spatial coordinates. We are able to pull the A out of the integral if it is a constant, which is usually the case or can be approximated as such.

Using machine learning tools to learn the relationship

$$\begin{array}{}\text{(4)}& \mathit{S}=\mathit{S}\left(\mathit{C},T,\text{RH},\mathrm{actinic}\mathrm{flux},\mathrm{stability},\mathrm{etc}.\right)\end{array}$$

has advantages over memorizing a concentration–concentration relationship because of the following:

• a.

  The formulation in Eq. (3) conserves mass.

• b.

  The R terms are simpler to memorize because they do not contain the complexity in A.

• c.

  There are fewer concentrations directly influencing S than C, so the machine learning algorithm should be simpler.

The difficulty resides in developing the training and testing sets needed to train and test the machine learning algorithm corresponding to Eq. (4). In principle, we can run a model many times, generate a data set, and then learn that data using machine learning techniques. That is, we can run many models that integrate Eq. (1) to find the relationship between concentrations at two time steps to develop our machine learning training set. But such models do not provide the value of S, and since the chemical system is stiff, the integrators make many complex calls to calculate the right-hand side of Eq. (1) and integrate it. Another way of saying this is that the ΔC is easily available from the models, but the S is not.

If we have many sets of ΔC values, in principle we can invert Eq. (3b) to obtain the corresponding S values. The difficulty with this approach is that there are more elements of S than ΔC, so a conventional inverse cannot be applied. Instead, we employ the generalized inverse of A to obtain S via the relationship

$$\begin{array}{}\text{(5)}& \mathit{S}={\mathbf{A}}^{\text{G}}\mathbf{\Delta }\mathit{C},\end{array}$$

where A^G is the generalized inverse of A. In the case that there are as many fluxes as quantities (A is a square matrix), and the quantities are coupled but linearly dependent (A is full rank), then A^G is the true inverse of A and readily calculable. If the system is overdetermined, where A is a rectangular matrix with more quantities than fluxes but has full column rank, then a left inverse can calculate A^G. An overdetermined system is typical in an aerosol module calculating condensation and evaporation, where fluxes depend on two quantities. However, if A is underdetermined, meaning there are more fluxes than known quantities, or A is otherwise rank-deficient from linear dependency, there is an infinite number of generalized inverses A^G. This means that given values for ΔC, Eq. (5) will not give reliable values for S.

Given sufficient constraints, A^G will be unique and provide the desired values of S that are needed to develop a machine learning training set. Ben-Israel and Greville (2003) show that the inverse can be unique if the solutions, S, are restricted to lie in a subspace that defines the “legal” solutions, and these restrictions are sufficiently constraining. The constrained generalized inverse of A produces solutions, S, that lie in the legal subspace defined by good examples of solutions. S is given by

$$\begin{array}{}\text{(6)}& {\mathbf{A}}_S^{\text{G}}={\mathbf{P}}_S{\left({\mathbf{AP}}_S\right)}^{\text{G}},\end{array}$$

where ${\mathbf{A}}_S^{\text{G}}$ is the generalized inverse of A restricted to the subspace of all possible solutions by the projection P_S, which in turn is defined by a set of basis vectors that define the subspace. Before we discuss obtaining the basis vectors, we first need to discuss how to obtain the projection, P_S.

Assume for the moment that we have the basis vectors S_k. We concatenate them (column-wise) to form the matrix U:

$$\begin{array}{}\text{(7)}& \mathbf{U}=〈{\mathit{S}}_{\mathrm{1}}\mathrm{|}{\mathit{S}}_{\mathrm{2}}\mathrm{|}\mathrm{\dots }\mathrm{|}{\mathit{S}}_k〉.\end{array}$$

The projection onto the subspace defined by these basis vectors S_k and the matrix U is then (Mukhopadhyay, 2014)

$$\begin{array}{}\text{(8)}& {\mathbf{P}}_S=\mathbf{U}{\left({\mathbf{U}}^{+}\mathbf{U}\right)}^{-\mathrm{1}}{\mathbf{U}}^{+},\end{array}$$

where U⁺ is the transpose of U.

Atmospheric-chemistry problems are stiff, so the U⁺U may be ill-conditioned. One source of this ill-conditioning, which can also hamper machine learning tools, is that the concentrations are often orders of magnitude apart. The modules use actual concentrations to make the mechanism easier to understand and debug. Normalizing the concentrations helps with both learning and stiffness and ill-conditioning. Ill-conditioned problems can hamper matrix inversion. Since the S vectors describe the subspace where the solutions must reside, their magnitude does not matter, just their direction. So we normalize the S vectors by dividing by the average of the nonzero values. Mathematically, we form a diagonal square matrix, N_S, with the averages on the diagonal and calculate the normalized S with

$$\begin{array}{}\text{(9)}& {\mathit{S}}_{\text{norm}}={\mathbf{N}}_S^{-\mathrm{1}}\mathit{S}.\end{array}$$

Since N_S is diagonal, the inverse is simply the reciprocal of each diagonal element. The ΔC values are recovered from the S_norm values via

$$\begin{array}{}\text{(10)}& \mathbf{\Delta }\mathit{C}={\mathbf{AN}}_S{\mathit{S}}_{\text{norm}}.\end{array}$$

Atmospheric-chemistry problems are also high-dimensional. Typical air quality models may have 100 to 200 chemical species, and since the vertical-column-mixing timescale is similar to the slower timescales of the chemistry, some models solve the vertical transport and chemistry simultaneously. Since typical air quality models have 10 to 20 vertical cells, the dimension of the problem is 1000 to 4000. Even though the inverse $({\mathbf{U}}^{+}\mathbf{U}{)}^{-\mathrm{1}}$ only has to be calculated once, this inversion may be intractable. Providing that the condition number of U⁺U is not too large, Gram–Schmidt orthonormalization can be performed on the S_i vectors before carrying out Eqs. (7) and (8), in which case they will describe the same subspace, but now the matrix U⁺U will be the identity matrix, which is its own inverse.

Now let us return to the question of how to find the basis vectors that define the “legal” subspace of S. These can be developed by solving Eq. (1) using Euler's method, in which case Eq. (3) becomes

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }{C}_i\sim A_{i,k}{\mathit{R}}_k\left(t\right)\mathrm{\Delta }t\sim A_{i,k}{\mathit{S}}_k,\end{array}$$

that is

$$\begin{array}{}\text{(12)}& {\mathit{S}}_k\sim {\mathit{R}}_k\left(t\right)\mathrm{\Delta }t.\end{array}$$

The value of Δt does not matter since it just changes the length of S_k, not its direction and therefore not its value in describing the subspace. The original module that calculates R_k can be run many times under many conditions to generate a set of S_k vectors that span the subspace. Then locality-preserving projections (LPPs), principal component analysis (PCA), or another similar algorithm can be used to find a minimum set of vectors that define the subspace.

3 Solution procedure for a photochemical module

The following overview aims to put into context the procedure outlined in this paper. The focus of this paper is on deriving and conducting the mass balancing framework and inverse problem detailed in steps 1–9. Steps 10–13 are provided for context: these include machine learning, operator replacement, and benchmarking. In principle, any machine learning algorithm can be used with the framework described here in steps 1–9.

1. Determine which species are active in the photochemical mechanism, that is, not the steady-state or buildup species.

2. From the mechanism, extract the A matrix for these species.

3. Using a representative set of atmospheric concentrations, T, RH, and actinic flux, use Eq. (10) and the photochemical module to generate data that match values of ΔC and S for many values of C, T, RH, and actinic flux for the models operator-splitting time step.

4. Normalize the S vectors by dividing each by the average of its nonzero elements. Use these averages to form the N_S matrix, which relates S to S_norm via Eq. (9).

5. Use the S_norm vectors and Eq. (7) to form the U matrix and then the U⁺U matrix. What is the condition number of the U⁺U matrix? If the system is large and not ill-conditioned, use Gram–Schmidt orthonormalization on the S vectors before calculating U and U⁺U, in which case U⁺U should be an identity matrix or a subset of one.

6. Use Eq. (8) to calculate P_S.

7. Use Eq. (6) to calculate the constrained generalized inverse ${\mathbf{A}}_S^{\text{G}}$.

8. Use Eq. (5) to calculate values of S from the values of ΔC.

9. Compare the values of S obtained from steps 3 and 7 to make sure they are very similar, using the dot product to calculate the angle between them. If they are, we have a good ${\mathbf{A}}_S^{\text{G}}$.

10. Use neural networks or another machine learning algorithm to memorize the S(C) relationship obtained using (a) Eq. (5) and (b) many runs of the mechanism for a wide range of C, T, RH, and actinic flux values.

11. Replace the mechanism with the neural network to calculate S(C) and Eq. (3b) to march forward.

12. Clock the speed improvement.

13. Calculate standard measures of performance such as mass balance, bias, and error compared to runs using the complete mechanism.

4 Photochemical mechanism

We tested the methods described above on the following very simplified set of photochemical reactions used by Michael Kleeman at the University of California, Davis, when teaching the modeling of atmospheric photochemistry. Although this mechanism is abbreviated, it contains the essential components of all atmospheric photochemical mechanisms related to ozone formation: NO_x chemistry, volatile organic compound (VOC) chemistry, and the formation of peroxy radicals from VOC chemistry that then react with NO to form NO₂ and OH, both of which may react to terminate.

The 10 reactions are given in Table 1. The oxygen atom and hydroxyl radical are assumed to be in a steady state, so there are six active species, which are listed in Table 2.

Table 1 Reaction mechanism.

Table 2 Active species.

The resulting A matrix represents the stoichiometry of the reactions, where the rows correspond to each species and the columns to each reaction:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \mathbf{A}=\\ & \left[\begin{array}{ccccccccccc}& R\mathrm{1}& R\mathrm{2}& R\mathrm{3}& R\mathrm{4}& R\mathrm{5}& R\mathrm{6}& R\mathrm{7}& R\mathrm{8}& R\mathrm{9}& R\mathrm{10}\\ {\mathrm{O}}_{\mathrm{3}}& \mathrm{0}& \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{NO}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\mathrm{NO}}_{\mathrm{2}}& -\mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{HCHO}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& -\mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\mathrm{HO}}_{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{2}& \mathrm{0}& \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ {\mathrm{HO}}_{\mathrm{2}}\mathrm{H}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& -\mathrm{1}\end{array}\right]\end{array}.\end{array}$$

As in prior efforts (Kelp et al., 2018; Keller and Evans, 2019), we employed a box model in Julia to generate 60 independent days of output for both ΔC and S, recording data every 6 min. We are interested in the set of S vectors that form a basis describing the subspace that contains the desired S vectors. First, the transformation in Eq. (9) is performed to normalize the sample S vectors. In this example, we use LPP (He and Niyogi, 2004), which is similar to PCA but more robust for this application. Here the LPP yields a basis set of seven vectors, which form the columns of the U matrix:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& \mathbf{U}=\\ & \left[\begin{array}{ccccccc}-\mathrm{0.6869}& \mathrm{0.1334}& -\mathrm{0.2068}& -\mathrm{0.1461}& \mathrm{0.0867}& -\mathrm{0.3715}& \mathrm{0.4761}\\ -\mathrm{0.6869}& \mathrm{0.1334}& -\mathrm{0.2068}& -\mathrm{0.1461}& \mathrm{0.0867}& -\mathrm{0.3715}& \mathrm{0.4761}\\ -\mathrm{0.1877}& -\mathrm{0.1444}& \mathrm{0.0260}& \mathrm{0.1967}& -\mathrm{0.0540}& \mathrm{0.5443}& -\mathrm{0.7353}\\ \mathrm{0.0406}& -\mathrm{0.5849}& -\mathrm{0.0080}& -\mathrm{0.2426}& \mathrm{0.1194}& \mathrm{0.1747}& \mathrm{0.0027}\\ \mathrm{0.0411}& -\mathrm{0.5911}& -\mathrm{0.0081}& -\mathrm{0.2452}& \mathrm{0.1207}& \mathrm{0.1765}& \mathrm{0.0027}\\ -\mathrm{0.0202}& -\mathrm{0.2414}& -\mathrm{0.0069}& -\mathrm{0.0875}& \mathrm{0.0844}& -\mathrm{0.3942}& -\mathrm{0.0570}\\ \mathrm{0.0149}& -\mathrm{0.2555}& \mathrm{0.0154}& \mathrm{0.0759}& -\mathrm{0.1568}& -\mathrm{0.2242}& -\mathrm{0.0519}\\ \mathrm{0.1063}& -\mathrm{0.1359}& \mathrm{0.0844}& \mathrm{0.4800}& -\mathrm{0.7829}& \mathrm{0.4002}& -\mathrm{0.0066}\\ -\mathrm{0.0670}& -\mathrm{0.0762}& \mathrm{0.9337}& \mathrm{0.0057}& \mathrm{0.0069}& -\mathrm{0.0136}& -\mathrm{0.0012}\\ -\mathrm{0.0354}& -\mathrm{0.3227}& -\mathrm{0.1856}& \mathrm{0.7455}& \mathrm{0.5554}& \mathrm{0.0107}& -\mathrm{0.0013}\end{array}\right]\end{array}.\end{array}$$

High condition numbers mean that the matrix inversion is problematic at best. The condition number of U⁺U is approximately 12. This is several orders of magnitude smaller than the same problem but without the N_S transformation of Eq. (9), where the condition number was 1888. This suggests that the N_S transformation has potential to reduce ill-conditioning arising from stiffness in S. The condition number of 12 after the N_S transformation ensures that the inversion needed to make the projection P_S in Eq. (8) is numerically tractable.

The resulting symmetric block diagonal projection P_S is equal to

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& {\mathbf{P}}_S=\\ & \left[\begin{array}{cccccccccc}\mathrm{0.500}& \mathrm{0.500}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.500}& \mathrm{0.500}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{1.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.495}& \mathrm{0.500}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.500}& \mathrm{0.505}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.587}& \mathrm{0.471}& -\mathrm{0.142}& \mathrm{0.001}& -\mathrm{0.005}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.471}& \mathrm{0.462}& \mathrm{0.163}& -\mathrm{0.001}& \mathrm{0.006}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& -\mathrm{0.142}& \mathrm{0.163}& \mathrm{0.951}& \mathrm{0.000}& -\mathrm{0.002}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.001}& -\mathrm{0.001}& \mathrm{0.000}& \mathrm{1.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& -\mathrm{0.005}& \mathrm{0.006}& -\mathrm{0.002}& \mathrm{0.000}& \mathrm{1.000}\end{array}\right]\end{array}.\end{array}$$

And Eq. (6) gives us

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& {\mathbf{A}}_S^{\text{G}}=\\ & \left[\begin{array}{cccccc}\mathrm{2.70}\mathrm{E}\mathrm{1}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}\\ \mathrm{2.70}\mathrm{E}\mathrm{1}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}\\ -\mathrm{4.18}\mathrm{E}\mathrm{1}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}\\ \mathrm{3.63}\mathrm{E}\mathrm{3}& -\mathrm{5.45}\mathrm{E}\mathrm{3}& -\mathrm{1.82}\mathrm{E}\mathrm{3}& \mathrm{3.63}\mathrm{E}\mathrm{3}& \mathrm{5.45}\mathrm{E}\mathrm{3}& \mathrm{3.63}\mathrm{E}\mathrm{3}\\ \mathrm{3.67}\mathrm{E}\mathrm{3}& -\mathrm{5.51}\mathrm{E}\mathrm{3}& -\mathrm{1.84}\mathrm{E}\mathrm{3}& \mathrm{3.67}\mathrm{E}\mathrm{3}& \mathrm{5.51}\mathrm{E}\mathrm{3}& \mathrm{3.67}\mathrm{E}\mathrm{3}\\ -\mathrm{2.84}\mathrm{E}\mathrm{3}& \mathrm{4.26}\mathrm{E}\mathrm{3}& \mathrm{1.42}\mathrm{E}\mathrm{3}& -\mathrm{3.45}\mathrm{E}\mathrm{3}& -\mathrm{4.26}\mathrm{E}\mathrm{3}& -\mathrm{2.84}\mathrm{E}\mathrm{3}\\ \mathrm{5.37}\mathrm{E}\mathrm{2}& -\mathrm{5.37}\mathrm{E}\mathrm{2}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}\\ \mathrm{0.0000}& -\mathrm{1.78}\mathrm{E}\mathrm{3}& -\mathrm{1.78}\mathrm{E}\mathrm{3}& \mathrm{0.0000}& \mathrm{0.0000}& \mathrm{0.0000}\\ -\mathrm{9.24}\mathrm{E}\mathrm{5}& \mathrm{1.14}\mathrm{E}\mathrm{6}& \mathrm{2.16}\mathrm{E}\mathrm{5}& -\mathrm{9.24}\mathrm{E}\mathrm{5}& -\mathrm{1.14}\mathrm{E}\mathrm{6}& -\mathrm{9.24}\mathrm{E}\mathrm{5}\\ \mathrm{9.81}\mathrm{E}\mathrm{4}& -\mathrm{1.21}\mathrm{E}\mathrm{5}& -\mathrm{2.29}\mathrm{E}\mathrm{4}& \mathrm{9.81}\mathrm{E}\mathrm{4}& \mathrm{1.21}\mathrm{E}\mathrm{5}& \mathrm{4.58}\mathrm{E}\mathrm{4}\end{array}\right]\end{array}.\end{array}$$

Since U and ${\mathbf{A}}_S^{\text{G}}$ have seven and six independent columns, respectively, but 10 rows, and the row rank is equal to the column rank, there must be linearly dependent rows. One manifestation of this is that the first two rows of U and ${\mathbf{A}}_S^{\text{G}}$ are identical or nearly so. The S values computed from ${\mathbf{A}}_S^{\text{G}}$ may not be identical to the original S corresponding to the ΔC values. However, all S values calculated from Eq. (5) using the above ${\mathbf{A}}_S^{\text{G}}$ are “legal”, in other words, within the subspace defined by the basis set U. Furthermore, the inverse ${\mathbf{A}}_S^{\text{G}}$ by definition satisfies ${\mathbf{AA}}_S^{\text{G}}=I$ so that even if a calculated S is not identical to the S from the original box model output, it can be used in Eq. (3b) to return a ΔC identical to that of the box model output.

5 Conclusions

Large models of the environment require the solution of large systems of equations over long periods of time. These models consume vast quantities of computational resources, so approximations are necessarily employed so that the models are computationally tractable. Machine learning tools can be used to dramatically improve the speed of these models, enabling a more faithful representation of the physics and chemistry while also improving runtime performance. But this field is in its infancy. To help facilitate the use of machine learning tools in these environmental models, we have developed a framework that (a) enables machine learning algorithms to learn flux terms, assuring that conservation principles dictated by the physics and chemistry are adhered to, and (b) allows parameters easily calculated by geophysical models to be used to back-calculate these flux terms that can then be used to train the machine learning algorithm of choice. Applications of this framework in environmental models include any process where conservation principles apply, such as conservation of atoms in chemical reactions; conservation of molecules during phase change; and conservation of energy in, say, radiative transfer calculations.

Appendix A: Glossary of symbols

Ci(t)	concentration at time t
Ci(t+Δt)	concentration at time t+Δt
$\mathrm{\Delta }{C}_i\equiv C_i(t+\mathrm{\Delta }t)-C_i\left(t\right)$
$i=\mathrm{1},n$	the number of molecular species
Δt	operator-splitting time step
Rj(t)	contribution to ΔCi from each reaction
$S_j\left(t\right)={\int }_t^{t+\mathrm{\Delta }t}{R}_j\left(\mathit{\tau }\right)\text{d}\mathit{\tau }$
$j=\mathrm{1},m$	the number of reactions, m>n
A	a sparse stoichiometry matrix relating ΔCi to Sj; most element values are 0, 1, or −1
AG	generalized inverse of A
${\mathbf{A}}_S^{\text{G}}$	constrained generalized inverse of A