Analytically tractable climate-carbon cycle feedbacks under 21st century anthropogenic forcing

. Changes to climate–carbon cycle feedbacks may signiﬁcantly affect the Earth system’s response to greenhouse gas emissions. These feedbacks are usually analysed from numerical output of complex and arguably opaque Earth system models. Here, we construct a stylised global climate–carbon cycle model, test its output against comprehensive Earth system models, and investigate the strengths of its climate–carbon cycle feedbacks analytically. The analytical expressions we obtain aid understanding of carbon cycle feedbacks and the operation of the carbon cycle. Speciﬁc results include that different feedback formalisms measure fundamentally the same climate–carbon cycle processes; temperature dependence of the solubility pump, biological pump, and

. Climate-carbon cycle feedbacks and state variables as represented in the stylised model introduced in this paper. Carbon stored on land in vegetation and soils is aggregated into a single stock c t . Ocean mixed layer carbon, c m , is the only explicitly modelled ocean stock of carbon; though to estimate carbon cycle feedbacks we also calculate total ocean carbon (Eq. 7).
bon uptakes even under identical atmospheric concentration or emission scenarios .
Here, we develop a stylised model of the global carbon cycle and its role in the climate system to explore the potential weakening of carbon cycle feedbacks on policy-relevant timescales (< 100 years) up to the year 2100. Whereas comprehensive Earth system models are generally used for projections of climate, models of the Earth system of low complexity are useful for improving mechanistic understanding of Earth system processes and for enabling learning (Randers et al., 2016;Raupach, 2013). Compared to comprehensive Earth system models, our model has far fewer parameters, can be computed much more rapidly, can be more rapidly understood by both researchers and policymakers, and is even sufficiently simple that analytical results about feedback strengths can be derived. Compared to previous stylised models (Gregory et al., 2009;Joos et al., 1996;Meinshausen et al., 2011a, c;Gasser et al., 2017a), our model features simple mechanistic representations, as opposed to parametric fits to Earth system model output, of aggregated carbon uptake both on land and in the ocean. Our stylised and mechanistically based climate-carbon cycle model also offers a workbench for investigating the influence of mechanisms that are at present too uncertain, poorly defined, or computationally intensive to include in current Earth system models. Such stylised models are valuable for exploring the uncertain but potentially highly impactful Earth system dynamics such as interactions between climatic and social tipping elements (Lenton et al., 2008;Kriegler et al., 2009;Schellnhuber et al., 2016) and the planetary boundaries (Rockström et al., 2009;Steffen et al., 2015).
Analyses of climate-carbon cycle feedbacks conventionally distinguish four different feedbacks ( Fig. 1) Ciais et al., 2013). (i) In the land concentrationcarbon feedback, higher atmospheric carbon concentration generally leads to increased carbon uptake due to the fertilisation effect, in which increased CO 2 stimulates primary productivity. (ii) In the ocean concentration-carbon feedback, physical, chemical, and biological processes interact to sink carbon. Atmospheric CO 2 dissolves and dissociates in the upper layer of the ocean to then be transported deeper by physical and biological processes. The concentration-carbon feedbacks are generally negative, dampening the effects of anthropogenic emissions. (iii) In the land climate-carbon feedback, higher temperatures, along with other associated changes in climate, generally lead to decreased storage on land at the global scale, for example due to the increase in respiration rates with temperature. (iv) In the ocean climatecarbon feedback, higher temperatures generally lead to reduced carbon uptake by the ocean, for example due to decreasing solubility of CO 2 . The climate-carbon feedbacks are generally positive, amplifying the effects of carbon emissions.
We begin by introducing our stylised carbon cycle model and testing its output against historical observations and future projections of Earth system models. Having thus established the model's performance, we introduce different formalisms used to quantify climate-carbon cycle feedbacks and describe how they can be computed both numerically and analytically from the model. We use our results to analytically study the relative strengths of different climate-carbon cycle feedbacks and how they may change in the future, as well as to compare different feedback formalisms. We conclude by speculating on how this stylised model could be used as a "workbench" for studying a range of complex Earth system processes, especially those related to the biosphere.

