Hysteresis of the Earth system under positive and negative CO₂ emissions

We use the Bern3D v2.0 s coupled to the LPX-Bern v1.4 dynamic global vegetation model. The Bern3D model includes a single layer energy-moisture balance atmosphere with a thermodynamic sea-ice component Ritz et al (2011), a 3D geostrophic-frictional balance ocean (Edwards et al 1998, Müller et al 2006) with an isopycnal diffusion scheme and Gent-McWilliams parameterization for eddy-induced transport (Griffies 1998), and a 10-layer ocean sediment module (Heinze et al 1999, Tschumi et al 2011, Jeltsch-Thömmes et al 2019). The horizontal resolution across Bern3D model components is 41 × 40 grid cells and 32 logarithmically spaced depth layers in the ocean. Wind stress is prescribed from the NCEP/NCAR monthly wind stress climatology (Kalnay et al 1996), and gas exchange at the ocean surface and calculation of carbonate chemistry follow OCMIP-2 protocols (Najjar and Orr 1999, Wanninkhof 2014, Orr and Epitalon 2015), with an adjusted gas transfer dependency on wind speed (Müller et al 2008). Marine productivity is restricted to the euphotic zone (75 m) and calculated as a function of light availability, temperature, and nutrient concentrations (P, Fe, Si) (see e.g. Parekh et al 2008, Tschumi et al 2011).

In the sediment, the transport, redissolution/remineralization, and bioturbation of solid material, the pore water chemistry, and diffusion are dynamically calculated in the top 10 cm (see Tschumi et al 2011, for more details). Burial (loss) of nutrients, carbon, and alkalinity from the sediment is balanced by a variable input flux to the coastal surface ocean during spin-up. These weathering input fluxes are set equal to burial fluxes at the end of the spin-up for transient simulations. The Bern3D model components have been evaluated comprehensively in several studies (e.g. Roth et al 2014, Battaglia et al 2016, Jeltsch-Thömmes et al 2019).

The Land surface Processes and eXchanges (LPX) model v1.4 (Lienert and Joos 2018) is a further development of the Lund-Potsdam-Jena model (Sitch et al 2003) and simulates the coupled nitrogen, water, and carbon cycles. Here, a resolution of 2.5° latitude × 3.75° longitude is applied. Vegetation is represented by plant functional types that compete within their bioclimatic limits for resources. A pattern scaling approach is used to pass climate information from the Bern3D on to LPX (Stocker et al 2013a). Spatial anomaly patterns in monthly temperature and precipitation, derived from a 21st century simulation with the Community Climate System Model (CCSM4), are scaled by the anomaly in global monthly mean surface air temperature as computed interactively by the Bern3D. These anomaly fields are added to the monthly baseline (1901 to 1931 CE) climatology of the Climate Research Unit (CRU) (Harris et al 2014). For further details on LPX-Bern the reader is referred to the literature (e.g. Keller et al 2017, Lienert and Joos 2018).

The Bern3D model is spun up under 1850 CE boundary conditions with constant CO₂ of 284.7 ppm for a total of 60 000 years. Similarly, LPX-Bern is spun up under 1850 CE boundary conditions with fixed land use for 1500 years, with an acceleration step towards equilibrium for soil pools, before the models are coupled and equilibrated for another 500 years. Simulations are started from 1850 CE conditions and run for 5 000 years, while all forcing except CO₂ is kept constant. Following the protocol of the carbon dioxide removal (CDR) experiment of the CDR model intercomparison project (Keller et al 2018), atmospheric CO₂ is increased at a rate of 1 %yr⁻¹ until 4 × CO₂ (simulation year 140). From there, prescribed CO₂ concentrations decrease at the following rates back to 284.7 ppm: 0.1, 0.3, 0.5, 0.7, 1.0, 2.0, 4.0, and 6.0 %yr⁻¹ (figure 1(a)). CO₂ is kept constant thereafter. Any remaining carbon flux from the ocean and land into the atmosphere are assumed to be removed by the sustained application of negative emission technologies. All simulations are conducted with three different model setups with equilibrium climate sensitivities (ECS) of 2 °C, 3 °C, and 5 °C. The ECS are diagnosed after 5000 years in instantaneous CO₂ doubling experiments. For each setup a control run with constant atmospheric CO₂ of 284.7 ppm is included.

Zoom In Zoom Out Reset image size

Cumulative carbon emissions to the atmosphere (figure 1(b)) are calculated for any period by summing the prescribed changes in the atmospheric carbon inventory (conversion: 2.12 GtC ppm⁻¹; figure 1(a), left axis) and the simulated changes in the cumulative atmosphere-ocean (figure 1(c)) and atmosphere-land biosphere (figure 1(d)) carbon fluxes.

