The oceanic origin of path-independent carbon budgets.

Virtually all Earth system models (ESM) show a near proportional relationship between cumulative emissions of CO2 and change in global mean temperature, a relationship which is independent of the emissions pathway taken to reach a cumulative emissions total. The relationship, which has been named the Transient Climate Response to Cumulative CO2 Emissions (TCRE), gives rise to the concept of a 'carbon budget'. That is, a finite amount of carbon that can be burnt whilst remaining below some chosen global temperature change threshold, such as the 2.0 °C target set by the Paris Agreement. Here we show that the path-independence of TCRE arises from the partitioning ratio of anthropogenic carbon between the ocean and the atmosphere being almost the same as the partitioning ratio of enhanced radiative forcing between the ocean and space. That these ratios are so close in value is a coincidence unique to CO2. The simple model used here is underlain by many assumptions and simplifications but does reproduce key aspects of the climate system relevant to the path-independence of carbon budgets. Our results place TCRE and carbon budgets on firm physical foundations and therefore help validate the use of these metrics for climate policy.

Introduction

Simulations with Earth system models consistently show a near-linear relationship between cumulative emissions of CO2 and change in global mean temperature[1-3], a relationship which has been shown to be consistent with observations[3]. The relationship is fundamental to the concept of a carbon budget, a finite quantity of cumulative CO2 emissions consistent with remaining below a chosen temperature change threshold[4], such as the 1.5 and 2 °C targets set by the Paris Agreement[5]. Carbon budgets have in turn been used to estimate the fraction of known fossil fuel reserves that can be safely burnt[6, 7] and to attribute historical responsibility for climate change[8]. Thus carbon budgets have proven a useful tool to conceptualize the challenges inherent in economic decarbonization[9].

The relationship between cumulative emission of CO2 and change in global mean temperature, the slope of which is the Transient Climate Response to Cumulative CO2 Emissions (TCRE), has two key features: (1) the relationship is close to linear, that is TCRE is nearly constant; and (2) the relationship is independent of the emissions pathway followed to get to a given quantity of cumulative CO2 emissions[1, 10]. Previous efforts to explain the physical origin of TCRE have largely focused on why the relationship is linear[11-13]. However, the second feature of TCRE, its path-independence, is arguably more important for climate policy as one can compute carbon budgets for non-linear temperature versus cumulative emission curves but only if the relationship remains path-independent. The physical origin of the path-independence of carbon budgets has not been fully explained[14], with only a single study examining the topic[15].

The near path-independence of carbon budgets is evident in the full spectrum of Earth System Models (ESMs) from simple box models[11, 15], to intermediate complexity climate models[10, 16] to full ESMs[17]. Therefore we hypothesize that the path-independence of carbon budgets must arise from some fundamental aspect of carbon cycle dynamics. To explore this hypothesis we will derive a simple model of the Earth system that exhibits path-independent carbon budgets and examine why the model exhibits path-independence. Thus by stripping an Earth system model down to its most fundamental components we hope to derive a model which can be more easily understood.

Results

The link between carbon-cycle dynamics and path independence

The first step in our analysis is to confirm that the path-independence of temperature versus cumulative emission curves arises from modulation of CO2 concentration by the simulated carbon cycle. We do this using model simulations with the University of Victoria Earth System Climate Model (UVic ESCM) (see methods) forced with CO2 that is decoupled from the model’s biogeochemical components (see Methods) (radiatively coupled – R-CO2)[18]. Experiments using radiatively coupled CO2 are part of the standard experiment-set conducted for carbon cycle model intercomparison projects[18]. Figure 1a,b shows temperature versus cumulative emission curves for two idealized climate model simulations, one experiment forced with fully coupled CO2 and the second experiment forced with R-CO2. The figure shows that the temperature versus cumulative emissions curves are nearly linear for R-CO2 but with large variation in slope of the curves for simulations with different rates of change of atmospheric R-CO2. Consequently the carbon budget for 2 °C of warming for R-CO2 increases by 17% when comparing an atmospheric CO2 change rate of 2 ppm a−1 to a change rate of 5 ppm a−1. The same metric for fully coupled CO2 changes by just 0.3%. Thus some aspect of the of carbon biogeochemistry is generating the near path-independence of carbon budgets.

Figure 1

