Technical Note: An efficient method for accelerating the spin-up process for process-based biogeochemistry models

To better understand the role of terrestrial ecosystems in the global carbon cycle and their feedbacks to the global climate system, process-based biogeochemistry models need to be improved with respect to model parameterization and model structure. To achieve these improvements, the spin-up time for those differential equation-based models needs to be shortened. Here, an algorithm for a fast spin-up was developed by finding the exact solution of a linearized system representing the cyclo-stationary state of a model and implemented in a biogeochemistry model, the Terrestrial Ecosystem Model (TEM). With the new spin-up algorithm, we showed that the model reached a steady state in less than 10 years of computing time, while the original method requires more than 200 years on average of model run. For the test sites with five different plant functional types, the new method saves over 90 % of the original spin-up time in site-level simulations. In North American simulations, average spin-up time savings for all grid cells is 85 % for either the daily or monthly version of TEM. The developed spin-up method shall be used for future quantification of carbon dynamics at fine spatial and temporal scales.

1 Introduction

Biogeochemistry models contain state variables representing various pools of carbon and nitrogen and a set of flux variables representing the element and material transfers among different state variables. Model spin-up is a step to get biogeochemistry models to a steady state for those state and flux variables (McGuire et al., 1992; King, 1995; Johns et al., 1997; Dickinson et al., 1998). Spin-up normally uses cyclic forcing data to force the model run and reach a steady state, which will be used as initial conditions for model transient simulations. The steady state is reached when modeled state variables show a cyclic pattern or a constant value and often requires a significant amount of computation time, which needs to be accelerated for regional and global simulations at fine spatial and temporal scales.

Spin-up is normally achieved by running the model repeatedly using one or several decades of meteorological or climatic data until a steady state is reached. The step could require that the model repeatedly run for more than 2000 annual cycles in some extreme cases. Specifically, the model will check the stability of the simulated carbon and nitrogen fluxes as well as state variables with specified threshold values. For instance, the model will check if the simulated annual net ecosystem production (NEP) is less than 1 g C m⁻² yr⁻¹ (McGuire et al., 1992). Another method to reach a steady state is to obtain the analytical solutions (King et al., 1995; Comins, 1997), which might also take a significantly long time.

For different biogeochemistry models, spin-up could take hundreds and thousands of years to reach a stability, normally longer than the model projection period (Thornton and Rosenbloom, 2005). Therefore, a more efficient method to reach the steady state will speed up the entire model simulation. Recently, a semi-analytical method (Xia et al., 2012) has been adapted to a carbon–nitrogen coupled model to speed up the spin-up process. The idea is to obtain an analytical solution very close to a steady condition, then start spin-up from the solution, which could significantly reduce spin-up time. This technique did not reach a cyclic pattern for state and flux variables and required an additional spin-up process to achieve the steady state. However, Lardy et al. (2011) and Martin et al. (2007) have implemented their spin-up methods for a linear problem of soil carbon dynamics including their seasonal cycles.

Here we developed a method to accelerate the spin-up process in a nonlinear model. We tested the method for representative PFTs and North America with both daily and monthly versions of the Terrestrial Ecosystem Model (TEM; Zhuang et al., 2003). In addition, we compared the performance of our algorithms with the semi-analytical version of Xia et al. (2012). The new algorithms will help us conduct very high spatial and temporal resolution simulations with process-based biogeochemistry models in the future.

2 Method

2.1 TEM description

