On the proportionality between global temperature change and cumulative CO₂ emissions during periods of net negative CO₂ emissions

We use version 2.9 of the University of Victoria Earth System Climate Model (UVic-ESCM) [19], a model of intermediate complexity with a horizontal grid resolution of ${1.8}^{\circ }\times {3.6}^{\circ }$. This version of the UVic-ESCM includes a 3-D ocean general circulation model with isopycnal mixing and a Gent-McWilliams parameterization of the effect of eddy-induced tracer transport [20]. For diapycnal mixing, a Brian and Lewis [21] profile of diffusivity is applied, with a value of 0.3 × 10⁻⁴ m² s⁻¹ in the pycnocline. The ocean model is coupled to a dynamic-thermodynamic sea-ice model and a single layer energy-moisture balance model of the atmosphere with dynamical feedbacks [22].

The land surface and vegetation are represented by a simplified version of the Hadley Centre's MOSES land-surface scheme coupled to the dynamic vegetation model TRIFFID. Land carbon fluxes are calculated within MOSES and are allocated to vegetation and soil carbon pools of the five plant functional types represented by the vegetation model [23]. Ocean carbon is simulated by means of a OCMIP-type inorganic carbon-cycle model [24] and a marine ecosystem/biogeochemistry model solving prognostic equations for nutrients, phytoplankton, zooplankton and detritus [25]. The marine ecosystem model also includes a representation of the nitrogen cycle. The version of the UVic ESCM used here includes a marine sediment component, which is based on [26].

Several sets of simulations were run with the UVic ESCM: prescribed atmospheric CO₂ simulations with CO₂ first increasing and then decreasing at 1% per year (referred to as 1% CO₂ simulations), zero CO₂ emissions simulations initialized from the point of peak CO₂ in the 1% CO₂ simulations, zero CO₂ emissions 'hiatus' simulations, whereby a 1% decrease in atmospheric CO₂ is prescribed from given points along a zero CO₂ emissions trajectory, and constant CO₂ concentration hiatus simulations, whereby a 1% decrease in atmospheric CO₂ is prescribed from given points along a constant CO₂ concentration trajectory.

The 1% CO₂ simulations were initialized from a pre-industrial spinup run. Atmospheric CO₂ was prescribed to increase at 1% per year until the time of doubling (2×CO₂), tripling (3×CO₂) and quadrupling (4×CO₂) of the pre-industrial atmospheric CO₂ concentration, and then decrease at the same rate until the pre-industrial CO₂ concentration was restored (figure 1(a), solid lines).

Zoom In Zoom Out Reset image size

The zero emissions (ZEs) simulations were initialized from the the point of peak CO₂ in the 2×CO₂ (2×CO₂-ZE), 3×CO₂ (3×CO₂-ZE) and 4×CO₂ (4×CO₂-ZE) simulations, with zero prescribed CO₂ emissions after that time (figure 1(a), dotted lines).

The zero CO₂ emissions hiatus simulations were initialized from given points (50, 100, 250, 500 and 1000 years after the time of peak CO₂) along the 4×CO₂-ZE simulation. From these points, a 1% decrease in atmospheric CO₂ was prescribed until the pre-industrial atmospheric CO₂ level was restored. These simulations are referred to as 4×CO₂-ZE-h50, 4×CO₂-ZE-h100, 4×CO₂-ZE-h250, 4×CO₂-ZE-h500, and 4×CO₂-ZE-h1000, respectively. Note that the atmospheric CO₂ level from which the 1% CO₂ decrease was prescribed differs among simulations and therefore the CO₂ change between the beginning and the end of the prescribed CO₂ decrease period is also different (figure 4(a)).

The constant CO₂ concentration hiatus simulations were initialized from given points (50, 100, 250, 500 and 1000 years after the time of peak CO₂) of a simulation with atmospheric CO₂ concentration held constant after the peak in the 4×CO₂ simulation. From these points, a 1% decrease in atmospheric CO₂ was prescribed, restoring atmospheric CO₂ to pre-industrial levels (supplementary figure 1). These simulations are referred to as 4×CO₂-CC-h50, 4×CO₂-CC-h100, 4×CO₂-CC-h250, 4×CO₂-CC-h500, and 4×CO₂-CC-h1000, respectively.