Model formulation
There is well-developed literature on stylised models used for gaining a deeper understanding of Earth system dynamics and even for successfully emulating the outputs of comprehensive coupled atmosphere-ocean and carbon cycle models (Anderies et al., 2013;Gregory et al., 2009;Joos et al., 1996;Meinshausen et al., 2011a, c;Gasser et al., 2017a). We developed a combination of existing models and new formulations to construct a global climate-carbon cycle model with the following characteristics. 3. All parameters have a direct biophysical or biogeochemical interpretation, although these parameters may be at an aggregated scale (for example, a parameter for the net global fertilisation effect, rather than leaf physiological parameters). We avoid models or model components constructed by purely parametric fits, such as impulse response functions, to historical data or projections of Earth system models (Kamiuto, 1994;Gasser et al., 2017b;Joos et al., 1996;Harman et al., 2011;Gregory et al., 2009;Meinshausen et al., 2011a). 4. The model is sufficiently simple that calculation of the model's feedback strengths is readily analytically tractable. This tractability may come at the expense of complexity, for example multiple terrestrial carbon compartments, or accuracy at millennial or longer timescales (Lenton, 2000;Randers et al., 2016).
Building on the work of Anderies et al. (2013), we constructed a simple model with globally aggregated stocks of atmospheric carbon in the form of carbon dioxide, c a ; terrestrial carbon, including vegetation and soil carbon, c t ; and dissolved inorganic carbon (DIC) in the ocean mixed layer, c m . The model's fourth state variable is global mean surface temperature relative to pre-industrial, T = T − T 0 . Compared to Anderies et al. (2013), our model includes more realistic representation of terrestrial and ocean processes but without increase in model complexity, as well as time lags for climate response to CO 2 .
We now describe the dynamics of the land carbon stock, the ocean carbon stock, and atmospheric carbon and temperature in our model.

Land
Net primary production (NPP) is the net uptake of carbon from the atmosphere by plants through photosynthesis. NPP is expected to increase with concentration of atmospheric carbon dioxide c a . A simple parameterisation of this socalled fertilisation effect is "Keeling's formula" for global NPP (Bacastow et al., 1973;Alexandrov et al., 2003): Throughout this article, the subscript "0" denotes the value of the quantity at a pre-industrial equilibrium, and "log" denotes natural logarithm. Keeling's formula incorporates all climate-change-related effects on global NPP occurring simultaneously with carbon dioxide changes, for example, precipitation and temperature effects, in addition to fertilisation effects. The curvature of the log function represents limitations to NPP such as changing carbon-use efficiency (Körner, 2003) or nutrient limitations (Zaehle et al., 2010). Constant climate sensitivity is also a key assumption, otherwise the relative weight of climate and CO 2 effects on NPP would change. At the same time, carbon loss from the world's soils through respiration, R, is expected to increase at higher global mean surface temperature, T . We approximate the net temperature response of global soil respiration using the Q 10 formalism R( T ) = R 0 Q T /10 R c t /c t0 (Xu and Shang, 2016), where Q R is the proportional increase in respiration for a 10 K temperature increase. We assume that preindustrial soil respiration is balanced by pre-industrial NPP, R 0 = NPP 0 . To avoid introducing multiple pools of carbon into the model, we also have to assume that global soil respiration is proportional to total land carbon (rather than soil carbon). Respiration in our model also implicitly includes other carbon-emitting processes such as wildfires or insect disturbances.
It follows that the change in global terrestrial carbon storage is In this expression we have also included loss of terrestrial carbon due to land use emissions LUC(t). We rearrange this expression to give where the terrestrial carbon carrying capacity is For model simplicity, we do not explicitly model factors affecting terrestrial carbon uptake such as seasonality, species interactions, species functionality, migration, and regional variability.