Surface air temperature (SAT) anomaly versus cumulative emissions of CO2 curves for: (a) simulations with the UVic ESCM forced with CO2 that is uncoupled from the model’s biogeochemical components (radiatively coupled CO2 – R-CO2). (b) Simulations with the UVic ESCM forced with fully coupled CO2. (c) Solutions for the ZD2OM not accounting for ocean carbonate chemistry or land carbon uptake. (d) Solutions for the full ZD2OM. All models are forced with scenarios where atmospheric CO2 concentration changes at a constant rate. Note that the horizontal scale is different in the right and left columns as R-CO2 emission have an airborne fraction of 1.

A simple Earth system model

The simple Earth system model devised here is an extension of the analytical approximation of TCRE developed by ref. 13. The function for TCRE derived by ref. 13 was intended to explore why the temperature versus cumulative CO2 emission curve is nearly linear, and thus why TCRE was near-constant with respect to cumulative CO2 emissions. The paper suggested that the diminishing efficiency of ocean heat uptake and the logarithmic relationship between CO2 concentration and radiative forcing compensate in such a way to create the near-constant TCRE with respect to changes in cumulative CO2 emissions. The relationship derived by ref. 13 also implied that for path-independence to exist the airborne fraction of carbon must be being modulated as a function of the rate of CO2 emissions. Such modulation is evident from the output of Earth System Models[14] and is believed to originates from ocean carbon uptake[13]. However, the physical mechanism underlying this modulation has remained poorly explained.

For this manuscript we have been able to incorporate a representation of ocean carbonate chemistry into the model using the carbonate-alkalinity approximation of ref. 19, opening the way to finding the physical origin of TCRE path-independence. We have named the model the Zero Dimensional Diffusive Ocean heat and carbon uptake Model (ZD2OM). The model is based on ref. 1 definition of TCRE, the forcing response equation from ref. 20, and a relationship for cumulative CO2 emissions (see Methods). Combining these relationships with an approximation of heat and carbon uptake as diffusion into a half-space, and the ref. 19 analytical approximation of ocean carbonate chemistry, we find an analytically solvable model approximating TCRE (see Methods). The ZD2OM has analytical solutions for linear and exponential growth in atmospheric CO2 concentration. Linear pathway solutions for the ZD2OM are shown in Fig. 1d and are compared with simulations with the UVic ESCM forced with the same CO2 concentration pathways. Also shown in Fig. 1c are linear pathway solutions for a version of the ZD2OM without carbon uptake by the ocean or land surface (see Methods). The figure shows that the full ZD2OM is able to capture the path-independence characteristic of the TCRE relationship, and that the modified ZD2OM is able to capture the path-dependent characteristic of simulations forced with R-CO2. Given the land carbon uptake in the ZD2OM is a fixed fraction of emissions the figure strongly implies that ocean carbonate chemistry is somehow generating the path-independence of TCRE.

The ZD2OM deviates from linearity at lower cumulative emissions than the UVic ESCM. The deviation from linearity is a common feature of models which possess oceans that are too diffusive relative to ocean models with higher fidelity and to the natural ocean[21]. The ZD2OM ocean is purely diffusive and therefore the deviation from linearity is expected. Here we are focused on finding the origin of path-independence of carbon budgets, not linearity in the TCRE relationship which is relatively well explored feature[11-13].

Figure 2a compares the rate of CO2 emissions computed by the ZD2OM for the historical era to data-based estimates of emission rates from the same era (see Methods). The ZD2OM produces an estimate of the CO2 emission rate that is generally consistent with the historical record despite the radical simplicity of the model. Notably the ZD2OM emission estimate has a higher inter-annual variability relative to the historical data. The source of this high variability in the ZD2OM output is likely the assumption that the land-borne fraction of carbon is constant in time, which contrasts with observational evidence that there is strong inter-annual variability in the uptake of carbon by the terrestrial biosphere[22]. The inter-annual variability in the land carbon sink has been linked to coupled ENSO-forest-fire dynamics and other complex feedbacks[23]. As the ZD2OM’s equation for emissions conserves mass, changes in the atmospheric CO2 change rate are interpreted as changes in emissions rate instead of changes in the land-borne fraction of carbon. Although the ZD2OM only has elementary form solutions for a small number of families of atmospheric CO2 pathway functions (linear or exponential growth), the equations which underlie the model can be integrated numerically for any CO2 pathway. This feature of the model is demonstrated in Fig. 2b and c, which show the change in global temperature and temperature versus cumulative emission curves produced by the ZD2OM when forced with the Representative Concentration Pathways (RCPs)[24]. Note that the ZD2OM produces relatively standard temperature versus cumulative emissions curves until the rate of CO2 slows near the end of each scenario where path-independence is lost.

