Accelerating the spin-up of the coupled carbon and nitrogen cycle model in CLM4

Abstract

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 pre-industrial 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), and a semi-analytical steady-state solution for soil organic C and N pools (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, 2014). Using our approach, we were able to implement the semi-analytical 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 spin-up 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: $$\begin{array}{}& & \frac{{\text{dC}}_{\text{Soil3}}}{\mathrm{d}t}=-k_{\text{S3}}{\mathrm{C}}_{\text{Soil3}}+k_{\text{S2}}{\text{C}}_{\text{Soil2}}+k_{\text{L3}}{\text{C}}_{\text{Litr3}},& & \frac{{\text{dC}}_{\text{Soil4}}}{\text{d}t}=-k_{\text{S4}}{\text{C}}_{\text{Soil4}}+k_{\text{S3}}{\text{C}}_{\text{Soil3}},\end{array}$$ where $k_{\text{L3}}$, $k_{\text{S2}}$, $k_{\text{S3}}$, and $k_{\text{S4}}$ are the turnover rates of the Litter3, Soil2, Soil3, and Soil4 C pools shown in Fig. 1, respectively. C${}_{\text{Litr3}}$, C${}_{\text{Soil2}}$, C${}_{\text{Soil3}}$ and C${}_{\text{Soil4}}$ are the amount of C in the Litter3, Soil2, Soil3 and Soil4 C pools, respectively. The first term on the right-hand sides of Eqs. (1) and (2) includes heterotrophic respiration. At the steady state, the left-hand sides of Eqs. (1) and (2) become 0; the amount of Soil3 C and Soil4 C can then be solved: $$\begin{array}{}& & {\text{C}}_{\text{Soil3}}=\frac{k_{\text{S2}}}{k_{\text{S3}}}{\text{C}}_{\text{Soil2}}+\frac{k_{\text{L3}}}{k_{\text{S3}}}{\text{C}}_{\text{Litr3}},& & {\text{C}}_{\text{Soil4}}=\frac{k_{\text{S3}}}{k_{\text{S4}}}{\text{C}}_{\text{Soil3}}.\end{array}$$ 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, $$\begin{array}{}& & {\text{C}}_{\text{Soil3,N}}=\frac{k_{\text{S2}}}{k_{\text{S3,N}}}{\text{C}}_{\text{Soil2}}+\frac{k_{\text{L3}}}{k_{\text{S3,N}}}{\text{C}}_{\text{Litr3}}=A_{\text{S3}}{\text{C}}_{\text{Soil3}},& & {\text{C}}_{\text{Soil4,N}}=\frac{k_{\text{S3}}}{k_{\text{S4,N}}}{\text{C}}_{\text{Soil3}}=A_{\text{S4}}{\text{C}}_{\text{Soil4}},\end{array}$$ where N denotes the native run, and $A$ is the accelerator.

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

[Figure omitted. See PDF]

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.

Methods

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_{\mathrm{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+\mathrm{1})m_{\mathrm{c}}$] years as a stop point for the accelerated decomposition (AD) run, where $n=300/m_{\mathrm{c}}$ is an integer. For example, if the number of years of forcing is $m_{\mathrm{c}}=\mathrm{7}$, the stop time will be at year 301. A stop point of $\sim$ 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: $$\text{C}\left(t_n\right)=\text{C}\left(t_{\mathrm{1}}\right)+\frac{\text{C}\left(t_{\mathrm{1}}\right)-C\left(t_{\mathrm{0}}\right)}{m_{\text{c}}}\left(t_n,-,t_{\mathrm{1}}\right),$$ where $t_{\mathrm{0}}$ is the beginning of the first cycle, $t_{\mathrm{1}}$ is the beginning of the next cycle, and $t_{\mathrm{1}}-t_{\mathrm{0}}=m_{\text{c}}$; $t_n-t_{\mathrm{1}}=\mathit{\tau }{m}_{\text{c}}$, $\mathit{\tau }$ 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^{\prime }\left(t\right)=ku\left(t\right)$, the stability requirement is $\mathrm{|}\mathrm{1}+kh\mathrm{|}\le \mathrm{1}$, in which $k$ is the rate constant and $h$ is the time step.

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