Ocean
In the upper-ocean mixed layer, mixing processes allow exchange of carbon dioxide with the atmosphere. The solubility and biological pumps then transport carbon from the mixed layer into the deep ocean. Since the residence time of deep ocean carbon is several centuries, we explicitly only model the dynamics of upper-ocean carbon while the deep ocean is treated merely as an extremely large carbon reservoir. We include the effects of ocean carbon chemistry, the solubility and biological pumps, and ocean-atmosphere diffusion on upper-ocean mixed layer carbon.
Ocean uptake of carbon dioxide from the atmosphere is chemically buffered by other species of DIC such as HCO − 3 and CO 2− 3 , which are produced when dissolved CO 2 reacts with water. The reaction of CO 2 with water, producing these other species, reduces the partial pressure of CO 2 in water allowing for more ocean CO 2 uptake before equilibrium with the atmosphere is achieved. The Revelle factor, r, is defined as the ratio of the proportional change in carbon dioxide content to the proportional change in total DIC (Sabine et al., 2004;Goodwin et al., 2007). For simplicity, we assume a constant Revelle factor, except for the temperature dependence, D T , of the solubility of CO 2 in seawater. Therefore CO 2 diffuses between the atmosphere and ocean mixed layer at a rate proportional to where since at pre-industrial equilibrium p(c m0 , 0) = c a0 . There are two main mechanisms by which carbon is transported out of the upper-ocean mixed layer into the deep ocean stocks: the solubility and biological pumps. In the solubility pump, overturning circulations exchange mixed layer and deep ocean water. We assume that the large size of the deep ocean means its carbon concentrations are negligibly changed over the 100-year timescales relevant for the model. The net transport of carbon to the lower ocean by the solubility pump can therefore be represented by where w 0 is the (proportional) rate at which mixed layer ocean water is exchanged with the deep ocean and w T parameterises weakening of the overturning circulation that is expected to occur with future climate change .
The biological pump refers to the sinking of biomass and organic carbon produced in the upper ocean to deeper ocean layers (Volk and Hoffert, 1985). In the models on which the IPCC reports are based, a weakening of the biological pump is predicted under climate change, mostly due to a decrease in primary production, in turn due to increases in thermal stratification of ocean waters . We represent this climate-induced weakening in a single approximately linear factor, so that the rate of carbon transported out of the upper-ocean mixed layer by the biological pump to lower deep sea layers is given by As on land, we assume a pre-industrial equilibrium at which the biological pump was balanced by transport of carbon back to the mixed layer by ocean circulation. We neglect deposition of organic carbon to the sea floor and the long timescale variations in the biological pump that may have contributed to glacial-interglacial cycles (Sigman and Boyle, 2000). We therefore add an additional term B( T ) − B(0) to the transport of carbon from the ocean mixed layer to the deep ocean. Organic carbon that does not sink to the deep ocean is rapidly respired back to forms of inorganic carbon; the ocean mixed layer stock of organic carbon is therefore small, around 3 Pg C (Ciais et al., 2013), and we do not count it in the model's carbon balance.
By combining the expressions for the solubility and biological pumps with ocean-atmosphere carbon dioxide diffusion, we obtain the rate of change of ocean mixed layer DIC, c m : The coefficient of the first term was chosen such that 1/D is the timescale on which carbon dioxide diffuses between the atmosphere and the ocean mixed layer (that is, the derivative of the first term with respect to c m , evaluated at the preindustrial equilibrium, is D).
The carbon content of the deep ocean does not explicitly enter Eq. (6). To evaluate ocean carbon feedbacks, however, we require the change in total ocean carbon content c M compared to pre-industrial conditions. We calculate this as ocean mixed layer carbon plus carbon transported to the deep ocean by the solubility and biological pumps: We do not explicitly model factors such as the thickness of ocean stratification layers, spatial variation of stratification, nutrient limitations to NPP, or changes in ocean circulation due to wind forcing, freshwater forcing, or sea-ice processes (Bernardello et al., 2014).

Atmosphere
We define c s to be the total carbon in our "system", comprised of carbon stocks in the ocean mixed layer, atmosphere, and terrestrial biosphere, that is The only processes that affect the total carbon are human emissions of fossil carbon into the atmosphere, e(t), and export of carbon into the deep ocean by the solubility and biological pumps, giving in which the initial value of c s is c a0 + c t0 + c m0 . To obtain the dynamics of atmosphere carbon stocks, we therefore solve the differential Eq. (9) and then use the carbon balance Eq. (8) to find c a . Increasing atmospheric carbon dioxide levels c a cause a change in global mean surface temperature, T , compared to its pre-industrial level. To model the response of T , we follow the formulation of Kellie-Smith and Cox (2011), which includes a logarithmic response as per the Arrhenius law and a delay of timescale τ . Physically, this time delay is primarily due to the heat capacity of the ocean.
The climate sensitivity λ specifies the increase in temperature in response to a doubling of atmospheric carbon dioxide levels. The climate sensitivity accounts for energy balance feedbacks such as from clouds and albedo. We use the transient climate sensitivity  as this specifies the response of the climate system over an approximately 100year timescale (see Sect. 3).

Model parameterisation and validation
Our climate-carbon cycle model has 12 parameters, four state variables, and three nontrivial initial conditions (by def-inition, the initial value of T is 0). We choose to parameterise each process with the best available knowledge about that process, rather than try to force the model to fit historical data. This is in line with our stated model purposes of understanding and learning, rather than prediction. Parameters for the response of climate to carbon dioxide (λ, τ ) and two parameters of the response of the ocean to changing temperature (B T and w T ) were set based on the output of atmosphere-ocean global circulation models. For the climate sensitivity λ, transient climate response was used. All other parameters are based on historical observations of the global carbon cycle (Table 1). Unless otherwise noted, we perform emissions-based model runs using harmonised historical data and future RCP scenarios on fossil fuel emissions [e(t)] and land use emissions [LUC(t)] (Meinshausen et al., 2011b). While the focus of our study is on future climate change, from the present day until 2100, we begin simulations in 1750 to compare our model against historical observations. Time series of the model output are displayed in Fig. 2. Model solutions were approximated in continuous time.

