Accelerating the Spin-up of the Coupled Carbon and Nitrogen Cycle Model in Clm4

The commonly adopted biogeochemistry spin-up process in an Earth system model (ESM) is to run the model for hundreds to thousands of years subject to periodic atmospheric forcing to reach dynamic steady state of the carbon– nitrogen (CN) models. A variety of approaches have been proposed to reduce the computation time of the spin-up process. Significant improvement in computational efficiency has been made recently. However, a long simulation time is still required to reach the common convergence criteria of the coupled carbon–nitrogen model. A gradient projection method was proposed and used to further reduce the computation time after examining the trend of the dominant carbon pools. The Community Land Model version 4 (CLM4) with a carbon and nitrogen component was used in this study. From point-scale simulations, we found that the method can reduce the computation time by 20–69 % compared to one of the fastest approaches in the literature. We also found that the cyclic stability of total carbon for some cases differs from that of the periodic atmospheric forcing, and some cases even showed instability. Close examination showed that one case has a carbon periodicity much longer than that of the atmospheric forcing due to the annual fire disturbance that is longer than half a year. The rest was caused by the instability of water table calculation in the hydrology model of CLM4. The instability issue is resolved after we replaced the hydrology scheme in CLM4 with a flow model for variably saturated porous media.


Introduction
The initial starting values of carbon-nitrogen (CN) models are not commonly available, especially for large-scale applications, but they have an important influence on the subsequent C-N states simulated by the models.Typically, Earth system model (ESM) simulations are initialized in the preindustrial period to allow sufficient time for the coupled system to respond to the various forcings.Initialization of the CN model is usually achieved by a spin-up run of the CN model given an arbitrary initial condition until an approximate C equilibrium is reached.This time marching of the model requires several hundreds to thousands of years of model simulations before a dynamic steady state is reached.The length of the transient state to reach a dynamic steady state is dependent on the initial conditions of the system.It has long been recognized that the spin-up process of CN models is time consuming due to the slow turnover rates of the soil carbon pools, which significantly affect the computational efficiency for global modeling.A number of approaches have been proposed in the past to improve upon the explicit forward time integration of ordinary differential equations in their native form and rate parameters for CN models.These approaches include the initialization of soil organic matter carbon pools with observations (Zhang et al., 2002), the accelerated decomposition method using a higher decomposition rate for litter and soil carbon pools (Thornton and Rosenbloom, 2005), decelerated bulk denitrification and leaching method (Shi et al., 2013) (Xia et al., 2012).Except for the semi-analytical approach, the other approaches mentioned above have been summarized and compared in Shi et al. (2013).The semi-analytical model needs initial spin-up values of net primary productivity (NPP), which still requires a long simulation time for stabilization, because C and N are coupled in CN models.We had previously restructured the CN model in Community Land Model version 4 (CLM4-CN) (Lawrence et al., 2011) and developed a steady-state solution directly using annually averaged rate parameters (Fang et al., 2013(Fang et al., , 2014)).Using our approach, we were able to implement the semianalytical method in Xia et al. (2012).Our numerical experiment showed that the semi-analytical method is not necessarily faster compared to the modified form of the "accelerated decomposition" approach in Koven et al. (2013).
Recently, Koven et al. (2013) used a modified form of the "accelerated decomposition" (hereafter referred to as the AD approach) by numerically increasing the decomposition rates for the two slowest soil carbon pools (named Soil3 C and Soil4 C) to a level so that their turnover rates are similar to the fast pools during the initialization.Numerical evaluation found that the approach significantly reduced the spinup time (Koven et al., 2013).Figure 1 shows the structure of the soil C pool represented in CLM4-CN.Note that heterotrophic respiration fractions are not shown.The reason that the AD approach can accelerate the spin-up is because these two slowest pools are essentially decoupled from the rest of the ordinary differential equations, in that all other pools do not need input from them.The approach, however, cannot be applied to the coarse woody debris (CWD) pool even though its turnover rate is on the same order of Soil3 C, because it is an input to the litter pools.Changing the rate of CWD will give a different solution of other pools during each integration step using the same initial condition, which will lead to a state far from equilibrium if the state is used in a restart simulation.
In the AD approach, once the solution is obtained from the accelerated run, the state of Soil3 C and Soil4 C can be analytically solved.From Fig. 1, the flux of Soil3 C and Soil4 C pools can be described by the following equations: where k L3 , k S2 , k S3 , and k S4 are the turnover rates of the Lit-ter3, Soil2, Soil3, and Soil4 C pools shown in Fig.
Equations ( 3) and ( 4) are applicable regardless of whether AD or a native run was used (the native run was defined here as the simulations without changing the decomposition rates of Soil3 and Soil4 C pools).Therefore, multiplying Eqs. ( 3) and ( 4) by their corresponding accelerator, the results should be close to the native runs.That is, where N denotes the native run, and A is the accelerator.
Even with this modified accelerator approach, long simulation times cannot be avoided in dry and cold places because the decomposition scaling factor is associated with soil water potential and temperature.Hence, new methods are needed to address the spin-up problem.In implicit time integration approaches, based on knowledge about the trajectory of the solution of the initial value problem, linear extrapolation from time integration was often used to find a good initial value for iterative multirate multidisciplinary processes (Birken et al., 2014, and references therein).A number of explicit Euler steps with small time steps followed by an explicit Euler step with a large time step when the change in components due to fast processes become negligible have been shown to efficiently solve stiff ordinary differential equations (Eriksson et al., 2003;Gear and Kevrekidis, 2003).We made use of these concepts, referred to as the gradient projection (GP) approach in this study, to further improve the computational efficiency of the biogeochemistry spin-up processes.