We first illustrate changes in the global carbon cycle budget (figure 1). Diagnosed cumulative emissions depend on ECS and the CDR rate. Generally, cumulative emissions are higher for low ECS than for high ECS and are in the range of 2 640 (ECS = 5 °C) to 2 885 GtC (ECS = 2 °C) in year 140 (figure 1(b)) when atmospheric CO₂ peaks at 4 × preindustrial (figure 1(a)). Cumulative emissions start to decline for all ECS with the onset of carbon dioxide removal (CDR), except for very slow CDR rates.

Anthropogenic carbon is redistributed in the Earth system (figure 1). The carbon uptake or release by the ocean depends on the difference in CO₂ partial pressure in the ocean and atmosphere, $\Delta$pCO₂ = pCO$_{2,\textrm{ocean}}-$pCO$_{2,\textrm{atm}}$ (McKinley et al 2020) and is governed by surface-to-deep mixing.. The ocean is the dominant net carbon sink. At peak CO₂ (year 140), 1 805 GtC have entered the atmosphere, 645 GtC the ocean, and 355 GtC the land biosphere for an ECS of 3 °C. The ocean continues to take up CO₂ for a few more centuries under slow CDR rates, as equilibrium is not yet reached due to the slow ocean overturning. For a CDR rate of 0.1 %yr⁻¹, for example, the oceanic sink peaks at 1 260 GtC (ECS = 5 °C) to 1 420 GtC (ECS = 2 °C) around simulation year 780.

The ocean turns from a sink to a source during CDR because CO₂ is removed from the atmosphere and $\Delta$pCO₂ becomes positive. In turn, the cumulative atmosphere-ocean C flux begins to decrease. The timing and magnitude of the sign change depends primarily on the CDR rate. The faster the CDR, the earlier the ocean turns from a sink into a source. In all cases, the ocean releases carbon to the atmosphere for centuries to millennia after CO₂ is stabilized at preindustrial concentration (c.f. figures 1(a) and (c)). The cumulative atmosphere-ocean C flux varies slightly among ECS; a smaller oceanic sink for higher ECS follows from lower CO₂ solubility under higher temperatures and larger reductions in overturning (c.f. figures 1(c) and 2(a), (c)).

Zoom In Zoom Out Reset image size

The cumulative atmosphere-land biosphere flux increases and declines largely in concert with CO₂ and SAT (figure 1(d)). CO₂ fertilization of plants stimulates carbon uptake under higher atmospheric CO₂ and respiration is enhanced under higher temperatures. The land biosphere sequesters about 250 GtC less in soils and litter (∼170 GtC), and vegetation (∼80 GtC) under high than low warming, mainly due to enhanced temperature-driven respiration (figure 1(d)). The cumulative atmosphere-land biosphere flux declines with CO₂ to approach zero for low and intermediate ECS. Under high ECS, the cumulative flux becomes temporarily negative. Thus, the land biosphere stores less carbon than at preindustrial. This highlights the risk that the land biosphere may not only loose excess 'anthropogenic' carbon, but even some of its 'natural' carbon inventory. The latter represents an irreversibility with potential long-term consequences for ecosystems and ecosystem services.

On millennial timescales, CaCO₃ compensation adds up to 200 GtC to the ocean. For most scenarios and by the end of the simulations, the cumulative contribution from sedimentation-weathering imbalances in the CaCO₃ cycle ranges between an addition of 40 GtC to a removal of 20 GtC, and no new equilibrium is reached. CaCO₃ compensation changes alkalinity twice as much as carbon, leading to changes in the carbon uptake capacity of the ocean, explaining the long-term perturbation in the cumulative atmosphere-ocean C flux (figure 1(c)). Imbalances in the sedimentation-weathering cycle of particulate organic carbon (POC) are small and mainly result from changes in POC preservation under lower O₂.

$\Delta$SAT rises in response to higher CO₂ forcing (figure 2(a)), and at maximum cumulative emissions amounts to 2.7 °C, 3.7 °C, and 5.4 °C for ECSs of 2 °C, 3 °C, and 5 °C, respectively. The spatial warming pattern at maximum emissions resembles the global warming pattern for RCP 8.5 scenarios (IPCC 2013), with largest warming in high latitudes and over continents. Continued warming after maximum emissions scales with the subsequent removal rate and ECS. After returning to initial CO₂ concentrations, it takes centuries to millennia before SAT returns to its preindustrial value. The slower the CDR and the higher the ECS, the longer SAT stays elevated (figure 2(a)). Slow CDR scenarios and a high ECS yields $\Delta$SAT of up to 1 °C at the end of the simulation, while for ECS = 2 °C remaining $\Delta$SAT is below 0.1 °C.