Feedback analysis
Our climate-carbon cycle model is sufficiently simple that the strengths of its feedbacks can be estimated analytically. Such computations are useful since the resulting symbolic expressions can be used to identify how parameters of interest affect feedback strengths and model dynamics. In this section we introduce definitions of feedback strengths, cal-  (Gregory et al., 2015). This result is consistent with simulations that indicate that maximum warming after a CO 2 pulse is reached after only a decade (Ricke and Caldeira, 2014) and with results from impulse response model experiments . Atmosphere-ocean mixed layer CO 2 equilibration rate D  Alexandrov et al. (2003) found that values between 0.3 and 0.4 are compatible with results from a process-based global NPP model. culate climate-carbon cycle feedbacks analytically and numerically, and estimate feedback non-linearities.

Definitions
There are multiple measures of carbon cycle feedbacks currently in use. We here review three of the most common measures.
Consider an emission of E Pg C over some time period to the atmosphere. In the absence of carbon cycle feedbacks, the atmospheric carbon content would increase by c off a ≡ E. With a feedback switched on, the atmospheric carbon content would actually change by c on a . The feedback factor is (Zickfeld et al., 2011) Out of the total atmospheric carbon change c on a , the carbon cycle feedback contributes (Hansen et al., 1984) Gain is the change in a feedback to atmospheric carbon content caused by changes in atmospheric carbon content: Gain and feedback factor are related by An alternative formalism, introduced by Friedlingstein et al. (2006), allows feedbacks to be characterised from carbon cycle model output. Climate models are not required, except as a forcing to the carbon cycle model. The formalism relates the changes in terrestrial and marine carbon stocks to changes in global mean temperature and atmospheric carbon dioxide as follows: Here the β L and β O feedback parameters are the land and ocean, respectively, carbon sensitivities to atmospheric carbon dioxide changes c a . Likewise, γ L and γ O are the land and ocean, respectively, carbon sensitivities to temperature changes T . Note that c M denotes the total marine carbon stock, both mixed layer and deep ocean. The differences c a , etc., are usually calculated over the duration of a simulation.
To isolate the β and γ feedback parameters, simulations are conducted with biogeochemical coupling only and with radiative coupling only (Gregory et al., 2009). In both the formalisms introduced thus far, the feedback measures are calculated by examining the changes in carbon stocks at the end point of model simulations. In contrast, Boer and Arora (2009) estimate sensitivities and B of the instantaneous carbon fluxes from atmosphere to land and ocean: These feedback parameters B and are usually computed for all time points during a simulation, again using biogeochemically coupled and radiatively coupled simulations. The two sets of parameters (B, ) and (β, γ ) are related by Accordingly, Boer and Arora (2013) refer to B and as direct feedback parameters and to β and γ as time-integrated feedback parameters.