Figure 2

(a) Comparison between the ZD2OM calculation of CO2 emissions rates and historical estimates of anthropogenic CO2 emissions rate for the industrial period. ZD2OM emission rates are computed from numerical integration of Equations (13), and (16) forced with the historical CO2 concentrations. Historical emission estimates are from ref. 46 updated to the present using data from ref. 23. Uncertainty bounds are from Table 6.1 of ref. 22 (b) ZD2OM estimate of global temperature change under the four Representative Concentration Pathway (RCP) scenarios. (c) Temperature versus cumulative CO2 emissions curves produced by the ZD2OM under the four RCPs.

A critical feature of the ZD2OM is that the model treats ocean heat and carbon uptake as diffusion into a half-space. Thus we must check that this assumption approximates heat and carbon uptake in full ESMs. Figure 3 shows the evolution of removal velocity (surface flux divided by the surface concentration anomaly of the fluxed quantity, see Methods) of heat and carbon computed from ESM output from the UVic ESCM and nine models from the Coupled Model Intercomparison Project phase 5 (CMIP5). Using removal velocity facilitates easier comparison of ocean heat and carbon uptake by examining a metric that has the same units for both quantities. If heat and carbon uptake by the ocean were governed by the exact same mechanism in ESMs then the heat and carbon removal velocity functions would be identical. For the models shown in Fig. 3 the velocities for heat and carbon are similar but not identical, consistent with spatial differences in uptake explored by ref. 25 and the biogeochemical processes that affect ocean carbon uptake[22] but not heat uptake. The figure also shows the fit between the removal velocities and the diffusive approximation of this quantity used to define the ZD2OM for heat and carbon, and for each ESM. The fits are generally close after the first 20 years of the simulations (see Supplementary Tables S2 and S3). The figure thus demonstrates that at the global scale ocean heat and carbon uptake in ESMs can be approximated well by diffusion, at least while carbon emissions continue. The values of the ocean diffusion parameter computed from fitting the diffusive approximation to the heat and carbon removal velocities are of a similar order of magnitude, but generally larger in value (0.20 to 8.48 cm2s−1, median 3.73 cm2s−1) than the enhanced ocean diffusion value of ~1 cm2s−1 needed to balance deep water formation, computed originally by ref. 26 (see Supplementary Tables S2 and S3). Enhanced ocean diffusivity schemes are commonly used to close meridional overturning circulation in global ocean models e.g. ref. 27, despite the mechanism having been shown to be very weak in the natural ocean[27].

Figure 3

Heat and carbon removal velocities at the ocean surface (see Methods) for nine CMIP5 models and the UVic ESCM. The nine CMIP5 models are those that stored all of the necessary variables to compute removal velocity. Some models stored all the necessary variables to compute heat removal velocity but not carbon removal velocity or vice-versa. See Figures S1 and S2 for all derivable heat and carbon removal velocities respectively. For all models the removal velocity is derived for the idealized 1% yearly change in atmospheric CO2 experiment[41]. Black lines show that best fit between the model derived removal velocity and removal velocity computed from the approximation of ocean uptake as diffusion into a half-space. Fits were created by altering the effective ocean diffusivity. If ocean heat and carbon uptake where governed by identical processes in ESMs then the heat and carbon removal velocity functions would be the same.

A simplified representation of the ZD2OM is:where Λ is TCRE, C is the size of the atmospheric carbon pool, C is the original size of the atmospheric carbon pool, C is the size of the oceanic carbon pool, κ is the ocean heat uptake parameter, λ is the climate feedback parameter, R is radiative forcing from an e-fold increase in atmospheric CO2, and l is the land-borne fraction of carbon. All variables and constants are defined in in Table 1.

Table 1

Variables and constants used in this manuscript.