$\Delta$OHC is, similar to $\Delta$SAT, largest for high ECS and slow CDR (figure 2(b)) and keeps increasing after maximum emissions. The surface ocean heat content equilibrates rather fast with the atmospheric temperature perturbation, while equilibration of the deep ocean is slow, with circulation as the rate limiting step. Thus, the slower the CDR, the longer the radiative imbalance lasts and the more of the perturbation will be transported to the deep ocean. Total OHC increases as long as the net downward heat flux at the top of the atmosphere is positive, resulting in an increase of $\Delta$OHC by a factor $\gt$2–3 after maximum CO₂, depending on CDR (see e.g. light colors in figure 2(b)).

$\Delta$AMOC is reduced by 4.8–8.8 Sv at maximum CO₂ forcing, depending on ECS and thus the corresponding temperature perturbation (figures 2(a) and (c)). AMOC starts to recover with declining forcing and overshoots to higher than initial values in all scenarios. The maximum overshoot is reached when CO₂ concentrations are back to initial. By the end of the simulations, AMOC states differ. Similar to SAT, slow CDR scenarios with an ECS of 5 °C show strongest deviations in AMOC strength at simulation year 5 000 and across the ensemble of runs. From the available simulations it is not clear if and when the AMOC would return to its initial strength.

Sea-ice reacts fast to changes in SAT and is substantially reduced in both hemispheres (figure 2(d)). The reduction is larger for higher ECS. Sea-ice recovers over most of the simulations. For ECS = 2 °C, however, sea-ice area seems to stabilize at ∼95% of its initial cover. The slow recovery and partial reversion into a new reduction of sea-ice area during and after CDR (see e.g. orange lines around year 600 in figure 2(d)) is driven by Southern Hemisphere (SH) sea-ice changes. Generally, the reduction in the SH is larger than in the Northern Hemisphere (NH) and longer lasting.

We quantify hysteresis as the difference in a variable between the up-path with increasing cumulative emissions and the down-path with decreasing cumulative emissions at 1000 GtC of cumulative carbon emissions. We thus consider large overshoot scenarios of up to about 2000 GtC (figure 1(b)).

We first examine differences in SAT between the up-path (dashed lines) and down-path (solid lines; figure 3(a)). Generally, the higher the ECS, the larger the hysteresis (see bars at the right of figure 3(a)). Non-linearities in SAT in the cumulative emission space do not scale linearly with ECS. For an ECS of 2 °C, hysteresis is small and largely negative. With ECS = 3 °C, SAT is in the mean 0.1 °C higher for the down- than up-path, while for ECS = 5 °C this difference increases to 1.5 °C. This means that at cumulative emissions of 1000 GtC and an ECS of 2 °C temperatures are mostly lower on the down-path compared to the up-path, while for ECS of 3 and 5 °C they are higher.

Zoom In Zoom Out Reset image size

Upper ocean (700 m) heat content (OHC₇₀₀, figure 2(b)) shows similar hysteresis as SAT. With decreasing cumulative emissions, OHC₇₀₀ first keeps increasing, until atmosphere and upper ocean are in equilibrium. For all three ECS, the mean hysteresis width is positive and ranges between 0.4 and 1.8 10²⁴ J.

Changes in O₂ result from changes in O₂ consumption due to remineralization, from changes in O₂ solubility, and to a lesser extent from changes in air-sea disequilibrium. Remineralization-driven changes in O₂ correlate well with changes in ideal age (not shown): the weaker the ventilation and thus the older the water masses, the more oxygen is consumed by remineralization of organic matter. In many regions, ocean overturning circulation and ventilation is reducing under global warming causing O₂ to decline. Higher temperatures at the surface lead to lower O₂ solubility and lower surface concentrations which are then transported to the ocean interior. Reductions in thermocline O₂ are larger for higher than lower ECS (figure 3(c)). The mean hysteresis width for an ECS of 5 °C is about –4 mmol m⁻³, while for ECS = 3 °C and ECS = 2 °C it is positive.