Analytical feedback strengths based on equilibrium changes
Analytical approximations to the strengths of carbon cycle feedbacks in our model require choosing a timescale on which the feedbacks will be calculated. Numerically estimated feedback factors (Eq. 11) and time-integrated feedback parameters (Eqs. 15-16) are conventionally calculated using carbon stock changes over 100 years or more. Responses on the longest timescales of our model are therefore most relevant if our analytical approximations are to approximate numerically calculated values. While recognising that the Earth's climate system is presently far from equilibrium, we use changes in the equilibrium state of the model to approximate model responses over long timescales. We analytically calculate the gains associated with each of the feedback loops in Fig. 1 as follows. We calculate the sensitivity (mathematically, partial derivative) of the equilibrium value of each quantity in the feedback loop with respect to the preceding quantity in the loop. We form the product of the derivatives (as per the chain rule of differentiation) to estimate the gain of that feedback loop. For example, to calculate the land climate-carbon gain we calculate the sensitivity of equilibrium temperature with respect to changes in atmospheric carbon content (∂T * /∂c a ), multiplied by the sensitivity of equilibrium terrestrial carbon with respect to changes in temperature (∂c * t /∂T ), multiplied by the sensitivity of equilibrium atmospheric carbon with respect to changes in terrestrial carbon (∂c * a /∂c t ). Land climate-carbon equilibrium gain. The subscript T denotes that the feedback involves temperature. Asterisks ( * ) denote equilibrium quantities. From these gains, the feedback factors F * TL , F * L , F * TO , and F * O can be calculated using Eq. (14). We label these gain and feedback factors g * and F * , respectively, to denote they are based on an equilibrium approximation, not directly from transient simulations as estimated by Zickfeld et al. (2011).
The derivatives of c * a are trivial to calculate: by carbon balance, ∂c * a ∂c t = ∂c * a ∂c M = −1. To calculate the derivatives of c * T , we set 0 = dc t dt , solve for c t , and calculate the necessary derivatives. A similar procedure provides ∂T * ∂c a . The remaining derivatives are ∂c M ∂T and ∂c M ∂c a . Carbon sunk into the deep ocean is substantial and cannot be neglected. Deep ocean carbon storage equilibrates on timescales of millennia or more, however, far longer than the timescales of interest in this model (we therefore write derivatives of c M rather than c * M ). We therefore cannot use the same equilibrium approach as for the other variables. Instead, we derive approximations to Eq. (7) as follows. First, we observe that in the SRES A2 scenario used below, both c m (t) and T (t) can be approximated as linear increases, starting at c m = c m0 and T = 0, respectively, over a time interval t lin . We estimate this time interval by t lin = (c m (t end ) − c m0 )/c m (t end ) us-  (Ciais et al., 2013) and temperatures in 2012 relative to temperatures in 1880 (Hartmann et al., 2013). Predicted future changes are carbon stocks in 2100 compared to 2012  and global mean surface temperature (GMST) averaged over 2081-2100 relative to 1986-2005 , under the range of RCP scenarios.

Ocean carbon changes (Pg C)
Land carbon changes (Pg C) GMST change, ing the value c m and gradient c m at the end of the simulation period. We obtain We use this equation to calculate the derivatives ∂c M ∂T and ∂c M ∂c a . Evaluating these derivatives will involve the derivatives ∂c m ∂T and ∂c m ∂c a . Since partial pressures across the air-sea interface equilibrate rapidly on the timescale of the model (D = 1 yr −1 , Table 1), we assume that c a ≈ p(c m , T ), rearrange for c m , and then calculate the appropriate derivatives from the resulting equation.
We analytically estimate equilibrium versions of the timeintegrated feedback parameters of Friedlingstein et al. (2006) using a similar approach: Since the ocean component of the model has multiple processes that respond to temperature, some analytical forms were too complicated for easy visual inspection (Table A1). We derived approximate analytical feedbacks by comparing the magnitudes of terms in the numerator and denominator of the feedback measures by expanding in power series of D T T and c a /c a0 .

Analytical feedback strengths based on carbon fluxes
We estimate the direct feedback parameters as follows: * Here dc t /dt and dc M /dt denote the atmosphere-land and atmosphere-ocean fluxes. The subscript T = 0 denotes a biogeochemically coupled (and radiatively decoupled) simulation and c a = c a0 denotes a radiatively coupled (and biogeochemically decoupled) simulation. The values of the feedback parameters are strongly scenario dependent . To calculate the direct feedback parameters, we assume a standard CO 2quadrupling concentration pathway in order to compare our results with Arora et al. (2013). This scenario has c a (t) = c a0 a t , where a = 1.01. In this scenario, 1 c a dc a dt = log a and, ignoring an initial exponential transient, dT dt = λ log a/ log 2. For the atmosphere-land carbon flux, the calculation is straightforward under the following assumptions. We assume that NPP 0 /c t0 log a so that c t tracks its carrying capacity c t ≈ K (Eq. 2). We also ignore land use change, so that dc t dt ≈ dK dt . Then we calculate dK dt | c a =c a0 = ∂K ∂T dT dt and dK dt | T =0 = ∂K ∂c a dc a dt . While the atmosphere-ocean flux could be read off directly from the first term of Eq. (6), this form is however not particularly useful. As it involves a small difference between two large quantities, c a and p(c m , T ), the size of the difference can only be estimated from numerical results and gives no immediate insight into how it depends on parameters. Furthermore, we seek to compare our analytical results to the results presented by Arora et al. (2013), in which the feedback parameters are presented as functions of c a or T only (not c m ).
We instead derive an approximation based on timescale separation as follows. The characteristic timescale of atmosphere-ocean diffusion is much faster than the solubility pump, biological pump, or human emissions into the atmosphere (D w 0 , B 0 /c m0 , log a). Since atmosphere-ocean Table 3. Feedback analysis. Gains (g), feedback factors (F ), time-integrated feedback parameters (γ and β), and direct feedback parameters ( and B) were calculated analytically and numerically. Analytical ocean feedbacks are approximations of the exact forms in Table A1 (see Sect. 4.2). Exact forms were used to calculate numerical values. In this table, p ≡ p(c m , T ). Units of the climate-carbon integrated feedback parameters are Pg C K −1 and concentration-carbon integrated feedback parameters are Pg C ppm −1 . Ranges for analytical results are written in the form (value at start of simulation) to (value at end of simulation). Emissions scenarios are as indicated; land use emissions were assumed to be zero. From the results of simulations using the SRES A2 scenario we use t lin ≈ 100 corresponding to a period between the years 2000 and 2100.
w 0 c m rc a -estimate from analytical form Fig. A1a Fig. A1b  Table A1. Taylor series expansions and L'Hôpital's rule were then used to derive the approximate forms in Table 3.