Variable	Description	Units
b	Change rate of atm. CO2	PgC a−1
C A	Atmospheric carbon pool	PgC
C O	Cumulative ocean carbon uptake	PgC
D s	Ocean suface DIC	mol m−3
E	Cumulative CO2 emissions	PgC
F	Radiative forcing	W m−2
L	Radiative response	W m−2
N	Ocean/planetary heat uptake	W m−2
T	Change in global temperature	K
T o	Change sea surface temperature	K
t	Time	a
q o	Ocean carbon uptake	PgC a−1
V q	Removal velocity	m a−1
β	Change rate of atm. CO2	% a−1
κ	Ocean heat uptake efficiency	W m−1 K−1
Λ	TCRE	K EgC−1
Constant
A c	Carbonate alkalinity: [HCO3 −] + 2[CO3 2−]	mol m−3
B o	Unit conversion constant (carbon)	m2 Pg mol−1
C Ao	Pre-Industrial atm. CO2	PgC
D so	Pre-Industrial ocean suface DIC	mol m−3
f o	Planetary ocean fraction	—
R	e-fold radiative forcing from CO2	W m−2
K 1	Equilibrium constant: [HCO3 −][H+]/[CO2*]	mol m−3
K 2	Equilibrium constant: [CO3 2−][H+]/[HCO3 −]	mol m−3
K o	Solubility of CO2	ppm−1 m−3 mol
l	Land-borne fraction of carbon	—
M r	CO2 atmospheric mixing ratio	PgC ppm−1
α	Ocean diffusivity	cm2s−1
Γ	Surface DIC change from e-fold change in atm. CO2	mol m−3
ε	Ratio of ocean to global temperature change	—
λ	Climate feedback parameter	W m−1K−1
ρC p	Heat capacity of water	J m−3K−1
τ	Unit conversion constant (heat)	s a−1

Values of the constants used in this manuscript are included in the Supplementary information Table S4. The abbreviation ‘atm.’ stands for atmosphere.

The ΔC term in Equation (1) is computed from an integral which includes the evolution of ocean surface Dissolved Inorganic Carbon (DIC) in time (see Methods). If ocean surface DIC evolution is approximated as a logarithmic function of atmospheric CO2 concentration then the integrand ΔC has elementary form function solutions for linear and exponential atmospheric CO2 concentration pathways (see Methods). Here we are interested in why TCRE does not change much given different pathways to get to the same atmospheric CO2 concentration. To examine this question within the context of linear changes in atmospheric CO2 concentration we select a value of atmospheric CO2 and explore the effect of rate of atmospheric CO2 change, b, on the TCRE relationship. A mathematically convenient value of atmospheric CO2 concentration is double the pre-industrial CO2 concentration (C  = 2 × C ), which gives the relationship:where B is a unit conversion constant, Γ change in ocean surface DIC from an e-fold increase in atmospheric CO2, μ is constant related to ocean diffusivity (see Methods), f is the fraction of the Earth covered by oceans, ρc is the heat capacity of water, τ is the conversion factor between seconds and years, ε is the ratio of sea surface to global surface temperature change. Two of the terms in the denominator of the Equation (2) are functions of b. Thus if TCRE is constant for different values of b then a change in b must cause an increase in the value of one term and a decrease in the value of the other term that cancel out almost exactly. Looking back to Equation (1) the two terms which contain b are equivalent to:

Figure 4a shows the evolution of these terms as a function of atmospheric CO2 change rate (b). That is, the x-axis represents different CO2 pathways. As b increases the ocean has less time to absorb carbon such that C is smaller than it would be for a pathway that takes more time to reach the same atmospheric CO2 concentration. Simultaneously the increase in b makes ocean heat uptake more efficient and thus κ larger. These effects nearly cancel to create stable value for the denominator of Equation (1), and thus a path-independent value of TCRE. The derivatives of the two terms cancel exactly where the two terms are of equal value and nearly cancel for a broad range of b values near the point of exact cancelation. This cancelation feature is a mathematical consequence of the one term being a function of and the other being a function of (see Supplementary information for a formal proof). This mathematical form arises in the ZD2OM from the ocean heat and carbon uptake both being driven by diffusion. In addition one quantity must be instantaneous (ocean heat uptake) and thus take on the inverse square root form, and the other quantity must be an time integral quantity (cumulative ocean carbon uptake) and thus take-on the square root function from, for the cancellation effect to emerge.

Figure 4

