Examining soil carbon uncertainty in a global model

. Reliable projections of future climate require land–atmosphere carbon (C) ﬂuxes to be represented real-istically in Earth system models (ESMs). There are sev-eral sources of uncertainty in how carbon is parameterised in these models. First, while interactions between the C, nitrogen (N) and phosphorus (P) cycles have been implemented in some models, these lead to diverse changes in land–atmosphere ﬂuxes. Second, while the ﬁrst-order parameterisation of soil organic matter decomposition is similar between models, formulations of the control of the soil physi-cal state on microbial activity vary widely. For the ﬁrst time, we address these sources of uncertainty simultaneously by implementing three soil moisture and three soil temperature respiration functions in an ESM that can be run with three degrees of biogeochemical nutrient limitation (C-only, C and N, and C and N and P). All 27 possible combinations of response functions and biogeochemical mode are equilibrated before transient historical (1850–2005) simulations are performed. As expected, implementing N and P limitation reduces the land carbon sink, transforming some regional sinks into net sources over the historical period. Meanwhile, regardless of which nutrient mode is used, various combinations of response functions imply a two-fold difference in the net ecosystem accumulation and a four-fold difference in equilibrated total soil C. We further show that regions with initially larger pools are more likely to become carbon sources, especially when nutrient availability limits the response of primary production to increasing atmospheric CO 2 . Simulating changes in soil C content therefore criti-cally depends on both nutrient limitation and the choice of respiration functions.


