Beyond Forcing Scenarios: Predicting Climate Change through Response Operators in a Coupled General Circulation Model

Global Climate Models are key tools for predicting the future response of the climate system to a variety of natural and anthropogenic forcings. Here we show how to use statistical mechanics to construct operators able to flexibly predict climate change for a variety of climatic variables of interest. We perform our study on a fully coupled model - MPI-ESM v.1.2 - and for the first time we prove the effectiveness of response theory in predicting future climate response to CO$_2$ increase on a vast range of temporal scales, from inter-annual to centennial, and for very diverse climatic quantities. We investigate within a unified perspective the transient climate response and the equilibrium climate sensitivity and assess the role of fast and slow processes. The prediction of the ocean heat uptake highlights the very slow relaxation to a newly established steady state. The change in the Atlantic Meridional Overturning Circulation (AMOC) and of the Antarctic Circumpolar Current (ACC) is accurately predicted. The AMOC strength is initially reduced and then undergoes a slow and only partial recovery. The ACC strength initially increases as a result of changes in the wind stress, then undergoes a slowdown, followed by a recovery leading to a overshoot with respect to the initial value. Finally, we are able to predict accurately the temperature change in the North Atlantic cold blob.


I. INTRODUCTION
Climate change is arguably one of the greatest contemporary challenges for humanity [1]. The provision of new and efficient ways to understand its mechanisms and predict its future development is one of the key goals of climate science. Global climate models (GCMs) are currently the most advanced tools for studying future climate change; their future projections are key ingredients of the reports of the Intergovernmental Panel on Climate Change (IPCC) and are key for climate negotiations [2]. For IPCC-class GCMs, future climate projections are usually constructed by defining a few climate forcing scenarios, given by changes in the composition of the atmosphere and in the land use, each corresponding to a different intensity and time modulation of the equivalent anthropogenic forcing. Typically, for each scenario an ensemble of simulations is performed, with each member differing in terms of initial conditions, applied forcing or chosen physical parametrizations. Subsequent phases of the Coupled Model Intercomparison Project (CMIP, now the sixth phase CMIP6 is active [3]), which is part of the Program for Climate Model Diagnosis and Intercomparison (PCMDI), allowed to define standardized experimental protocols ("scenarios") for numerical simulations with several GCMs and for the evaluation of the model runs [4,5].
A bottleneck of this approach is that the consideration of an additional forcing scenario requires running a new ensemble of simulations. Additionally, for each forcing scenario, it is hard to disentangle the impact of each component contributing to the scenario, e.g. different greenhouse gases with their concentration pathways and land surface alterations in geographically distinct regions. Finally, no rigorous prescription exists for translating the climate change projections in the case one wants to consider different time modulations of a given forcing scenario, e.g. a faster or slower increase of CO 2 .