All simulations were forced with prescribed changes in atmospheric CO₂ or CO₂ emissions only. All other other forcings were held constant at their pre-industrial level. In simulations with prescribed atmospheric CO₂ concentration, CO₂ emissions were diagnosed from the rate of increase in atmospheric CO₂ and carbon fluxes to the land and the ocean as the residual term in the carbon mass balance.

Negative CO₂ emissions are applied generically, without specifying any particular technology. We have assumed that the CO₂ captured from the atmosphere is removed permanently from the climate system (e.g. via underground storage).

First, we examine the relationship between global mean surface air temperature change and cumulative CO₂ emissions for the 1% CO₂ ramp-up, ramp-down simulations (figure 1, solid lines). On the upward limb of the 1% CO₂ simulations, surface air temperature and cumulative CO₂ emissions increase nearly linearly in response to exponentially rising atmospheric CO₂ concentrations. Surface air temperature peaks a few decades after the peak in atmospheric CO₂, with a longer lag for simulations with higher peak atmospheric CO₂. After the peak, surface air temperature declines, albeit at a slower rate than the rate at which it increased on the upward limb of the 1% CO₂ simulations. The lagged response of temperature to atmospheric CO₂ decrease is due to the slow thermal equilibration of the deep ocean. Cumulative CO₂ emissions start to decrease (i.e the rate of CO₂ emissions starts to become negative) right after the peak in atmospheric CO₂, indicating that the prescribed CO₂ decline cannot be achieved by CO₂ uptake by natural carbon sinks alone, but requires artificial removal of CO₂ from the atmosphere. Similarly to temperature, the decrease in cumulative emissions is slower than the increase on the upward limb of the 1% CO₂ simulations. This 'carbon inertia' is largely due to the slow biogeochemical equilibration of the deep ocean, with a small contribution from the lagged response of the terrestrial biosphere, as will be discussed below.

The relationship between surface air temperature change (ΔT) and cumulative CO₂ emissions (CE) exhibits hysteresis behaviour: for a given amount of cumulative CO₂ emissions, ΔT (relative to year 1) is larger on the upward than on the downward limb of the 1% CO₂ simulations (figure 2). Consistently with previous studies [1, 3, 4, 11], the ΔT versus CE relationship is approximately linear during the CO₂ ramp-up phase. During the CO₂ ramp-down phase, however, the ΔT versus CE curve is nonlinear, with larger nonlinearity for simulations with higher peak atmospheric CO₂ concentration. This nonlinearity is due to the inertia in physical and biogeochemical climate processes mentioned in the previous paragraph. The lagged response of temperature to CO₂ decline is due to the slow thermal equilibration of the deep ocean, which continues to take up heat for several decades after the peak in atmospheric CO₂ (figure 3(c)), and is consistent with the results of earlier studies [27–29]. Biogeochemical inertia is largely due to the ocean, with a smaller contribution from the terrestrial biosphere. The ocean continues to take up carbon even after atmospheric CO₂ levels start to decline, as the deep ocean is still equilibrating with past (positive) CO₂ emissions (figure 3(f)). On land, vegetation also continues to take up carbon for several years after atmospheric CO₂ starts to decline (figures 3(d) and (e)). This is mostly because of forest expansion at northern high latitudes, which lags atmospheric CO₂ due to the long timescales associated with vegetation shifts (not shown). The lag in the biogeochemical response to CO₂ decline is consistent with the results of [29, 30].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Ocean thermal and biogeochemical inertia affect the ΔT versus CE relationship on the downward limb of the 1% CO₂ simulations in opposite ways. This can be seen by separating the ratio ${\rm{\Delta }}T/\mathrm{CE}$ into two terms [1]:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}T}{\mathrm{CE}}=\displaystyle \frac{{\rm{\Delta }}{C}_{A}}{\mathrm{CE}}\displaystyle \frac{{\rm{\Delta }}T}{{\rm{\Delta }}{C}_{A}},\end{eqnarray} \tag{ 1 }$$

