A semi-analytical solution to accelerate spin-up of a coupled carbon and nitrogen land model to steady state

A semi-analytical solution to accelerate spin-up of a coupled carbon and nitrogen land model to steady state J. Xia, Y. Luo, Y.-P. Wang, E. Weng, and O. Hararuk Department of Botany and Microbiology, University of Oklahoma, OK, USA CSIRO Marine and Atmospheric Research, Centre for Australian Weather and Climate Research, Aspendale, Victoria, Australia Department of Ecology and Evolutionary Biology, Princeton University, NJ, USA


Introduction
Modeling ecosystem biogeochemical cycles is highly dependent on initial values because of long-term persistence of ecosystem state properties.It requires setting up initial values of all state variables (e.g., carbon and nitrogen pool sizes) in any biogeochemical models before scientists can use the models for any analysis.The initial values are either estimated from observations (D'Odorico et al., 2004;Luo and Reynolds, 1999) or assumed to be at steady state.The latter is usually achieved by traditional spin-up methods that perform long model simulations until no trend of change in pool sizes over many periods of the repeated climate forcing, even though the pool sizes vary seasonally and inter-annually within one period of the repeated forcing (Johns et al., 1997;McGuire et al., 1997;Thornton et al., 2002;Yang et al., 1995).Spinning biogeochemical models to steady state is computationally expensive, especially so when simulations are performed with global biogeochemical models of coupled C and N cycles (Thornton and Rosenbloom, 2005).In general, a fully coupled earth system model with Published by Copernicus Publications on behalf of the European Geosciences Union.
biogeochemical cycles needs to be spun up in several sequential steps, and each step can take several thousands of simulation years (Doney et al., 2006).Spin-up of a fully coupled earth system model with a relatively coarse resolution (≈ 3.75 • ) atmospheric model is estimated to take hundreds of real-world days of computation using the present National Center for Atmosphere Research (NCAR) supercomputers (Jochum et al., 2010).As a consequence, spin-up has become a serious constraint on global modeling analysis of biogeochemical cycles.Thornton and Rosenbloom (2005) have explored a few spin-up methods for achieving steady states of a coupled carbon-nitrogen ecosystem model.The spin-up of their model begins with initial values of no soil organic matter (SOM) and very small plant pools (Thornton et al., 2002).During the spin-up, the accumulation rate of SOM strongly depends on the nitrogen addition rate in their model.They periodically increase the mineral nitrogen supply during the early stage of the spin-up to acclerate the spin-up.Each nitrogen addition period is followed by a period with reducing nitrogen input linearly to the normal level (Thornton and Rosenbloom, 2005;Thornton et al., 2002).The efficiency of this punctuated nitrogen addition method is low without an optimal combination of the nitrogen addition and reduction periods.A prior analysis for searching such optimal combinations is needed when this method is applied to a new model.An accelerated decomposition method is based on an assumption of linear scaling between decomposition rates and litter/soil pool sizes.A linearity test between decomposition rates and pool sizes is needed before performing this method for each new model (Thornton and Rosenbloom, 2005).The accelerated decomposition method can distort interactions between carbon and nitrogen cycles in the a global biogeochemical model of terrestrial carbon, nitrogen, and phosphorus (CASACNP) (Y.-P.Wang, unpublished data) and ocean physics for spinning up an ocean model (Bryan and Lewis, 1979;Danabasoglu et al., 1996).Although the above methods have been used in some modeling studies (Randerson et al., 2009), most models still use the traditional spin-up method with long-term iterative simulations to achieve the steady states of variables.
Carbon processes in terrestrial ecosystems can be represented by first-order, linear differential equations (Bolker et al., 1998;Luo and Weng, 2011;Luo et al., 2012).This property renders a possibility to obtain an analytical solution of steady states for terrestrial carbon cycle models (Bolker et al., 1998;Comins, 1994;Govind et al., 2011;King, 1995;Luo et al., 2001).Ludwig et al. (1978), for example, have used a two-step approximation.They first calculated the steadystate pool sizes of the fast variables by holding the slow variables fixed, and then analyzed the slow variables with the fast variables held at corresponding steady-state pool sizes.Most of the analytical solutions (Comins, 1994;Govind et al., 2011;King, 1995;Luo et al., 2001) are obtained by using constant net primary productivity (NPP).However, temporal climate fluctuations and seasonal plant growth make it difficult to get an analytical solution for state variables in the models.Some studies (Lardy et al., 2011;Martin et al., 2007) have attempted to obtain the analytical solution of steadystate pool sizes via managing the climate fluctuations.These methods have used matrix-based analysis but still need to solve several relatively complicated equations.For these reasons, no effective analytical method has been developed to save the spin-up time for global land models.
An analytical solution is still possible to obtain steadystate pool sizes if we can overcome two obstacles.First, we need to get time-averaged approximations of the timevarying variables, such as environmental scalars and NPP.Since most spin-up uses repeated climate forcing variables to estimate steady states of pool sizes, we can estimate averages of those time-varying variables within one period of the repeated forcing variables.Second, for most carbonnitrogen coupled models, nitrogen regulates carbon cycle via its influences on photosynthesis and decomposition.Nitrogen pool sizes are usually related with carbon cycle via carbon/nitrogen (C/N) ratios (Gerber et al., 2010;Wang et al., 2010).Nitrogen influences on photosynthesis are fast processes and can be accounted for by short, initial spin-up to reach a stabilization of NPP.The C/N ratios in the end loop of initial spin-up for stable NPP can be used to estimate steadystate nitrogen pool sizes.Using the above approximations of those time-varying variables will generate errors for estimating steady states of carbon and nitrogen pools.Thus, some additional spin-up may be necessary to achieve steady states of all pools.
This study was intended to develop a semi-analytical solution (SAS) to accelerate spin-up of global carbon-nitrogen coupled models.We first discussed biogeochemical principles underlying SAS.SAS becomes permissible because a set of first-order ordinary differential equations can adequately describe carbon transfers within ecosystems over time and be analytically solved to obtain steady-state pool sizes.We applied SAS to the Australian Community Atmosphere Biosphere Land Exchange (CABLE) model and developed a general procedure of SAS for spin-up.There are three key steps of SAS.The first step is a short spin-up to obtain steady-state NPP, averaged environmental factors within one period of repeated forcing variables, and C/N ratios.The second step is to analytically solve the differential equation to calculate steady-state carbon and nitrogen pool sizes.The last step is to make an additional spin-up to meet the steady-state criterion for all pools.We evaluated the computational efficiency of SAS for the carbon-only and the coupled carbon-nitrogen models against the traditional spin-up, applications of SAS to other biogeochemical models, and possible model analyses enabled by SAS.

Leaf
Canopy Photosynthesis CWD (X 6 ) Fig. 1.Diagram of the carbon processes of the CASACNP model on which Eq. ( 1) is based.SOM stands for soil organic matter.

Biogeochemical principles for the semi-analytical solution (SAS)
The semi-analytical solution (SAS) we introduced in this study is built upon principles of biogeochemical cycles in terrestrial ecosystems.The biogeochemical cycle of carbon in an ecosystem is usually initiated with plant photosynthesis, which fixes CO 2 from the atmosphere into an ecosystem.The photosynthetic carbon is partitioned into leaf, root, and woody biomass.Dead biomass becomes litter to metabolic, structural, and coarse woody debris (CWD) litter pools.The litter carbon is partially released to the atmosphere as respired CO 2 and partially converted to soil organic matter (SOM) in fast, slow, and passive pools (Fig. 1).The mean carbon residence time varies greatly among different pools, from several months in leaves and roots to hundreds or thousands of years in woody tissues and SOM (Torn et al., 1997).
The carbon processes in a terrestrial ecosystem can be mathematically expressed by the following first-order ordinary differential matrix equation (Luo et al., 2003): where X(t) = (X 1 (t), X 2 (t), ..., X 9 (t)) T is a 9×1 vector describing nine carbon pool sizes in leaf, wood, root, metabolic litter, structural litter, CWD, fast SOM, slow SOM, and passive SOM, respectively, in the Community Atmosphere Biosphere Land Exchange (CABLE) model (Wang et al., 2011).
The diag(c) is a 9 × 9 diagonal matrix with diagonal entries given by vector c = (c 1 , c 2 , . . ., c 9 ) T ; components c j (j = 1, 2, . . ., 9) quantify the fraction of carbon left from pool X j (j = 1, 2, . . ., 9) after each time step.ξ is an environmental scalar accounting for effects of soil type, temperature and moisture on carbon decomposition.B = (b 1 , b 2 , b 3 , 0, . . ., 0) T represents the partitioning coefficients of the photosynthetically fixed carbon into different pools.U is the input of fixed carbon via plant photosynthesis.In general, Eq. ( 1) can adequately summarize C cycle processes in most land models (Cramer et al., 2001;Parton et al., 1987).Equation ( 1) cannot be directly solved to obtain the steadystate values of carbon pools because matrices A and B, the environmental scalar ξ , and ecosystem carbon influx U (t) vary with time and driving variables.Since carbon influx involves fast processes, its steady-state value U ss can be obtained from short spin-up.Most model spin-ups use repeated driving variables in long-term, iterative simulations.Thus it is possible to calculate averaged values of the environmental scalar ( ξ ), the carbon transfer ( Ā) and partitioning ( B) coefficients within one period of repeated driving variables.With U ss and the mean values of the time-varying variables ( ξ, Ā, and B), we can analytically calculate the steady-state carbon pool sizes X ss by letting Eq. (1) equal zero as We divided X ss by C/N ratio in each pool to obtain steadystate nitrogen pool sizes N ss .The C/N ratios were the temporal average values of the last loop of the initial spin-up for U ss .Using temporal averages of those time-varying variables ( ξ , Ā, B, and mean C/N ratios) will yield approximation errors to estimate the steady states of carbon and nitrogen pools.Equations (1)-(4) assume the terrestrial carbon cycle as a linear system, while biogeochemical models simulate some carbon processes nonlinearly.For example, allocations of NPP to new leaf and stem growths are determined by a dynamic function of NPP in the Community Land Model (CLM-CN) (Oleson et al., 2010).As a result, using the linear system will generate some additional approximation errors to estimate the steady states of carbon and nitrogen pools.Thus, the steady-state carbon and nitrogen pool sizes that are analytically calculated by Eq. ( 4) need to be further adjusted with additional spin-up to meet the criterion of steady states for all carbon and nitrogen processes.
Overall, our semi-analytic solution of spin-up consists of three steps: (1) an initial spin-up to obtain steady-state carbon influx U ss , temporally averaged values of the time-varying variables in Eq. (3) ( ξ , Ā, and B), and C/N values; (2) calculation of the steady-state carbon pool sizes X ss using Eq. ( 3) and the steady-state N pool sizes N ss from dividing X ss by C/N ratios; and (3) additional spin-up to meet the steadystate criteria for all carbon and nitrogen processes.

Model description
We applied the semi-analytic solution of spin-up to the CA-BLE model, which is one of the land surface models for simulation of biophysical and biochemical processes.Kowalczyk (2006) and Wang et al. (2011) have described the CA-BLE model in detail.The model includes 5 submodels: radiation, canopy micrometeorology, surface flux, soil and snow, and biogeochemical cycles.The CABLE model calls the radiation submodel first to compute absorption and transmission of both diffuse and direct beam radiation in the two-bigleaf canopy and at soil surface (see the details in Wang and Leuning, 1998).The canopy micrometeorology submodel estimates the canopy roughness length, zero-plane displacement height and aerodynamic transfer resistance based on the theory developed by Raupach (Raupach, 1989a(Raupach, , b, 1994) ) and Raupach et al. (1997).The surface flux submodel uses the absorbed radiation to estimate the water extraction and ground heat flux, which are required in the soil and snow submodel.The biogeochemical-cycles submodel is called last to compute the respiration of non-leaf plant tissues, the soil respiration, and the net ecosystem CO 2 exchange.
The biogeochemical submodel of CABLE evolves from the CASACNP model, which was developed by Wang et al. (2010).It adopted the model structure of carbon processes from the CASA' model (Randerson et al., 1997) and contains coupled carbon, nitrogen, and phosphorus cycles.In this study, phosphorus cycle and its coupling with carbon and nitrogen cycles were not activated.The CASACNP model has 9 pools, which include three plant pools (leaf, wood, and root), three litter pools (metabolic litter, structural litter, and coarse woody debris), and three soil pools (microbial biomass, slow and passive soil organic matter) (Fig. 1).There is an additional pool of inorganic nitrogen (NO − 3 + NH + 4 ) when nitrogen cycle is coupled with carbon cycle.The equations that describe changes in pool size with time have been presented by Wang et al. (2010) and can be represented by Eq. ( 1).In Eq. ( 1), parameter C is set to be a constant in the CABLE model, while matrices A and B, the environmental scalar ξ , and ecosystem carbon influx U (t) vary with time and driving variables.In the A matrix, the carbon transfer coefficients are determined by lignin/nitrogen ratio from plant to litter pools, lignin fraction from litter to soil pools, and soil texture among soil pools.In the B matrix, the carbon partitioning coefficients of the photosynthetically fixed carbon into plant pools are determined by availabilities of light, water and nitrogen as the carbon allocation scheme described by Friedlingstein et al. (1999).The environmental scalar (ξ ) regulates the leaf turnover rates by cold and drought stresses on leaf senescence rate, the turnover rates of litter carbon pools via limitations of soil temperature, moisture, and nitrogen availability, and SOM turnover rates by soil temperature, moisture, and texture.The soil texture is spatially fixed in the CABLE model.The soil nitrogen will limit litter decomposition if the gross mineralization is less than the immobilization (Wang et al., 2010).We spun up the model for about one hundred of simulation years to obtain the steady-state carbon influx U ss .
In the CABLE model, the optimal carbon decay rates (the C matrix in Eq. 1) are preset and vary with vegetation types.The vegetation types for each 1 • × 1 • grid cell in the model were derived from the 0.5 • × 0.5 • International Geosphere-Biosphere Programme (IGBP) classification (Loveland et al., 2000).During the spin-up of the coupled carbon-nitrogen model, the carbon influx (U ) and litter decomposition rate are regulated by the soil nitrogen availability (Wang et al., 2010).The nitrogen regulation may periodically occur until all the nitrogen processes of the model reach steady states.In the CABLE model, the nitrogen inputs of deposition, fertilizer application and fixation are explicitly estimated.The nitrogen deposition in 1990 was estimated from Dentener (2006) and nitrogen fixation from Wang and Houlton (2009).The global fertilizer application of nitrogen is taken as 0.86 Gt N yr −1 from Galloway et al. ( 2004) and is distributed uniformly within the cropland biome (Wang et al., 2010).
The forcing variables required for the CABLE model include incoming short-and long-wave radiation, air temperature, specific humidity, air pressure, wind speed, precipitation and ambient CO 2 concentration.The CABLE model first generated daily meteorological forcing (surface air temperature, soil temperature and moisture).Then the daily forcings were used to integrate the full model with a time step of one day.In this study, the meteorological forcings of 1990 were used to run the global version of the CABLE model to steady states.A site version of CABLE (v2.01), which has been calibrated by datasets from Harvard Forest, was used for the site-level analysis in this study.We used forcing data of Harvard Forest from 1992 to 1999 with the time step of half an hour for the site-level simulation.A detailed description of the data sources was provided by Urbanski et al. (2007).

The procedure of semi-analytic solution to spin-up
For modeling analyses of biogeochemical responses to global change, models often have to be first spun up to steady

Initial spin-up.
• Read in the initial parameters and spin-up NPP (or plant carbon pools) to steady state • Meanwhile, save all the values of the variables in the equations in step 4 7. Final spin-up.
• Read in the analytical solved carbon and nitrogen pools, and spin up all pools to steady states state for all pools and fluxes.Traditionally, the biogeochemical models first read in all meteorological input parameter values and initial pool sizes.Then the models continuously run with recycled meteorological forcing variables for thousands of simulation years until steady states are reached for all pools and fluxes.
To implement SAS with the CABLE model, we did the following things (Fig. (3) creating equations to calculate the analytical solutions of each pool according to the structures of matrix A, C and vector B; and (4) setting up a criterion for the steady state of the slowest pool for the final spin-up.More details about the recoding can be found in Text S1.
5. Making an initial spin-up by running the model using repeated meteorological forcing until NPP (or all plant pools) reached stabilization (U ss ).In this study, we ran the model until the mean change in NPP over each loop (8 yr) of site simulation at Harvard Forest was smaller than 10 −4 g C m −2 .For the global simulation, we ran the model until the mean changes in plant carbon pools over each loop (1 yr) were smaller than 0.01 % per year compared to the previous cycle.Meanwhile, the mean values of all the time-varying parameters in Eq. ( 3) were written to those newly created variables.Those parameters are the stable NPP (U ss ), the mean environmental scalar ( ξ ) and matrices of carbon transfer ( Ā) and partitioning ( B) coefficients within one period of repeated forcing variables, as well as C/N ratios at the end of the initial spin-up.
6. Calculating the analytical solution of the steady states of carbon and nitrogen pools.The steady-state carbon pools were solved by letting carbon influx equal efflux for each pool (Eq.4).Nitrogen pools are obtained by dividing the steady-state carbon pools by the mean C/N ratios of the end loop of the initial spin-up.
7. Making the final spin-up by using the analytically solved carbon and nitrogen pools as initial values until the steady-state criterion for the soil carbon pools was met.The steady-state criterion set in this study was that the change in any soil carbon pool ( C soil ) within each simulation cycle was smaller than 0.5 g C m −2 yr −1 (as one criterion in Thornton and Rosenbloom, 2005).According to the difference in turnover rate, a slower pool needs a longer time to reach steady state during the spinup.The final spin-up is determined by the dynamic of the slowest carbon pool when the criterion of steadystate soil carbon pools is small enough.

Performances of SAS at Harvard Forest
For Harvard Forest, the traditional spin-up method took 9768 and 6768 yr (1220 and 846 loops, respectively) to get the steady states of the carbon-only and coupled carbon-nitrogen simulations ( C soil < 0.5 g C m −2 yr −1 ), respectively.The first step of the SAS was to spin up the model to reach stable NPP.It took 64 and 336 yr (8 and 42 loops, respectively) for the carbon-only and coupled carbon-nitrogen simulations, respectively (Fig. 3).After the semi-analytical solution of steady-state values was obtained for all carbon and nitrogen pools, it took another 45 and 89 loops of the carbon-only and coupled carbon-nitrogen simulations, respectively, for the change in any soil carbon pool in each loop of simulation ( C soil ) less than 0.5 g C m −2 yr −1 .In comparison with the traditional spin-up method, the SAS method saved about 95.7 % and 84.5 % of computational time for getting steady states of the carbon-only and coupled carbon-nitrogen simulations, respectively.The differences in steady-state carbon pools between the SAS and traditional spin-up methods were small, being 0.85 % and 0.24 % of total ecosystem carbon content for carbon-only and coupled carbon-nitrogen simulations, respectively (Table 1).The stable NPP in the coupled carbon-nitrogen simulation (1.45 g C m −2 d −1 ) was 32.24 % less than that in the carbononly simulation (2.14 g C m −2 d −1 ).The steady-state value of the total ecosystem carbon pool, which represents the ecosystem carbon storage capacity, in the coupled carbonnitrogen simulation (28.83 kg C m −2 ) was 31.44 % lower than that in the carbon-only simulation (42.05 kg C m −2 ; Table 1).Although the passive SOM pool determined the spinup time of CABLE, the slow SOM pool had the largest pool size (Figs. 4 and 5) because the steady-state pool size (i.e., carbon storage capacity) is jointly determined by carbon influx and residence time.

Application of SAS to global simulations
The traditional spin-up method spent 2780 and 5099 simulation years for carbon-only and coupled carbon-nitrogen simulation, respectively, before the change in the slowest carbon pool met the steady-state criterion ( C soil < 0.5 g C m −2 yr −1 ; Fig. 6).For SAS, the initial spin-up took 200 simulation years for obtaining steady states of plant carbon pools in the global carbon-only model and 201 yr for the coupled carbon-nitrogen model.With the SAS spinup method, all carbon pools in the carbon-only model reached steady states ( C soil < 0.5 g C m −2 yr −1 ) after analytical calculation without any final spin-up (as shown in the black arrow in Fig. 6).In the coupled carbon-nitrogen model, the SAS needed another 483 simulation years to obtain the steady states of all pools (as shown in the gray arrow in Fig. 6) after the analytical calculation.The SAS method saves about 92.4 % and 86.6 % of the computational time for spin-up of the global carbon-only and coupled carbonnitrogen models to steady states, respectively (Fig. 6).
With the traditional spin-up method, the SOM pool continued to decrease after it reached the steady-state criterion ( C soil < 0.5 g C m −2 yr −1 ) (Fig. 6).Additional thousands  of simulation years were needed for the traditional method to reach steady states of all SOM pools, which were analytically obtained by the SAS with no time (Fig. 6).

Steady-state pools and fluxes as regulated by N
The capacity of an ecosystem to store carbon is determined by ecosystem carbon influx and residence times of different pools (Luo et al., 2003).The global mean steady-state NPP was greater in the carbon-only (0.37 kg C m −2 yr −1 ) than the coupled carbon-nitrogen (0.35 kg C m −2 yr −1 ) simulation (Fig. 7a and b).A larger proportion of photosynthetically fixed carbon was partitioned to pools with long residence time (e.g., plant wood, slow and passive SOM; Figure 7).The global mean of the ecosystem carbon pool sizes at steady state decreased from 14.1 in the carbononly model to 12.1 kg C m −2 in the coupled carbon-nitrogen model (Fig. 7).The mean residence time (as dividing steadystate ecosystem carbon pool size by NPP) of the ecosystem carbon pool at steady state in the coupled carbon-nitrogen model (34.0 yr) was 10.5 % shorter than that in the carbononly model (38.0 yr).In the CABLE model, large fractions of photosynthate went to plant pools, the CWD pool, and the slow SOM pool, but a very small fraction (∼ 0.1 %) to the passive SOM pool (Fig. 7).Globally, nitrogen processes down-regulated soil carbon storage much more substantially at the high (e.g., temperate conifer and mixed forests) than the low latitudes (e.g., arid and semi-arid deserts and tropical forest) (Fig. 8).

Computational efficiency
The SAS method saved 92.4 % of computational time for spin-up of the global carbon-only model and 86.6 % for the global coupled carbon-nitrogen model in comparison with the traditional spin-up method.At the site level in Harvard Forest, SAS saved 95.7 % and 84.5 % of computational time for the carbon-only model and coupled carbon-nitrogen model, respectively.That means the spin-up with the SAS method can be up to 20 times as fast as the traditional method.The computational efficiency with the SAS method is higher than the best method (the accelerated decomposition method) explored by Thornton and Rosenbloom (2005) for site-level spin-up in an evergreen needle-leaf forest.Lardy et al. (2011) have recently developed an iterative matrix method to accelerate spin-up with the Pasture Simulation Model (PaSim).They reported that their method speeds up the spin-up by up to 20 times as well.
The SAS method can be easily implemented for biogeochemical models at site, regional, and global scales.As described in the Methods section, implementation of the SAS     method involves some light recoding of original models to enable analytical calculation of steady states of carbon and nitrogen pools together with initial and final spin-ups (Text S1).The accelerated decomposition method examined by Thornton and Rosenbloom (2005) has been found difficult to be applied to an age-structured model (Lardy et al., 2011) and requires re-parameterization with a linear scaling factor for each new model.The SAS method described in this study is easily programmed into an existing model, for example in less than 150 lines for the CABLE model.The computational cost for spinning up models strongly depends on the criterion used for steady state.The more precise the criterion (i.e., the smaller value), the longer time the spin-up needs for a model to reach the steady state.With the criterion in this study ( C soil < 0.5 g C m −2 yr −1 ), the CA-BLE model was spun up for thousands of simulation years with the traditional method.Additional thousands of simulation years were needed beyond the traditional spin-up to reach a steady state of any SOM pool, which was analytically solved by the SAS (Fig. 6).This suggests the SAS method is more efficient to estimate the steady states with high precision than the traditional spin-up method.

Applications of SAS to various types of biogeochemical models
The developed SAS is generally applicable to most of the terrestrial biogeochemical models that share a similar structure with the CABLE model in this study (Fig. 1).The CENTURY model (Parton et al., 1993), for example, has 4 plant pools for grassland/crop (shoot, root, grain, standing dead) and 8 for forests (leaf, fine roots, fine branches, dead branches, large wood, dead large wood, coarse roots, dead coarse roots).The dead plant materials go into 4 litter pools (surface structure and metabolic litter, soil structure and metabolic litter) and then 4 soil organic matter (SOM) pools (surface microbes, soil microbes, slow SOM, and passive SOM).Similar to the CABLE model, the decomposition rate of each carbon pool (C matrix in Eq. 1) in the CEN-TURY model is preset with an optimal value and modified by soil texture, soil temperature and moisture.The transfer coefficients (A matrix in Eq. 1) in the CENTURY model are determined by lignin to nitrogen ratio from plant to litter pools, lignin fraction from litter to soil pools, and soil texture from soil to soil pools (Bolker et al., 1998;Parton et al., 1987).
If SAS is applied to CENTURY, a flow diagram as in Fig. 1 would be needed separately for grassland/crop and forest systems due to different numbers of plant pools.The rest of the SAS procedure as described in Fig. 2 can be exactly applied to CENTURY for spin-up.The Rothamsted carbon (RothC) model has four active SOM pools (decomposable plant material, resistant plant material, microbial biomass, and humified organic matter) and one inert organic matter pool (Coleman and Jenkinson, 1999).Although the biogeochemical models differ on the mechanisms of ecological feedbacks (Hurtt et al., 1998), carbon transfers among these models all follow first-order decay equations as in Eq. ( 1).Thus, the SAS procedure can be applied to spin-up of the RachC model.An application of SAS to spin up other coupled carbonnitrogen models may require additional steps, as different models couple the nitrogen cycle to the carbon cycle in different ways.For example, the Princeton Geophysical Fluid Dynamics Laboratory (GFDL) LM3V model uses one arbitrary nitrogen storage (or buffering) pool in plants to avoid short-term switches between N sufficiency and limitation in plants (Gerber et al., 2010).The optimum size of the nitrogen buffering pool equals the total nitrogen losses from living pools (leaf, root, and sapwood) over one year.The rest of the model structure of LM3V for carbon and nitrogen transfers within the ecosystem is similar to CABLE.Thus, the SAS can be applied to the LM3V model with one additional step to calculate the steady state of the optimal nitrogen buffering pool from the analytically solved steady-state live biomass carbon pools.Similarly, the CLM-CN and O-CN models have nitrogen buffering pools (Thornton and Zimmermann, 2007;Zaehle et al., 2010).Therefore, it is necessary to determine how each model couples the nitrogen cycle with the carbon cycle before we apply the SAS procedure to spin-up of the coupled carbon-nitrogen models.
One assumption of the SAS method is the linearity of terrestrial carbon cycle (Eqs.1-4).However, most biogeochemical models simulate some processes in nonlinear ways.For example, allocation of NPP to different plant pools is a nonlinear function of NPP in the CLM-CN model (Oleson et al., 2010).The more nonlinear processes a model includes, the larger the approximation errors that will be generated by the analytical solution, and a longer time is needed for the final spin-up step to adjust all variables to steady states.Another assumption of the SAS method is NPP stabilizes faster than soil carbon pool size.However, some models have very complex vegetation submodels and their NPP cannot stabilize quickly.For example, the CLM-CN model has 20 and 19 pools for vegetation carbon and nitrogen, respectively (Oleson et al., 2010).Thus, the SAS method cannot save spin-up time of these models as much as the CABLE model in this study.NPP of some other models, e.g., the PnET model (Aber and Federer, 1992), is a function of plant nitrogen, and therefore NPP will be different between the initial and final spin-up steps.For these models, an iteration of the analytical solution at the end of each recycle of meteorological forcing would be needed.That means in the fifth step of the SAS method (Fig. 1), the values of all variables in Eq. ( 4) will be saved at the end of each recycling of meteorological forcing instead of after NPP reaches stabilization.Then the steps 5-7 would be iterated until all pools reach steady states.Such an iterative procedure has been successfully applied to the Pasim model (Lardy et al., 2011).
In principle, The SAS method can be applicable to spin up ocean biogeochemical models.Most of the ocean biogeochemical models use the traditional method for spin-up for more than 10 000 yr to reach steady state (Schmittner et al., 2008).A key issue of applying the SAS method to spin up ocean biogeochemical models is whether or not the ocean carbon and nitrogen cycles can be mathematically described by a matrix form similarly as in Eq. ( 1).If yes, this SAS method can facilitate high resolution analysis of or ensemble simulations with ocean models in the future.

SAS-facilitated model analyses
Land models have been developed by the modeling community in the past two decades to predict future states of ecosystems and climate.Model intercomparison has recently become a popular method to improve our understanding of the land model performances (Friedlingstein et al., 2006;Johns et al., 2011).A common protocol for many of the model intercomparison projects is to spin up all the involved models to steady states before global change scenarios are applied to project future changes (Schwalm et al., 2010).Presently, individual models use their own methods for spin-up, often with different criteria of steady states.Standardized spinup with a fast, easily implemented method can help reduce www.geosci-model-dev.net/5/1259/2012/uncertainties in model-model and model-data intercomparison studies.The SAS method has the potential to serve those projects in such a capacity.
Accelerated spin-up reduces computational overhead for modeling analyses and makes some computationally costly analyses feasible.For example, model parameters usually represent the average physiological properties of plant functional types or mean soil attributes.Most of these parameters in the model are assigned values based on relatively few field and/or laboratory observations (Stockli et al., 2008).More and more databases have been developed to indicate that key plant physiological properties, such as leaf traits (GLOP-NET; Reich et al., 2007), carboxylation capacity (Vcmax; Kattge et al., 2009), and biomass allocation (Poorter et al., 2012), greatly vary among plants of different species at different locations.Similarly, properties of soil processes, such as temperature sensitivity of soil respiration (Peng et al., 2009), also greatly vary over time and space.The natural variations in key plant and soil properties can be adequately represented only by probability distributions of parameters, which would propagate in the land models to generate uncertainties in model projections (Weng and Luo, 2011;Xu et al., 2006).The model projection uncertainties can be quantified through ensemble analysis.However, such ensemble analysis of land models against parameter variations is computationally not feasible, because each ensemble element requires spin-up at least once up to a thousand and even million times for one ensemble analysis.Without the ensemble analysis against parameter variations at regional or global scales, uncertainties in model projections cannot be fully assessed.Fast spin-up methods, including SAS, could reduce computational cost and enable the ensemble analysis that is impossible with traditional methods.
The SAS method not only accelerates spin-up with high computational efficiency but also possibly offers an analytical framework for comparative study of structural differences in modeled ecosystem carbon storage capacity.The sum of steady-state carbon pool sizes within one ecosystem obtained from an analytical solution of Eq. ( 4) represents the ecosystem carbon storage capacity (Luo et al., 2003).At the site of Harvard Forest, the carbon storage capacity was simulated to be 42.05 kg C m −2 with the carbon-only model.The capacity was down-regulated by nitrogen processes to be 28.83 kg C m −2 in the coupled carbon-nitrogen model.As indicated by Eq. ( 4), the carbon storage capacity is determined by carbon influx (U ss ) and residence time [( ξ ĀC) Indeed, the steady-state solution of Eq. ( 4) can be used to analyze the determinants of ecosystem carbon storage capacity.Such a linear analytical solution has been successfully used to estimate the steady-state soil organic carbon pools of the CENTURY model (Bolker et al., 1998).The steady-state carbon pool in the biogeochemical model is determined by the carbon influx and residence times (reverses of turnover times).The temporal and spatial variations of the residence times are controlled by environmental scalars to describe effects of temperature, moisture, and soil types on carbon transfer processes.The environmental scalars usually are substantially different among models.For example, the temperature scalar on carbon decay rates is a generalized Poisson function in the CENTURY model (Parton et al., 1994) while it is a simple exponential equation in the Terrestrial Ecosystem Model (TEM) (McGuire et al., 1997).Analyses of carbon cycles in an analytical framework, such as Eq. ( 4), can help compare structural differences in determining ecosystem carbon storage capacity among biogeochemical models.

Conclusions
We developed a new method -the semi-analytical solution (SAS) -to accelerate the spin-up of a process-oriented biogeochemical model to steady states.The SAS described in this study mainly contains 3 steps: (1) making an initial spinup to get steady-state values of photosynthetic carbon input and plant pools; (2) calculating the semi-analytical solution of the steady-state pool sizes; (3) having a final spin-up to meet the criterion of steady states for slowest pools (Fig. 2).For spin-up of the carbon-only model, the SAS method can save about 95.7 % and 92.4 % computational time for sitelevel and global simulations, respectively.The efficiency of the SAS method decreased for coupled carbon-nitrogen simulations, but it still resulted in 84.5 % and 86.6 % reduction in computational cost for site-level and global simulation, respectively.For those models with complex vegetation dynamics, iterations of the SAS method would be needed.We suggest that the SAS described in this study would be a candidate for solving the "spin-up problem" in global models.This method enables many modeling analyses, such as ensemble analysis with regard to parameter variability, which otherwise are impossible with computationally costly spinup methods.Edited by: D. Lawrence

Supplementary material related to this article is available online at
32 Figure 1.

Fig. 2 .
Fig.2.The spin-up strategies of the spin-up method with the semianalytical solution (SAS) used in this study.
2): 1. Developing a flow diagram to link carbon pools and fluxes within ecosystems as in Fig. 1. 2. Organizing the linkage between carbon pools and fluxes into carbon transfer matrices A and C, and plant carbon partitioning coefficients into vector B. The preset values of the optimal carbon decay rates were organized into the C matrix. 3. Figuring out how each element of the time-varying variables (A, B and ξ ) in the Eq.(1) was determined in the model.4. Recoding a section in the model (e.g., the biogeochemical cycle submodel in the CABLE).The recoding of the model includes 4 steps: (1) setting up a criterion for the stable NPP for the initial spin-up; (2) creating new variables to store the mean values of the time-varying parameters;

Fig. 7 .
Fig. 7. Structure of CASACNP model for (a) carbon-only and (b) coupled carbon-nitrogen simulations on which Eq. (1) is based.The fraction of carbon that flows through differential pathways (the numbers near the arrows) is partitioned to the 9 pools.The numbers in the boxes are steady states of ecosystem carbon influx (kg C m −2 yr −1 ) and pool size (kg C m −2 ).The fraction to plant pools is determined by partitioning coefficients in the vector B in Eq. (1).The fractions to litter and soil pools are determined by the transfer coefficient matrix A in Eq. (1).The values of B and A are the global mean values at steady states.SOM stands for soil organic matter.Upper and bottom values near the arrows represent fractions from structural litter and coarse wood detritus (CWD), respectively.

Fig. 8 .
Fig. 8. Global distributions of soil carbon density (kg C m −2 ) at steady state simulated by the carbon-only (upper panel) or coupled carbon-nitrogen (bottom panel) models.
The 45.86 % decrease in carbon storage largely resulted from the NPP decrease in the coupled carbon-nitrogen model NPP (1.45 g C m −2 d −1 ) in comparison with that in the carbon-only model (2.14 g C m −2 d −1 ).
0743778, DEB 0840964, and EF 1137293.Parts of the model runs were performed at the Supercomputing Center for Education & Research (OSCER), University of Oklahoma.

Table 1 .
Mean steady-state values (kg C m −2 ) of all pools and their total value (C tot ) from spin-up with traditional and the SAS methods, and the corresponding relative differences ( C; %) for multiple carbon pools.M-litter, metabolic litter pool; S-litter, structural litter pool; FSOM, fast SOM; SSOM, slow SOM; PSOM, passive SOM.