(a) Rate dependent terms from Equation (1) evaluated at the time of doubled atmospheric CO2. Terms are a function of the rate of CO2 change b At high atmospheric CO2 change rates cumulative ocean carbon uptake is lower while ocean heat uptake efficiency (κ) is higher. These two factors compensate exactly where the two terms have the same value and nearly compensate over a wide range of rates of atmospheric CO2 change centred about the point of exact cancelation. (b) Grey envelope encompasses region where the ZD2OM predicts TCRE should be nearly path-independent, defined as where change in the rate dependent terms compensate by greater than 95%. Thick black line is the atmospheric CO2 change rate at which the terms cancel exactly. Plotted atop this envelope are the historical rates of atmospheric CO2 change and the rates from the four RCP scenarios. (c) Fifth to ninety-fifth percentile ranges for the carbon and heat ratios as calculated from observation derived parameters. Red square is the median estimate for each. Black line is the value of the ratios that the ZD2OM predicts exact rate compensation would occur. Grey lines are labeled by the fraction of rate compensation that is predicted to occur with a heat ratio at that value.

The point where the two rate-dependent terms are equal is where their derivatives are of equal magnitude but of opposite sign and this point is therefore at the core of the pathway-independence of TCRE. By examining the exact cancelation point we can illuminate physically why the terms are close in value under ambient conditions. Setting the terms equal and collecting all of the carbon parameters to one side and the heat parameters to the other we find:where L is the enhanced longwave radiation to space (radiative response) from the Earth’s surface (L = λT), and N is ocean heat uptake (N = κT). Thus, these equations state that for a certain rate of atmospheric CO2 change the cumulative partitioning of anthropogenic carbon between the ocean and the atmosphere is the same as the instantaneous partitioning of anthropogenically enhanced radiative forcing between the ocean and space. The equivalency depends on the particular values for the parameters that govern ocean carbonate chemistry (summarized in Γ). As far as we can discern there is no fundamental reason that these ratios are so close in value, and no equivalent relationship is expected for non-CO2 greenhouse gases at their present day emission rates. For example N2O and CH4 are both held largely in the atmosphere, with atmospheric fractions at ~80%[22, 28] and 99%[29] of the atmosphere-ocean reservoir respectively. Therefore ocean uptake of these gases is expected to be small[22] and the equivalent of the left-hand ratio of Equation (5) should be much larger than the right-hand ratio for these gases (the right-hand side is independent of radiative forcing agent). However, very slow emissions of N2O may allow a match between the ΔN 2 O /ΔN 2 O ratio and the L/N ratio.

Figure 4b shows the envelope where TCRE should be within 95% of its theoretical peak value, as a function of atmospheric CO2 concentration and rate of atmospheric CO2 change. Plotted on top of this envelope are the pathways of historical CO2 change and the four RCP scenarios used in the last Intergovernmental Panel on Climate Change (IPCC) assessment report[30]. The figure shows that the historical CO2 pathway falls within the envelope expected to maintain path-independence of TCRE. The RCP scenarios fall within the envelope until atmospheric CO2 concentration is stabilized near the end of each scenario. Therefore the ZD2OM predicts that carbon budgets should be near path-independent for all plausible future emission scenarios.

Observation based estimates of the heat and carbon ratios

Values for the ‘heat ratio’ and the ‘carbon ratio’ can be computed from observation based values of atmospheric CO2 change, ocean carbon uptake, total radiative forcing, and ocean heat uptake (see Methods). These observation based estimates give a contemporary value of the carbon ratio of 1.5 [1.3 to 1.9] and a value for the heat ratio of 2.5 [0.45 to 10] where the values in brackets are the 5th to 95th percentile range. Figure 4c shows the range of these observation based estimates in addition to the theoretical bounds where the ratios should be close to path-independent, derived from Equation (15). The figure shows that the central estimates of the heat and carbon ratios are within the bounds that should create near path-independence. However, the large uncertainty associated with radiative forcing[31] and ocean heat uptake[32] mean that we cannot exclude the possibility that the ratios in the natural world are beyond the bounds that would create path-independence TCRE and carbon budgets.

Discussion

Our analysis suggests that one of the key features that generates the path-independence of TCRE in the ZD2OM is that heat and carbon uptake by the ocean are approximated by diffusion. However, simple ESMs and intermediate complexity ESMs generally have more diffusive oceans than full ESMs[21]. Additionally, because full ESMs are less diffusive they exhibit more linear temperature versus cumulative emission curves[33]. Thus the two properties that make carbon budgets useful for climate policy, linearity and path-independence, may play off one-another. More diffusive models are less linear and more path-independent, and less diffusive models are more linear and less path-independent. Such an effect may explain the path-dependent temperature versus cumulative emission curves shown by ref. 34 generated by the not very diffusive GFDL-ESMG2 (see Fig. 3). Full ESMs are likely to better represent the natural world than intermediate complexity models and thus carbon budgets for policy purposes may have to account for some path-dependence.