where ${\rm{\Delta }}{C}_{A}$ is the change in atmospheric CO₂. Note that for the CO₂ ramp-up phase (positive emissions) changes are defined relative to year 1, whereas for the CO₂ ramp-down phase (negative emissions) changes are defined relative to the time of peak CO₂ concentration. As is evident in figure 3(a), biogeochemical inertia results in less cumulative negative emissions required to achieve the prescribed drop in atmospheric CO₂ (i.e. larger ${\rm{\Delta }}{C}_{A}/\mathrm{CE}$), acting to increase the ratio ${\rm{\Delta }}T/\mathrm{CE}$. On the other hand, thermal inertia results in a smaller temperature drop for the prescribed atmospheric CO₂ decline (i.e. smaller ${\rm{\Delta }}T/{\rm{\Delta }}{C}_{A}$) (figure 3(b)), which acts to decrease ${\rm{\Delta }}T/\mathrm{CE}$. As will be discussed later, the TCRE is smaller on the downward than on the upward limb, indicating that thermal inertia is slightly larger than biogeochemical inertia. This difference between thermal and biogeochemical inertia is the result of the different vertical profiles of heat and carbon in the ocean. Ocean heat and carbon uptake are governed by CO₂ and temperature gradients between the ocean mixed layer and the atmosphere, and the mixing with the deep ocean. Since the deep ocean is colder than the mixed layer throughout the 1% CO₂ simulations, mixing with the deep ocean cools the mixed layer, enhancing the temperature gradient at the atmosphere-ocean interface. This delays the reversal in air-sea heat flux to nearly the end of the CO₂ decline phase (figure 3(f)). On the other hand, since the deep ocean is enriched in carbon relative to the mixed layer, upward mixing of deeper waters acts to decrease the CO₂ gradient at the sea surface, causing the ocean to turn from a sink to a source of carbon only a few decades after the start of the CO₂ decline phase (figure 3(c)).

The balance between ocean thermal and biogeochemical inertia, which determines the width of the hysteresis in the global mean temperature change versus cumulative CO₂ emissions curve, is set by the timescale of mixing of heat and carbon into the deep ocean and could therefore be dependent on a model's ocean mixing parameterization. Simulations with the UVic ESCM using different mixing parameterizations and a range of mixing parameters indicate that the balance between ocean heat and carbon uptake in simulations with a 1% per year increase in atmospheric CO₂ concentration to 4× pre-industrial levels and constant CO₂ concentration thereafter is largely independent of the mixing parameterization [31]. Therefore we expect this balance to be maintained also for zero end negative CO₂ emissions scenarios, and the results discussed above to be robust against the choice of mixing parameterization. This inference is further supported by simulations with a complex Earth system model with a symmetric 1% per year increase and subsequent decrease in atmospheric CO₂, which also exhibits larger thermal than biogeochemical inertia [29].

The system's response on the downward limb of the 1% CO₂ simulations results from the combination of the lagged response to positive CO₂ emissions prior to the decline in atmospheric CO₂, and the response to negative CO₂ emissions. The continued warming and carbon uptake in response to past positive emissions can be quantified in simulations with prescribed zero emissions after peak atmospheric CO₂ (figure 1, dotted lines). Over the time of the CO₂ ramp down in the 1% CO₂ simulations (70, 110 and 140 years, respectively), atmospheric CO₂ levels drop by 67 ppm, 106 ppm, 138 ppm and surface air temperature increases by 0.1 °C, 0.4 °C, 0.7 °C in the 2×CO₂-ZE, 3×CO₂-ZE and 4×CO₂-ZE simulations, respectively. If the temperature and carbon cycle response on the downward limb of the 1% CO₂ simulations is corrected for the response in the ZE simulations, the ΔT versus CE relationship is closer to linear (figure 2, dashed lines): the 'bulge' at the beginning of the return path is reduced, and the curves for the 2×CO₂, 3×CO₂ and 4×CO₂ simulations converge to a similar value. This suggests that the relationship between global mean temperature change and cumulative CO₂ emissions is approximately linear also during periods of net negative CO₂ emissions, provided that the negative emissions are applied from a state where the system is at equilibrium with past CO₂ emissions.