Model description
Community land model CLM4 is the land component of the Community Earth System Model (CESM) (Lawrence et al., 2011).Processes simulated in CLM4 include biogeophysics (solar and longwave radiation, momentum, heat transfer in soil and snow, hydrology of canopy, soil, and snow, and stomatal physiology and photosynthesis) and biogeochemistry (phenology, autotrophic respiration, heterotrophic respiration, carbon and nitrogen allocation, and nitrogen source/sink).The vegetation structures (leaf area index, stem area index and height) in CLM4-CN are represented through the predictive state variables of leaf and stem carbon, which are coupled to simulate fluxes of carbon and nitrogen state variables in vegetation, litter, and soil organic matter (Lawrence et al., 2011;Thornton and Zimmermann, 2007).The tree, shrub and grass plant functional types (PFTs) are divided into tropical, temperate and boreal climate groupings using the PFT physiology and climate rules of Nemani and Running (1996) and C3/C4 photosynthetic pathways in the case of grasses (Lawrence and Chase, 2007).For this study, we used CLM4-CN in offline mode, which is not coupled to an atmosphere model.

Gradient projection method
If m c is the number of years (one cycle) of atmospheric forcing that will be used repeatedly in the spin-up run, we use a spin-up time of [(n + 1)m c ] years as a stop point for the accelerated decomposition (AD) run, where n = 300/m c is an integer.For example, if the number of years of forcing is m c = 7, the stop time will be at year 301.A stop point of ∼ 300 years for the modified AD approach was selected based on the model results in Koven et al. (2013), but it is not an absolute requirement.The best approach is to stop when NPP reaches a dynamic steady state.
At the end of the accelerated run, a dynamic steady-state water table should be reached in the soil column.Due to the slow turnover rates, the total soil carbon gradually approaches steady state from one cycle to the next (Fig. 2a).We can approximate C at a future time t n (Fig. 2b) using the C gradient between two consecutive cycles expressed in the following equation: where t 0 is the beginning of the first cycle, t 1 is the beginning of the next cycle, and t 1 − t 0 = m c ; t n − t 1 = τ m c , τ is an integer close to the turnover years (reciprocal of turnover rate) of the Soil4 C pool to satisfy the stability requirement of forward or explicit time integration that is used in CLM4-CN to solve the time-dependent ordinary differential equations.The explicit method can be numerically unstable (convergence of solution is not guaranteed) if the time step is too big (LeVeque, 2007).For the first-order kinetic type problem, i.e., u (t) = ku(t), the stability requirement is |1 + kh| ≤ 1, in which k is the rate constant and h is the time step.
We call the method shown in Eq. ( 7) the gradient projection (GP) method.This method is analogous to that described in Eriksson et al. (2003), which uses a large time step that satisfies the stability requirement for integrating the slowest processes once the contributions from fast processes become negligible.We allow j p to be chosen based on the time period needed to stabilize the components from fast processes between cycles after perturbation, or set as an integer equaling m c × (100/m c + 1) years of simulation after restart from the accelerated run before using this approach, and also to perform j p years of simulation followed by each projection until the solution meets the common convergence criteria of 0.5 g m −2 for total soil C during two consecutive cycles (Shi et al., 2013;Thornton and Rosenbloom, 2005).During each projection, the balance check for C and N is turned off.The GP method is only applied to the dominant C and N pools, i.e., coarse wood debris, dead stems, dead coarse roots and the Soil4 pool.