The study of ref. 15 used a similar approach to the approach used here to examine the origin of path-independence in temperature versus cumulative emission curves. The study examined a linear two-box ocean model and showed that for path-independence to exist the timescale for increase in cumulative emissions must be much greater than the timescale for decrease in airborne fraction of CO2 and much smaller than the damping-timescale of the ocean[15]. Thus path-independence should exist for all but very slow emission rates, consistent with the lower and upper bounds of the shaded area from Fig. 4b. The study of ref. 15 implicitly accounts for carbon biogeochemistry using impulse-response functions derived form ESM simulations[35]. Thus the study was unable to explicitly account for the role in ocean carbonate chemistry in generating path-independence.

The relationship between cumulative emissions of CO2 and change in global temperature has three features that make derivation of carbon budgets useful. (1) The relationship is nearly linear, thus a gram of carbon emitted 100 years ago has the same impact on climate as a gram of carbon emitted today or 100 years from now. (2) The relationship is nearly path-independent such that the pathway followed to get to an emission total is irrelevant. (3) Warming following cessation of CO2 emissions is either close to zero[36], or will peak and stabilize[2, 25]. Here we have provided a simple hypothesis for the existence of the second feature, path-independence. However, the model used here is unsuited for exploring the stability of warming following cessation of emissions. The close-to-zero post-cessation warming intuitively should arise from compensation between ocean heat and carbon uptake[25, 36] but more work is needed to fully explain the mechanism behind this feature of ESM output.

The analysis presented here suggests that the path-independence of TCRE will breakdown as emissions approach zero (Fig. 4c). As most emission scenarios envision a gradual reduction of emissions to zero[30] this effect may be important for computing the final value of the carbon budget for a given temperature target. Thus investigating the effect of long-tail low-rate emissions on the carbon budget may be an opportunity for research with more complex Earth system models.

Here we have shown that the path-independence of TCRE and therefore carbon budgets arises from: (1) ocean heat and carbon uptake being governed by a similar mechanism which at the global average is consistent with diffusion into a half-space. (2) The partitioning ratio of anthropogenic carbon between the ocean and the atmosphere being close in value the partitioning ratio of enhanced radiative forcing between the ocean and space. The first factor should apply to any long-lived greenhouse gas that is soluble in sea-water. The second factor, however, is a deep coincidence likely applicable only to CO2 and is the core property that generates the path-independence of carbon budgets. The heat and carbon ratios do appear to be close in the natural world, however confirming whether at the global average the natural ocean heat and carbon uptake can be approximated by diffusion, or whether the fit to diffusion is a property of ESMs that is an artifact of using refs 26, 27 based closure schemes in ESMs, is critical for validating the path-independence of carbon budgets.

Methods

UVic ESCM

The University of Victoria Earth System Climate model 2.9 (UVic ESCM 2.9) is a climate model of intermediate complexity composed of a full three dimensional ocean general circulation model coupled to a simplified energy and moisture balance atmosphere[37]. The model contains a full realization of the oceanic[38] and terrestrial carbon cycles[39, 40]. The simulations with the UVic ESCM shown in Fig. 1 were forced with atmospheric CO2 concentration pathways where CO2 changes at a constant rate of 1 to 5 ppm a−1 in intervals of 1 ppm a−1. Simulations with radiatively coupled CO2 were conducted using the protocol of ref. 18. Output from the UVic ESCM used to compute removal velocities for heat and carbon are from the 1% annual change in atmospheric CO2 concentration experiment, a standard benchmarking experiment used for model intercomparison[41]. For all simulations atmospheric CO2 concentration is prescribed and emissions of CO2 are diagnosed as the residual term of the carbon-cycle mass balance.

Derivation of the ZD2OM

The zero dimensional diffusive ocean model (ZD2OM) is a model of energy and carbon uptake by the Earth’s ocean that is simple enough to be amenable to analysis. The model is derived from the ref. 1 definition of the TCRE, the forcing-response equation[20], the ocean carbonate chemistry approximation of ref. 19, and the assumption that the ocean heat and carbon uptake are governed by diffusion into a half-space. Here we present a summarized derivation of the model. A complete derivation is included in the Supplementary materials.

