(a-mathematical-formulation)=
# Mathematical formulation
The Climate Module starts from the global emissions of CO2, CH4, and N2O, as generated by the TIMES global model, and proceeds to compute successively:
- the changes in CO2, CH4, and N2O concentrations via three separate sets of equations,
- the total change (over pre-industrial times) in atmospheric radiative forcing resulting from the three gases plus an exogenously specified additional forcing resulting from other causes (other anthropogenic and/or natural causes, as defined by the user), and
- the temperature changes (over pre-industrial times) in two reservoirs (surface and deep ocean).
The Climate Equations used to perform these calculations were initially adapted from Nordhaus and Boyer (1999), who proposed linear recursive equations for calculating concentrations and temperature changes based on the CO2 life cycle. These linear equations give results that are good approximations of those obtained from more complex climate models (Drouet *et al.*, 2004; Nordhaus and Boyer, 1999). The non-linear radiative forcing equation used by these authors and in TIMES is the same as the one used in most models. The choice of the Nordhaus and Boyer's climate equations is motivated by the simplicity of their approach and by the fact that their climate module is well-documented and acceptably accurate. In our implementation, the forcing equation has been replaced by a linear approximation whose values closely approach the exact ones as long as the useful range is carefully selected. This was done in order to keep the entire model linear, and therefore to allow the user to set constraints on forcing and on temperature as well as on concentrations and on emissions.
Rigorously, the concentration and forcing equations used in the climate module are applicable only to CO2 emissions, since the concentration equations simulate the carbon cycle. In order to model other GHGs, one way is to use these same equations, while replacing CO2 emissions by CO2-equivalent emissions of any number of gases endogenous to the model. However, a more detailed and generally preferable approach is to model separately the life cycle of each endogenous emission separately, and this is the approach used in TIMES. The additional forcing due to the remaining (non endogenous) emissions, is accounted for via an exogenous forcing quantity directly defined by the user.
We now describe the mathematical equations used at each of the three steps of the climate module.
## Concentrations (accumulation of CO2, CH4, N2O)
a\) CO2 accumulation is represented as the linear three-reservoir model below[^46]: the atmosphere, the quickly mixing upper ocean + biosphere, and the deep ocean. CO2 flows in both directions between adjacent reservoirs. The 3-reservoir model is represented by the following 3 equations when the step of the recursion is equal to one year:
$$M_{atm}(y) = E(y) + (1 - \varphi_{atm-up}) M_{atm}(y-1) + \varphi_{up-atm} M_{up}(y-1)$$ (2-7-1)
$$M_{up}(y) = (1 - \varphi_{up-atm} - \varphi_{up-lo}) M_{up}(y-1) + \varphi_{atm-up} M_{atm}(y-1) + \varphi_{lo-up} M_{lo}(y-1)$$ (2-7-2)
$$M_{lo}(y) = (1 - \varphi_{lo-up}) M_{lo}(y-1) + \varphi_{up-lo} M_{up}(y-1)$$ (2-7-3)
with
- $M_{atm}(y)$, $M_{up}(y)$, $M_{lo}(y)$: masses of CO2 in atmosphere, in a quickly mixing reservoir representing the upper level of the ocean and the biosphere, and in deep oceans (GtC), respectively, in year $y$ (GtC)
- $E(y-1)$: CO2 emissions in previous year (GtC)
- $\varphi_{ij}$: transport rate from reservoir $i$ to reservoir $j$ ($i, j = atm, up, lo$) from year $y-1$ to $y$
b\) CH4 accumulation is represented by a so-called single-box model in which the atmospheric methane concentration obeys the following equations assuming a constant annual decay rate of the anthropogenic concentrations $\Phi_{CH4}$ (whereas the natural concentration is assumed in equilibrium):
$$CH4_{atm}(y) = (1 - \Phi_{CH4}) \cdot CH4_{atm}(y - 1) + EA_{CH4}(y)$$ (2-7-1a)
$$CH4_{up}(y) = CH4_{up}(y - 1)$$ (2-7-1b)
$$CH4_{tot}(y) = CH4_{atm}(y) + CH4_{up}(y)$$ (2-7-1c)
where
- $CH4_{atm}$, $CH4_{up}$, and $EA_{CH4}$ are respectively: the atmospheric concentration, the natural concentration[^47] (both expressed in Mt), and the anthropogenic emission of CH4 (expressed in Mt/yr). $EA_{CH4}$ is generated within the model, but $CH4_{up}$ is fully exogenous (see values for CH4-UP and CH4-ATM in {numref}`cli-parameters`). All quantities are indexed by year.
- $d_{CH4}$ = 2.84 (the density of CH4, expressed in *Mt/ppbv*) is then used to convert concentration in Mt into ppbv.
- $1 - \Phi_{CH4}$ is the one-year retention rate of CH4 in the atmosphere, see {numref}`cli-user-input-parameters`.
c\) N2O accumulation is also represented by a single-box model in which the atmospheric N2O concentration obeys the following equations:
$$N2O_{atm}(y) = (1 - \Phi_{N2O}) \cdot N2O_{atm}(y - 1) + EA_{N2O}(y)$$ (2-7-2a)
$$N2O_{up}(y) = N2O_{up}(y - 1)$$ (2-7-2b)
$$N2O_{tot}(y) = N2O_{atm}(y) + N2O_{up}(y)$$ (2-7-2c)
where
- $N2O_{atm}$, $N2O_{up}$, and $EA_{N2O}$ are respectively: the atmospheric concentration, the natural concentration (both expressed in Mt), and the anthropogenic emission of N2O (expressed in Mt/yr). $EA_{N2O}$ is generated within the model, but $N2O_{up}$ is fully exogenous (see values for N2O-UP and N2O-ATM in {numref}`cli-parameters`). All quantities are indexed by year.
- $d_{N2O}$ = 7.81 (the density of *N2O*, expressed in *Mt/ppbv*) is then used to convert concentration in Mt to ppbv units.
- $1 - \Phi_{N2O}$ is the one-year retention rate of N2O in the atmosphere, see {numref}`cli-user-input-parameters`.
*Note*: For both CH4 and N2O, the total atmospheric concentrations (UP+ATM) are used in the forcing expressions (see below) and are reported in the results.
(a-radiative-forcing)=
## Radiative forcing
We assume, as is routinely done in atmospheric science, that the atmospheric radiative forcing caused by the various gases are additive (IPCC, 2007). Thus:
$$\Delta F(y) = \Delta F_{CO2}(y) + \Delta F_{CH4}(y) + \Delta F_{N2O}(y) + EXOFOR(y)$$ (2-7-4)
We now explain these four terms.
a\) The relationship between CO2 accumulation and increased radiative forcing, $\Delta F_{CO2}(y)$, is derived from empirical measurements and climate models (IPCC 2007).
$$\Delta F_{CO2}(y) = \gamma \times \frac{\ln(M_{atm}(y)/{M_{0}})}{\ln2}$$ (2-7-4a)
where:
- $M_0$ (i.e. CO2ATM_PRE_IND) is the pre-industrial (circa 1750) reference atmospheric concentration of CO2 = 596.4 GtC
- $\gamma$ is the radiative forcing sensitivity to atmospheric CO2 concentration doubling = 3.7 W/m2
b) The radiative forcing due to atmospheric CH4 is given by the following expression (IPCC, 2001)
$$\Delta F_{CH4}(y) = 0.036 \cdot \left( \sqrt{CH4_{y}} - \sqrt{CH4_{0}} \right) - \left\lbrack f(CH4_{y},N2O_{0}) - f(CH4_{0},N2O_{0}) \right\rbrack$$ (2-7-4b)
c) The radiative forcing due to atmospheric N2O is given by the following expression (IPCC, 2001)
$$\Delta F_{N2O}(y) = 0.12 \cdot \left( \sqrt{N2O_{y}} - \sqrt{N2O_{0}} \right) - \left\lbrack f(CH4_{0},N2O_{y}) - f(CH4_{0},N2O_{0}) \right\rbrack$$ (2-7-4c)
where:
$$f(x,y) = 0.47 \cdot \ln\left\lbrack 1 + 2.01 \cdot 10^{- 5} \cdot (xy)^{0.75} + 5.31 \cdot 10^{- 15} \cdot x(xy)^{1.52} \right\rbrack$$ (2-7-4d)
Note that the $f(x,y)$ function, which quantifies the cross-effects on forcing of the presence in the atmosphere of both gases (CH4 and N2O), is not quite symmetrical in the two gases. As usual, the 0 subscript indicates the pre-industrial times (1750)
d) $EXOFOR(y)$ is the increase in total radiative forcing at period *t* relative to pre-industrial level due to GHGs that are not represented explicitly in the model. Units = W/m2. In Nordhaus and Boyer (1999), only emissions of CO2 were explicitly modeled, and therefore O(y) accounted for all other GHG's. In TIMES, N2O and CH4 are fully accounted for, but some other substances are not (e.g. CFC's, aerosols, ozone, etc.). Therefore, our values for $EXOFOR(y)$ will differ from those in Nordhaus and Boyer. It is the modeler's responsibility to include in the calculation of $EXOFOR(y)$ only the forcings from those gases and other causes that are not modeled. {numref}`tiam-world-exoforcing-example` shows a possible trajectory for $EXOFOR$.
The parameterization of the three forcing equations (4a, 4b, 4c) is not controversial and relies on the results reported by Working Group I in the IPCC. IPCC (2001, Table 6.2, p.358) provides a value of 3.7 for $\gamma$, smaller than the one used by Nordhaus and Boyer ($\gamma = 4.1$). We have adopted this lower value of 3.7 W/m2 as default in TIMES. Users are free to experiment with other values of the $\gamma$ parameter. The same reference provides the entire expressions for all three forcing equations.
## Linear approximations
In TIMES, each of the three forcing expressions is replaced by a linear approximation, in order to preserve linearity of the entire model. All three forcing expressions (4a, 4b, 4c) happen to be concave functions. Therefore, two linear approximations are obvious candidates. The first one is an approximation from below, consisting of the chord of the graph between two selected points. The second one has the same slope as the chord and is tangent to the graph, thus approximating the function from above. The final approximation is taken to be the arithmetic average of the two approximations. These linear expressions are easily derived once a range of interest is defined by the user.
As an example, we derive below the linear approximation for the CO2 forcing expression. The other approximations are obtained in a similar manner, and the parameters of the linear approximations are shown in the next section.
:::{admonition} Linear approximation for the CO2 forcing expression
First, an interval of interest for the concentration M must be selected by the user. The interval should be wide enough to accommodate the anticipated values of the concentrations, but not so wide as to make the approximation inaccurate. We denote the interval $(M_1,M_2)$.
Next, the linear forcing equation is taken as the half sum of two linear expressions, which respectively underestimate and overestimate the exact forcing value. The underestimate consists of the chord of the logarithmic curve, whereas the overestimate consists of the tangent to the logarithmic curve that is parallel to the chord.
By denoting the pre-industrial concentration level as $M_0$, the general formulas for the two estimates are as follows:
*Overestimate:*
$$F_{1}(M) = \frac{\gamma}{\ln 2} \cdot \left\lbrack \ln(\frac{\gamma}{slope \cdot \ln(2) \cdot M_{0}}) - 1 \right\rbrack + slope \cdot M$$ (2-7-5)
*Underestimate*:
$$F_{2}(M) = \gamma \cdot \ln(M_{1}/M_{0})/\ln 2 + slope \cdot (M - M_{1})$$ (2-7-6)
*Final approximation*:
$$F_{3}(M) = \frac{F_{1}(M) + F_{2}(M)}{2}$$ (2-7-7)
where:
$$slope = \gamma \cdot \frac{\ln(M_{2}/M_{1})/\ln 2}{(M_{2} - M_{1})}$$
:::
The linearized forcing expression implemented in TIMES is the average of the two linear estimates.
## Temperature increase
In the TIMES Climate Module as in many other integrated models, climate change is represented by the global mean surface temperature. The idea behind the two-reservoir model is that a higher radiative forcing warms the atmospheric layer, which then quickly warms the upper ocean. In this model, the atmosphere and upper ocean form a single layer, which slowly warms the second layer consisting of the deep ocean.
$$\Delta T_{up}(y) = \Delta T_{up}(y-1) + \sigma_1\{F(y) - \lambda \Delta T_{up}(y-1) - \sigma_2 [\Delta T_{up}(y-1) - \Delta T_{low}(y-1)]\}$$ (2-7-8)
$$\Delta T_{low}(y) = \Delta T_{low}(y-1) + \sigma_3[\Delta T_{up}(y-1) - \Delta T_{low}(y-1)]$$ (2-7-9)
with
- $\Delta T_{up}$ = globally averaged surface temperature increase above pre-industrial level,
- $\Delta T_{low}$ = deep-ocean temperature increase above pre-industrial level,
- $\sigma_1$ = 1-year speed of adjustment parameter for atmospheric temperature (also known as the *lag* parameter),
- $\sigma_2$ = coefficient of heat loss from atmosphere to deep oceans,
- $\sigma_3$ = 1-year coefficient of heat gain by deep oceans,
- $\lambda$ = feedback parameter (climatic retroaction). It is customary to write $\lambda$ as $\lambda=\gamma/C_s$, $C_s$ being the *climate sensitivity* parameter, defined as the change in equilibrium atmospheric temperature induced by a doubling of CO2 concentration.
:::{admonition} Remark
In contrast with most other parameters, the value of $C_s$ is highly uncertain, with a possible range of values from 1°C to 10°C. This parameter is therefore a prime candidate for sensitivity analysis, or for treatment by probabilistic methods such as stochastic programming. In {numref}`cli-parameters`, a best estimate value of 2.9 °C is shown, as per IPCC (2001, 2007).
:::
In the next section we describe all the input parameters required to define the climate equations and those needed to define climate constraints. With few exceptions (such as the densities of the gases), all parameters are modifiable by the user, should the need arise. We also provide {numref}`cli-parameters` summarizing the default values of the parameters.
[^46]: There exists another well-known representation of CO2 accumulation equations, using a five-box model.
[^47]: Note that the subscripts *atm* and *up*, which for the CO2 equations referred to the atmosphere and upper reservoirs, have been reused for the CH4 and N2O equations to stand for anthropogenic and natural concentrations.