We used a process-based biogeochemistry model, TEM (Zhuang et al., 2003), as a test bed to demonstrate the performance of the new algorithms of spin-up. TEM simulates the dynamics of ecosystem carbon and nitrogen fluxes and pools (McGuire et al., 1992; Zhuang et al., 2010, 2003). It contains five state variables: carbon in living vegetation (C_v), nitrogen in living vegetation (N_v), organic carbon in detritus and soils (C_s), organic nitrogen in detritus and soils (N_s), and available inorganic soil nitrogen (N_av). Carbon and nitrogen dynamics in TEM are governed by the following equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{v}}}{\mathrm{d}t}}=\text{GPP}-R_{\mathrm{A}}-L_{\mathrm{C}},\text{(2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{v}}}{\mathrm{d}t}}=\text{NUPTAKE}-L_{\mathrm{N}},\text{(3)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{s}}}{\mathrm{d}t}}=L_{\mathrm{c}}-R_{\mathrm{H}},\text{(4)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{s}}}{\mathrm{d}t}}=L_{\mathrm{N}}-\text{NETNMIN},\text{(5)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{av}}}{\mathrm{d}t}}=\text{NINPUT}+\text{NETNMIN}-\text{NLOST}-\text{NUPTAKE},\end{array}$$

where GPP is gross primary production, R_A is autotrophic respiration, L_Cis carbon in litterfall, NUPTAKE is nitrogen uptake by vegetation, L_N is nitrogen in litterfall, R_H is heterotrophic respiration, NETNMIN is net rate of mineralization of soil nitrogen, NINPUT is nitrogen input from the outside ecosystem, and NLOST is nitrogen loss from the ecosystem. Key carbon fluxes are defined as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{GPP}=C_{max}f\left(\text{PAR}\right)f\left(\text{PHENOLOGY}\right)f\left(\text{FOLIAGE}\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }f\left(T\right)f\left(C_{\mathrm{a}},G_{\mathrm{v}}\right)f\left(\text{NA}\right)f\left(\text{FT}\right),\text{(7)}& {\displaystyle }& {\displaystyle }\text{NPP}=\text{GPP}-R_{\mathrm{A}},\text{(8)}& {\displaystyle }& {\displaystyle }\text{NEP}=\text{GPP}-R_{\mathrm{A}}-R_{\mathrm{H}}.\end{array}$$

For detailed GPP definition, see Zhuang et al. (2003). NEP will be near zero when the ecosystem reaches a steady state. Therefore, the spin-up goal is to keep running the model driven with repeated climate forcing data until NEP is close to zero with a certain tolerance value (e.g., 0.1 g C m⁻² yr⁻¹).

2.2 Spin-up acceleration method

TEM can be reformulated as

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}\mathit{x}}{\mathrm{d}t}}=g\left(\mathit{x},t\right)+\mathit{h}\left(t\right),\end{array}$$

where x is a vector of state variables (e.g., C_V); h(t) is the vector of carbon–nitrogen input from the atmosphere (such as nitrogen input), independent of x; and g(x,t) is the process rate function of element pools (e.g., GPP).

By linearizing the model in terms of pools, we could obtain

$$\begin{array}{}\text{(10)}& {\displaystyle }g\left(\mathit{x},t\right)=g\left({\mathit{x}}_{\mathrm{0}},t\right)+\mathbf{J}\left(\mathit{x}-{\mathit{x}}_{\mathrm{0}}\right),\end{array}$$

where x₀ is initial pool sizes and J is the Jacobian matrix of the process rate:

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathbf{J}={\displaystyle \frac{\mathrm{d}g}{\mathrm{d}x}}=\left(\begin{array}{ccc}\frac{\partial g_{\mathrm{1}}}{\partial x_{\mathrm{1}}}& \mathrm{\cdots }& \frac{\partial g_{\mathrm{1}}}{\partial x_n}\\ \mathrm{⋮}& \mathrm{⋮}\\ \frac{\partial g_n}{\partial x_{\mathrm{1}}}& \mathrm{\cdots }& \frac{\partial g_n}{\partial x_n}\end{array}\right),\end{array}$$