A. Elements of response theory
A possible strategy to deliver flexible and accurate climate change projections is the construction of response operators able to transform inputs given by forcing scenarios into outputs in the form of climate change signal. In this regard, it is tempting to use the fluctuation-dissipation theorem (FDT) [6], which provides a dictionary for translating the statistics of free fluctuations of a system into its response to external forcings. The idea of using the FDT to predict climate change from climate variability has been first proposed by [7] and used by several authors thereafter [8][9][10]. The FDT has recently been key to inspiring emergent constraints, which are tools for reducing the uncertainties on climate change by looking at empirical relations between climate response and variability of some given observables [11,12] Unfortunately, the applicability of the FDT to systems outside equilibrium faces some serious theoretical and practical challenges; see discussion in [13][14][15][16].
The climate is a non-equilibrium system whose dynamics is primarily driven by the heterogeneous absorption of solar radiation. The motions of the geophysical fluids with the associated transports of mass and energy, as well as the exchanges of infrared radiation, tend to re-equilibrate the system and allow it to reach a steady state [16][17][18]. In the case of nonequilibrium systems, forced fluctuations contain features that are absent from the free fluctuations of the unperturbed system [19]. In the specific case of the climate system, this means that climate change projects only partially on the natural modes of climate variability, whilst climate surprises -unprecedented events -are indeed possible when forcings are applied [14]. Nonetheless, the unperturbed state does contain useful information for predicting its response to perturbations. Response theory is a generalisation of the FDT that allows one to exploit this information via explicit computation of response operators able to predict how general systems -near or far from equilibrium, deterministic or stochastic -change as a result of forcings. After the pioneering work by [6], response theory has been firmly grounded in mathematical terms for stochastic [20] and deterministic [19,21,22] systems.
A detailed treatment including a discussion on the conditions for the applicability of response theory in the context of climate modelling can be found in [23] and [15] and in the Materials and Methods section. Here we introduce the minimal ingredients that will be necessary in the following. Let us consider a dynamical system described by the state vector x , evolving according toẋ = F (x) and an observable of interest Φ = Φ(x), for example the global surface temperature.
We apply to the system an additional forcing to the right-hand side of the equations of motion as Ψ (x , t) = X (x )f (t), characterized by a pattern in the phase space X and a time modulation f . Response theory gives a way to compute the response of the expectation value Φ f (t) of the observable, where · indicates the expectation value over an ensemble of independent realizations of the system -in a deterministic dynamical system, a sampling starting from many different initial conditions [24,25].
The linear part of the response (see the Materials and Methods section) is given by the convo- where G (1) Φ is the linear Green function of Φ for the forcing field X. As discussed in the Materials and Methods section, G Φ can be conveniently computed for any observable Φ with a particular choice [15,23] of the time modulation f : a step function, as in the case of a classic CO 2 doubling experiment. Once G (1) Φ is known, equation 1 allows to compute the response to any other time modulation of the forcing, without the need to repeat the numerical simulations.
Response theory is thus able to treat in a unified and comprehensive way forcings with any temporal modulation, ranging from instantaneous to adiabatic changes, and to perform climate change projections that are relevant also at local scale. Furthermore, different sources of forcings can be treated independently. In the linear response regime, they can be treated in an additive way.
One can then construct flexible response operators able to ease some of the limitations of classical climate change scenarios modelling exercises, and allow for a more efficient use of climate change data [15].

B. Response theory and climate change
The applicability of response theory in the context of climate science has been successfully tested in systems of various degrees of complexity, ranging from the Lorenz 1963 model [13], to the Lorenz 1996 model [26,27], to simple quasi-geostrophic models [10], up to intermediate complexity climate models like PLASIM [15,23]. Through response theory one can put the concept of equilibrium climate sensitivity (ECS) on solid theoretical grounds, able to clarify how to interpret state-dependent sensitivity. Additionally, response theory allows one to generalize at all time scales the concept of transient climate response (TCR) [28], which is defined as the change in the globally averaged surface temperature at the end of the ramp of increase of CO 2 concentration growing at 1% per year up to doubling (cfr. [23]). Some properties of the climate response operators can be associated with the features of the variability of the unperturbed system; this is especially clear in the case of resonant response [16,29,30]. The climate models used so far to test response theory [15,23] lacked an active and dynamic ocean, so that the multiscale nature of climate processes was only partially represented. Capturing the slow oceanic processes is essential for a correct representation of the multidecadal and long-time climatic response. Encouragingly, response theory has recently been shown to have a great potential for predicting climate change in multi-model ensembles of CMIP5 atmosphere-ocean coupled GCMs outputs [31]. Blending together data coming from different models is outside response theory theoretical framework, yet heuristically justifiable. However, a proper treatment within the boundaries of the theory, based on ensemble of simulations with the same model featuring an active ocean component, is still lacking.
In this paper we construct response operators from an ensemble of climate change experiments using a coupled atmosphere-ocean climate model, the MPI-ESM v.1.2 [32]. We perform two sets of experiments. One is a classic CO 2 doubling experiment (2xCO2), from which we compute the Green functions of several observables of interest, as described in the Materials and Methods section. We then perform an ensemble of experiments with a different forcing scenario: the 1% per year increase in the CO 2 concentration until doubling (1pctCO2). We use the Green function computed with 2xCO2 to reconstruct the response of 1pctCO2, comparing the results of the prediction with the direct numerical simulations. The goal is to prove the effectiveness of response theory for performing climate change projections using complex climate models.
First, we analyse the response of the globally averaged near-surface temperature (T 2m ), investigating TCR and ECS and linking these quantities to the response on all time scales. We then focus on two classical indicators of the large-scale ocean circulation, namely the strength of the Atlantic Meridional Overturning Circulation (AMOC) and the strength of the Antarctic Circumpolar Current (ACC), and give evidence of the ability of response theory to predict the slow modes of the climate response. We also address the coupling between the ocean and the atmosphere, analysing the net energy exchange at surface between the two fluids, usually referred to as global ocean heat uptake (OHU). This is by definition the change in global ocean heat content [33]. A non-vanishing OHU indicates the presence of a global net energy imbalance. Indeed, in current conditions, the ocean is well-known [34] to absorb a large fraction of the Earth's energy imbalance due to global warming and to store it through its large thermal inertia, up to time scales defined by the deep ocean circulation. Finally, we prove the validity of the approach for a local climatic variable, examining the surface temperature change in the North Atlantic, where the ocean deep water formation takes place. This region features a complex interplay between radiative exchanges and meridional energy transport, thus being particularly sensitive to the strength of the forcing and the changes in the large-scale circulation.
In what follows we discuss the climate change projections obtained using linear response theory. Details on the computation of the Green functions and on the procedure used to derive them are reported in the Materials and Methods section.