Results
A total of 38 single point tower sites from the FLUXNET (Baldocchi et al., 2001) were selected to assess the gradient projection method.These sites include temperate, boreal, tropical, and subtropical climatic environments and four ecosystem types (tropical forests, temperate forests, boreal forests, grasslands, and Mediterranean-type ecosystems) (Table 1).
The meteorological forcing, site information such as soil texture, vegetation cover, and satellite-derived phenology at each site are provided by the North American Carbon Program (NACP) site synthesis team for the sites located in North America and by the Large Scale Biosphere-Atmosphere Experiment in Amazônia Model Intercomparison Project (LBA-MIP) for the sites located in South America.The NACP site synthesis and LBA-MIP data sets are detailed in Schwalm et al. (2010) and at http://www.climatemodeling.org/lba-mip.Each site has two runs, one using the AD method and the other using the GP method.The a Approximate height of the wind/temperature and flux measurements above the surface.b Abbreviated PFTs are BDBorl -broadleaf deciduous boreal tree; BDTmp -broadleaf deciduous temperate tree; BDTrop -broadleaf deciduous tropical tree; BEShr -broadleaf evergreen shrub; BETrop -broadleaf evergreen tropical tree; crop -C3 crop; C3NAGrs -C3 non-arctic grass; C4Grs -C4 grass; NEBorl -needleleaf evergreen boreal tree; and NETmp -needleleaf evergreen temperate tree.c The site information and meteorological forcing are from the LBA-MIP data set.d mc is the number of years of atmospheric forcing.
available forcing (Table 1) is applied repeatedly during the simulation for each site.
Table 2 shows the comparison of total simulation years till a certain convergence criterion is met.Three convergence threshold values in C TOC (3.0, 1.0, and 0.5 g m −2 yr −1 ) were compared.The quality of total soil C is better when the threshold value is smaller (Thornton and Rosenbloom, 2005).Compared to the modified AD approach, the reduction in computation cost is shown in Fig. 3. Figure 3 shows that when a high-quality solution ( C TOC ≤ 0.5 g m −2 yr −1 ) is required, the average total reduction in computation cost is 40 %.On average, 23 % of computation time is reduced in achieving the low-quality solution ( C TOC ≤ 3 g m −2 yr −1 ).
Note that the computation cost reduction for sites US-Me2, RJA, US-IB1 and US-SO2 is not shown in Fig. 3. Site US-Me2 met the convergence criteria before the GP method is applied.Sites RJA and US-IB1 show oscillation of the annual average total C from one full length (multiple years) of forcing cycle to the next, and site US-SO2 shows a carbon periodicity much longer (81 years) than that of the atmo-Table 2. Comparison between the gradient projection method (in bold) and the accelerated spin-up method.