where g represents g(x,t). x_n represents each state variable in the TEM (e.g., V_C). The numerical discretization of Eq. (9) is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{x}_{i,k}-x_{i,k-\mathrm{1}}=\mathit{\tau }\cdot {\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot x_{i,k-\mathrm{1}}\\ \text{(12)}& {\displaystyle }& {\displaystyle }+\mathit{\tau }\left(g\left(x_{\mathrm{0}},k-\mathrm{1}\right)-{\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot x_{\mathrm{0},k-\mathrm{1}}+h_{k-\mathrm{1}}\right),\end{array}$$

where τ is time step (month), x_i, k is the pool x_i size at time k, and J_k−½ is a Jacobian matrix at time step $k-\frac{\mathrm{1}}{\mathrm{2}}$. Here $\frac{\mathrm{1}}{\mathrm{2}}$ refers to the half time step in the middle of a month, at which values of J are calculated as the mean value at time steps k and k−1. x_i,0 refers to the initial pool x_i size.

We introduce

$$\begin{array}{}\text{(13)}& {\displaystyle }{f}_{k-\mathrm{1}}=g\left(x_{\mathrm{0}},k-\mathrm{1}\right)-{\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot x_{\mathrm{0},k-\mathrm{1}}+h_{k-\mathrm{1}}.\end{array}$$

The Eq. (12) can then be written as

$$\begin{array}{}\text{(14)}& {\displaystyle }{x}_{i,k}-x_{i,k-\mathrm{1}}=\mathit{\tau }\cdot {\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot x_{i,k-\mathrm{1}}+\mathit{\tau }\cdot f_{k-\mathrm{1}},\end{array}$$

where J_k−½ is a Jacobian matrix at the time step $k-\frac{\mathrm{1}}{\mathrm{2}}$. After running a large number of annual cycles, the model approaches a cyclo-stationary state, which can be expressed by condition $x_{T+i}=x_i$, where T is the number of time steps in one cycle, for example, when spin-up is made at a monthly time step using monthly climatology of temperature, precipitation and other forcing data, T=12, and $\overrightarrow{x_{\mathrm{1}}}$ is the size of carbon pools on 1 January, while J^1.5 is the matrix of mean process rate constants for January.

By introducing

$${\displaystyle }{A}_k=\mathit{\tau }\cdot {\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}},y_k=\mathit{\tau }{f}_{k-\mathrm{1}},B_k=\mathbf{I},C_k=\mathbf{I}+A,$$

where I is an identity matrix, Eq. (12) can be further written as

$$\begin{array}{}\text{(15)}& {\displaystyle }-C_k\cdot x_{i,k-\mathrm{1}}+B_k\cdot x_{i,k}=y_k.\end{array}$$

The cyclic boundary condition is $x_{\mathrm{1}}=x_{T+\mathrm{1}}$.

Then Eq. (13) will become

$$\begin{array}{}\text{(15a)}& {\displaystyle }-C_{\mathrm{1}}\cdot x_{i,T}+B_{\mathrm{1}}\cdot x_{i,\mathrm{1}}=y_{\mathrm{1}}.\end{array}$$

Thus Eq. (15) and (15a) become a formulation of a linear problem with T unknown vectors $\overrightarrow{x_T}$, which can be solved using LU (lower and upper) decomposition or Gaussian elimination (Martin et al., 2007). Xia et al. (2012; see Eq. 4) and Kwon and Primeau (2006) also had linear equations for a steady state, but conducted the model simulation at an annual time step. Going for annual average form reduces the size of the problem and prevents Xia et al. (2012) from obtaining the exact solution of the problem including seasonal cycle (see their Eqs. 15, 15a). While our new approach runs the model at a monthly time step with the cyclic boundary conditions for state variables x, it still targets a steady state for the ecosystem at an annual time step instead of a monthly time step.

2.3 Numerical implementation

Equation (15a) is explicitly expressed as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left(\begin{array}{cccccccc}B& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& -C\\ -C& B& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -C& B& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& -C& B& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& -C& B& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& -C& B\end{array}\right)\\ \text{(16)}& {\displaystyle }& {\displaystyle }\times \left(\begin{array}{c}{x}_{\mathrm{1}}\\ x_{\mathrm{2}}\\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ x_T\end{array}\right)=\left(\begin{array}{c}{y}_{\mathrm{1}}\\ y_{\mathrm{2}}\\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ y_T\end{array}\right).\end{array}$$

Equation (16) can be shown in the form Mx=Y.

We apply the Gaussian elimination to the upper block that reduces M to a lower triangular form (Fig. 1). The resulting matrix is lower diagonal:

$$\begin{array}{}\text{(17)}& {\displaystyle }{M}^{\prime }=\left(\begin{array}{cccccccc}B{}^{\prime }& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -C& B& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -C& B& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& -C& B& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& -C& B& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& -C& B\end{array}\right).\end{array}$$

Equation (16) is thus reduced to the form $M^{\prime }x=Y^{\prime }$, where M^′ is the lower diagonal, and solution of Eqs. (15a) and (16) will be readily obtained for x.

Figure 1 Algorithms and procedures of the new spin-up method.

2.4 Algorithm implementation for TEM

In the original TEM, carbon fluxes can be defined as

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle }\text{NPP}=\text{GPP}-\text{MR}-\text{GR},\text{(19)}& {\displaystyle }& {\displaystyle }\text{MR}=V_{\mathrm{C}}\cdot K_T,\text{(20)}& {\displaystyle }& {\displaystyle }\text{GR}=\left\{\begin{array}{ll}\mathrm{0.25}\cdot \left(\text{GPP}-\text{MR}\right),& \text{if}\text{GPP}>\text{MR}\\ \mathrm{0},& \text{otherwise}\end{array}\right.,\end{array}$$

where NPP is defined as the difference of GPP and plant maintenance respiration (MR) and growth respiration (GR). MR is assumed as a function of C_V and temperature (K_T). Here we revised MR calculation:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{MR}\\ \text{(21)}& {\displaystyle }& {\displaystyle }=\left\{\begin{array}{ll}{V}_{\mathrm{C}}\cdot K_T,& \mathrm{\text{if}}\text{GPP}>V_{\mathrm{C}}\cdot K_T\\ \mathrm{0.75}\cdot V_{\mathrm{C}}\cdot K_T+\mathrm{0.25}\cdot \text{GPP},& \text{otherwise}.\end{array}\right.\end{array}$$

Table 1 Test sites for new spin-up algorithms.

Figure 2 The time for NEP (g C yr⁻¹ m⁻²) to reach a steady state with the original spin-up method at the Harvard Forest site. x represents model simulation years.

The NEP is defined as the difference between NPP and heterotrophic respiration (R_H).

The basic workflow to implement the method is (1) linearizing TEM first to obtain a sparse matrix with n variable (for TEM, n = 5) system, (2) performing Gaussian elimination for the linear system and (3) solving the sparse matrix to acquire the state variable values (Fig. 1). To adapt this method to a daily version of TEM, we changed the cyclic condition T from 12 to 365. The other steps are the same as for the monthly version. We tested the new method for the carbon-only version and carbon–nitrogen coupled version of TEM for different PFTs (Table 1). Specifically, for the carbon-only version, we only solved the differential equations that govern the carbon dynamics, while for the carbon–nitrogen coupled version, we solved the differential equations that govern both carbon and nitrogen dynamics in the system. For both versions, the spin-up process strives to reach a steady state for carbon pools and fluxes.

Figure 3 Simulated NEP (g C m⁻² yr⁻¹) with the original spin-up method after different spin-up times of (a) 50, (b) 100, (c) 150 and (d) 200 years. After these spin-up times, 63, 89, 93 and 98 % of grids reach their steady states, respectively.

3 Results and discussion

At the Harvard Forest site, the traditional spin-up method took 564 years to reach the steady state for both the carbon-only and coupled carbon–nitrogen simulations with an annual NEP of less than 0.1 g C m⁻² yr⁻¹ (Fig. 2). In contrast, the improved method took 72 years for the carbon only and 122 for the coupled carbon–nitrogen simulations. For carbon and nitrogen pools, it took another 45 years (equivalent cyclic time) to reach a steady state with a NEP of less than 0.1 g C m⁻² yr⁻¹. In comparison with the traditional spin-up method (Zhuang et al., 2003), the new method saved 65 % of computational time to reach the steady state in the carbon-only simulations (Table 2). The differences in steady-state carbon pools between using the new method and traditional spin-up methods were small (less than 0.85 %). Similarly, for the coupled carbon–nitrogen simulations, the new method saves a similar amount of time to reach the steady state.

For all seven test sites, the original spin-up method in TEM takes 204–564 years (1.1–2.5 s of computing time) to reach the steady state at different sites. In contrast, our new method only takes 0.3–0.6 s, while the semi-analytical method (Xia et al., 2012) will need 0.5–0.9 s to reach the steady state at different sites (Table 2). Compared to the original spin-up method, the new method is not only faster but also computationally stable.

The time of spin-up to reach a steady state of NEP varied for different PFT grids using the original method (Fig. 2). In general, to allow 98 % of grid cells to reach their steady states of NEP, it takes 250 annual model runs, while the new method will only need on average of 0.6 s (equivalent to 60-year annual model runs with the original method) (Fig. 3). For regional tests in North America, we found that the average saving time with the new method with monthly TEM is 25, 32 and 22 % for Alaska, Canada and the conterminous US, respectively.

Table 2 Spin-up time comparison for different methods for seven study sites; seconds represent real computation time, and years refer to the annual spin-up cycles.

To compare the performance of the new method with other existing methods, we adapted the semi-analytical method (Xia et al., 2012) to the TEM model. To do that, we first revised the TEM model structure to

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{\mathrm{d}P\left(t\right)}{\mathrm{d}t}}=\mathit{\epsilon }\mathbf{AC}P\left(t\right),\end{array}$$

where P(t) is a vector of pools in TEM (e.g., C_V and C_S). ε is a scalar. A is a pool transfer matrix (in which A_ij represents the fraction of carbon transfer from pool j to i). C is a diagonal matrix with pool components (where diagonal components quantify the fraction of carbon left from the state variables after each time step). With this method, we obtained an analytical solution for the intermediate state. We then kept running TEM with the traditional spin-up process. Specifically, we started TEM simulation to estimate the state variable values. Based on these values, the spin-up runs were conducted to reach the final steady state. We found that the semi-analytical solution is better than the original spin-up method but slower than the new method proposed in this study (Table 2).

The TEM model has a relatively small set of state variables for carbon and nitrogen. The version we used is TEM 5.0, which has only five state variables (Zhuang et al., 2003). Thus, the linearization process is relatively easy and the matrix size is relatively small; consequently, the computing is not a burden. To accelerate the spin-up for multiple soil carbon pool models with relatively simple and linear decomposition processes, implementing our method will still be relatively easy but will take a great amount of computing time to equilibrate. For models such as CLM, multiple methods have been tested to accelerate their spin-up process (e.g., Fang et al., 2015), the direct analytical solution method introduced in this study might be time consuming to achieve.

4 Summary

We developed a new method to speed up the spin-up process in process-based biogeochemistry models. We found that the new method shortened 90 % of the spin-up time using the traditional method. For regional simulations in North America, average spin-up time saving is 85 % for either daily or monthly version of TEM. We consider our method is a general approach to accelerate the spin-up process for process-based biogeochemistry models. As long as the governing equations of the models can be formulated as the form in Eq. (9), the algorithm could be adopted accordingly. Our method will significantly help future carbon dynamics quantification with biogeochemistry models at fine spatial and temporal scales.