In contrast to SAT or O₂, changes in $\Omega_{\textrm{arag}}$ depend weakly on the ECS, but vary with the rate of CDR (figure 3(d)). $\Omega_{\textrm{arag}}$ decreases with increasing dissolved CO₂. Therefore, $\Omega_{\textrm{arag}}$ follows the increase and decline in atmospheric CO₂ without much delay in near surface waters and with delay in the deep ocean. At atmospheric CO₂ of about 565 ppm and cumulative emissions close to 1 000 GtC, there are almost no waters left with $\Omega_{\textrm{arag}}\ \gt$ 3. The decrease in cumulative emissions is delayed relative to the decrease in atmospheric CO₂, and $\Omega_{\textrm{arag}}$ starts to recover at much higher cumulative emissions during the down-path relative to the up-path. This hysteresis is particularly large for low CDR rates (figure 3(d)). Recovery of $\Omega_{\textrm{arag}}$ is generally faster for low ECS in the cumulative emissions space, resulting from larger contributions of changes in ocean and land carbon to the emission budget (see figure 1). CDR is potentially highly effective in mitigating near surface ocean acidification quickly, but deep ocean acidification is, as the perturbation in ocean carbon (figure 1(c)), long-lived.

Next, we focus on the spatial patterns of hysteresis in SAT and O₂ and how they vary with ECS. We compare spatial anomalies for the CO₂ down-path minus the up-path at cumulative emissions of 1 000 GtC and a CDR rate of 1 %yr⁻¹ (figure 4). 1 000 GtC are reached in simulation year 67, 68, and 69 on the up-path, and in simulation year 257, 249, and 231 on the down-path for ECS of 2 °C, 3 °C, and 5 °C, respectively. We again consider large levels of overshoot above the baseline of 1 000 GtC.

Zoom In Zoom Out Reset image size

The hysteresis of SAT is largely negative for an ECS of 2 °C and becomes positive for higher ECS (figures 3(a) and 4(a), (c) and (e)). A clear polar amplification is visible across all three ECSs. Highest positive temperature differences are located over Alaska and Eastern Siberia and in the Bering and Chukchi Seas, and in the Atlantic sector of the Southern Ocean. These regions show reduced sea-ice area and thus a lower albedo for the down-path compared to the up-path (see figure 2(d)).

Thermocline O₂ shows spatially varying positive and negative hysteresis. Hysteresis of O₂ is negative in the Arctic, the Atlantic south of the equator, and most of the Pacific. The Atlantic north of the equator and the northern Pacific show large positive hysteresis in thermocline O₂ (figures 4(b), (d) and (f)). The pattern is generally comparable between different ECSs. In parts of the North Pacific and in the Southern Ocean, however, thermocline O₂ shows positive hysteresis for ECSs of 2 °C and 3 °C, while for ECS = 5 °C the hysteresis is negative. The hysteresis patterns result from the interplay of consumption of O₂ during remineralization, ventilation of water masses, and contributions from temperature anomalies to the solubility component of oxygen, O$_{2,\textrm{sol}}$. Hysteresis in thermocline O$_{2,\textrm{sol}}$ is generally negative, except in the Southern Ocean for ECS of 2 and 3 °C with slightly positive values (not shown). The most negative O$_{2,\textrm{sol}}$ differences coincide with the O₂ hysteresis minima in the Arctic and the Atlantic south of the equator.

The Indian Ocean, the tropical Western and Northeastern Pacific, the Atlantic north of the equator and the Atlantic sector of the Southern Ocean show minima in ideal age hysteresis (not shown), indicating higher ventilation of waters in the thermocline. Directly around Antarctica, water mass age increases for ECS = 5 °C. On top, export production exhibits most negative hysteresis in the North Atlantic in response to a weaker AMOC and thus reduced nutrient supply, and most positive hysteresis around Antarctica as a result of reduced sea-ice cover. Together, reduced water mass age and less export production lead to less oxygen consumption (i.e. Indian Ocean, Atlantic north of the equator, parts of the Pacific), while less ventilation and increased export production as seen around Antarctica yield lower O₂ on the down-path than on the up-path for ECS of 2 °C and 3 °C (figures 4(b) and (d)). For ECS = 5 °C, and to a lesser extent ECS = 3 °C, negative hysteresis in O$_{2,\textrm{sol}}$ and positive hysteresis in ideal age around Antarctica lead to the sign change in O₂ hysteresis there (figures 4(d) and (f)).

Finally, in the case of $\Omega_{\textrm{arag}}$, differences between the three ECS are small (figure 3(d)). Ω_arag is generally lower during rising CO₂ than during CDR for the same cumulative emissions (figures 3(d) and 5). The recovery is strongest in the sub-tropics, where $\Omega_{\textrm{arag}}$ reaches values $\gt$3 in most areas (figure 5(b)). The tropics show slightly higher $\Omega_{\textrm{arag}}$ values for CDR compared to emissions, while the high latitudes are similar for the two time slices (figure 5).

Zoom In Zoom Out Reset image size