To further test this idea, we examine the system's response in a set of simulations whereby the 1% atmospheric CO₂ decline is prescribed from different points along the 4×CO₂ ZE simulation (figure 4). The hypothesis is that the later the CO₂ ramp down is applied, i.e. the closer the system is to equilibrium with past emissions, the more linear the ΔT versus CE relationship. The linearity of the ΔT versus CE relationship is quantified as the relative deviation of ${\rm{\Delta }}T/\mathrm{CE}$ from the slope of the line calculated by linear regression. Results indicate that the later the CO₂ decline is applied, the smaller the deviation from linearity (figure 5). The maximum deviation from nonlinearity is about 10% on the downward limb of the 4×CO₂ simulation, where the CO₂ ramp down is applied instantly after peak CO₂ levels are reached, but only 2% in the simulation where 1000 years elapse before the CO₂ decline is applied (4×CO₂-ZE-h1000).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given that a near-proportional relationship exists between global mean temperature change and cumulative CO₂ emissions during periods of net negative CO₂ emissions, a proportionality constant can be defined, similarly to the TCRE for positive emissions. This proportionality constant will in the following be referred to as TCRE− to distinguish it from the TCRE during periods of positive emissions. Both TCRE and TCRE− are calculated by linear regression of the ΔT versus CE curve (table 1; see caption for details about the calculation). On the upward limb of the 1% CO₂ simulations, TCRE is lowest for the simulation with the largest peak atmospheric CO₂ concentration (4×CO₂), consistent with previous studies that show a lower TCRE at higher cumulative emissions (e.g. [13]). On the downward limb of the 1% CO₂ simulations, TCRE is also lowest for the simulation with the largest peak atmospheric CO₂ concentration, suggesting that negative emissions are less effective at cooling if applied at higher atmospheric CO₂ concentration. In the 4×CO₂ ZE hiatus simulations, TCRE− is lowest for the simulation in which the CO₂ decline is applied from the latest point on the zero CO₂ emissions trajectory, i.e. when the system is closer to equilibrium with past emissions (4×CO₂-ZE-h1000). As discussed in the previous section, this is also the simulation for which the ΔT versus CE relationship is closest to linear. Therefore, it can be expected that the TCRE− for the 4×CO₂-ZE-h1000 simulation (1.4 °C/TtC) is closest to the 'true' TCRE− value (i.e. the value that would result if the CO₂ decline were applied when the system is fully equilibrated with past CO₂ emissions). This value is somewhat lower than the TCRE on the upward limb of the 4×CO₂ simulation (1.6 °C/TtC), suggesting that negative CO₂ emissions are less effective at cooling than prior positive emissions are at warming. It should be noted that this difference does likely not arise from an asymmetry in the response of the climate-carbon cycle system to positive and negative CO₂ emissions, but rather from differences in the background atmospheric CO₂ concentration from which the positive and negative CO₂ emissions are applied. This state dependence can be explained by the stronger ocean stratification at 4×CO₂, from which the CO₂ decline is applied, compared to pre-industrial CO₂ levels, from which the CO₂ increase is applied. Due to this stronger stratification, it takes longer to mix heat upward into the ocean mixed layer than it took to mix heat downward during the period of CO₂ increase.

In the following, we examine ΔT and CE for the period of atmospheric CO₂ increase (between year 1 and the year of peak CO₂) and decrease (between the start of the CO₂ decline and the year of restoration of the pre-industrial CO₂ concentration) for the various simulations (table 1). The magnitude of both ΔT and CE is larger on the upward than on the downward limb of the 1% CO₂ simulations. For instance, the climate warms by 4.4 °C if atmospheric CO₂ is ramped up to 4×CO₂, but cools only by 3.3 °C if atmospheric CO₂ is ramped down to pre-industrial levels at the same rate. Also, cumulative CO₂ emissions consistent with ramping up atmospheric CO₂ to 4×CO₂ are 2890 GtC, whereas 2430 GtC need to be removed to restore atmospheric CO₂ from 4×CO₂ to pre-industrial levels. The lower ΔT and CE on the downward limb of the 1% CO₂ simulations are associated with the continued warming and CO₂ uptake in response to past positive CO₂ emissions discussed in the previous section. If ΔT and CE are corrected for the warming and CO₂ uptake in the ZE simulations, both the amount of cooling and the required negative CO₂ emissions are larger during the CO₂ ramp-down phase (3.9 °C and 2710 GtC for the 4×CO₂ simulation). It should be noted, however, that additional negative CO₂ emissions are required to maintain atmospheric CO₂ at the pre-industrial level.

