Journal topic
Biogeosciences, 15, 3967–3973, 2018
https://doi.org/10.5194/bg-15-3967-2018
Biogeosciences, 15, 3967–3973, 2018
https://doi.org/10.5194/bg-15-3967-2018

Technical note 03 Jul 2018

Technical note | 03 Jul 2018

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

Technical Note: An efficient method for accelerating the spin-up process for process-based biogeochemistry models
Yang Qu1, Shamil Maksyutov2, and Qianlai Zhuang1,3 Yang Qu et al.
• 1Department of Earth, Atmospheric, and Planetary Sciences, Department of Agronomy, Purdue University, West Lafayette, IN 47907, USA
• 2National Institute for Environmental Studies, 16-2 Onogawa, Tsukuba, Ibaraki, 305-8506, Japan
• 3Department of Agronomy, Purdue University, West Lafayette, IN 47907, USA

Correspondence: Qianlai Zhuang (qzhuang@purdue.edu)

Abstract

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−2 yr−1 (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 (Cv), nitrogen in living vegetation (Nv), organic carbon in detritus and soils (Cs), organic nitrogen in detritus and soils (Ns), and available inorganic soil nitrogen (Nav). Carbon and nitrogen dynamics in TEM are governed by the following equations:

$\begin{array}{}\text{(1)}& & \frac{\mathrm{d}{C}_{\mathrm{v}}}{\mathrm{d}t}=\text{GPP}-{R}_{\mathrm{A}}-{L}_{\mathrm{C}},\text{(2)}& & \frac{\mathrm{d}{N}_{\mathrm{v}}}{\mathrm{d}t}=\text{NUPTAKE}-{L}_{\mathrm{N}},\text{(3)}& & \frac{\mathrm{d}{C}_{\mathrm{s}}}{\mathrm{d}t}={L}_{\mathrm{c}}-{R}_{\mathrm{H}},\text{(4)}& & \frac{\mathrm{d}{N}_{\mathrm{s}}}{\mathrm{d}t}={L}_{\mathrm{N}}-\text{NETNMIN},\text{(5)}& & \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, RA is autotrophic respiration, LCis carbon in litterfall, NUPTAKE is nitrogen uptake by vegetation, LN is nitrogen in litterfall, RH 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}& \text{GPP}={C}_{max}f\left(\text{PAR}\right)f\left(\text{PHENOLOGY}\right)f\left(\text{FOLIAGE}\right)\\ \text{(6)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}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)}& & \text{NPP}=\text{GPP}-{R}_{\mathrm{A}},\text{(8)}& & \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−2 yr−1).

## 2.2 Spin-up acceleration method

TEM can be reformulated as

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

where x is a vector of state variables (e.g., CV); 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)}& g\left(\mathbit{x},t\right)=g\left({\mathbit{x}}_{\mathrm{0}},t\right)+\mathbf{J}\left(\mathbit{x}-{\mathbit{x}}_{\mathrm{0}}\right),\end{array}$

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

$\begin{array}{}\text{(11)}& \mathbf{J}=\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). xn represents each state variable in the TEM (e.g., VC). The numerical discretization of Eq. (9) is

$\begin{array}{ll}& {x}_{i,\phantom{\rule{0.125em}{0ex}}k}-{x}_{i,\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}=\mathit{\tau }\cdot {\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot {x}_{i,\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}\\ \text{(12)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}+\mathit{\tau }\left(g\left({x}_{\mathrm{0}},\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}\right)-{\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot {x}_{\mathrm{0},\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}+{h}_{k-\mathrm{1}}\right),\end{array}$

where τ is time step (month), xi, k is the pool xi size at time k, and Jk−½ 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. xi,0 refers to the initial pool xi size.

We introduce

$\begin{array}{}\text{(13)}& {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},\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}+{h}_{k-\mathrm{1}}.\end{array}$

The Eq. (12) can then be written as

$\begin{array}{}\text{(14)}& {x}_{i,\phantom{\rule{0.125em}{0ex}}k}-{x}_{i,\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}=\mathit{\tau }\cdot {\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot {x}_{i,\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}+\mathit{\tau }\cdot {f}_{k-\mathrm{1}},\end{array}$

where Jk−½ 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 $\stackrel{\mathrm{\to }}{{x}_{\mathrm{1}}}$ is the size of carbon pools on 1 January, while J1.5 is the matrix of mean process rate constants for January.

By introducing

${A}_{k}=\mathit{\tau }\cdot {\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}},\phantom{\rule{1em}{0ex}}{y}_{k}=\mathit{\tau }{f}_{k-\mathrm{1}},\phantom{\rule{1em}{0ex}}{B}_{k}=\mathbf{I},\phantom{\rule{1em}{0ex}}{C}_{k}=\mathbf{I}+A,$

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

