Reversibility in an Earth System model in response to CO2 concentration changes

We use the HadGEM2-ES Earth System model to examine the degree of reversibility of a wide range of components of the Earth System under idealized climate change scenarios where the atmospheric CO2 concentration is gradually increased to four times the pre-industrial level and then reduced at a similar rate from several points along this trajectory. While some modelled quantities respond almost immediately to the atmospheric CO2 concentrations, others exhibit a time lag relative to the change in CO2. Most quantities also exhibit a lag relative to the global-mean surface temperature change, which can be described as a hysteresis behaviour. The most surprising responses are from low-level clouds and ocean stratification in the Southern Ocean, which both exhibit hysteresis on timescales longer than expected. We see no evidence of critical thresholds in these simulations, although some of the hysteresis phenomena become more apparent above 2 × CO2 or 3 × CO2. Our findings have implications for the parametrization of climate impacts in integrated assessment and simple climate models and for future climate studies of geoengineering scenarios.


Introduction
There is considerable uncertainty in future emission and concentration pathways of greenhouse gases. Mitigation scenarios aiming to avoid potentially dangerous climate change are likely to require significant reductions in greenhouse gas emissions within the next two decades (Mignone et al 2008, Rogelj et al 2011. However, it is conceivable that early CO 2 emissions reduction may not be achieved and the current rapid growth in greenhouse gas concentrations (Friedlingstein et al 2010) may continue for at least several decades. This creates an urgent need to understand whether the undesirable elements of climate 1 Present address: Laboratoire de Météorologie Dynamique, IPSL, CNRS/UPMC, Paris, France. 2 Present address: EUMETSAT, Darmstadt, Germany. change that occur during the global warming period can be reversed in a usefully short time if ever it becomes possible in the future to lower the atmospheric concentration of greenhouse gases through economically viable CO 2 air capture (Matthews 2010), possibly associated with a rapid decarbonization of the economy.
Earth's climate is a dynamic system whose components respond on different timescales to external forcings. While the atmosphere adjusts quickly to its boundary conditions, the response of the system as a whole is largely controlled by the ocean, the cryosphere, the land surface and some biogeochemical cycles, all of which respond over a much wider range of timescales. The longer timescales are associated with processes such as mixing into the deep ocean and ocean circulation, diffusion of heat and moisture in the deeper soil, vegetation growth and carbon cycling, and ice sheet formation and melting. These components introduce inertia into the system, leading to significant time lags between a forcing and the response it induces. For instance it is well known that in stabilization scenarios the global-mean surface temperature (GMST) response lags the radiative forcing because of the thermal inertia of the global ocean, a process known as constant concentration warming commitment (Wigley 2005, Hare andMeinshausen 2007). Scenarios where atmospheric concentrations of greenhouse gases exceed a particular target before decreasing towards that target are known as 'overshoot' or 'peak-and-decline' scenarios (Huntingford and Lowe 2007). Such scenarios may be characterized by hysteresis effects and irreversible changes. Irreversibility means that the system cannot be restored to its initial state or only does so on a timescale far longer than those normally considered practical from a human perspective. Hysteresis refers here to the dependence of the climate system not only on its current environment (measured either by the atmospheric CO 2 concentration or the GMST) but also on its past environment. In many cases this arises because there is a lag between the input variable (e.g., atmospheric CO 2 concentration) and the output variable (e.g., global-mean precipitation). In this study irreversibility implies hysteresis, but hysteresis does not necessarily imply irreversibility. Finally it should be noted that hysteresis effects do not necessarily lead to the existence of multiple equilibria in the system as we only consider transient simulations.
Reversibility in the climate system has been looked at for the carbon cycle in the context of zero-emission commitment scenarios either with an Earth System Model of Intermediate Complexity (Solomon et al 2009) or coupled carbon-climate models (Lowe et al 2009, Frölicher andJoos 2010). Reversibility has also been investigated in three-dimensional climate models for physical climate variables such as precipitation, sea-ice, sea-level rise or cloud feedbacks (Yoshida et al 2005, Frölicher and Joos 2010, Samanta et al 2010, Wu et al 2010, Cao et al 2011. The closest experimental setup to ours is that of Samanta et al (2010) who increased the CO 2 atmospheric concentration from 350 to 700 ppmv and decreased it back to 350 ppmv in an atmospheric model coupled to a slab ocean model. The recent advent of second generation of Earth System models such as the Met Office Hadley Centre HadGEM2-ES model (Collins et al 2011, Jones et al 2011 allows reversibility and hysteresis effects to be examined in a self-consistent framework and for a much wider range of Earth System components than has previously been possible. The objective of this study is to provide an overview of these effects in the HadGEM2-ES model and guide future research on specific components and processes of the Earth system.

Model description and experimental design
The Hadley Centre Global Environment Model version 2 Earth System model (HadGEM2-ES) is a state-of-the-art model that can be used to investigate how couplings and feedbacks may affect the response of the Earth System to climate forcing. It comprises a physical coupled atmosphere-ocean-sea-ice model with the addition of component models for atmospheric aerosols, tropospheric chemistry and the terrestrial and oceanic carbon cycle to characterize various aspects and couplings within the Earth system. Some processes, such as permafrost processes, ice sheets and glacier dynamics, and ground water recharge, are not or only partly included in the model, although it is recognized that they could play a role in reversibility issues. A detailed description of the Earth System model can be found in Collins et al (2011). The model includes a representation of the carbon cycle for terrestrial and oceanic ecosystems with interactive vegetation (five different plant functional types) and plankton distributions (diatoms and non-diatoms phytoplankton and zooplankton). The impact of dust fertilization on plankton growth is represented, and DMS production by plankton is a source of aerosols. The ozone and methane concentrations as well as the sulfate oxidant fields vary with the meteorology and the climate. Seven aerosols species are represented in the model, namely sulfate, biomass burning aerosol, black and organic carbon from fossil fuel burning, sea-salt, mineral dust, and biogenic aerosols. They interact with the atmosphere through direct and indirect effects and are coupled with other model components.
A traditional benchmark CMIP experiment consists of increasing the CO 2 atmospheric concentration by 1% yr −1 starting from an initial concentration of 286 ppmv which corresponds to that of year 1850. This experiment is generally used to investigate the transient climate response to the forcing by greenhouse gases. We build on this 1% CO 2 experiment by constructing five different scenarios. In the first scenario, we ramp-up the CO 2 atmospheric concentration by 1% yr −1 starting from the initial concentration (286 ppmv) and finishing at four times that concentration (1144 ppmv). Four ramp-down scenarios have been branched at different points along this increasing concentration. The CO 2 concentration is decreased at 1% yr −1 starting at four different concentration levels (4 × CO 2 , 3 × CO 2 , 2 × CO 2 and 1.5 × CO 2 ), giving a range of radiative forcing increases of up to around 7.5 Wm −2 covering much of the range spanned by the new policy-relevant simulations for the next IPCC assessment (Moss et al 2010). Several techniques have been proposed to remove carbon dioxide from the atmosphere, either through large-scale biomass growth coupled with carbon sequestration, modification of the ocean chemical composition, acceleration of the carbonate cycle or direct air capture. The potential for successful carbon dioxide removal using these techniques is uncertain, and the current proposed methods cannot achieve the level required to sustain a 1% decrease per year in atmospheric CO 2 concentration. However it is not uncommon to investigate unusual climate scenarios with more simple climate models. Performing such peak-and-decline scenarios with a full Earth System model is a good test of whether the simple models have the right parametrizations to operate under conditions of decreasing concentrations of greenhouse gases.
It should be noted that because the CO 2 atmospheric concentration is prescribed, the C budget is not closed. The CO 2 fluxes between the atmosphere and the ocean and between the atmosphere and the land surface are only diagnostics in these simulations. Nonetheless they can be used to understand how the carbon-climate feedback operates as described in Hibbard et al (2007) and Jones et al (2011). It is also recognized that the ramp-down part of the scenario is unrealistic as it would necessitate very large negative C emissions as discussed further below.

Results and discussion
3.1. Atmospheric component GMST increases at a near-constant rate when CO 2 concentration increases exponentially and is about 5.5 • C above the pre-industrial value at the time of CO 2 quadrupling (figure 1). When the CO 2 concentration starts to decrease, the temperature change reverses immediately over land but with a time lag over the ocean, as found previously for lower rates of greenhouse gas reduction. We find that the larger the CO 2 peak concentration, the longer is this time lag (it is up to 20 yr when CO 2 concentration peaks at four times the pre-industrial value). For the 2 × CO 2 ramp-up and ramp-down, the time lag is about 10 yr and is consistent with that found by Samanta et al (2010). The GMST change remains 0.5-1.2 • C above its initial value by the time CO 2 returns to its pre-industrial value (figure 1(b)). Since many components of the climate system are driven by surface temperature, we present further changes in the system as a function of the GMST but we recognize that some variables such as precipitation and continental ecosystems also respond to the CO 2 concentration (Cao et al 2011, Andrews et al 2010. We try to identify critical CO 2 concentrations and GMST changes beyond which the changes in some components either become irreversible or relax to their pre-industrial state on a much longer timescale. The variations in global-mean precipitation are dominated by changes over the oceans, while land precipitation changes are noisy and less systematic (figures 2(a) and (b)). The global-mean precipitation shows a clear hysteresis behaviour in all four experiments, confirming earlier findings using the HadCM3 model (Wu et al 2010, Cao et al 2011. For a given GMST change, precipitation is also greater on the ramp-down CO 2 pathway than on the ramp-up path. Two explanations for this phenomenon have been provided in previous studies. It has been shown that the precipitation response to an increase in CO 2 atmospheric concentration can be decomposed into a (positive) slow response proportional to surface temperature change and a (negative) fast adjustment to the CO 2 radiative forcing (Andrews et al 2010, Cao et al 2011. The hysteresis is due to the decoupling between the two terms that occurs during the ramp-down phase. It is also possible to explain the hysteresis as a response to a changing imbalance in the ocean surface energy budget (Wu et al 2010). Preference for one or the other explanation depends on whether one considers the hydrological cycle to respond to the atmospheric or surface heat budget. It is interesting that precipitation extremes (measured here as the land average of the local daily maximum precipitation rate in a year) show no hysteresis, pointing to a potentially higher degree of reversibility for hydrological climate impacts (figure 2(c)). During the ramp-up experiment (black), CO 2 concentration increases at 1% yr −1 starting at the pre-industrial level of 286 ppm. During the ramp-down experiments, CO 2 is prescribed to decrease at 1% yr −1 starting from four different levels: 4 × CO 2 (red), 3 × CO 2 (green), 2 × CO 2 (blue), and 1.5 × CO 2 (pink). The dashed lines show the approximate trends the global-mean surface temperature would have followed in the absence of any lag in the response to the CO 2 concentration, with the slope of the dashed lines being taken equal to the inverse of the slope of a linear fit to the black crosses.
Clouds also respond to the CO 2 forcing, with high-level cloud increasing and low-and mid-level cloud decreasing during the ramp-up phase, and a general tendency to opposite behaviour as the CO 2 is reduced (figures 2(d)-(f)). Such changes in the global-mean cloud amounts are well established (Mitchell and Ingram 1992) and have also been confirmed in more recent studies (Ringer et al 2006). Examination of the radiative impact of the individual cloud types (not shown) indicates a high degree of compensation between the net effects of the changes in high-and middle-level clouds in both the ramp-up and ramp-down phases. The radiative impact of the changes in low-level cloud is small because it primarily results from changes in optically thin clouds and smaller changes of opposite sign in optically thicker clouds. An interesting feature of the ramp-down phase is the lagged response of low-level clouds which show a Figure 2. Globally-and annually averaged model variables as a function of the CO 2 concentration (ppmv) and GMST (K) for (a) ocean precipitation (mm day −1 ), (b) land precipitation (mm day −1 ), (c) maximum daily precipitation rate over land (mm day −1 ), (d) high-level cloud fraction, (e) mid-level cloud fraction, (f) low-level cloud fraction, (g) change in global ocean heat content (10 24 J), (h) ocean carbon store (GtC), (i) surface ocean nutrient (mmol N m −3 ), (j) sea-ice cover (%), (k) average drought fraction (%), (l) land permafrost area (10 6 km 2 ), (m) vegetation carbon store (GtC), and (n) soil carbon store (GtC). Same colour code as in figure 1. small hysteresis effect. Further examination reveals that this is primarily an oceanic effect, depending, like the precipitation, on the balance between fast and slow responses and also on the geographical pattern of the SST change. In contrast, over land changes in all cloud types reverse almost immediately with surface temperature when CO 2 is decreased.

Oceanic component
Ocean heat uptake relates closely to the thermal expansion component of sea-level rise, which may dominate 21st century sea level increases . By the time the CO 2 concentration has risen to four times pre-industrial levels the model projects the ocean heat content (figure 2(g)) to have increased by 3 × 10 24 J, with a thermal expansion of about 40 cm. After the forcing declines from 4 × CO 2 the heat content continues to increase by a further 60% over the next century. The time between the peaks in forcing and ocean heat uptake varies strongly with the magnitude of the peak forcing, ranging from approximately 25 yr (1.5 × CO 2 ) to a century (4 × CO 2 ). This shows that only a rapid reduction in radiative forcing might be able to limit and reduce the thermal expansion component of sea-level rise before some of the most severe coastal impacts are realized. The annual mean Arctic and Antarctic sea-ice covers are entirely reversible up to 4 × CO 2 ; the recovery has a 50 yr lag over CO 2 but no lag over the GMST in line with earlier findings ( j)). This is in contrast with the results from a simple model used by Abbot et al (2011) who found that cloud and heat transport feedbacks could lead to bifurcations in the state of Arctic sea-ice, which would lead to possible hysteresis effects in the ramp-up and ramp-down scenario we consider. The Atlantic meridional overturning circulation (AMOC) seems to be reversible in all four cases, although there is a lag of 2-3 decades in the 3× and 4 × CO 2 ramp-down experiments. For these two cases, the AMOC weakening trend continues to undershoot further after CO 2 ramp-down before the recovery starts. In all cases, the AMOC recovers to its normal strength during the later stabilization period after CO 2 returns to the pre-industrial level (not shown).
Unlike the ocean heat content, the ocean carbon content shows a decrease soon after the forcing reverses (figure 2(h)). To first order, both the ocean heat and carbon contents are governed by atmosphere-ocean gradients and the exchange between the mixed layer and the deep ocean. Whilst throughout these simulations the surface ocean is colder than the atmosphere by the mixing up of cold water from depth, until atmospheric CO 2 concentrations become very high, mixing with the deep ocean is always enriching the surface waters in CO 2 . Deep mixing is therefore acting to increase the ocean-atmosphere temperature gradient, delaying the turnover in air-sea heat flux until towards the end of our simulations, but decreases the ocean-atmosphere CO 2 gradient, causing the ocean to become a source of CO 2 to the atmosphere only a few decades after the ramp-down starts. This results in a timescale for equilibrating ocean heat in the ramp-down that is much longer than that for equilibrating the CO 2 , at least in the ocean mixed layer. The ocean surface pH responds very quickly to atmospheric concentrations and is almost totally reversible in the global mean. However some of the impact of ocean acidification on marine ecosystems, for example on calcification (Riebesell et al 2000, Iglesias-Rodriguez et al 2008, may not be reversible in the same way, and may even feed back to limit the pH reversibility through changes in alkalinity not considered in our simulations. Deeper ocean pH and changes in the depth of the lysocline do not experience the reversibility experienced in the surface ocean. Major nutrient concentrations simulated in the surface ocean waters decreases almost linearly under the increasing greenhouse gas forcing. Nutrients available for photosynthesis (i.e., those at the surface ocean where light can penetrate) are expected to decrease under global warming (Behrenfeld 2006), and this occurs in the model as the ocean warms and stratifies, slowing the replenishment of nutrients from depth. An unexpected feature of these experiments is that, despite the mixed layer thickness in most of the ocean recovering in line with the removal of heat from the surface ocean, the globally averaged recovery of surface nutrient concentrations is very slow (figure 2(i)). The limited nutrient recovery results from a significant (>century) lag in the recovery of freshwater-driven stratification over much of the Southern Ocean, preventing upwelling of nutrient-rich waters from depth. This effect could have long-lasting consequences for marine productivity. The persistence of a low-salinity lid across much of the Southern Ocean can be understood as a simple function of the time lag between supply of freshwater through precipitation and removal of that water by mixing. Even if the high-latitude precipitation response to our forcing were to be perfectly reversible and simply scale with that forcing, since the salinity at any point in time integrates the changing freshwater flux over some interval prior to that point in time, and since the preceding precipitation is different to the present-day one during the ramp-up and ramp-down periods, a lag in surface salinity, independent of the complexities of the ocean physical circulation, is to be expected.

Terrestrial component
Soil moisture changes differ by depth. Moisture changes in the top (10 cm) soil layer responds linearly to changes and are completely reversible. Moisture changes in the bottom (1-3 m) soil layer shows much more of a lag when responding back to a decrease in CO 2 concentration because of the inertia associated with the deep soil hydrology. We calculated a drought diagnostic from the anomaly of soil moisture in the top 1 m of the soil. A soil moisture anomaly threshold was created for each grid cell based on the 20th percentile distribution of monthly soil moisture anomalies for the first 10 yr of the simulation. Any time the soil moisture anomaly falls below this threshold the grid cell is assumed to be in drought for that month. Figure 2(k) shows the proportion of the land surface that is in drought at any particular year (it is 20% on average during the first 10 yr of the simulation). The hysteresis arises because it takes time for the changes in the forcing meteorology to impact the soil moisture down to a depth of 1 m. The permafrost fraction over land (diagnosed from monthly soil temperatures) shows a hysteresis effect as a function of the GMST, with approximately 20% less permafrost area in the ramp-down than in the ramp-up experiment figure 2(l). Although permafrost extent is eventually fully reversible, any loss of organic carbon during the time thawed would be effectively irreversible as this has accumulated over many millennia. The various vegetation types in the model exhibit different behaviour when the atmospheric CO 2 concentration changes. Changes in vegetation fraction are not linear with temperature because they are driven by their carbon balance which depends on both climate and CO 2 fertilization (Gregory et al 2009, Denman et al 2007 and competition between vegetation types. When the CO 2 concentration increases, the bare soil fraction decreases while trees, shrub and grass fractions increase globally. Gross primary production (GPP) over land responds linearly with the logarithm of CO 2 in the ramping up phase. However both vegetation and soils continue to take up carbon many years after the atmospheric CO 2 concentration is diminished (figures 2(m) and (n)). For vegetation carbon, this is because of the long response time associated with vegetation competition. Trees continue to increase in the NH high latitudes after the peak in CO 2 and temperature and this overwhelms the decrease in the fertilization effect. For soil carbon, this is because of the long lag between the increase in soil carbon pools and the increase in heterotrophic respiration. Figure 3 shows the time-latitude Hovmöller diagrams of surface pH, surface nutrient concentration, vegetation carbon and soil carbon. There is a small hysteresis in pH ( figure 3(a)). In the Arctic for CO 2 concentrations of 300-400 ppmv (roughly the first and last 40 yr), pH is on average 0.06 smaller in the ramp-down phase. While the global-mean surface ocean pH is almost totally reversible, regional variations are to be expected because of changes in sea-ice and freshwater budget as underlined by Steinacher et al (2009). Much of the global-mean reduction in surface nutrient concentration seen in figure 2(i) is due to large reductions seen at high latitudes ( figure 3(b)) caused by fresher, more stable, surface waters. At lower latitudes the zonal signal is dominated by reductions in upwelling regions. Figures 3(c) and (d) show the terrestrial carbon stores. In the tropics the terrestrial carbon stores are fully reversible within the timescale of the simulation. The timescale of reversibility is longer at higher latitudes and these regions are mostly responsible for the hysteresis seen in figures 2(m) and (n). Carbon store changes at lower latitudes are caused by increases/decreases in the carbon density of forests, while at higher latitudes the increase is caused by an increase in forest area. Soil carbon peaks before vegetation carbon, because soil carbon varies with NPP which is closely linked to the atmospheric CO 2 concentration, while the peak in vegetation carbon is delayed by the slow expansion in the area of forest. Figure 4 shows the positive and negative emissions (GtC yr −1 ) necessary to close the CO 2 atmospheric budget in the ramp-up and ramp-down phases. It can be observed that the peak decrease is only about half the peak increase. This is because there is still some natural uptake after the peak so the negative emissions are offset by some continued terrestrial and oceanic increase. The ramp-down negative emissions are almost flat (especially when ramping down from 4 × CO 2 ). This is because the offsetting from natural uptake lessens and then reverses while the effort to bring atmospheric concentrations down keeps decreasing as the 1% profile levels off (i.e., the absolute rate of ramp-down decreases in time). The cumulative emissions amount to 3250 GtC to get to the 4 × CO 2 peak. The atmospheric increase is about 1820 GtC, which corresponds to an average airborne fraction of 56%. On the way down, cumulative negative emissions are about −2400 GtC so a smaller equivalent amount of CO 2 has to be taken out than was initially introduced. A fact relevant to CDR technology is that any negative emission will be subject to an airborne fraction in the same way as positive emissions. During the ramp-down phase, the average airborne fraction is 76%, which is significantly greater than during the ramp-up phase (as the natural sinks are still responding at the beginning of the ramp-down).

Summary and conclusions
We have investigated hysteresis and irreversibility in the Earth system using a highly idealized experiment in which the atmospheric CO 2 concentration is first ramped up and then ramped back down to its initial level. Although the scenario is unrealistic because of the requirement for large negative C emissions in the ramp-down phase, it is useful to investigate hysteresis and irreversibility in the Earth system and can provide insight into a hypothetical scenario of carbon dioxide removal. In summary only the ocean surface pH responds almost immediately to the atmospheric CO 2 concentrations in our experiments. Other quantities exhibit some time lag relative to the change in CO 2 and most also exhibit a lag relative to the GMST, which we describe as a hysteresis behaviour in CO 2 or GMST space. For certain quantities, such as precipitation, the hysteresis in GMST space occurs as a response to asynchronous variations in CO 2 and GMST. For ocean heat content or permafrost area, the hysteresis arises because the quantity is not in equilibrium with GMST and depends on processes operating on longer timescales. Longer experiments would be required to discriminate hysteresis phenomena due to inertia from those reflecting a bi-stability or multiple equilibria in the system. The most surprising responses are those in low-level clouds and ocean stratification in the Southern Ocean, which both exhibit hysteresis on timescales longer than expected. We did not see evidence of critical thresholds in these simulations, although some of the hysteresis phenomena become more apparent above 2 × CO 2 or 3 × CO 2 . We summarize in table 1 potential implications of our results for climate impact studies. Overall it appears that the most significant hysteresis in the HadGEM2-ES model is for variables relating to terrestrial and marine ecosystems. In contrast physical changes in the model show little hysteresis (except for sea-level rise). The latter is consistent with the findings from Samanta et al (2010) who found that physical climate parameters return to their initial values soon after the CO 2 concentration does in a similar 2 × CO 2 ramp-up and ramp-down experiment.
It is common to parameterize climate impacts in integrated assessment models and simple climate models as a function of the global-mean surface temperature change. While such a parametrization may be justified by the large uncertainties in damage cost functions and the monotonous and phased increase in CO 2 concentration, surface temperature and climate impacts in most scenarios, our results bring evidence that climate impacts cannot be parameterized as a function of GMST alone in unconventional scenarios such as those with CO 2 overshooting. A geoengineering scenario where the CO 2 atmospheric concentration peaks at levels much higher than today's before decreasing rapidly may experience delayed climate impacts. Determining more precisely which impacts may be irreversible will require even more complex representation of the Earth system with a focus on components which are potential sources of irreversibility, such as ice sheets, soil carbon in permafrost regions and terrestrial and oceanic ecosystem biodiversity.