The ref. 1 definition of TCRE is:where Λ is TCRE, T is the change in global mean temperature, and E is cumulative emissions of CO2. For the ZD2OM both T and E are subject to approximations. Global mean temperature change is derived from the forcing-response equation[20]:where F is the radiative forcing, λ is the climate feedback parameter, and N is planetary heat uptake. We assume that all radiative forcing is from changes in the atmospheric concentration of CO2 and that all heat uptake goes into the ocean. Thus we approximate F using the classical approximation of radiative forcing from CO2 [42]:where R is the radiative forcing from an e-fold increase in CO2, C is the atmospheric CO2 content, and C is the pre-industrial atmospheric CO2 content. N is approximated as heat diffusion into a half-space:where f is the fraction of the planet covered by ocean, ρC is the volumetric heat capacity of water, T is the global mean sea surface temperature anomaly, τ is the conversion factor between seconds and years, α is the diffusivity of the ocean, and t is time. To solve for temperature, N is represented by N = κ(t)T, where κ is the classical ocean heat uptake efficiency e.g. ref. 43, here taken to be a time-evolving function instead of a constant e.g. ref. 13. Substituting the above relation into Equation (9) we find:where ε is the ratio of sea surface temperature change to global temperature change, taken to be a constant. Analysis of National Oceanic and Atmospheric Administration (NOAA) ocean and global temperature change products shows that the decadal average values of ε vary between 0.79 to 0.85 since the 1950s without a clear secular trend[44] (see Table S1). Analysis of the output of climate model simulated sea surface temperature anomalies and 2 m air temperature anomaly output also suggest that ε is relatively constant (See Figure S1).

To reduce the number of constants we define and γ = f ρC τε. Solving for temperature we find:

Cumulative emissions are the sum of the carbon added to the atmosphere, ocean, and terrestrial biosphere from anthropogenic sources. For the ZD2OM we assume that the uptake by land is a constant fraction of anthropogenic emissions (l) and that the change in the atmospheric carbon pool is known. The estimated value of the land-borne fraction of carbon exhibits high uncertainty and strong inter-annual and decadal variability[22]. With no clear pattern in the decadal value of the parameter assuming a constant value in time appears to be the most parsimonious assumption. Therefore the evolution of emissions is given by the equation:where q is the flux of carbon into the ocean. We approximate q using a similar approximation as that used for κ:where D is the surface Dissolved Inorganic Carbon (DIC) of the ocean, D is the pre-industrial surface DIC of the ocean and B is a conversion factor equal to , and A is the area of the ocean, M is the molar mass of carbon, and 1 × 1015 is the conversion factor between g and Pg. Substituting Equation (13) into Equation (12) we find:

With functions for temperature and cumulative emissions of CO2 we can derive the function for TCRE:

To further evaluate Equation (15) we must link the change in atmospheric CO2 to the change in the surface concentration of DIC (C and D ). To investigate this link we incorporate an approximation of ocean carbonate chemistry into our model.

For our purposes we wish to solve the carbonate system analytically to derive a relationship between atmospheric CO2 concentration and ocean surface DIC. To make this possible we will make three simplifying assumptions: 1) we will ignore the temperature and salinity dependence of the carbonate equilibrium constants K1, K2, and CO2 solubility constant K. 2). We will approximate alkalinity as carbonate alkalinity following refs 3, 19) we will assume that at the surface the partial pressure of CO2 in the ocean is equal to the partial pressure of CO2 in the atmosphere. For the full derivation of Equation (16) see the Supplementary Information. From these assumptions we derive the approximation:where A is the carbonate alkalinity A  = [HCO3 −] + 2[CO3 2−] and M is the conversion factor between the mass of CO2 in the atmosphere and the partial pressure of CO2 in the atmosphere. The above equation can be substituted into the equation for TCRE (Equation (15)) to solve to the evolution of TCRE in time. However, the integral in the equation, , can only be solved using numerical integration. As we wish to derive as a system that can be analyzed analytically we therefore approximate D using a log-space Taylor series approximation of Equation (16) centred about C  = C :

