Response of microbial decomposition to spin-up explains CMIP5 soil carbon range until 2100

Abstract. Soil carbon storage simulated by the Coupled Model Intercomparison Project (CMIP5) models varies 6-fold for the present day. Here, we confirm earlier work showing that this range already exists at the beginning of the CMIP5 historical simulations. We additionally show that this range is largely determined by the response of microbial decomposition during each model's spin-up procedure from initialization to equilibration. The 6-fold range in soil carbon, once established prior to the beginning of the historical period (and prior to the beginning of a CMIP5 simulation), is then maintained through the present and to 2100 almost unchanged even under a strong business-as-usual emissions scenario. We therefore highlight that a commonly ignored part of CMIP5 analyses – the land surface state achieved through the spin-up procedure – can be important for determining future carbon storage and land surface fluxes. We identify the need to better constrain the outcome of the spin-up procedure as an important step in reducing uncertainty in both projected soil carbon and land surface fluxes in CMIP5 transient simulations.


Introduction
The land surface currently absorbs about a third of anthropogenic emissions of CO 2 (Canadell et al., 2007;Le Quéré et al., 2009) and so helps to offset global warming. Future global warming may enhance microbial decomposition and emissions of CO 2 from respired soil organic carbon (SOC), the largest carbon pool in the terrestrial biosphere (Jobbágy and Jackson, 2000). Higher emissions from SOC could accelerate increases in atmospheric CO 2 concentrations even if plant carbon uptake by photosynthesis increased under higher atmospheric CO 2 (Ahlström et al., 2013;Friedlingstein et al., 2014;Nishina et al., 2014). Conversely, if the soil remains a carbon sink (Le Quéré et al., 2009;Lund et al., 2010) the negative feedback on rising atmospheric CO 2 (Davidson and Janssens, 2006) would help limit rates of increase. How soil carbon is represented in models and how it responds to climate is critical to resolving whether the land will remain a sink or become a source of CO 2 .
Recent model intercomparisons, such as the fifth phase of the Coupled Model Intercomparison Project (CMIP5; Taylor et al., 2012) and the Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP; Warszawski et al., 2014), have highlighted a lack of consensus among models on whether the soil carbon sink will be sustained during the 21st century (Friedlingstein et al., 2014;Nishina et al., 2014). These models also exhibit large discrepancies in stores of SOC they simulate. For example, Todd-Brown et al. (2013) report that total SOC simulated by CMIP5 models for the present day represents a 6-fold variation ranging from ∼ 510 to ∼ 3040 Pg C. Another large range (∼ 1090 to ∼ 2645 Pg C) exists in the present day SOC simulated by ISI-MIP models despite being driven by a harmonized weather data set . These latter results indicate that a significant fraction of the uncertainty in estimates of total SOC arises from the representation of land processes rather than differences in climate drivers. Soil carbon pools of widely different sizes have the potential to react differently to future climate change. We therefore examine the likely reasons for the large differences between CMIP5 models in their simulation of SOC. This work is founded in the recognition that the SOC varies among the CMIP5 models for the present day over a 6-fold range (Todd-Brown et al., 2013) and this range contributes to model-tomodel variations in SOC change in the future (Todd-Brown et al., 2014). We explore why this 6-fold range exists and ultimately show that individual model responses to the spinup procedure, particularly the dominant role of turnover time relative to SOC input, are the key reason for this range. Explaining why the amount of carbon mobilized in the active cycle varies greatly between models is critical but has been largely ignored in the literature to date. As noted by Knutti and Sedláček (2013), there may be multiple sources of disagreement between models such as a lack of process understanding, or the reduced availability of relevant observational data sets to constrain models. Technical aspects of climate modelling, such as how different state variables are initialized or spun up to an equilibrated state prior to an experiment being conducted, and how equilibration is defined in this context, can also lead to major differences between model simulations. Discriminating between these sources of uncertainty to understand why CMIP5 models differ so significantly in the amount of SOC in the present day, and subsequently in the total amount of C mobilized in the global cycle under a future climate, enables an improvement in model projections. Increasing the consistency between models is required to improve our confidence in the sign of the soil carbon feedback in the future.
To avoid misconceptions, we define and differentiate between two states that are commonly called "initial" states in land modelling. Our definition of "initial state", which is not known or reported in CMIP5 models, is the state at the beginning of a climate model integration. This "initial state" may come from a previous simulation, from off-line simulations, from observations or via expert judgement. In the case of SOC, it may be initialized as a "cold start" or in a state equilibrated with an atmosphere that reflects the period prior to the beginning of a simulation. This model state is then commonly integrated forward in time until those model states that are considered important are in equilibrium with the atmospheric model over some period of time and to a degree that is defined by the modeller (but not reported). This generates what we define as an "equilibrated state". In CMIP5, simulations are then reported from the beginning of the historical period (say 1850), initialized with this "equilibrated state" and integrated forward in time to the present day under observed forcings, and then into the future using a representative concentration pathway (Taylor et al., 2012). The values of a climate model's state variables at 1850 are commonly thought of as the "initial state" but they are not; it is the model-specific equilibrated state under pre-industrial forcing and this reflects the ability of the climate model to represent global and regional temperatures, rainfall and so forth. We therefore call this the "equilibrated state" and note that this differs from the "initial state" due to the earth system model's simulated climate, the definition of "equilibrium" over time and space and crucially how the state variables are parameterized. Here we show that a great deal of the 6-fold range in SOC in the CMIP5 models at the "equilibrated state" assumed representative of 1850 (and consequently in the present day reported by Todd-Brown et al., 2013) is a consequence of the procedures used to evolve the model from the "initial state" to the "equilibrated state". These procedures may influence how SOC changes through to 2100 (Todd-Brown et al., 2014) due to the current stateof-the-art representation of SOC decomposition.