A. Global Mean Surface Temperature
In terms of globally averaged surface temperatures, the 2xCO2abrupt and 1pctCO2 (cfr. not captured at all by models featuring a non-dynamic ocean [15,23]. The importance of the slow modes of climate response, associated with the oceanic thermal inertia, can be quantified considering the ratio between T CR and ECS. Here we have T CR/ECS ≈ 0.5 (ECS ≈ 3.5K), which is much smaller than what found (≈ 0.85) in [23], indicating a very prominent role of the slow modes of variability.
As discussed in [23], the structure of the Green function ( Figure 5) significantly departs from an exponential relaxation behavior, differently from what assumed in simplistic one-time scale models of the response of the system [35]. Here, the presence of slow oceanic time scales makes the simple one-time-scale scenario even more insufficient. After a sharp decrease, the Green function decreases at a much slower pace in the range 70 y -400 y. As said, this range is responsible alone for about half of the total final warming. Inspecting the time series of the response, one realises that it is hard to model the complex properties of the Green function with a simple sum of two exponentials; see also [30].
The prediction of the 1pctCO2 T 2m change is able to capture accurately both the fast response

B. Atlantic Meridional Overturning Circulation
The AMOC is strongly influenced by buoyancy perturbations in the Atlantic basin [36]. It is relevant at climatic level because it encompasses a large portion (25%) of the total (atmo-spheric+oceanic) meridional heat transport [37]. On longer time scales (>100 y), a negative feedback acts as a a restoring mechanism, associated with a positive sign in the Green function. The presence of fast (meaning here decadal) response associated with the GHG forcing has already been found in other models [38,39], and is most likely related to the timescales of the sea-ice melting, consistently with paleoclimate simulations of the last interglacial climate with prescribed freshwater influx from reconstructed sea-ice melting [40]. The slow recovery of the AMOC might be understood in terms of the heat advection feedback [41][42][43].
In the 1001-2000 years period, teh response theory shows that a steady state is progressively reached over multi-centennial scales. The newly established AMOC is significantly weaker than the unperturbed AMOC, although a large ensemble spread is found. The results that have been predicted via response theory are remarkably consistent with simulations obtained from higher resolution versions of the same model [44], intermediate-complexity models [39] and other fully coupled models with comparable resolution but also including an interactive carbon cycle [45], indicating that the AMOC is simulated to never recover to the unperturbed value in a warmer climate.

C. The Antarctic Circumpolar Current
The ACC is by far the strongest large-scale oceanic current and its role in the general circulation is two-fold. On one hand, it isolates Antarctica from the rest of the system, being associated with a very strong zonal circulation in the Southern Ocean and the overlying atmospheric synopticscale eddies. On the other hand, although eminently wind-driven, it marks the area of outcropping of deep water occurring at the southern flank of the subtropical gyre, as part of the global-scale overturning circulation.
The Green function is shown in the inset of Figure 6b. We found that the initial strengthening of the ACC can be associated with an increase in surface zonal wind stress (not shown here). Such a surface forcing determines an enhanced Eulerian mean ACC transport, consistently with previous low resolution simulations [46]. On decadal scales, we have a loss in the correlation between wind stress and ACC, corresponding to the Green function turning negative after about 30 years.
Beyond these time scales, we have an agreement in the response of the AMOC and ACC, which indicates a coherent response of the global ocean circulation. Other models [47,48] also feature such a behavior on intermediate time scales, consistently with the idea that the two circulations are related via the thermal wind balance [49].
The prediction of the ACC strength evolution in the 1pctCO2 scenario is rather accurate for the first 1000 years, except for an underestimation of the positive short-term response, which is smoothed out. This points to an insufficient ability of response theory in representing such scales, given the complex coupling between surface wind stress and downward momentum transfer. We observe the presence of a strong variability (on decadal time scales) of the predicted signal. This might result from either the use of a too small ensemble size or, more interestingly, could be the signature of the natural variability encoded by a so-called Ruelle-Pollicott resonance, according to recently described mathematical properties [16,29,30], as mentioned in the introduction.
Note that in the 1001-2000 yrs period the ACC reaches an approximate steady state a bit later than the AMOC, possibly as a result of having a larger inertia, consistently with the different depth scales of the two currents. The AMOC maximum overturning depth scale is indeed located at about 1 km, whereas the outcropping in the Southern Ocean is related to isopycnal surfaces reaching much deeper [49,50]. This has profound implications for setting the time scales of the ACC and AMOC response. The propagation of deep water formation anomalies in the Northern Hemisphere is in fact mediated by Kelvin waves in the Northern Atlantic, whereas much slower interior adjustment through Rossby waves communicates the anomaly to the Southern Ocean, as analysed in detail in [51].

E. The North Atlantic cold blob
Finally, we study the surface temperature response for a domain covering the North Atlantic, thus including the areas where the deep water formation occurs, see Figure 4. The domain is comprised between 26W and 53W degrees of longitudes, and 53N and 69N degrees of latitudes.
The region is identified as a peculiar spot for the effects of the GHG forcing, since the sea-ice melting has been hypothesized to delay the surface warming of this area compared to surrounding regions, through a weakening of the overturning circulation [53]. Indeed, the surface warming over the North Atlantic region is remarkably different from the behavior of the rest of the extratropics, which features a time dependent response (not shown here) similar in shape but somewhat amplified with respect to the global mean depicted in Figure 1. Indeed, a long-lasting plateau -a hiatus in the temperature increase -is observed around the end of the CO 2 increase ramp in the North Atlantic. The plateau is well captured by response theory, and comes in agreement with the AMOC weakening [53] predicted in Figure 2a. This result is non-trivial, given that such local response results from an interplay of local factors (sea-ice melt mechanisms around Iceland [54]) and, as mentioned, the response of the large-scale oceanic circulation. This hints at the potential of using response theory in order to identify suitable global-scale quantities that can be used as predictors for the response of local observables, see the mathematical framework in [30]. Figure 4. Comparison between 1pctCO2 T 2m anomalies in the North Atlantic cold blob window (in K) as simulated (in thick red, with orange shadings denoting the ensemble error) and reconstructed via response theory prediction (in K).

III. DISCUSSION AND CONCLUSIONS
We have shown here that response theory is a valuable tool for predicting crucial aspects of the response in the climate system to a prescribed forcing. Predicting the response relies on obtaining the observable-dependent Green function from model simulations with a forcing having the same spatial pattern as the one of interest and a conveniently defined temporal modulation. The Green function allows one to deal with a continuum of time-dependent forcings, beyond the standard use of reference scenarios. Our findings provide guidance to the climate modellers' community on how to set up the experimental protocols aimed at defining climate change scenarios, minimising the need for computational resources. Note that, as an example, one can define exact relationships between quantities like equilibrium climate sensitivity and a generalised measure of the transient climate response, see [23], [15] and [16].
It is important to remark that the qualitative properties of the response of the investigated climatic observables are vastly different. In some cases the response is monotonic, in other cases not.
In some cases the overall long-term sensitivity is positive, in others negative. In all considered cases response theory successfully predicts the time-dependent change.
We have made here use of a fully coupled model, highlighting the slow components of the response associated with the oceanic modes of variability. We argue that these are responsible for about half of the global mean temperature response, pointing towards the crucial role of the thermal inertia by the oceans. The presence of a vast range of active time scales in the system makes the prediction of the response theoretically challenging and practically extremely relevant.
The predicted changes in the AMOC and ACC feature clearly distinguishable fast and slow regimes of response. The former is essentially different in the two circulations, being the ACC subject to the effect of surface wind stress anomalies (substantially underestimated, compared to actual simulations). The latter was found to be well correlated between AMOC and ACC, as a signature of the forced response of the global ocean circulation circulation, which they are both part of. Consistently with previous findings [49,50], the ACC reaches a steady state much later than the AMOC. The plateau in near-surface warming over the North Atlantic is also related to the slow response of the buoyancy-driven circulation. Indeed, the AMOC and ACC initial reduction is likely the consequence of the sea-ice albedo feedback and the related anomalies in North Atlantic freshwater influx.
Note that, for illustrative purpose and consistency with the scenarios used in the CMIP phases, we dealt with forcings characterised by changes in CO 2 only. This means that the pattern of the forcing is determined by the spatially homogeneous CO 2 mixing ratio. Beyond that, response theory has been proposed [55] as a tool for framing geoengineering strategies and understanding its limitations. The next step is to investigate other kinds of forcings, also exhibiting nontrivial spatial patterns (e.g. aerosol forcings, land-use change, land glaciers location and extension).
Response theory has been shown to provide rigorous tools for finding functional relations between the response of different observables of system to forcings, in the spirit of the theory of emergent constraints. This allows to treat comprehensively the concept of feedback across different time scales and define causal links between different variables [30]. This is another promising application of response theory to climate change studies.
We have only addressed here the first-order (linear) response, but the theory is also applicable to higher orders [13]. Some hints on the non-linear component of the response could be also driven by comparing the response of the systems to abrupt forcings with different magnitude.
Finally, we remark that the possible failure of the prediction with response theory with other observables and/or forcings has itself fundamental implications for the knowledge of the dynamics of the system one is studying. It has been highlighted that in the vicinity of critical transitions the linear response operator diverges as a result of the increase in the time correlations of the system [29], signalling the crisis of the chaotic attractor [56]. The experimental design here provided is thus also a clear and mathematically sound strategy for the study of tipping points and their role for the climate response [16,57] in state-of-the-art climate models.

Retrieval of AMOC and ACC
Typically, the large-scale circulation in the ocean is measured in terms of the mass transport across a suitably chosen section of a basin. The strength of the AMOC is computed as the vertically integrated mass weighted meridional mass streamfunction across the 26.5 • of latitude [65]. This is provided as a standard diagnostic in the outputs of the MPI-ESM model. The ensemble mean unforced strength of the AMOC amounts to 17.3±1.1 Sv, which is consistent with recent available measurements from the RAPID monitoring array [66]. The ACC is roughly zonally symmetric, and its location is closely related to the isopycnal slopes in the Southern Ocean. Traditionally, it has been measured in terms of the strength of the mass transport across the Drake passage.
Similarly, we take the vertically integrated barotropic streamfunction difference between the 2 • × 2 • boxes centered around the 68 • W, 54 • S and 60 • W, 65 • S locations [47]. The ACC at the beginning of the simulations amounts to roughly 138 Sv, which is consistent with the multi-model mean estimate found in [47], amounting to 155±51 Sv. It is also not far from the value commonly used as benchmark for the assessment of climate models (173 Sv [67]).

Linear response theory
Response theories allow one to predict how the statistical properties of a system changes as a result of acting modulations in its external or internal parameters. The validity of the corresponding response formulas is heavily dependent on the hypothesis that the unperturbed system is structurally stable, i.e., roughly speaking, far from bifurcations, or, in terms of geophysical systems, from tipping points (see related discussions in a climate context [15,23,27]). Rigorous derivations of response theories have been provided for the case of deterministic [19,21,22] and stochastic [20] dynamics. We only remark here that statistical mechanical arguments encoded by the chaotic hypothesis [68] (a non-equilibrium analogue of the ergodic hypothesis) indicate the feasibility of the methodology proposed here.
In this paper we follow to a large extent the approach presented in [15] and [23] (see also [13,27]) for the study of a large ensemble of intermediate-complexity atmospheric model runs and follow the deterministic route for response theory [19,21,22]. Let us consider a dynamical system described by the state vector x , whose dynamics is described by the set of differential equationsẋ = F (x ). We add a perturbation vector field of the form Ψ (x , t) = X (x )f (t), where X is the structure of the forcing in the phase space and f its time modulation. The expectation value of any observable Φ = Φ(x) can be written as: where Φ 0 is the expectation value in the unperturbed state, and the term Φ Φ is the linear Green function of the generic observable Φ. For ease of notation we have not indicated in equation A1 the dependence of the response on X, as in the applications considered in this paper X is fixed and only the time modulation f is varied. Note that for a time modulation f such that lim t→0 f (t) = f 0 , |f 0 | finite, and f (t) = 0 if t < 0, as in the case of where H is the Heaviside function (H(t) = 0 for t ≤ 0 and H(t) = 1 for t > 0), one typically has that Φ 1) f (0) = 0, as observed in this paper for all observables except the OHU. In this latter case, one has lim t→0 Φ f (t) = 0 because the Green function has a singularity (in the form of a Dirac's δ contribution) for t = 0 [30].

Procedure for the retrieval of the Green functions
The strategy for testing the prediction of the mentioned key variables with the coupled model ensembles is as follows. First, we compute the Green function from the 2xCO2 experiment. Conveniently, the time modulation of the forcing is given in this case by f (t) = f 2xCO 2 H(t), where H(t) is the Heaviside function, and f 2xCO 2 is a constant depending on the amplitude of the forcing.
The time array is set in such a way that the instantaneous doubling occurs at t = 0. Eq. A2 can be thus rewritten as: The outputs of the 2xCO2 experiments and the corresponding Green functions for the observables described above are presented in Figures 5-8. In particular, the time evolution of OHU for this scenario is shown in Figure 7. The positive forcing due to the instantaneous CO 2 doubling leads to an instantaneous jump in the OHU, leading to an annual average value of more than 1 PW in the first year. Given eq. A3, this indicates that the Green function G OHU has a singular behaviour at t = 0 (cfr. [30]), while a regular behaviour is found for t > 0, corresponding to the negative radiative Planck feedback.
For all observables, we then use the Green functions above to perform predictions for the 1pctCO2 scenario using Eq. A2. Since the radiative forcing is approximately proportional to the logarithm of the CO2 concentration, the time modulation of such forcing can be expressed as a ramp function (cfr. [15] and [23]): where the time scale τ ≈ 70 years denotes the time needed to reach the doubling in the CO2 concentration. Again, note that the forcing is such that X, i.e. the heating pattern induced by