Numerical estimation of feedback strengths
In addition to the analytical approximations to carbon cycle feedbacks derived from our model, we also estimate feed-back factors from direct numerical simulations of our model. To compare the results of our model to previous studies, we use the following scenarios. To compare our results with the time-integrated feedback parameters reported by Friedlingstein et al. (2006) and the feedback factors and gains of Zickfeld et al. (2011), we employ the SRES A2 emissions scenario used in those articles. To compare our results with the direct feedback parameters of Arora et al. (2013), we use the same scenario used in that article in which CO 2 concentration increases 1 % yr −1 until it quadruples (approximately 140 years). For each scenario, we perform four simulations: 1. Fully coupled simulation.
2. Completely uncoupled simulation, giving c off a (t) = c a0 + t E(t)dt for the emissions-driven scenario and the specified concentration pathway for the concentrationdriven scenario.
3. Biogeochemically coupled simulation. We switch off feedbacks involving temperature response to atmospheric CO 2 by setting λ = 0. Since our model contains no radiative forcing other than changes in CO 2 , temperature T = 0 in this simulation. From this simulation we estimate the carbon-concentration feedback

Feedback non-linearity
Zickfeld et al. (2011) found, in emissions-driven scenarios, that the fully coupled simulation sunk more terrestrial and marine carbon than the sum of the biogeochemically and radiatively coupled scenarios. They refer to this difference as the non-linearity of the feedback, with the land sink contributing 80 % of the non-linearity and the ocean sink 20 %. Our analytical expressions for the feedbacks can be used to obtain an alternative measure of feedback non-linearity. We evaluate the non-linearity in the ocean and land climate-carbon feedbacks using F * TO (c a , c m , c t , T ) − F * TO (c a0 , c m0 , c t0 , T ) and F * TL (c a , c m , c t , T ) − F * TL (c a0 , c m0 , c t0 , T ), respectively, where the F * (c a0 , c m0 , c t0 , T ) are analytical approximations of feedback factors from a radiatively coupled simulation (all carbon stocks are fixed at pre-industrial levels). We evaluate the non-linearities in the ocean and land concentrationcarbon feedbacks using F * O (c a , c m , c t , T )−F * O (c a , c m , c t , 0) and F * L (c a , c m , c t , T )−F * L (c a , c m , c t , 0), respectively, where the F * (c a , c m , c t , 0) are analytical approximations of feedback factors from a biogeochemically coupled simulation (temperature is fixed at its pre-industrial level). These expressions indicate the effect of temperature and atmospheric carbon on the concentration-carbon and climate-carbon feedbacks, respectively, We used the SRES A2 scenario.

Model evaluation
Most predictions of our model are within the range of model predictions produced for the IPCC's Fifth Assessment Re-port (Table 2). Our model estimates around 55 Pg C more historical land carbon uptake than the IPCC multi-model mean, possibly due to our simplification to a single land carbon pool. Because it omits radiative forcing due to greenhouse gases other than CO 2 , our model consistently underestimates future temperature changes, although in all except the RCP8.5 scenario the projections are within the IPCC model range. The purpose of our model is not to precisely predict future climate change, but to serve as an approximate, mechanistically based emulator of the carbon cycle component of Earth system models (see Sect. 2). If we choose parameters to fit historical observations rather than based on the best available knowledge about each process (see Sect. 3), then our results remain mostly within IPCC model range although ocean and land uptake are consistently above and below the IPCC multi-model mean, respectively (Table A2a). We conclude that the model emulates historical observations and future projections of Earth system models with sufficient accuracy for this purpose.