SOC in earth system models
In all global terrestrial models participating in recent intercomparison projects such as CMIP5 and ISI-MIP, the SOC balance and its change ( SOC) are represented in a similar way. First, inputs of carbon into the soil are derived from plant pools. Plant carbon uptake and turnover times respond to climate change, climate variability and atmospheric CO 2 independent of the size of the SOC pools. Meanwhile, modelled microbial decomposition releases carbon by heterotrophic respiration (R h ). The balance can be summarized by where SOC in is the input to the SOC pools from plant and litter pools. Microbial decomposition is commonly represented as a first-order process and applied to a succession of pools. In each pool, a parameter k reflects the specific baseline decomposition rate (Xia et al., 2013;Exbrayat et al., 2013a, b) at a reference soil temperature and non-limiting moisture conditions. Then, the decay rate is adjusted at each time step by an environmental scalar (Todd-Brown et al., 2013;Xia et al., 2013;Exbrayat et al., 2013a, b;Nishina et al., 2014) that describes the instantaneous response of microbial activity to the soil physical state as the product of a soil temperature (f T ) and a soil moisture respiration function (f W ). Various formulations of f T and f W have been implemented in model codes (Lloyd and Taylor, 1994;Falloon et al., 2011;Todd-Brown et al., 2013;Exbrayat et al., 2013a, b;Nishina et al., 2014), usually assuming a space-and time-invariant response to the same conditions. Their effect on decay rate varies according to local soil conditions and therefore climate.
The actual decay rate (k × f T × f W ) is applied to the amount of substrate available, SOC, to determine the amount of microbial decomposition D m at each model time step: where k × f T × f W is equivalent to the fraction of respired substrate, the inverse of the turnover time SOC/R h . A part of the decomposed organic matter is routed to pools with longer turnover time and the rest is emitted as CO 2 . There may be variations between models in the number of pools they represent (Todd-Brown et al., 2013;Nishina et al., 2014) and the formulations of the environmental response functions (Falloon et al., 2011;Exbrayat et al., 2013a) but at the ecosystem scale, R h is proportional to the amount of substrate, i.e. SOC, available in the soil. This parameterization may be inconsistent with our current understanding of microbial decomposition (Allison et al., 2010;Schmidt et al., 2011;Wieder et al., 2013) because it lacks the representation of processes such as microbial activity and priming effect (e.g. Xenakis and Williams, 2014). However, the first-order dependency of R h on SOC, soil temperature and moisture is able to explain complex phenomena such as the apparent acclimation of decomposers to warming by quick depletion of the most labile substrate pools (Luo et al., 2001;Kirschbaum, 2004;Knorr et al., 2005).