$\begin{array}{}\text{(15)}& -{C}_{k}\cdot {x}_{i,\phantom{\rule{0.125em}{0ex}}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)}& -{C}_{\mathrm{1}}\cdot {x}_{i,\phantom{\rule{0.125em}{0ex}}T}+{B}_{\mathrm{1}}\cdot {x}_{i,\phantom{\rule{0.125em}{0ex}}\mathrm{1}}={y}_{\mathrm{1}}.\end{array}$

Thus Eq. (15) and (15a) become a formulation of a linear problem with T unknown vectors $\stackrel{\mathrm{\to }}{{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}& \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)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}×\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)}& {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 1Algorithms 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)}& & \text{NPP}=\text{GPP}-\text{MR}-\text{GR},\text{(19)}& & \text{MR}={V}_{\mathrm{C}}\cdot {K}_{T},\text{(20)}& & \text{GR}=\left\{\begin{array}{ll}\mathrm{0.25}\cdot \left(\text{GPP}-\text{MR}\right),& \text{if}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\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 CV and temperature (KT). Here we revised MR calculation:

$\begin{array}{ll}& \text{MR}\\ \text{(21)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}=\left\{\begin{array}{ll}{V}_{\mathrm{C}}\cdot {K}_{T},& \mathrm{\text{if}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\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 1Test sites for new spin-up algorithms.

Figure 2The time for NEP (g C yr−1 m−2) 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 (RH).

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 3Simulated NEP (g C m−2 yr−1) 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−2 yr−1 (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−2 yr−1. 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 2Spin-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)}& \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., CV and CS). ε is a scalar. A is a pool transfer matrix (in which Aij 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.

Data availability
Data availability.

All data used in this study are available from the authors upon request.

Author contributions
Author contributions.

QZ and SM designed and supervised the research. YQ performed model simulations and data analysis. All authors contributed to the paper writing.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

This study is supported through projects funding to Qianlai Zhuang from the Department of Energy (DESC0008092 and DE-SC0007007) and the NSF Division of Information and Intelligent Systems (NSF-1028291). The supercomputing resource is provided by the Rosen Center for Advanced Computing at Purdue University.

Edited by: Christopher A. Williams
Reviewed by: two anonymous referees

References

Baldocchi, D. D., Xu, L., and Kiang, N.: How plant functional-type, weather, seasonal drought, and soil physical properties alter water and energy fluxes of an oak–grass savanna and an annual grassland, Agr. Forest Meteorol., 123, 13–39, 2004.

Comins, H. N.: Analysis of nutrient-cycling dynamics, for predicting sustainability and CO2-response of nutrient-limited forest ecosystems, Ecol. Model., 99, 51–69, 1997.

Dickinson, R. E., Shaikh, M., Bryant, R., and Graumlich, L.: Interactive canopies for a climate model, J. Climate, 11, 2823–2836, 1998.

Fang, Y., Liu, C., and Leung, L. R.: Accelerating the spin-up of the coupled carbon and nitrogen cycle model in CLM4, Geosci. Model Dev., 8, 781–789, https://doi.org/10.5194/gmd-8-781-2015, 2015.

Gilmanov, T. G., Demment, M. W., Wylie, B. K., Laca, E. A., Akshalov, K., Baldocchi, D. D., and Emmerich, W. E.: Quantification of the CO2 exchange in grassland ecosystems of the world using tower measurements, modeling and remote sensing, in: XX International Grassland Congress: offered papers, 26, p. 587, 2005

Goulden, M. L., Winston, G. C., McMillan, A. M. S., Litvak, M. E., Read, E. L., Rocha, A. V., and Rob Elliot, J.: An eddy covariance mesonet to measure the effect of forest age on land–atmosphere exchange, Glob. Change Biol., 12, 2146–2162, 2006.

Gu, L., Meyers, T., Pallardy, S. G., Hanson, P. J., Yang, B., Heuer, M., and Wullschleger, S. D.: Influences of biomass heat and biochemical energy storages on the land surface fluxes and radiative temperature, J. Geophys. Res.-Atmos., 112, https://doi.org/10.1029/2006JD007425, 2007.

Gu, L., Massman, W. J., Leuning, R., Pallardy, S. G., Meyers, T., Hanson, P. J., and Yang, B.: The fundamental equation of eddy covariance and its application in flux measurements, 152, 135-148, Agr. Forest Meteorol., 2012.

Johns, T. C., Carnell, R. E., Crossley, J. F., Gregory, J. M., Mitchell, J. F., Senior, C. A., Tett, S. F. B., and Wood, R. A.: The second Hadley Centre coupled ocean-atmosphere GCM: model description, spinup and validation, Clim. Dynam., 13, 103–134, 1997.

King, D. A.: Equilibrium analysis of a decomposition and yield model applied to Pinus radiata plantations on sites of contrasting fertility, Ecol. Model., 83, 349–358, 1995.

Kwon, E. Y. and Primeau, F.: Optimization and sensitivity study of a biogeochemistry ocean model using an implicit solver and in situ phosphate data, Global Biogeochem. Cy., 20, GB4009, https://doi.org/10.1029/2005GB002631, 2006.

Lardy, R., Bellocchi, G., and Soussana, J. F.: A new method to determine soil organic carbon equilibrium, Environ. Modell. Softw., 26, 1759–1763, 2011.

Martin, M. P., Cordier, S., Balesdent, J., and Arrouays, D.: Periodic solutions for soil carbon dynamics equilibriums with time-varying forcing variables, Ecol. Model., 204, 523–530, 2007.

McGuire, A. D., Melillo, J. M., Joyce, L. A., Kicklighter, D. W., Grace, A. L., Moore, B. I. I. I., and Vorosmarty, C. J.: Interactions between carbon and nitrogen dynamics in estimating net primary productivity for potential vegetation in North America, Global Biogeochem. Cy., 6, 101–124, 1992.

NVIDIA CUDA Compute Unified Device Architecture Programming Guide, NVIDIA Corporation, 2007.

Ollinger, S. V. and Smith, M. L.: Net primary production and canopy nitrogen in a temperate forest landscape: an analysis using imaging spectroscopy, modeling and field data, Ecosystems, 8, 760–778, 2005.

Thornton, P. E. and Rosenbloom, N. A.: Ecosystem model spin-up: Estimating steady state conditions in a coupled terrestrial carbon and nitrogen cycle model, Ecol. Model., 189, 25–48, 2005.

Turnipseed, A. A., Anderson, D. E., Blanken, P. D., Baugh, W. M., and Monson, R. K.: Airflows and turbulent flux measurements in mountainous terrain: Part 1. Canopy and local effects, Agr. Forest Meteorol., 119, 1–21, 2003.

Turnipseed, A. A., Anderson, D. E., Burns, S., Blanken, P. D., and Monson, R. K.: Airflows and turbulent flux measurements in mountainous terrain: Part 2: Mesoscale effects, Agr. Forest Meteorol., 125, 187–205, 2004.

Van Gorsel, E., Delpierre, N., Leuning, R., Black, A., Munger, J. W., Wofsy, S., and Chen, B.: Estimating nocturnal ecosystem respiration from the vertical turbulent flux and change in storage of CO2, Agr. Forest Meteorol., 149(11), 1919–1930, 2009.

Xia, J. Y., Luo, Y. Q., Wang, Y.-P., Weng, E. S., and Hararuk, O.: A semi-analytical solution to accelerate spin-up of a coupled carbon and nitrogen land model to steady state, Geosci. Model Dev., 5, 1259–1271, https://doi.org/10.5194/gmd-5-1259-2012, 2012.

Zhuang, Q., McGuire, A. D., Melillo, J. M., Clein, J. S., Dargaville, R. J., Kicklighter, D. W., Myneni, R. B., Dong, J., Romanovsky, V. E., Harden, J., and Hobbie, J. E.: Carbon cycling in extratropical terrestrial ecosystems of the Northern Hemisphere during the 20th century: a modeling analysis of the influences of soil thermal dynamics, Tellus B, 55, 751–776, 2003.

Zhuang, Q., He, J., Lu, Y., Ji, L., Xiao, J., and Luo, T.: Carbon dynamics of terrestrial ecosystems on the Tibetan Plateau during the 20th century: an analysis with a process-based biogeochemical model, Global Ecol. Biogeogr., 19, 649–662, 2010.