[Figure omitted. See PDF]

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_{\text{c}}\times (100/m_{\text{c}}+\mathrm{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${}^{-\mathrm{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.

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

ID	Site	Longitude	Latitude	Elev.	Height${}^{\mathrm{a}}$	CLM4	Clay	Sand	Silt	$m_{\mathrm{c}}^{\mathrm{d}}$
		(${}^{\circ }$ E)	(${}^{\circ }$ N)	(m)	(m)	PFT${}^{\mathrm{b}}$	(%)	(%)	(%)	Years
1	US-Ha1 (Goulden et al., 1996)	$-$72.1715	42.5378	343	30	BDTmp	6	66	29	1991–2006
2	US-WCr (Yi et al., 2004)	$-$90.0799	45.8059	520	30	BDTmp	20.17	42.52	37.32	1998–2006
3	US-Syv (Desai et al., 2005)	$-$89.3477	46.242	450	37	BDTmp	16.43	46.56	37.01	2001–2006
4	US-PFa (Davis et al., 2003)	$-$90.2723	45.9459	470	122	BDTmp	20.17	42.52	37.32	1995–2005
5	US-UMB (Curtis et al., 2005)	$-$84.7138	45.5598	234	50	BDTmp	0.6	92.6	6.8	1998–2006
6	US-MOz (Gu et al., 2007)	$-$92.2	38.7441	219	30	BDTmp	24.68	46.38	28.94	2004–2007
7	US-Dk2 (Katul et al., 2003)	$-$79.1004	35.9736	163	42	BDTmp	21.62	54.43	23.95	2003–2005
8	US-MMS (Sims et al., 2005)	$-$86.4131	39.3231	275	48	BDTmp	63	34	3	1999–2006
9	US-Ton (Baldocchi et al., 2010)	$-$120.966	38.4316	169	23	BDTmp and C3NAGrs	15	41	44	2001–2007
10	BAN${}^{\mathrm{c}}$	$-$50.1591	$-$9.82442	120	40	BDTrop	37	24	39	2004–2006
11	CA-Oas (Griffis et al., 2003)	$-$106.198	53.6289	530	30	BDBorl	18.8	50.32	30.87	1997–2006
12	CA-Gro (McCaughey et al., 2006)	$-$82.1556	48.2167	300	30	BDBorl	20	65	25	2004–2006
13	US-Ho1 (Hollinger et al., 1999)	$-$68.7403	45.2041	79	29	NETmp	15.9	50.35	33.75	1996–2004
14	CA-Ca1 (Humphreys et al., 2006)	$-$125.334	49.8672	300	45	NETmp	2.63	84.42	12.94	1998–2006
15	CA-TP4 (Arain and Restrepo-Coupe, 2005)	$-$80.3574	42.7098	219	30	NETmp	0	98	2	2002–2007
16	US-NR1 (Turnipseed et al., 2002)	$-$105.546	40.0329	3050	26	NETmp	21.43	43.13	35.45	1998–2007
17	US-Dk3 (Katul et al., 2003)	$-$79.0942	35.9782	163	21	NETmp	13.66	51.59	34.81	1998–2005
18	US-Me2 (Hudiburg et al., 2013)	$-$121.557	44.4524	1310	30	NETmp	7	67	26	2002–2007
19	CA-Obs (Griffis et al., 2003)	$-$105.118	53.9872	629	30	NEBorl	4.12	80.89	14.97	2000–2006
20	CA-Qfo (Chen et al., 2006)	$-$74.3421	49.6925	382	25	NEBorl	4	51.5	29	2004–2006
21	CA-Ojp (Griffis et al., 2003)	$-$104.692	53.9163	579	30	NEBorl	2.5	94.47	3.02	2000–2006
22	K67${}^{\mathrm{c}}$	$-$54.9589	$-$2.85667	130	63	BETrop	90	2	8	2002–2004
23	K83${}^{\mathrm{c}}$	$-$54.9714	$-$3.01803	130	64	BETrop	80	18	2	2001–2003
24	RJA${}^{\mathrm{c}}$	$-$61.9309	$-$10.0832	191	60	BETrop	10	80	10	2000–2002
25	K77${}^{\mathrm{c}}$	$-$54.8944	$-$3.01983	130	18	Crop	80	18	2	2001–2005
26	FNS${}^{\mathrm{c}}$	$-$62.3572	$-$10.7618	306	8.5	Crop	10	80	10	1999–2001
27	US-Ne2 (Suyker and Verma, 2012)	$-$96.4701	41.1649	362	6	Crop	31.68	30.7	37.62	2001–2006
28	US-Ne1 (Suyker and Verma, 2012)	$-$96.4766	41.1651	361	6	Crop	31.68	30.7	37.62	2001–2006
29	US-IB1 (Matamala et al., 2008)	$-$88.2227	41.8593	225	4	Crop	37.2	7.8	55.4	2005–2007
30	US-Ne3 (Suyker and Verma, 2012)	$-$96.4397	41.1649	363	6	Crop	31.68	30.7	37.62	2001–2006
31	US-ARM (Fischer et al., 2007)	$-$97.4884	36.605	311	65	Crop	43.1	27.98	28.92	2000–2007
32	PDG${}^{\mathrm{c}}$	$-$47.6499	$-$21.6195	690	21	C3NAGrs and C4Grs	3	85	12	2001–2003
33	US-IB2 (Matamala et al., 2008)	$-$88.241	41.8406	226	4	C3NAGrs and C4Grs	34.8	12.18	53	2004–2007
34	CA-Let (Flanagan et al., 2002)	$-$112.94	49.7093	960	4	C3NAGrs	35.6	28.1	34.8	1997–2006
35	US-Var (Baldocchi et al., 2004)	$-$120.951	38.4133	129	2	C3NAGrs	12.5	29.5	58	2001–2007
36	US-Shd (Suyker et al., 2003)	$-$96.6833	36.9333	350	5	C3NAGrs	38.4	5.1	56	1997–2000
37	US-Los (Yi et al., 2004)	$-$89.9792	46.0827	480	10	BEShr	16.43	46.56	37.01	2000–2006
38	US-SO2 (Lipson et al., 2005)	$-$116.623	33.3739	1406	5	BEShr	21.31	43.94	34.75	1998–2006

${}^{\mathrm{a}}$ Approximate height of the wind/temperature and flux measurements above the surface. ${}^{\mathrm{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. ${}^{\mathrm{c}}$ The site information and meteorological forcing are from the LBA-MIP data set. ${}^d$ $m_{\text{c}}$ is the number of years of atmospheric forcing.

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 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 $\mathrm{\Delta }{\text{C}}_{\text{TOC}}$ (3.0, 1.0, and 0.5 g m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{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 ($\mathrm{\Delta }{\text{C}}_{\text{TOC}}\le 0.5$ g m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{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 ($\mathrm{\Delta }{\text{C}}_{\text{TOC}}\le \mathrm{3}$ g m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}})$.

Comparison between the gradient projection method (in bold) and the accelerated spin-up method.

ID	Site	Number of simulation years to reach
		$\mathrm{\Delta }{\text{C}}_{\text{TOC}}\le \mathrm{3}$	$\mathrm{\Delta }{\text{C}}_{\text{TOC}}\le \mathrm{1}$	$\mathrm{\Delta }{\text{C}}_{\text{TOC}}\le 0.5$
		(g m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$)	(g m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$)	(g m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$)
1	US-Ha1	416/416	816/480	1024/624
2	US-WCr	828/558	1395/729	1809/828
3	US-Syv	930/600	1314/744	1680/924
4	US-PFa	1375/617	2057/891	2255/946
5	US-UMB	387/387	855/567	1116/567
6	US-MOz	772/596	1400/924	1776/1096
7	US-Dk2	453/408	888/696	1215/849
8	US-MMS	872/552	1136/704	1400/744
9	US-Ton	413/402	749/574	959/672
10	BAN	1281/1050	1677/1284	1959/1452
11	CA-Oas	1870/860	2170/1090	2470/1240
12	CA-Gro	351/351	966/774	1455/1023
13	US-Ho1	972/585	1503/756	1845/864
14	CA-Ca1	597/468	972/549	1215/558
15	CA-TP4	798/528	1170/636	1224/828
16	US-NR1	1310/740	1910/870	2470/1000
17	US-Dk3	640/432	848/472	992/488
18	US-Me2	312/312	318/318	354/354
19	CA-Obs	509/441	1029/700	1351/770
20	CA-Qfo	384/384	819/612	1143/816
21	CA-Ojp	520/441	812/539	1143/651
22	K67	543/447	720/510	831/612
23	K83	354/354	477/411	633/510
24	RJA	NA${}^{*}$/348	NA/510	NA/606
25	K77	310/310	440/425	585/445
26	FNS	468/432	1008/750	1350/954
27	US-Ne2	954/726	1368/942	1620/1098
28	US-Ne1	732/534	1200/714	1506/858
29	US-IB1	NA/309	NA/333	NA/459
30	US-Ne3	312/312	648/474	948/612
31	US-ARM	864/552	1192/616	1400/672
32	PDG	309/309	663/540	1077/807
33	US-IB2	536/440	804/588	900/628
34	CA-Let	630/450	960/490	1160/560
35	US-Var	608/427	881/651	1568/679
36	US-Shd	1784/1168	2128/1368	2320/1472
37	US-Los	490/427	903/539	1169/651
38	US-SO2	NA/NA

${}^{*}$ NA – not evaluated

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 atmospheric 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.

Stacked bar chart of percent reduction in computation cost for three convergence threshold values.

[Figure omitted. See PDF]

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 layers (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. (2005, 2007) parameterizations are then used to simulate groundwater–soil 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.

Annual average total soil C with respect to time at sites RJA, US-IB1 and US-SO2.

[Figure omitted. See PDF]

The STOMP simulator was developed to predict non-isothermal 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 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.

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 omitted. See PDF]

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. 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 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).

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 omitted. See PDF]

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.

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.

[Figure omitted. See PDF]

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.

Code availability

The source code of CLM4.0 and STOMP can be requested through http://www.cesm.ucar.edu/models/cesm1.0/ and http://stomp.pnnl.gov/licensing.stm, respectively. The method implemented in this study can be obtained upon request. Contact: [email protected].

Acknowledgements

This research has been accomplished through funding support from Pacific Northwest National Laboratory's Laboratory Directed Research and Development Program. We thank the North American Carbon Program Site-Level Interim Synthesis team, the Large Scale Biosphere-Atmosphere Experiment in Amazônia Model Intercomparison Project team, and the site investigators for collecting, organizing, and distributing the data required for this analysis. We thank the anonymous referee and Yiqi Luo for comments that improved the manuscript. A portion of this research was performed using PNNL Institutional Computing at Pacific Northwest National Laboratory. PNNL is operated by Battelle for the US Department of Energy under contract DE-AC05-76RL01830. Edited by: S. Arndt