CMIP5 data
From the CMIP5 archive we downloaded monthly soil carbon density (cSoil in metadata), litter carbon density (cLitter) and heterotrophic respiration (rh) for 15 CMIP5 models from 10 international institutions. A list of models can be found in Table 1 while further details about models and land components have been summarized in Table 2. We note that four of these models, namely BCC-CSM1.1 (model A), CCSM4 (model C), NorESM1-M and NorESM1-ME (grouped as model J), represent nitrogen limitation on plant productivity while the others do not. We selected data for the historical (1850-2005) and the most intensive Representative Concentration Pathway 8.5 (RCP 8.5, 2006(RCP 8.5, -2100 experiments. A total of 79 simulations for the historical experiment, including 34 simulations continuing for RCP 8.5 (Table 1) were available. When cLitter was reported, we added it to cSoil as both pools are parameterized to generate R h following first-order kinetics.
To calculate stock sizes we first multiplied spatially explicit data of cSoil and cLitter in kg C m −2 by corresponding grid-cell areas (areacella in metadata) and integrated their values globally. Similarly, we calculated global fluxes of R h by multiplying monthly fluxes in kg C m −2 by grid-cell areas and integrating them globally. Fluxes were summed to obtain annual averages. Annual soil carbon input (SOC in ) from above-ground biomass was not available from the database. Therefore, we calculated it by inverting the SOC balance: As models did not start their historical simulations at the same time, we focus our analyses on the overlapping period of 1861-2100. We also averaged all simulations from the same model or institution in an attempt to account for model dependence (see Bishop and Abramowitz, 2013, for a discussion on the topic).
In the following, we report values of stocks and fluxes averaged for three periods of time, the pre-industrial (1861-1870), modern (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) and future (2091-2100) periods. While the period 1861-1870 is not part of the preindustrial control runs sensu stricto, the minor increase in atmospheric CO 2 between pre-industrial times (i.e. before 1850) and 1870 is unlikely to have led models to simulate  (Chylek et al., 2011) CTEM ( Table 3. Then, we multiply the organic carbon content of the dominant soil units (in % weight) by the bulk density (provided in kg dm −3 ) to obtain the carbon density (in kg C m −2 ) in each 0.5 • × 0.5 • grid cell. We multiply the density by the surface area of each grid cell and sum results to obtain a total soil carbon content of ∼ 1170 Pg C. Following Todd-Brown et al. (2013), a confidence interval of 29 % below the mean (i.e. ∼ 830 Pg C) to 32 % above the mean (i.e. ∼ 1550 Pg C) was considered to take variations in soil carbon content and the mapping processes into account. The range we obtain is slightly smaller than reported by Todd-Brown et al. (2013) (890-1660 Pg C) because we use an updated version of the HWSD and did not replace bulk density values for Andisols and Histosols.

Results
We first compare total SOC for pre-industrial (1861-1870), modern (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) and future (2091-2100) periods. Figure 1 compares the total SOC range in CMIP5 models for 1861-1870 (563-2938 Pg C), 1996-2005 (576-3047 Pg C), and 2091-2100 (582-3266 Pg C, derived using the RCP 8.5 scenario). All three periods show very similar distributions of SOC among the models and the present day and future ranges already exist at the beginning of the historical simulations. Figure 1 highlights that the size of SOC pools of individual CMIP5 models remains largely consistent over the three time periods. Indeed, pre-industrial SOC predicts modern SOC, modern SOC predicts future SOC and pre-industrial SOC predicts future stocks with a high degree of precision (Fig. 1). Also represented in Fig. 1 is the 95 % confidence interval of total SOC estimated from HWSD that we use as a reference for modern total SOC (i.e. in 1996SOC (i.e. in -2005. We note that only three models fall within this range: BCC-CSM1.1 (model A), CanESM2 (model B) and HadGEM2 (model F). Models based on the CLM4 land surface model (i.e. models C and J) underestimate modern SOC while all remaining models overestimate it. Note that these models C and J include nitrogen limitation of the vegetation response to increasing CO 2 .
We next investigate the likely reasons for the existence of this pre-industrial CMIP5 range in total SOC. The first obvious step is to check whether models are at equilibrium prior to climate change experiments. Models may not agree on total SOC simply because some of them, and especially those at the extremes of the CMIP5 spectrum, are still drifting towards their own steady-state and therefore do not comply with our experiment protocol. In Fig. 2 we show the relationship between pre-industrial SOC in and R h . This relationship is highly significant (R 2 = 1; p < 0.001) and strongly suggests that all models were equilibrated under pre-industrial boundary conditions. This removes the possibility that models were not in equilibrium and means that the 6-fold CMIP5 range is likely linked with the internal terrestrial processes represented in these models.
Two major internal terrestrial processes are involved: SOC in , the amount of SOC that enters the soil pools, and the turnover time of organic matter that corresponds to the amount of SOC that is released from soil pools. The relationship between SOC in and total SOC during the preindustrial period is shown in Fig. 3. Overall, the relationship is not significant (R 2 = 0.04; p = 0.604). Further, the models that equilibrate with the largest total SOC stock (models E, H, I) are not the models with the largest SOC input. Similarly, the small equilibrated SOC pool size of models C and J seems unrelated to SOC in despite these models including N limitations on plant productivity and SOC in . In short, the amount of SOC in cannot explain the size of the equilibrated pools. In Fig. 4, we therefore present the relationship between the pre-industrial SOC turnover time (i.e. the inverse of the decay rate expressed as SOC/R h ) and total SOC. This relationship is highly significant (p < 0.001) and linear (R 2 = 0.84) and models with a longer turnover time, i.e. a low decay rate, require larger pools to offset the same SOC input, and vice versa. Further, turnover times are not affected by the number of SOC pools represented. Models with the longest turnover time have alternatively nine (model E) or two pools (models H and I), while models with the shortest turnover time have eight (model A), six (models C and J) or four pools (model F).

Discussion
Despite the change imposed on boundary conditions during global warming experiments (Anav et al., 2013;Friedlingstein et al., 2014), CMIP5 present day and projected SOC stocks are largely determined by their equilibrated pool size (Fig. 1) in 1860. This was not unexpected due to the slow response of SOC pools but it clearly shows that modern and future stocks are mostly defined by the equilibrated pool size while changes can be explained by a combination of changes in the input and output fluxes (see Todd-Brown et al., 2014, for a detailed account of these mechanisms). Further, as SOC in 1860 is unknown from observations, CMIP5 models use a spin-up procedure from an initial state assuming steady preindustrial boundary conditions (Xia et al., 2012) to obtain an equilibrated state for pre-industrial SOC. In order to reach equilibrium, iterative or semi-analytical methods (e.g. Xia et al., 2012) are employed to reach the pool sizes required to balance input (SOC in ) and output fluxes (R h ). Steady-state is assumed when the trend in SOC becomes negligible. Hence, it is not the actual value of SOC that defines the equilibrium but its lack of variation in time (Xia et al., 2013;Exbrayat et al., 2013b). It is worrisome that these procedures are not clearly documented and therefore how a model is evolved from its true "initial state" to its "equilibrated state" is not known.
However, we have verified that all CMIP5 models were close to equilibrium prior to the initiation of climate change experiments. Following Eqs. (1) and (2), the model-specific value of SOC obtained by a model via spin-up depends on two factors. First, if SOC in is large, a larger SOC pool is required to offset it through microbial decomposition and R h , for a given decay rate, k × f T × f W . Conversely, low values of SOC in lead SOC pools to equilibrate to lower values for a particular decay rate. Second, if the decay rate is high (short turnover time) during spin-up, SOC pools will remain small, for a given SOC in . Conversely, low decay rates, or long turnover time, will require large pools of substrate to offset the same input SOC in . Both factors are modelspecific: SOC in is derived from plant primary productivity fluxes (Davidson and Janssens, 2006) while the baseline decomposition rate k and the shape of the response functions f T and f W are highly model-dependent (Falloon et al., 2011;Exbrayat et al., 2013a, b;Todd-Brown et al., 2013).
Here we have shown that the large range exhibited by CMIP5 SOC is principally due to the response of microbial decomposition during the spin-up process. This is a long process that corresponds to multiple centuries of steady climate conditions but as noted is not reported as part of CMIP5 and might represent a short period if the "initial state" is already well equilibrated or may represent many centuries if not. Throughout this period, however, for each CMIP5 model, model-specific parameter k and environmentalresponse functions f T and f W drive SOC pools to the size required by the Future total SOC [Pg C] y = 1.04x + 34.1 R 2 =0.96; p < 0.001 Figure 1. Relationship between total SOC in CMIP5 models at two different times: modern stocks as a function of pre-industrial stocks (upper panel), future stocks as a function of modern stocks (middle panel) and future stocks as a function of pre-industrial stocks (lower panel). Letters correspond to models as in Table 1 and models in green (i.e. C and J) integrate nitrogen limitation. The grey area is the 95 % confidence interval of modern total SOC derived from the HWSD. Equation, R 2 and p values correspond to the linear relationship between stocks built using data from all models (solid line). The dotted line is the 1 : 1 line.  Table 1 and models in green (i.e. C and J) integrate nitrogen limitation. The solid line is a linear relationship constructed using all models with equation, R 2 and p values indicated in the top left corner. The dotted line represents the 1 : 1 relationship. turnover time they simulate to compensate for SOC in . This observation corroborates the predominance of turnover time in the uncertainty of ecosystem response to climate change  and Fig. 4 shows that it is independent of the number of pools considered in each model. The resulting equilibrated state obtained prior to the initiation of CMIP5 transient simulations propagates through the present and into the future even when one is using RCP 8.5.
Our results raise a critical problem linked to model initialization and then equilibration by spin-up. According to our analysis of the CMIP5 models, a simple solution to reduce the uncertainty in simulated SOC stocks would be to modify model parameters, especially those related to SOC turnover, to obtain a steady-state consistent from model to model with SOC values representative of pre-industrial conditions. Alternatively, because of the millennial timescales of soil genesis, as well as land use changes, steady-state of global SOC stocks is not guaranteed to have existed at the end of the preindustrial era. Therefore, one could choose to consider only model parameters that achieve modern stocks in accordance with observations in response to past changes (e.g. Exbrayat et al., 2014). However, this would require multiple realizations of computationally expensive models, or the use of emulators. Furthermore, it would be necessary to represent site history, and especially disturbances, with a high degree of confidence during simulations to avoid over-fitting parameters and this may not be realistic at global scale. Therefore, assuming an equilibrated pre-industrial state is a more readily available option that is supported by the lack of variations in simulated SOC during historical experiments despite changing boundary conditions. Pre-industrial total SOC [Pg C] y = 11.31x + 914.0 R 2 =0.04; p = 0.604 Figure 3. Relationship between pre-industrial SOC input and preindustrial total SOC stocks at the beginning of the historical experiment. Letters correspond to the same models as in Table 1 and models in green (i.e. C and J) integrate nitrogen limitation. The solid line is a linear relationship constructed using all models with equation, R 2 and p values indicated in the top left corner.
Thus, we suggest that one could use available estimates and confidence interval of modern SOC stocks to constrain the pre-industrial equilibrated state. These estimates include global data sets such as HWSD and other (Shangguan et al., 2014) but also regional data that may better represent high latitude stocks and permafrost (e.g. Northern Circumpolar Soil Carbon Database; Hugelius et al., 2013). Of course, while changing parameter values corresponding to SOC turnover time is relatively straightforward, it would be important to ensure that these pools are sustained by an input representative of carbon uptake. At equilibrium SOC in equals net primary productivity (NPP) because plant pools do not vary in size. Here all models predict SOC in within two standard deviations of the uncertainty range of modern, high confidence, NPP estimates (56.4 ± 8-9 Pg C yr −1 ; Ito, 2011). Although not directly comparable with pre-industrial values, this global estimate indicates that models simulate acceptable values of global carbon uptake.
As decomposition processes are represented following first-order kinetics, simulating more realistic SOC stocks from an initial condition, and through spin-up to an equilibrated state in response to adequate uptake fluxes would likely lead models to represent more correct modern stocks. Nevertheless, as each model relies on its own formulation of the response functions f T and f W , the ensemble would still exhibit different sensitivities of SOC stocks to climate change. However, by removing a degree of freedom associated with spin-up procedures, we believe that these observational data sets are a valuable tool for increasing the consistency between models and making them more comparable. It would improve the confidence we can have in projections of SOC fluxes and feedbacks on future climate change.

Conclusions
We have demonstrated that the 6-fold range in SOC stocks simulated by CMIP5 models can be explained by the modelspecific response of microbial decomposition to spin-up under pre-industrial conditions. Model-dependent parameter and response functions drive the size of the pools to the amount required by decay rates to offset SOC in under the steady-state assumption. Once established, the resulting pool sizes remain similar through to the present and into the future even under the high-emission RCP8.5 scenario that generates future conditions the least similar to current ones. We therefore identify the spin-up procedure, and especially the response of microbial decomposition during this very long model integration, as a key source of uncertainty in the simulation of SOC in CMIP5 models. Critically, this involves the interaction of a technical and a process-linked uncertainty in CMIP5 models' experimental framework. The technical methods used for spin-up are model specific and not commonly reported. Interlinked with the technical uncertainty is the parameterization of processes within the spin-up period.
A model that equilibrates to a soil carbon store well outside the observed range should be examined with care. A very large amount of stored carbon increases the potential for the land surface to become a source as even a tiny relative change in decay rate can strongly enhance R h and possibly reach a tipping point where it offsets increases in SOC in . Conversely, a very small SOC store increases the likelihood that it will remain a sink. Such results are likely to be artefacts of model implementation when SOC values are largely inconsistent with observed ranges.
In conclusion, we recommend that future intercomparisons should constrain model parameters so that each model achieves an equilibrated state similar to observations as the outcome of the spin-up procedure. This would remove a degree of freedom associated with the process linking initialization to equilibration via a poorly constrained spin-up procedure when comparing differences in projected changes.