Consistent with the TCRE discussed in the previous paragraph, the ratio ${\rm{\Delta }}T/\mathrm{CE}$ on both upward and downward limbs of the 1% CO₂ simulations decreases with higher peak CO₂ levels. The value for ${\rm{\Delta }}T/\mathrm{CE}$ on the upward limb is very similar to the TCRE calculated by linear regression, reflecting the near-linear relationship between ΔT and CE. On the other hand, the value for ${\rm{\Delta }}T/\mathrm{CE}$ on the downward limb is smaller than TCRE− reflecting the nonlinear ΔT versus CE curve (the 'bulge' in the curve steepens the slope calculated by linear regression).

In the 4×CO₂ ZE hiatus simulations, both ΔT and CE are smaller for simulations which are more closely equilibrated with past CO₂ emissions. This is a result of the smaller CO₂ decrease during the CO₂ ramp-down phase (figure 4) for simulations which follow the 4×CO₂ ZEs trajectory for longer (for instance, atmospheric CO₂ decreases by 230 ppm less in 4×CO₂-ZE-h1000 than in the 4×CO₂-ZE-50 during the CO₂ ramp-down phase). The ratio ${\rm{\Delta }}T/\mathrm{CE}$ is larger the longer the system follows the 4×CO₂-ZE trajectory, which is inconsistent is with the decrease in TCRE−, a result of the different ΔCO₂ during the CO₂ ramp-down phase in the different simulations.

To remediate the problem of different ΔCO₂ in 4×CO₂-ZE-h simulations, we examine ΔT and CE for 4×CO₂ constant concentration hiatus simulations (4×CO₂-CC-h), whereby a 1% decrease in atmospheric CO₂ is prescribed from given points along a constant 4×CO₂ concentration simulation (which results in the same CO₂ decline for all simulations). The cooling over the CO₂ decline phase exhibits non-monotonic behaviour as a function of the time CO₂ is held constant before the 1% CO₂ decline is applied (referred to as 'hiatus'): it increases up to a hiatus of 250 years and then decreases (table 1). In contrast, the amount of negative cumulative emissions required to restore atmospheric CO₂ from 4×CO₂ to pre-industrial levels increases with the length of the hiatus. The ratio ${\rm{\Delta }}T/\mathrm{CE}$ decreases with the length of the hiatus, consistently with the decrease in TCRE− for these simulations, with the value of ${\rm{\Delta }}T/\mathrm{CE}$ being similar to the TCRE− value for simulations with longer hiatus. The ΔT and CE behaviour in 4×CO₂-CC-h simulations can be understood in terms of two competing processes: the equilibration with past positive CO₂ emissions, which, based on the discussion above, can be expected to lead to a larger cooling and larger cumulative negative CO₂ emissions in simulations with longer hiatus; and the equilibration with 4×CO₂ levels, which leads to reduced cooling and larger required negative CO₂ emissions in simulations with longer hiatus. The longer the system is exposed to constant 4×CO₂ levels, the warmer and more carbon rich the deep ocean, such that an imposed decline in atmospheric CO₂ is less effective at cooling (as the water that mixes to the surface from the deep ocean from the deep ocean is warmer) and requires larger amounts of negative emissions to attain a desired (lower) CO₂ level.

So far, we have examined ΔT and CE only over the period of atmospheric CO₂ decline. If we consider a longer time period (400 years from the start of the CO₂ decline, with atmospheric CO₂ restored to pre-industrial concentrations for several centuries), we find that in the 4×CO₂ ZE hiatus simulations the amount of cooling is smallest and the temperature is furthest away from pre-industrial the longer the system is left to equilibrate with past CO₂ emissions (figure 4). The amount of negative emissions is also largest for simulations with the longest ZE hiatus. These results indicate that the later along the zero CO₂ emissions trajectory the CO₂ decline is prescribed, the longer it takes to restore global mean surface air temperature to pre-industrial levels.