This approximation of D is within 0.5% accurate of the full function up to an atmospheric CO2 concentration of 600 ppm (see Figure S5 in Supplementary information). As a check on the logarithmic approximation of DIC evolution and on the other simplifying assumptions outlined above, the surface DIC anomaly computed from Equation (17) is compared to output from the UVic ESCM in Figure S4. The figure shows that the surface DIC anomaly produced by Equation (17) closely matches that from the model output. At the point of doubled atmospheric CO2 concentration Equation (17) underestimates the DIC anomaly by 2.6%. With the D approximation in-hand we can derive an analytically amenable function for TCRE.

For simplicity we define a constant Γ to represent the factor in front of the logarithm:

This factor is equivalent to the increase in ocean surface DIC given an e-fold increase in atmospheric CO2 concentration. Substituting this relationship in the equation for ocean carbon flux (Equation (13)) we find:and the relationship for TCRE becomes:

The approximate form of the ZD2OM shown in Equation (20) can be solved analytically for certain CO2 concentration pathways. Two pathways relevant to the climate problem are a linear increase in atmospheric CO2 concentration C  = bt + C , and a exponential increase in atmospheric CO2 concentration C  = C e , where b is the rate of change in atmospheric CO2 content (Pg a−1) and β is the growth rate of atmospheric CO2 content (% a−1). The linear solution is:and the solution for an exponential increase in atmospheric CO2 concentration is:

ZD2OM without carbonate chemistry

To compare the ZD2OM with ESM simulations forced with R-CO2 a version of the ZD2OM with no carbonate chemistry and no uptake of carbon by land was derived. The temperature function is unaffected by carbon uptake and thus remains Equation (11). The function for emissions will simply be:

The equation for TCRE of R-CO2 therefore is:

Removal velocity

The removal velocity describes how efficient the ocean is at redistributing quantities at the ocean surface into the ocean interior. Removal velocity is defined as the ratio of the flux into the ocean over the anomaly of the surface concentration of the quantity being fluxed. That is, if we are given a flux the removal velocity is the instantaneous rate at which the substance must be removed from the surface to maintain the present surface concentration of the substance. The relationship has units of velocity (ms −1). The relationship for heat is:where V is the removal velocity. The relationship for carbon is:

The diffusive approximation of ocean heat and carbon uptake predicts that the removal velocity to be:

Removal velocity can be computed from model output saved in the CMIP5 archive and therefore can be used to test whether the diffusive approximation of heat and carbon uptake is accurate, to estimate values for μ and to test whether heat and carbon uptake are linked in models. For heat removal velocity both N and T are subject to large internal variability, therefore both variables are filtered with a 7 year moving average for plotting in Fig. 3. Removal velocity is analogous to the ‘piston velocity’ parameter used in diffusive gas exchange problems[19].

Heat and carbon ratios from observations

The carbon ratio can be estimated from the change in atmospheric CO2 and change in ocean carbon content since the beginning of industrialization. Observation derived estimates of both of these quantities were produced for the fifth assessment report of the IPCC (IPCC AR5)[22]. The estimate for change in atmospheric CO2 is 240 ± 10 PgC and for cumulative ocean carbon uptake is 155 ± 30 PgC. To compute an estimate of the carbon ratio we take the IPCC ranges to be the 90% uncertainty range and assume a normal distribution of uncertainty. A simple Monte-Carlo error propagation is performed to estimate the uncertainty on the carbon ratio.

The heat ratio can be estimated from the total radiative forcing and planetary heat uptake. Total radiative forcing was estimated by IPCC AR5[31] at 2.3 [1.1 to 3.3] W m−2. The estimate of planetary heat uptake was taken from ref. 32 and is 0.62 ± 0.43 W m−2. The heat ratio is , which is equivalent to . We assume that the uncertainty in F and N are normally distributed and use Monte-Carlo error propagation to compute the uncertainty bound for the heat ratio.

Historical emissions comparison

ZD2OM emission rates are computed from numerical integration of Equations (13), and (16) forced with the historical CO2 concentrations. The historical CO2 concentrations are taken from the RCP archive to year 2005[30] and extended to 2016 using data from Mauna Loa ref. 45. Historical CO2 emission estimates are from ref. 46 updated to the present using data from ref. 23. Uncertainty bounds are from Table 6.1 of ref. 22 for decades since the 1980s. Before the 1980s uncertainty bounds for yearly CO2 emissions are unavailable.

Supplementary material