Introduction
A major step in the transition from climate system models to Earth system models (ESMs) is the inclusion of additional biogeochemical processes. If the carbon (C) cycle was in equilibrium this would be an academic exercise, but how terrestrial C stores respond to warming resulting from human emissions of atmospheric carbon dioxide (CO 2 ) is of critical importance. Vegetation stores around 450-650 Pg C (Prentice et al., 2001), while soils store 1500-2400 Pg C (Batjes, 1996), with additional carbon stored in peatland and wetland soils (∼ 530 Pg C; Bridgham et al., 2006) and in permafrost (∼ 1670 Pg C; Tarnocai et al., 2009). If vegetation and soil processes respond to global warming by increasing the terrestrial C sink, this could help offset human emissions. Conversely, any decrease in the magnitude of the terrestrial sink, or any progressive loss of stored terrestrial C would provide a positive feedback on global warming.
Our current understanding is that human-induced increases in atmospheric CO 2 likely enhanced the terrestrial C uptake during the 20th century (Sarmiento et al.,  than global warming has enhanced microbial decomposition and corresponding release by heterotrophic respiration (R h ). It is, however, uncertain whether this increase can be sustained into the future (McCarthy et al., 2010;Norby et al., 2010;Zak et al., 2011). Indeed, some ecosystems appear to lack any significant response to increasing atmospheric CO 2 (Adair et al., 2009;Bader et al., 2009;Norby et al., 2010). According to previous modelling studies, any additional terrestrial carbon uptake linked to CO 2 fertilisation is also likely to be more than offset in the future by the increase in R h following warming (Cox et al., 2000). This extra C released into the atmosphere would further accelerate global warming (Kirschbaum, 2000), and a climate-change-driven acceleration of soil organic C decomposition rates would therefore represent a positive feedback on climate (Kirschbaum, 2004). However, there is a lack of agreement between modelbased estimates of when and at what rate soil C storage might begin to decline (Friedlingstein et al., 2006). Further, net primary production (NPP) and microbial decomposition are controlled by the availability of nitrogen (N) and phosphorus (P), with N limitation tending to dominate in temperate and boreal ecosystems, and P limitation tending to dominate in the tropics (Luo et al., 2004;Vitousek et al., 2010;Wang et al., 2010;Goll et al., 2012). Recently, there has been an extensive effort to implement these processes in terrestrial ecosystem models (e.g. Thornton et al., 2007;Sokolov et al., 2008;Yang et al., 2009;Wang et al., 2010;Zaehle and Friend, 2010;Esser et al, 2011;Menge et al., 2012). Generally, adding N limitation reduces the simulated global land C uptake during the 20th century relative to non-nutrientlimited simulations. Early results suggest P limitation makes a negligible difference to the global terrestrial carbon uptake, but can introduce very large regional differences particularly in the tropics (Zhang et al., 2011). However, despite the recognition of the importance of interactions between these biogeochemical cycles, interactions between terrestrial C and N cycles are represented in just three of the ESMs used in the Coupled Model Intercomparison Project, Phase 5 (CMIP5; Taylor et al., 2012), among which two models use the Community Land Model as their terrestrial component: CCSM4 and NorESM (Todd-Brown et al., 2013). Meanwhile, the terrestrial P cycle is omitted in all CMIP5 simulations (Todd-Brown et al., 2013). This introduces critical uncertainties in projections as nutrient limitation prevents vegetation growth at the rate allowed for by CO 2 fertilisation in different ways between models.
An additional uncertainty resides in the current parameterisation of microbial decomposition and corresponding R h in ESMs. So far, all CMIP5 models represent decomposition as a first-order process (Todd-Brown et al., 2013) in which instantaneous soil moisture and soil temperature are used to adjust a time-invariant decay rate that is applied to the amount of substrate available (i.e. C pool size). Put in a mathematical way, at each time step, the actual amount of microbial decomposition D m in a specific C pool is calculated as with k the reference decay rate that is scaled by f T , a function of soil temperature T s , and by f W , a function of soil moisture θ s (usually expressed as a fraction of water saturation; Moyano et al., 2012), and C s is the amount of C in the pool. The product f W (θ s )×f T (T s ) is sometimes referred to as an environmental scalar. Part of the decomposition is emitted as CO 2 , and this flux corresponds to R h , while the rest is typically assigned to different soil C pools. The lack of physiological control has been recently identified as being inconsistent with our current understanding of decomposition process (e.g. Allison et al., 2010;Schmidt et al., 2011). Nevertheless, first-order kinetics applied to a succession of C pools with different residence time are able to explain complex processes including the apparent thermal acclimation of decomposers to warming (Luo et al., 2001) with a quick depletion in the most labile pools (Kirschbaum, 2004;Knorr et al., 2005). However, in current models, simple changes in the formulation of f W , the soil moisture-respiration function, and f T , the soil temperature-respiration function can have a major influence on R h (Falloon et al., 2011;Exbrayat et al., 2013a). These impacts on R h can determine whether soil carbon stores increase or decrease for the same NPP, meaning they control whether the soil will remain a sink or convert to a source of CO 2 in the future. The various representations of R h are also responsible, at least in part, for the six-fold range in soil C achieved by CMIP5 simulations at the end of the 20th century in response to a three-fold range in NPP (Todd-Brown et al., 2013).
In this paper, we address two questions arising from the current parameterisation of the land carbon cycle. First, how do N and P limitations on plant productivity affect the response of soil C to different combinations of f W and f T over the 20th century? Second, through which mechanism is R h sensitive to the formulation of its response to changes in soil moisture and soil temperature? We explore these two sources of uncertainty in combination and quantify their influence on both equilibrated states and the response of the terrestrial component of a global ESM to the historical increase in atmospheric CO 2 and associated warming. Therefore, we examine at global and regional scales how the simulated historical carbon cycle is affected by the way soil moisture and soil temperature control R h using three formulations of f W , three formulations of f T , and N and P limitation.

Modelling system
We use the Carnegie-Ames-Stanford Approach with a Carbon-Nitrogen-Phosphorus (CASA-CNP) land biogeochemical model (Wang et al., 2010) coupled with the Table 1. Formulations of f W implemented in the CASA-CNP model (θ s : soil moisture; θ wilt : moisture at wilting point; θ fc : moisture at field capacity; θ opt : optimum moisture; θ lopt : lower optimal moisture -all expressed relative to moisture at saturation).

Function
Equation Community Atmosphere Biosphere Land Exchange (CA-BLE) land surface model . CASA-CNP simulates the turnover of terrestrial carbon based on three vegetation, three litter and three soil pools. Soil R h sums the CO 2 fluxes from the decomposition of litter and soil carbon. In each pool, R h is represented as a first-order process that depends on substrate availability, soil moisture and soil temperature, and these two latter terms are calculated in CABLE in response to meteorological forcing. CASA-CNP can be run in a carbon-only (C-only), C with nitrogen limitation (CN), and CN with phosphorus limitations (CNP) mode. Effectively, NPP is limited by the concentration of N (in CN mode) as well as P (in CNP mode) in leaves. The uptake of mineral N and labile P depends on their availability in soils, while mineralisation rates are tightly linked to C decomposition rates (Wang et al., 2010). We use parameter values for CASA-CNP that were previously reported by Wang et al. (2010). The CABLE + CASA-CNP terrestrial system has been coupled to the CSIRO Mk3L climate system model (Phipps et al., 2011;Zhang et al., 2011). The relatively coarse resolution of the model (5.6 • lat. 3.2 • long.) makes it a computationally efficient candidate of choice to create multiple simulations for sensitivity analyses, while simulated climate is still representative of the historical period (Phipps et al., 2011). Since we address the terrestrial carbon balance, our setup uses prescribed sea surface temperatures (SSTs) from the CSIRO Mk3.6 model (Rotstayn et al., 2012) using CMIP5 historical forcing data from 1850 to 2005 (Taylor et al., 2012) that were re-gridded at the resolution of Mk3L.

Model versions
To examine the uncertainty linked to the choice of biophysical response functions, three variations of f W and three variations of f T were implemented in the CASA-CNP model ( Fig. 1). These represent the key features of a larger suite of functions used in previous offline site-scale studies with the CABLE + CASA-CNP modelling system (Exbrayat et al., 2013a), and their exact formulation can be found in Table 1 and Table 2 for f W and f T , respectively. As shown in Fig. 1, the general consensus is that soil respiration is enhanced by intermediate moisture associated with warm temperatures. For example, while the bell-shaped f W used in CASA-CNP simulates a smooth response of soil respiration to drying conditions, SOILN (Jansson and Berg, 1985) and TRIFFID (Cox, 2001) predict a constantly null or lowmoisture adjustment below wilting point, respectively. Further, SOILN considers a whole range of optimal moisture conditions, while the two other formulations of f W both have a single, though different, optimal moisture. In saturated conditions, TRIFFID allows a higher respiration rate than the other f W . Comparing the temperature functions f T , CASA-CNP allows higher respiration rate for temperatures below +10 • C, while K1995 (Kirschbaum, 1995) is higher than the others between +10 • C and +40 • C. Finally PnET (Aber et al., 1997) displays the highest temperature-based adjustment of R h for soil temperatures above +40 • C. Interestingly, while CASA-CNP and PnET continue to increase, K1995 starts decreasing above +37 • C.

Experiments
Simulations were performed using the modelling system de-  Table 2. Formulations of f T implemented in the CASA-CNP model (T s : soil temperature in • C).
Function Equation ×12.64 −1 * Last terms in the equations are used to scale the original functions to the CASA-CNP model as explained by Exbrayat et al. (2013a).
Mk3L-CABLE-CASA-CNP simulation. Once equilibrated offline, total C storage for each of the 27 equilibration runs was used to reinitialise the coupled climate model and a further spin-up was undertaken until soil C storage achieved steady state. Finally, historical transient runs including increasing atmospheric CO 2 , based on CMIP5 specifications, and driven by corresponding CSIRO Mk3.6 SSTs were performed for each of the 27 model versions for 1850-2005. By prescribing atmospheric CO 2 we recognise that we limit the land-atmosphere coupling to energy and water exchanges between the land and the atmosphere. This is not a full coupling (e.g. Friedlingstein et al., 2006) where atmospheric carbon is also affected by terrestrial primary production and respiration fluxes. Further, we use a fixed land use map for each model version with dominant plant functional types set following the land cover data for 2005 from Hurtt et al. (2006). However, our experiments permit an assessment of how the choice of f W , f T and nutrients affect terrestrial systems and R h more simply than if we allowed for these feedbacks, or if we allowed ocean-atmosphere exchanges or land use change to affect atmospheric CO 2 . We note, of course, that these fluxes are included implicitly in the prescribed atmospheric CO 2 data. Finally, by using the same radiative forcing in all simulations, we isolate the effect of the different R h parameterisations on the terrestrial carbon cycle more simply than if variations in atmospheric CO 2 occurred in our simulations. For simplicity, we could have driven CASA-CNP with prescribed historical weather observations but using a climate model provides the opportunity to perform 21st century projections in future analyses.

Global land carbon balance
Equilibrated total soil carbon is presented in Table 3 for each model version. Globally, for the same combination of f T and f W , soil carbon in CN mode equilibrates at a level between 72 and 86 % of the C-only mode, while differences between CN and CNP are negligible. However, regardless of the biogeochemical mode adopted, differences in f T and f W introduce a difference of 4.5 times between the version that simulates the largest soil carbon pool (SOILN f W with PnET f T ) and the version that simulates the smallest (the original CASA-CNP version) in response to similar NPP during spin-up: 48.0 ± 0.1 Pg C a −1 , 46.2 ± 0.7 Pg C a −1 and 45.1 ± 0.7 Pg C a −1 in C-only, CN and CNP mode, respectively (±1 standard deviation).
The cumulative global net ecosystem accumulation (NEA) of terrestrial carbon since 1850 is shown in Fig. 2 (a positive accumulation corresponds to a net terrestrial sink). Each panel in Fig. 2 shows results for all nine combinations of f W and f T for a given C-only (Fig. 2a), CN (Fig. 2b) or CNP (Fig. 2c) mode using thin lines, and their style represents which f T was used in each model as indicated. The shaded area represents the total simulated range for a given nutrient limitation mode. All simulations show a net accumulation of carbon over the 20th century at the global scale as NPP increases on average due to a combination of CO 2 fertilisation and warmer temperatures driven by the observed CO 2 increases. However, there are major differences between the results from the C-only, CN and CNP modes and between the various moisture and temperature functions.
For the C-only mode, the range in the simulated NEA introduced by different f W and f T is very large, ranging from 207 to 438 Pg C (Fig. 2a). To illustrate the magnitude of the terrestrial sink, this represents ∼ 43 to ∼ 92 % of the ∼ 475 Pg C of accumulated emissions from fossil-fuel and land-use change represented in each panel of Fig. 2 from data by the Carbon Dioxide Information Analysis Center (CDIAC - Houghton, 2008;Boden et al., 2010). Adding N limitation reduces the terrestrial sink to between 61 and 175 Pg C for the 20th century, or 13 and 37 % of anthropogenic fossil-fuel emissions. This is more in accordance with Global Carbon Project estimates of uptake that vary around 30 % (Le Quéré et al., 2009). The range of NEA simulated in the CN mode, resulting from the choice of f W and f T , also decreases by about a factor of two relative to the C-only mode. The results from the CNP mode demonstrate a further reduction in both the magnitude and variability in NEA to 41-134 Pg C or ∼ 9 to ∼ 28 % of anthropogenic emissions. Note that while the uncertainty ranges of CN and CNP modes overlap, the ∼ 40 Pg C more than the highest member of the CN and CNP simulations (Fig. 2b, c). It is interesting to note that in each mode there is a great interaction between f T and f W as shown by the position of model versions using the same temperature function alternatively at the higher or lower end of the simulated range. We compare these simulations with previous estimates of global terrestrial NEA. Figure 3 compares our simulation results with estimates using the time period overlapping our simulations from the Global Carbon Project (1959-2005Canadell et al., 2007, data accessible at http://www. globalcarbonproject.org/carbontrends) and an intercomparison study of dynamic global vegetation models (1958-2002Sitch et al., 2008). Every simulation in the C-only mode, for all combinations of f T and f W , overestimates NEA as compared to previous studies. Generally, results for the mean terrestrial uptake from the CN and CNP models are more consistent with the estimates by Canadell et al. (2007) and Sitch et al. (2008). CN mode is actually the most similar to previous estimates with almost all combinations of f T and f W within the ∼ 0.75 Pg C a −1 uncertainty range provided by Sitch et al. (2008) and matching the estimate of Canadell et al. (2007). Over the same period all CNP simulations are below the range suggested by Sitch et al. (2008). The lower panel of Fig. 3 shows the variability of the land sink as illustrated by the standard deviation of annual NEA. C-only mode simulates excessive variability compared to Canadell et al. (2007). Both the CN and CNP modes are closer to Canadell et al. (2007) though slightly high.
Although N and P limitations reduce the absolute range in NEA, differences in f W and f T lead some model versions with equivalent N and P limitations to simulate twice as much NEA as other versions. This is true of C-only, CN and CNP modes (Figs. 2 and 3). To illustrate this, we analyse the annual NEA normalised by the annual NPP. We chose NPP because while it is affected by NP limitations, it is very similar between versions within the same nutrient mode and hence appears sensitive to the choice of f T and f W (Figs. S1, S2 and S3 in Supplement). The lower panel in Fig. 3 indicates that, by reducing all C fluxes and turnover processes, NP limitations reduce the inter-annual variability of land-atmosphere carbon fluxes. As a result, CN and CNP modes (Fig. 4b, c) do not exhibit the post-1960 step change in NEA that corresponds to a greater carbon sink in the C-only model (Fig. 4a) in response to the sudden increase in the growth rate of atmospheric CO 2 concentrations (Supplement Fig. S4). As a result, uncertainty ranges in Fig. 4b and c show that the ensembles of NP-limited simulations more often contain both net sources and net sinks during a same year as they overlap the dotted line that represents zero NEA. However, differences in f T and f W still introduce a two-fold uncertainty in cumulative historical NEA as shown in Fig. 2 even though the ensemble spread appears small relative to the effect of introducing NP limitations. To compare the spread generated by the different f W and f T relative to variability in NEA/NPP, we calculate a measure analogous to a signal-to-noise ratio for each C-only, CN and CNP mode. We define the "signal" as the temporal variability in NEA/NPP, calculated as the standard deviation of the annual mean NEA/NPP (the annual mean is represented by the black line in Fig. 4). The "noise" is calculated as the intra-annual variability between combinations of f W and f T . It corresponds to the standard deviation of the distance between all models and the mean NEA/NPP for all years (full ranges with maximum distances are in grey in Fig. 4). This signal-to-noise ratio decreases from 3.8 in C-only mode to 1.7 and 1.4 in CN and CNP modes, respectively. This indicates that the uncertainty due to f W and f T relative to the variability in NEA increases when NP limitations are added. From Fig. 4, it is clear that this results from the lack of signal in response to post-1960 increase in atmospheric CO 2 . We next investigate the regional implications of the choice of f W and f T to explain the four-fold range in soil carbon and two-fold range in NEA simulated between simulations with the same nutrient limitation.

Regional variations
We are aware that data by Canadell et al. (2007) and Sitch et al. (2008) are based on model simulations that do not integrate NP limitation. However, they integrate processes that were not represented in our modelling system for the purpose of simplifying the understanding of our results. Further, the good agreement of these previous studies makes us more confident on the reliability of CN and CNP simulations in terms of response to the historical increase in atmospheric CO 2 . Therefore, we will use CN simulations as reference in the remainder of the manuscript, as they were the most similar to these independent estimates of global NEA for 1959-2005 (Sect. 3.1).
The soil C density at equilibration for all CN simulations is shown in Fig. 4. Each panel in Fig. 5 represents a combination of a f W (rows) with a f T (columns). Large differences are observed in pool sizes as a function of f W and f T (similar patterns exist in C-only and CNP simulations). For example, the K1995 and PnET f T both equilibrate at much higher carbon density than CASA-CNP functions in the mid-to high latitudes in the Northern Hemisphere. PnET also has a higher soil C density in warmer regions. Differences implied by f W are more localised and do not seem to depend on a latitudinal gradient. SOILN equilibrates at a higher level of soil C in dry regions of south-west Australia, southern Africa and the   western edge of South America, while the two other f W provide relatively similar results when used with the same f T . Besides leading model versions to equilibrate at various levels of soil carbon both locally and globally, the different functions do not have the same sensitivity to a similar change in the soil temperature or moisture as shown by their respective shapes (Fig. 1). Figure 6 shows the change in the average value of the environmental scalar f W (θ s )×f T (T s ) between the 10 first and 10 last years of the CN simulations. There are large variations between the different model versions depending on both f W and f T . First, the decay rate (i.e. R h per unit of C s ) does not increase everywhere and there are significant decreases in dry regions (e.g. Arabian Peninsula, western Sahara and western Australia) with the SOILN f W used with the K1995 or PnET f T as it is the most constraining function in dry conditions (Fig. 1). All model versions using the K1995 f T have the highest relative increase in decay rate in northern Eurasia, while the other f T do not imply an increase of more than 20 % in R h except when used with SOILN.
Differences in the regional response of f T and f W are likely to generate regional differences in the land carbon balance. Therefore, we now investigate spatial variations in NEA to understand the roughly two-fold difference in global NEA simulated by all model versions in each nutrient mode (Fig. 2). Figure 7 shows NEA between the periods of 1996-2005 and 1850-1859 for each CN simulation. As it is significantly correlated (p < 0.001) with changes in soil carbon (Fig. S5), we can link NEA to how soil conditions affect the response of R h depending on the choice of f W and f T . Although regional differences appear depending on the choice of a f W when keeping the same f T , most of the uncertainty in NEA is related to the choice of a f T . a large region of negative NEA (i.e. a net source of CO 2 ) is simulated in the northern latitudes of eastern Eurasia. This negative NEA occurs irrespective of the choice of f W . The K1995 and PnET f T also simulate a negative NEA over high latitudes of North America which we attribute to soil warming that triggers higher R h that offset any increases in NPP, but only if the CASA-CNP and SOILN f W are used. The moisture functions nevertheless lead to large regional differences, especially over central Europe. Figure 7 shows a gradual increase in NEA from the CASA-CNP f W , through the SOILN f W to the TRIFFID f W . This additional sink is most apparent when the PnET f T is used as TRIFFID accumulates soil carbon in places where the CASA-CNP and SOILN f W lose carbon (e.g. Americas). Although we see regional differences in NEA depending on f T and f W , we note that places where the environmental scalar and hence the decay rate (or R h per unit of C s ) increased the most (red colours on Fig. 6) are not where most soil C is lost (blue colours on Fig. 7 and S5), except perhaps when using the SOILN f W with the K1995 f T . Based on linear regressions (not shown) the change in environmental scalar values does not explain more than 12 % of the variability in the change of soil carbon and NEA. The spatial influence of N limitation on the C-only simulations can be examined by comparing total NEA between C-only and CN simulations (Fig. 8) as well as differences in soil carbon changes (Fig. S6). The most obvious differences are in the mid-and high latitudes of the Northern Hemisphere (in red on Fig. 8 and S2). There, C-only simulations store up to 5 kg C m −2 (Fig. 8) more than CN simulations as a result of no nutrient limitation applied to NPP (also Fig. 2). Differences of the same magnitude are observed in soil carbon changes (Fig. S6). In large regions across the Northern Hemisphere, the difference is large enough to change the sign of NEA over the historical period from net sources in CN simulations to net sinks in C-only simulations (stippling on Fig. 6), especially with the K1995 f T . Elsewhere, differences in NEA and soil carbon change between the C-only and CN simulations are commonly small (< 1 kg C m −2 ). Local responses are, however, dependent on the choice of both f W and f T . For example, for each f T , there is a relatively large variation in soil carbon accumulation between the CASA-CNP f W and the TRIFFID f W (Supplement Fig. S6).
In comparison with the effect of including N limitation, the consequences of including P limitation are minor ( Fig. S7  and S8). In terms of NEA, there is a further reduction in comparison to the CN-limited simulations but the magnitude of the reduction is at most ∼ 1 kg C m −2 and only a few places see a change in the sign of NEA between CNP and CN simulations. However, we note that some regions have higher NEA in CNP simulations than in CN simulations despite the additional P limitation on NPP. Similarly, the effect on soil carbon from the addition of P limitation is generally small in comparison to the addition of N limitation and broad differences between the temperature and moisture functions are even smaller than for NEA (Supplement Fig. S8). Overall tiny regional differences between CN and CNP simulations sum up to a similar global signal.

Equilibration of soil carbon
Using different formulations of f T and f W leads to significant differences in the outcome of the spin-up procedure (Fig. 5). In our simulations, all model versions were brought to equilibrium until C pools achieved a steady state. This is a standard procedure (e.g. Wang et al., 2010;Xia et al., 2012) that would most likely have been used in all CMIP5 simulations that incorporated carbon. Spinning up a model means integrating the model with steady boundary conditions until the trend in carbon pool is negligible or R h ≈ NPP. According to Eq. (1), in ESMs, the amount of decomposition, and therefore R h , is controlled by a time-invariant reference k parameter, the f W (θ s ) × f T (T s )product and the amount C s of carbon available in soil for decomposition. Model equili-bration consists of achieving the carbon pool size needed to simulate R h at a level that compensates for NPP, while integrating the model under steady boundary conditions. Given that NPP is similar between simulations within the same nutrient mode (Supplement Figs. S1, S2 and S3), it is only our modifications to f W and f T that have led to total soil C in our CN model to range from 765 Pg C to 3495 Pg C. As shown in Fig. 5, differences in total soil carbon at equilibrium are due to large regional differences, especially at high latitudes. This can be explained by the relative position of these functions for cold temperatures (Fig. 1): the CASA-CNP f T is systematically above the two other f T for soil temperatures below 10 • C. It therefore requires less substrate to simulate R h at a level that compensates for the same NPP than K1995 and PnET. Conversely, PnET causes the model to equilibrate at a higher soil C density in warmer regions as it is well below the two other functions for soil temperature corresponding to Africa and South America. As there are dry and wet regions at any latitude, differences implied by f W are more localised. However, the same relationship between the relative positions of the curves can be seen. The most noticeable feature is that SOILN, a very limiting f W in dry conditions, equilibrates at a higher level than the two other f W in south-west Australia, southern Africa and the western edge of South America, where it requires more substrate to achieve the same R h to compensate NPP.
In Fig. 5 we adopted a colour scale similar to Fig. 3 of Todd-Brown et al. (2013) to present soil C in different CMIP5 ESMs. The regional differences implied by the different f W and f T map particularly well onto the diversity shown by the CMIP5 models and total soil carbon approximates the six-fold range found in CMIP5 models (Todd-Brown et al., 2013), although these include models well outside of the 95 % confidence interval of 890-1660 Pg C in current soil (Todd-Brown et al., 2013). For example, five of our CN simulations do not fall within the 890 to 1660 Pg C 95 % confidence interval of present soil carbon reported by Todd-Brown  et al. (2013). Dismissing them leads our projected range of cumulative historical NEA to shrink by about a third from 61-175 Pg C to 86-161 Pg C. We do not explore this more in detail here but we suspect that these similarities between our simulations and CMIP5 results nevertheless strongly indicate that the formulation of the time and space invariant f W and f T is a key source of uncertainty in the equilibrium of these models. We note, however, that CMIP5 models also vary three-fold in their NPP, with much lower values in N-limited models, as well as in the number of pools they employ. They are also likely to use different values of k as shown by Todd-Brown et al. (2013) with their reduced complexity models. Equation (1) highlights that in the first-order parameterisation of microbial decomposition, substrate availability C s is also a regulating factor. Therefore, we explore hereafter the influence of initial conditions on the soil carbon response, especially since the change in the environmental scalar through time cannot explain NEA.

Response to transient conditions
Despite variations in total soil carbon content, all model versions simulate a positive NEA linked with increasing CO 2 emissions (Fig. 2) over the 20th century irrespective of nutrient limitation mode or the choice of f W or f T . This means that each model version simulates the terrestrial biosphere as a net carbon sink associated with the fertilisation effect of increasing atmospheric CO 2 and temperature over the 20th century. This has been described previously in CASA-CNP and other coupled models whether or not they included NP interactions with C (Cox et al., 2000;Friedlingstein et al., 2006;Sitch et al., 2008;Zhang et al., 2011). The similarity of NPP for simulations with the same nutrient limitation indicate that there is only poor interaction between response functions and mineralisation of nutrients for these historical simulations, although this may be a non-negligible process for NPP in a warmer future (Zaehle et al., 2010a).
However, comparing our range of results due to different formulations of f W and f T across the three nutrient modes with estimates of NEA from Canadell et al. (2007) and Sitch et al. (2008) suggests that the C-only mode simulates an excessively high terrestrial uptake. In general, the CN mode appears most consistent with other estimates. Earlier studies, each using a single f W or f T , predicted a reduction in the CO 2 fertilisation effect by up to 72 % when considering N limitation (Thornton et al., 2007;Sokolov et al., 2008;Jain et al., 2009;Zaehle et al., 2010b;Bonan and Levis, 2010;Zhang et al., 2011), and our average reduction of 64 % in global NEA (Fig. 2) due to N limitations is consistent with these previous estimates. Our results suggest that reductions in CO 2 fertilisation simulated by CASA-CNP in response to N and P limitations (Zhang et al., 2011) are robust to the choice of f T or f W .
Of course, it would be straightforward to calibrate the model in C-only mode to reduce the overestimation of NEA, but then the CN and CNP modes would grossly underestimate observations. Further, parameter values optimised to reproduce observed data would likely compensate for the lack of representation of key biogeochemical processes (N and P), introducing a high risk of obtaining acceptable simulations for the wrong reasons. It has been demonstrated that over-fitted parameters that provided acceptable calibration results were not able to capture the response of a system to changes if some processes were missing in the model structure (e.g. Exbrayat et al., 2013b). Since the availability of N and P has a key influence on NPP that supplies substrate for decomposition and affects NEA (Vitousek and Howarth, 1991;Luo et al., 2004;Vitousek et al., 2010;Goll et al., 2012), results from CN and CNP are likely more robust than C-only. Rather than calibrating a C-only version, adding N and P to more correctly reflect the response of the biogeochemical system to increase in atmospheric CO 2 and temperature is preferable.
The two-fold range in NEA introduced by f W and f T remains under all nutrient limitation modes, and this can be attributed to large regional differences in the Northern Hemisphere (Fig. 7). There, NEA changes from a carbon sink to a carbon source in response to very similar NPP within the same nutrient mode. Further, differences in NEA and soil carbon change cannot be explained by the change imposed on the environmental scalar ( Fig. 6 and Sect. 3.2). Therefore, only the amount of available substrate can explain these differences in R h . Figure 9 presents the change in soil C as a function of initial conditions (i.e. in 1850 after spin-up) in all simulations. Grid boxes with low initial values can gain or lose C for all combinations of a f W with a f T . The Conly simulations accumulate soil C almost everywhere because the lack of nutrient limitation allows for a stronger response of NPP which dominates the response and offsets R h . In contrast, in the nutrient-limited CN and CNP, soil carbon losses are observed where substrate availability is initially high. This is particularly true for K1995 and PnET f T (Fig. 9, central and rightmost columns). In these simulations where NPP is limited by nutrient availability, even a small relative increase in the environmental scalar f W (θ s ) × f T (T s ) applied to initially large pools enhances R h enough to transform some local carbon sinks into sources. Further, CN simulations generally equilibrate at higher soil carbon content than CNP simulations (Table 3). Hence, in regions where large initial soil carbon pools trigger losses during simulations, CN runs lose a bit more carbon than the corresponding CNP runs. Nevertheless, it is worth noting that CN and CNP simulations are quite comparable under historical forcing despite the added complexity of P limitation. This characteristic may still involve heterogeneous responses of CN and CNP modes to future climate change, something that remains to be explored beyond the work presented here (see e.g. Zhang et al., 2011;Goll et al., 2012).
We see here an analogy with the model-specific nature of soil moisture described by Koster et al. (2009). That is, the amount of soil C as simulated in ESMs is not something that can be directly compared with a quantity that might be measured in the field. Rather, soil C in each model is the value required by the model to reach steady state, and through which variations trigger an acceptable response of land-atmosphere exchanges to historical changes. Since observations are not available to constrain the model in the future, a large uncertainty arises that can be exemplified by the lack of consensus between model projections in a previous inter-comparison project despite a rather good agreement for historical simulations (Friedlingstein et al., 2006). This adds to recently stated concerns that the current parameterisation of decomposition is not representative of our understanding of this process (Allison et al., 2010;Schmidt et al., 2011;Todd-Brown et al., 2012). Our results point to a critical need to refine the initialisation of ESMs by spin-up as it controls the sign of change in soil C and NEA, and hence the carbon-climate feedback from the land on the atmosphere. This is especially true when the dominant CO 2 -fertilisation effect is reduced by more realistic nutrient limitation. Although we do not explore this here, we suggest that gridded data sets of current soil carbon content such as those presented by Tarnocai et al. (2009) for high latitudes, or Todd-Brown et al. (2013 globally, could possibly be used as guidelines to constrain soil carbon. Since methods now exist that greatly speed up the spin-up procedure (Xia et al., 2012) -the most expensive part of such simulations -trial and error procedures are feasible.

Conclusions
We have used 27 combinations of f T , f W and nutrient limitations in an ESM to explore how the land carbon balance responds to changing atmospheric CO 2 over the period . Various formulations of f T and f W generate a range of equilibrated soil carbon stores very similar to the six-fold range of global soil C achieved by CMIP5 models regardless of whether nutrient limitation is implemented. That is, the range in soil carbon in CMIP5 is likely the result of equilibration methods.
Implementing N and P limitations on plant productivity in the CASA-CNP ecosystem model better constrains the simulation of the historical response of the terrestrial C cycle irrespective of the f T or f W used because of the lack of response to post-1960 rapid increase in atmospheric CO 2 . However, in these simulations the initial carbon pool size is the main driver of the response of soil carbon to global warming. As the magnitude of the available substrate controls the sensitivity of R h to changes in temperature and moisture, larger pools are more likely to deplete under global warming. Due to the size of soil carbon pools even small changes in forcing can lead to large C losses and drive the whole land C balance response to warming. Therefore, the wide range of responses in CMIP5 in terms of soil carbon may well be an artefact of the initialisation procedure used.
Based on our experiments, we recommend representing at least CN interactions in ESMs in order to capture the correct magnitude of historical land-atmosphere carbon fluxes and the response of the system to increasing atmospheric CO 2 . The other clear implication of our results is that a more concerted effort in how microbial decomposition processes are represented in ESMs is required. We need to address how equilibrium should be defined or constrained to match some estimates, how nutrients should be represented and how we develop these efforts with limited global databases of soil carbon.