Feedback analysis
All feedback measures calculated directly from our stylised model simulations, as well as most of our analytically estimated feedback measures, are within a factor of 2 of the mean output from Earth system models reported by  Arora et al. (2013) for direct feedback parameters). This is a remarkably good agreement considering the highly stylised nature of our model. Furthermore, all of the numerically time-integrated feedback parameters from our stylised model are within the multi-model range reported by Friedlingstein et al. (2006). The agreement observed here serves as additional validation of our model as well as validation of the approximations used to calculate analytical feedback factors.
An exception to the close agreement is the ocean concentration-carbon feedback, with the analytically estimated feedback factor and time-integrated feedback parameter indicating a weaker negative feedback than the numerical estimates from our stylised model or Earth system models. This mismatch is primarily due to two approximations: one in the numerical simulation and one in the analytical approximation. The numerical approximation is that disconnecting climate feedbacks in the biogeochemically coupled simulation leaves less carbon available to be sunk into the ocean, placing the ocean feedback at a different point in the highly non-linear (as parameterised by the Revelle factor) ocean carbon uptake dynamics. The analytical approximation is that Eq. (21) underestimates carbon sunk into the deep ocean.
We used parameters (Table 1) based on the best available data about each process (see Sect. 3). With a set of parameters based instead on fit to historical changes (Table A2), the numerically estimated feedbacks became slightly stronger: that is, the already positive climate-carbon feedbacks be-came more positive and the already negative concentrationcarbon feedbacks more negative. The numerical feedback estimates retained, however, good agreement with analytical estimates as well as with previous numerical estimates by Friedlingstein et al. (2006) and Zickfeld et al. (2011). One exception was the ocean concentration-carbon feedback, for which the analytical estimate remained outside Friedlingstein et al.'s range as noted above, but the direct numerical estimate moved to also be outside their range. We conclude that our results are relatively insensitive to parameter values, though mechanistically based parameter values perform slightly better than fitted parameter values.
Focusing on the analytical expressions, we observe that the approximate analytical expressions for the three different feedback measures all have similar dependences on state variables and parameters. All measures of the land climate-carbon feedback have dependence of the form c t0 log Q R /Q T /10 R . The ocean climate-carbon feedbacks all have terms of the form B 0 B T and w 0 D T c m /r. The land concentration-carbon feedback has the form c t0 K C /c a and the ocean concentration-carbon feedbacks have the form w 0 c m /rc a . We conclude that for each type of carbon cycle feedback, all three feedback formalisms detect similar features of the climate-carbon system.
The analytical expressions provide rapid insight into how feedback strengths depend on state variable and parameter values that could otherwise only be studied numerically or by qualitative reasoning. The analytical forms show that increasing Revelle factor r, as is likely to occur in an increasingly acidic ocean (Sabine et al., 2004), will decrease the strengths of ocean climate-carbon and concentration-carbon feedbacks. Weakening overturning circulation, via w 0 , would also decrease the strength of the ocean carbon cycle feedbacks. On land, the parameters Q R and K C control the terrestrial carbon cycle feedbacks.
We can compare likely trends in feedback strengths from the analytical expressions for the direct feedback parameters. According to our model, the destabilising ocean climate-carbon feedback is almost constant, while the ocean concentration-carbon feedback weakens with c m (since c m /c a ∼ c 1−r m ). Similarly, according to our model the destabilising land climate-carbon feedback will weaken less than the stabilising concentration-carbon feedback (under CO 2 doubling, ∼ Q − T /10 R weakens by 9 % at the new temperature equilibrium while ∼ 1/c a weakens by 50 %). This difference between the land climate-carbon and concentrationcarbon feedbacks stems from the differing curvatures of K(c a , T ) as a function of T (close to linear) and c a (concave). We conclude that under continued carbon emissions, according to our model, both land and ocean feedbacks will overall become more positive.
Where multiple processes contribute in parallel to a feedback, inspection of analytical forms can indicate the relative contributions of the different processes to the feedback. In the ocean component of the model, CO 2 solubility, the biological pump, and the solubility pump are all temperature dependent and therefore contribute to the ocean climate-carbon feedback. Remarkably, all three processes contribute temperature dependences of a similar magnitude; we therefore list all three in the approximate analytical gain and timeintegrated feedback parameter in Table 3. The three terms represent temperature dependence of the biological pump, CO 2 solubility, and the solubility pump.