ID Site
Number of simulation years to reach spheric forcing (9 years) (Fig. 4).The oscillation noted in the simulations at RJA and US-IB1 differs from the variability within the forcing cycle, which happens when soil C has a fast turnover rate such that soil C dynamics are primarily controlled by variability of the forcing (Luo et al., 2014).Due to the aforementioned reasons, the GP method failed at those three sites.
We first checked whether the oscillation and longer periodicity were caused by fire disturbance.However, this can only explain the oscillation at site US-SO2.The annual fire disturbance at site US-SO2 is longer than half a year, while it is less than a month at the other two sites.In the original CLM4, soil water is calculated first for the top ten soil lay- ers (3.8 m below the ground surface) and one aquifer layer using a water-content-based formulation for water mass conservation and a groundwater table as the bottom boundary condition (Oleson et al., 2010); the Niu et al. (2005Niu et al. ( , 2007) ) parameterizations are then used to simulate groundwatersoil water interaction and update the water table depth.If the water table is below 3.8 m, groundwater does not contribute to the soil moisture in the overlaying soil layers.We found that, after 100 years, the water table calculation scheme in CLM4 has resulted in a significantly different evolution of water table depth from one cycle to the next when repeatedly forcing the model with atmospheric data at sites RJA and US-IB1.The issue has also been found previously and an effort has been made to eliminate the oscillations (Oleson et al., 2010), but such oscillations can still occur under specific conditions such as at RJA and US-IB1.When we turned off the groundwater component, i.e., applying a zero flux boundary condition at the bottom of the soil column, we did not see oscillations in SOC at RJA and US-IB1.In the Niu et al. (2007) groundwater model, after solving the mass conservation equations (Richards' equation) in the top ten layers, water is then moved around to account for recharge and subsurface runoff and in the meantime to satisfy two conditions for water content in each layer; i.e., the water content has to be greater than the minimum content and smaller than the effective porosity of the layer.By moving water mass around after the Richards' equation is solved, the Richards' equation at each node is no longer satisfied if its moisture deviates from its previous solution.We have confirmed the local mass conservation error of water in the original model of CLM4.The error is large when recharge or subsurface runoff is high.The water content formulation itself has been previously shown to cause solution instability for soils near saturation (Hills et al., 1989).Instead of solving the soil water and groundwater separately, we use a flow model for variably saturated porous media, STOMP (Subsurface Transport Over Multiple Phases) (White and Oostrom, 2000), to see if it can resolve the oscillation in the total soil C.
The  of water mass crossing the control volume surface.The nonlinear equations describing mass conservation are discretized spatially on structured orthogonal grids using the integral finite difference approach of Patankar (1980), which is locally and globally mass conserving.The equations are discretized temporally using first-order backward Euler differencing or implicit time stepping that is suitable for the solution of the equations that are numerically unstable (LeVeque, 2007).Newton-Raphson iteration is used to resolve the nonlinearities from the constitutive equations that relate the primary and secondary variables.Detailed information regarding STOMP, such as the user's guide, theory guide and code availability, can be found at http://stomp.pnnl.gov.For each soil column, the number of vertical grids used for STOMP is 15 and it is the same as that in CLM4.In CLM4, the top ten grids (3.8 m below the ground) are used in the soil water scheme.The same initial saturation condition as that in CLM4 is prescribed.For the grid at the top, the Neumann boundary condition is used.For the bottom (42 m below the ground), a zero flux boundary condition is used.Because the aquifer is unconfined, we use the bottom node pressure to calculate water table depth.replacing the soil water and groundwater-soil water interaction scheme with STOMP at sites RJA and US-IB1.Using STOMP, mass conservation is improved, and the moisture content calculated is more accurate, resulting in wetter and cooler soil (Figs.5b, c and 6b, c).
Figure 7a and b shows the oscillations of water table depth resolved at both RJA and US-IB1; i.e., the oscillations between forcing cycles noted in the original hydrology scheme in CLM4 are caused by the local water mass balance error.Each cycle of atmospheric forcing at both sites has a length of 3 years.The GP method is successful at those two sites.The little jumps in Fig. 7c and d are where the GP method is applied.Both sites show higher total soil C predicted compared to Fig. 4a and b because of the new flow model.The issue of the original groundwater model in CLM4 might explain why it took so long (> 4000 years) for some of the grids in the global simulation to converge as shown in Shi et al. (2013).In addition, the gradient projection model is not recommended for sites where the length of the fire season is too long.For those sites, the overall time it takes for the spin-up run to steady state is much shorter compared to others; therefore, no improvement in the spin-up time is necessary.

Conclusions
We described a gradient projection method to further speed up the spin-up process based on the slow nature of soil organic C decomposition.Comparison between our approach and the modified accelerator approach showed that 20-69 % of simulation years can be reduced with our approach.While the approach was specifically evaluated using CLM4-CN, it can also be readily applied to other CN models in Earth system models.No matter what modification is made to improve the spin-up efficiency, a final spin-up is always needed to reach a converged solution due to disequilibrium caused by the modification.Our approach is especially useful when a new model formulation is proposed and a high-quality solution (small convergence threshold) is needed for a fair comparison.
In addition, we also found that the original numerical hydrology scheme, especially the water table calculation in CLM4, creates numerical oscillations in simulated water tables, leading to a challenge in achieving the common convergence criteria for soil C. To resolve the issue, we replaced the hydrological model using a flow model for variably saturated porous media.The new flow model caused an increase of about 10 % in computation time, but gives more accurate results that corrected the oscillation behavior of the original hydrological model.Comparing the total C predicted by the old and new flow models, we also see more C being predicted using the new flow model.Whether the prediction of more C is realistic depends on other factors besides the hydrology, so we have not attempted to evaluate the simulated C using observations.Nevertheless, a correct implementation of numerical schemes is always desirable for reducing uncertainty in model prediction.
, and a semianalytical steady-state solution for soil organic C and N pools Published by Copernicus Publications on behalf of the European Geosciences Union.Y. Fang et al.: Accelerating the spin-up of the coupled C and N cycle model in CLM4

Figure 1 .
Figure 1.Soil carbon pool structure of CLM4-CN.The arrows represent the decomposition pathways, and k is the turnover rate of each pool.

Figure 2 .
Figure 2. Annual average total soil carbon change with respect to time (a) and the gradient projection over a shorter time interval (b).

Figure 3 .
Figure 3. Stacked bar chart of percent reduction in computation cost for three convergence threshold values.
STOMP simulator was developed to predict nonisothermal hydrological flow and reactive transport in variably saturated subsurface environments.In STOMP, the water mass conservation equation balances the time rate of change of water mass within a control volume with the flux www.geosci-model-dev.net/8

Figure 4 .
Figure 4. Annual average total soil C with respect to time at sites RJA, US-IB1 and US-SO2.
Figures 5 and 6 show the model comparison at the beginning of the first 3 years between the simulations using the original soil hydrology scheme in CLM4 and the simulation after

Figure 5 .
Figure 5.Comparison of water table depth (a), average soil moisture content (b), and average soil temperature (c) using the original soil hydrology model and STOMP in CLM4 at site RJA.

Figure 6 .
Figure 6.Comparison of water table depth (a), average soil moisture content (b), and average soil temperature (c) using the original soil hydrology model and STOMP in CLM4 at site US-IB1.

Figure 7 .
Figure 7.Comparison of water table depth simulated by the original soil hydrology scheme in CLM4 (solid line) and STOMP (dashed line) at site RJA (a) and site US-IB1 (b) in the last 42 years of the simulation; (c) and (d) are annual average total soil carbon at sites RJA and US-IB1 using STOMP in CLM4 and the GP method.

Table 1 .
Location, PFT, soil type and number of years of atmospheric forcing for each site.