Feedback non-linearity
As shown in Sect. 4.5, our analytical feedback expressions enable a new way of estimating feedback non-linearities that is not possible from direct numerical simulation. Since the sum of the four non-linearities is negative (Table 3), we conclude that summing feedbacks found by individual decoupled simulations will overestimate the atmospheric carbon levels, that is, underestimate land and ocean sinks. This result matches the findings of Zickfeld et al. (2011) and Matthews (2007). Terrestrial feedbacks contributed 83 % of the total non-linearity in our model, compared to 80 % reported by Zickfeld et al. (2011). Furthermore, we can distinguish the non-linearities in the climate-carbon and concentrationcarbon feedbacks. We found that the non-linearity in the terrestrial carbon-climate feedback was almost 4 times larger than any other (Table 3). By inspecting the analytical derivation of the gains we conclude that this dominance is likely due to a combination of three reasons: first, due to the sensitivity of temperature to carbon dioxide, ∂T /∂c a = λ/c a log 2, the carbon-climate feedbacks are much more sensitive to c a than the concentration-carbon feedbacks are to T . Second, the non-linearity in the land climate-carbon feedback is larger than the ocean climate-carbon feedback because its feedback factor is larger and therefore more sensitive to changes in gain (see Eq. 13). Third, the century timescale of the simulation prevented ocean carbon dynamics, which generally take place on longer timescales, from being exhibited. We conclude that care must be taken when summing results for feedbacks from decoupled simulations, especially for simulations involving land feedbacks.

Conclusions
Earth system models span a wide range of complexity. Here, we constructed a highly stylised, globally aggregated climate-carbon cycle model. Despite the model's simplicity -just four state variables -the model emulated globally aggregated historical trends and future projections of Earth system models. The model's simple form allowed climatecarbon cycle feedbacks to be estimated analytically, providing mechanistic insight into these processes. For example, we showed that carbon-climate feedbacks are less sensitive than carbon-concentration feedbacks; on land, this is due to the shape of K(c a , T ). The simple but accurate climate-carbon cycle model could be a starting point for model-based investigations of Earth system processes that are too poorly understood to be incorporated in more comprehensive models.
Stylised models such as ours have significant value in policy contexts. When confronted with difficult policy decisions involving long time periods and significant uncertainty, collaborative work with scientists allows policymakers to identify and clarify the impacts of various policy actions. In this context, the utility of a model is to increase stakeholders' understanding of a system and its dynamics under various conditions (Voinov and Bousquet, 2010;Anderies, 2005). This is in stark contrast to the use of more comprehensive models to predict impacts of policies in which mechanisms underlying dynamics and trade-offs are not transparent and quick explorations with stakeholders are not practical. The utility of a stylised model is in facilitating a learning process rather than in "accurately" predicting outcomes.
We foresee at least two strands of valuable future research based on the climate-carbon cycle model developed in this paper. First, our climate-carbon cycle model could be extended by including further processes relevant on different timescales of interest for Earth system analysis. This would enable a more in-depth analytical analysis of the feedback strengths and gains relating to other aspects of Earth system dynamics, such as the Earth's energy balance, carbon storage in the tropics compared to extratropics, albedo changes, the cryosphere, nutrient cycles, and even societal feedbacks. The task of characterising the Anthropocene as an epoch could thus move beyond qualitative comparison of human-impact trends to an improved characterisation of the feedbacks that maintain different Earth system "regimes". The effects on feedback strengths of different functional forms, such as the fertilisation effect K C , and how to constrain these functional forms from data could also be investigated and could yield insight into the continued divergence of Earth system model projections.
Second, the model could comprise a workbench for the systemic understanding of planetary boundary interactions and hence generate crucial insights into the structure of the safe operating space for humanity delineated by the planetary boundaries (Rockström et al., 2009;Steffen et al., 2015). Such extensions should focus on linking the core abiotic and biotic dimensions of the planetary boundary framework. The present lack of well-developed model representations of the dynamics and ecosystem structure of biosphere diversity, heterogeneity, and resilience, despite ongoing efforts in this direction (Purves et al., 2013;Bartlett et al., 2016;Sakschewski et al., 2016), emphasises the need for a more conceptual understanding of biosphere integrity, its vulnerability to anthropogenic perturbation, and its role for Earth system resilience.
Data availability. Input parameters are listed in Table 2 and time series inputs are from publicly available data as listed in Sect. 3.    Table A2. Testing parameters fitted to historical data. The following changes to parameter values were made to those listed in Table 1: K C = 0.25, Q R = 2.45, λ = 1.91 K, w 0 = 0.185 yr −1 . (a) Historical and projected changes of carbon stocks. See Table 2 for further information on how the figures were calculated and sources for model comparison. (b) Feedback analysis. See Table 3