Multidecadal variability in the climate system: phenomena and mechanisms

An open and fundamental issue in climate dynamics is the origin of multidecadal variability in the climate system. Resolving this issue is essential for adequate attribution of human-induced climate change. The purpose of this paper is to provide a perspective on multidecadal variability from the analysis of observations and results from model simulations. Data from the instrumental record indicate the existence of large-scale coherent patterns of multidecadal variability in sea surface temperature. Combined with long time series of proxy data, these results provide ample evidence for the existence of multidecadal sea surface temperature variations. Results of a hierarchy of climate models have provided several mechanisms of this variability, ranging from pure atmospheric forcing, via internal ocean processes to coupled ocean-atmosphere interactions. An important problem is that current state-of-the-art climate models underestimate multidecadal variability. We argue that these models miss important processes in their representation of ocean eddies and focus on a robust mechanism of multidecadal variability which is found in multi-century simulations with climate models having a strongly eddying ocean component.


Introduction
The climate system, with its atmosphere, oceans, cryosphere, land surface and biosphere, displays variability on a multitude of time scales [26]. Some of these time scales, say the daily and seasonal variability, are directly related to processes external to the climate system, such as variations in solar insolation. These external processes are therefore considered as a forcing to the climate system resulting in forced variability. However, variability also arises through processes internal to the climate system, mainly because of the existence of feedbacks; such variability is called internal or intrinsic climate variability. Both internal and forced variability are referred to as natural variability, as these would exist without the presence of humans on Earth. Climate variability is also caused by emissions of greenhouse gases and aerosols due to human activities; such variability is referred to as anthropogenic climate variability, or simply climate change. a e-mail: h.a.dijkstra@uu.nl (corresponding author) One of the fundamental issues in climate research is to unravel the physical mechanisms of intrinsic climate variability on 10-100 years time scales [51], which we indicate in this paper as MV (Multidecadal Variability; a list of acronyms is provided in the "Appendix"). The most important reason to study MV is that its longer time scale makes it difficult to distinguish it from climate change, hampering the attribution of present-day climate and weather phenomena (e.g., exceptional regional droughts) to atmospheric greenhouse gas increases. Indeed, in the absence of MV, global temperature records (e.g., as shown in Sect. 2) would have increased more smoothly with increasing greenhouse gas concentrations. The clearer causal relationship would likely have led to faster mitigation measures (e.g., emission reductions) much more easily. A prominent example of the complications in attribution was the recent so-called 'hiatus' in global mean surface temperature rise following the decade after the 1997/1998 El Niño-Southern Oscillation (ENSO) event, where MV played a substantial role [47].
The instrumental sea surface temperature (SST) record, spanning about 150 years, is too short to adequately detect and understand the processes causing MV, but pronounced variability on these time scales has been found. However, as these datasets show a strong non-stationarity, in particular during the last 40 years, detrending techniques are necessary to apply traditional statistical analysis methods, such as principle component (PC) analysis [53], and results depend on the particular detrending method chosen [19,62]. The significance of the particular time scale and pattern can therefore not be adequately demonstrated. Proxy data, such as records from ice cores or sediment cores, are often long enough to demonstrate significance of a particular time scale, but lack the resolution to determine spatial patterns, or to investigate cause-effect relations [51].
To determine the physical processes causing MV, models are a necessary tool. A hierarchy of models has been applied to study MV, from Conceptual Climate Models (CCMs) described by, for example, low-dimensional ordinary differential equations to General Circulation Models (GCMs) described by systems of partial differential equations. MV has been found in many of these models, and a range of mechanisms have been suggested to explain this variability. CCMs usually identify a specific transparent mechanism which is subsequently often hard to link to GCM output and observations. Analysis of GCM output usually leads to several different and less transparent storylines of mechanisms, which are subsequently often hard to connect to basic geophysical fluid dynamics principles. One of the current problems is that MV is clearly underestimated in state-of-the-art GCMs (Coupled Model Intercomparison Project phase 5; CMIP5) compared with observations [3,44]. This underestimation might point to the fact that GCMs are highly susceptible to errors associated with respect to their representation of small-scale processes [36].
Although an enormous amount of work has been done to characterize MV, there is no consensus on the physics behind it. The aims of this paper are to (i) create some 'beacons in the haze' and (ii) to propose a new robust mechanism of MV due to internal ocean dynamical processes. In Sect. 2, a brief summary of analyses of MV from observations is presented. In Sect. 3, suggested mechanisms of MV based on modeling results are categorized. In Sect. 4, we describe a new robust mechanism of MV generation caused by the presence of mesoscale eddies in the ocean. A summary, discussion and outlook follow in Sect. 5.

Observations
Prominent regions of MV in SST are the North Atlantic, the North Pacific and the Southern Ocean. We will first discuss results on global quantities and then turn to these specific regions.

Global variability
The global mean surface temperature (GMST) is an important variable in the climate system, as it integrates many large-scale processes. It is the chosen measure to communicate climate change to the general public and policy makers. We here consider GMST from the HadCRUT4 [50] dataset and SST from the HadISST dataset, which is provided on a 1 • × 1 • grid [54]. The SST data are deseasonalised by removing the mean annual cycle (at each location) and detrended with the two-factor detrending method where an estimate of the contributions due to natural climate forcing and due to the anthropogenic forcing is removed [18]. We estimate these contributions from the CMIP5 ensemble average of the 'natural forcings only' experiment (orange dashed line in Fig. 1a) and the difference between the 'all forcings' and the 'natural forcings only experiments' (green dashed), respectively. The 'all forcings' estimate was extended from 2006 to 2018 with the RCP8.5 scenario and with the absence of major volcanic outbreaks the 'natural forcings' were assumed constant.
As can be seen in Fig. 1b, clear signatures of MV can be found, with several periods where the trend in GMST was negative. Again, the short instrumental record limits accurate determination of any dominant periods. However, MV in the GMST has been detected over the last two millennia in a large number of proxy data [51]. Regressing the detrended SST anomalies in Fig. 1b onto the multidecadal GMST signal results in significant positive correlations throughout most of the tropics and extending along the Eastern Pacific (Fig. 1c). Indeed, this pattern is reminiscent of the first PC of the global SST identified as a hyperclimate (SST-) mode [14] which is forced by local air-sea interaction and exists because of the large oceanic heat capacity and atmospheric teleconnections. Another global variability mode has been suggested with the stadium wave [36].

North Atlantic
The first analyses of MV in North Atlantic SSTs [17,37,58] indicated the existence of variability on a time scale of 50-70 years, which we will refer to here as Atlantic Multidecadal Variability (AMV). Following Enfield et al. [16], an AMV index is calculated as the area average SST in the region [0 • N,60 • N] × [80 • W,0 • E] which is subsequently 13-year low-pass filtered with a Butterworth filter. Figure 2a shows both the deseasonalized monthly values and the low-pass filtered index. Recent warm periods were from the 1930s to the 1960s and from 1995 up to present day, whereas during the 1970s to the mid-1990s the North Atlantic surface ocean was relatively cold.
Because of the relatively short observational time series of the instrumental record (∼ 150 year of SST), it is difficult to extract a dominant pattern of multidecadal variability with much confidence. Kushnir [37], Kaplan et al. [33] and Delworth and Greatbatch [7] present a reconstruction of a pattern with a ∼ 50-year period which shows a near standing signal in SST. This SST pattern is basin wide with the largest anomalies appearing south of Greenland. Figure 2b shows the global regression pattern of the detrended SST fields on the AMV index from the HadISST data set, indicating the characteristic horseshoe pattern with smaller anomalies near the subtropical North American coast. In regional patterns, there is a more pronounced negative SST anomaly near the coast of North America and a positive SST anomaly over the rest of the North Atlantic basin [37].
Low-frequency variability in the North Atlantic SST has been determined from proxy data stretching back at least 300 years [9], and within these data there is a statistically significant peak (above a red noise background) at about 50 years. Significant multidecadal peaks in the (a) (b) (c) Fig. 1 a The HadCRUT4 Global Mean Surface Temperature (blue) and the CMIP5 multi-model ensemble means of the 'natural forcings only' experiment (dashed orange) and the estimate of the anthropogenic forcing GMST signal derived as the difference between the 'all forcings' and 'natural forcings only' CMIP5 experiments. b The annual two-factor detrended GMST (grey) with a 13 year low-pass filtered curve to isolate MV where 7 years are removed at each end to avoid filter edge effects (blue). c Regression map of the detrended annual SST anomalies on the low-pass filtered GMST anomaly signal. Purple dashed lines demarcate the areas of significant correlations at the 98% level. Significance in the regression maps was established as in Jüling et al. [30] spectrum were, for example, found in Greenland ice-core data [4], where five overlapping records between the years 1303 and 1961, were available with annual resolution. Mantua et al. [45] analyzed monthly mean SST data for the North Pacific (20-70 • N) and showed that the PC of the leading Empirical Orthogonal Function (EOF) displays MV. Subsequent analyses have indicated that the SST field has robust MV statistical modes, where the first EOF is usually referred to the Pacific Decadal Oscillation or the Pacific Multidecadal Variability (PMV). The second EOF is usually called the North Pacific Gyre Oscillation [13] and the PC also displays MV. The PMV index used here is the low-pass filtered first principal component of the Pacific SSTs north of 20 • N [45] and a plot of the this index from the detrended HadISST dataset [54] is provided in Fig. 3a.
The global regression pattern in Fig. 3b shows the characteristic cold anomaly extending from the Kuroshio region into the mid-Pacific and the warm anomalies surrounding it, also seen in the global data (Fig. 1c). Over much of the 1950s to 1970s, the western (eastern) North Pacific was warmer (colder) than normal. The situation reversed in the late 1970s but returned at the end of the 1990s. The change in the sign of the PMV in the late 1970s is often referred to as the Pacific Climate shift [64] while the negative tropical SST anomalies after the year 2000 contributed to the GMST trend slowdown [35]

Southern Ocean
Data are sparse in the Southern Ocean, and it is therefore difficult to identify MV in this region. Indications for MV come from satellite-derived sea ice data [59], which covers only a 34-year period (1979-2012). Long-term trends in Antarctic sea-ice season length are found in both the Ross Sea and Bellingshausen Sea. Another signature of MV may be the occurrences of Maud Rise Polynyas, i.e., open ocean areas in the sea ice of the Weddell Sea [23].
One can define an SMV index over the region [80 • S, 50 , the choice being motivated by high-resolution model results [31]. Although the detrended HadISST data shows little variability prior to 1940, probably due to sparse data, the SMV index displays MV (Fig. 4a). This indicates that the Weddell Sea was warmer than normal over the period . This result depends on the detrending method used but even without detrending MV is found in the Southern Ocean SSTs [73]. For the SMV index, the global regression pattern (Fig. 4b) is mostly confined to the Southern Ocean with largest positive regression coefficients in the Weddell and the Ross Seas and negative values in the Amundsen and Bellingshausen Seas.

Interpretation framework
Note that although MV is shown through the filtered time series in the figures above, and large-scale patterns seem to be associated with MV, it is a priori unclear whether there are large-scale processes generating MV as a quasi-periodic signal [48]. Such a view has been advocated for many climate variability phenomena, for example ENSO at interannual time scales, based on the clear physical oscillation mechanisms proposed [12]. Alternatively, interactions involving small-scale processes can be responsible for MV leading through an (inverse) energy cascade to regimes of (power-law) scaling between the amplitude of the variability and the time measure [42].
Often the spectral power, S(ω), of a certain time series X (t) displays a power-law distribution where ω indicates the frequency. In this case, there is a specific relation between fluctuations at different frequencies, given by the power-law exponent η. A special case for which this happens is a self-similar process [21], for which X (λt) = λ H X (t) in the distribution sense, and H is the Hurst exponent; in this case η = 2H + 1.
To determine scaling behavior in the time series, one considers the behavior of fluctuations X ( t) versus t, t being a multiple of the minimum (i.e., sampling) time scale. The fluctuations are said to display scaling behavior when the fluctuation function S 2 behaves as where the brackets indicate ensemble mean. When H < 0 and hence η = 2H + 1 < 1, fluctuations will decrease with scale t, and when H > 0, and η > 1 fluctuations will increase with time scale t [41]. An analysis of a wide range of observational time series has led to the identification of different scaling regimes. Up to about 5-10 days, a power-law spectrum with η ≈ 2 is found, giving H ≈ 0.5, which is the Hurst exponent of Brownian motion. Up to about 50 years, the macroweather regime is characterized by η ≈ 0.2, giving H ≈ −0.4. In the climate regime, up to 50 kyr, a positive Hurst exponent H ≈ 0.2 is found, associated with η ≈ 1.4 [43]. Other regimes at even larger time scales have also been identified [41].
One interpretation of these results is that apart from some specific forced responses, such as the annual and daily cycle, and those associated with volcanic eruptions, all internal climate variability is contained within these scaling regimes [42]. The variability then may originate from multifractal cascading processes, such as those in turbulent flows [41]. Whereas there is no dispute that these processes lead to fluctuation variations in the weather regime, other processes may be responsible for the scaling behavior at larger time scales [20]. Moreover, variability associated with specific phenomena, such as ENSO, may not be part of this scaling regime [56].
Rypdal and Rypdal [56] show, for example, that when the scaling analysis is applied to interstadial and stadial ice core data separately, the scaling exponent η ≈ 1, H ∼ 0. They suggest that the so-called 1/ f or pink noise (η = 1) scaling holds over both the macroweather and climate regimes. As another example, Moon et al. [49] analyze the GISS monthly averaged SST data from 1901 to 2012 and proxy data from 32 paleoclimate data sets using multifractal temporally weighted detrended fluctuation analysis. They also find scaling characteristics of pink noise at multidecadal time scales. The associated spatial pattern shows a connection between the PMV region, the Eastern Pacific and the Southern Ocean.

Models and mechanisms
In Sect. 3.1, a short overview is given of the model hierarchy that has been used to study MV. In Sect. 3.2, mechanisms derived from studies using conceptual to intermediate complexity climate models will be summarized. In Sect. 3.3, we will discuss the MV-relevant results of the analysis of simulations of the latest large-scale global climate model inter-comparison project (CMIP5).

The climate model hierarchy
Scales and processes are important properties of climate models, motivating classification of models [11] using these two traits (Fig. 5). Here the trait 'scales' refers to both spatial and temporal scales as there exists a relation between them: Smaller spatial scales are usually associated with faster processes. 'Processes' refers to either physical, chemical or biological processes that occur in different climate subsystems.
Conceptual climate models (CCMs; lower left hand corner of Fig. 5) represent only a limited number of processes and scales. In these models, only very specific interactions in There are, of course, overlaps between the different model types, but for simplicity they are sketched here as non-overlapping the climate system are described. For example, models of glacial-interglacial cycles can be formulated by small-dimensional systems of ordinary differential equations [6]. Here, only the interactions between ice sheet volume, atmospheric CO 2 concentration and global mean ocean temperature are included.
Scales can be added to CCMs by spatially discretizing the governing partial differential equations up to three dimensions, even with only a limited number of processes. The atmosphere, ocean, ice and land domains are divided into grid boxes. Over each grid box the budgets of momentum, mass and heat are formulated as differential equations. In a class of climate models that are known as Earth System Models of Intermediate Complexity (EMICs, [5]; centre right in Fig. 5) a larger number of scales are represented than in CCMs. For example, the atmospheric model may consist of a quasi-geostrophic or shallow-water model and the ocean component may be a zonally averaged model. The advantage of EMICs is that longer time scale processes, such as land-ice and carbon-cycle processes, can be included. Individual component models of an EMIC may also be used to study the interaction of a limited number of processes. Such models are usually referred to as Intermediate Complexity Models (ICMs; centre left in Fig. 5). A prominent example is the Zebiak-Cane model of the ENSO phenomenon [71].
GCMs having a higher spatial resolution and with more processes represented are located in the upper part in Fig. 5, in particular those used in CMIP5 and CMIP6 simulations. When they represent other features than the ocean/atmosphere circulation, such as biogeochemistry, they are usually referred to as Earth System Models (ESMs). In time, the ESMs of today will be the EMICs of the future and the state-of-the-art ESMs will shift toward the upper right  [25] or advective ocean [57], b atmospheric forced view [15], c noise excited normal mode view such as in Te Raa and Dijkstra [63] and Weng and Neelin [70]. d MV through coupled ocean-atmosphere interactions [39] corner in Fig. 5. The early ESMs (before CMIP5, so including the models used in CMIP3 simulations) are considered here as the present-day EMICs.

CCS and EMIC results
There have been numerous studies into MV with reduced complexity models. We choose here to focus on examples of model analyses which have led to key hypotheses on the origin of MV in the climate system. The first conceptual model addressing the origin of low-frequency variability of SST due to atmospheric high-frequency heat flux forcing was presented almost 50 years ago [25]. An ocean mixed layer has a slow response to an approximately white noise atmospheric forcing and can be integrated to give a so-called red noise response in SST with elevated power in the MV frequency range; there is no preferred frequency. As an extension of the Hasselmann [25] model, Saravanan and McWilliams [57] have considered the effect of a fixed advective ocean and show that increased MV power can be obtained through a process of spatial resonance (Fig. 6a). For the specific case of Atlantic MV, Eden and Jung [15] used an ocean-only model forced by observed atmospheric winds and claim that the AMV is forced by the atmospheric MV (Fig. 6b).
A class of ocean models in which AMV was found to be purely due to intrinsic oceanic processes is that of the idealized single hemispheric basin models [2,28,63]. Here an ocean flow is only forced by an atmospheric heat flux which creates a meridional temperature gradient in the surface ocean. This drives an Atlantic Meridional Overturning Circulation (AMOC) due to meridional density differences. The AMOC is not stable under the prescribed atmospheric heat flux and leads to an oscillating flow on multidecadal time scales. The instability is clearly related to an internal mode of variability of the ocean which is usually referred to as a Thermal Rossby Mode (Fig. 6b). In Delworth et al. [8], MV was found in the GFDL-R15 model and connected to variability in the AMOC. Later, Delworth and Mann [10] showed that the noise in the atmospheric freshwater forcing was crucial to allow for MV to occur (Fig. 6b) and that the multidecadal time scale resulted from AMOC variability. Knight et al. [34] demonstrated the connection between AMOC variability and AMV in GCMs and an extensive review on this issue can be found in Zhang et al. [74].
Weng and Neelin [70] study the variability arising through coupling of an ocean basin and atmospheric noise and find that this coupling can lead to increased power at MV time scales (Fig. 6c). The underlying pattern is an ocean basin mode, which is quite different from a Thermal Rossby Mode. Weijer and van Sebille [69] perform simulations with the CCSM4, with focus on the PMV and suggest that it arises through an excitation of Rossby basin modes by atmospheric noise (Fig. 6c). Work with earlier GCMs [38,39] tried to convey the idea that MV was caused by midlatitude ocean-atmosphere coupling. Here, feedback processes between the atmosphere and ocean are thought to be crucial for MV to occur (Fig. 6d).

CMIP5 models
The Coupled Model Intercomparison Project phase five (CMIP5) archive contains results from GCMs which have been used in the 5th assessment report of the Intergovernmental Panel of Climate Change (IPCC). The archive contains two types of model simulations which are relevant here. The first type are long (multi-centennial) control simulations under constant (pre-industrial) forcing (solar, greenhouse gases, aerosols) and the second type are simulations under historical forcing conditions, usually over the period 1850-2010. An overview of both types of simulations can be found in Cheung et al. [3], where the number of historical ensembles and length of the control simulation is also shown. This simulation length is in most cases over 500 years, and hence, MV due to intrinsic processes can be simulated.
From several CMIP5 models, there is clear evidence that the AMV index displays MV [24]. The six models in Han et al. [24] that show the 'best' agreement (in terms of correlation, amplitude and spatial pattern) indicate a dominant period of 20-70 years. Of these models, the GFDL-CM3 model is able to simulate the negative SST anomalies near Newfoundland, a feature found in observations but difficult to capture in models. This model actually displays two significant periods of variability, one around 25 years and one around 50 years.
Spectra for the PMV in the CMIP5 control and historical simulations were summarized in the review by Newman et al. [52]. There is no sign that the PMV index displays variability beyond a red noise background. Oceanic Rossby waves, atmospheric noise, shifts in the Kuroshio Current and reemergence of SST anomalies due to ENSO are all considered to play a role in the PMV [52].
Reintges et al. [55] analyzed a variety of CMIP5 GCMs and found that several models show deep convection events in the Southern Ocean at multidecadal timescales, for example the GFDL-CM [72]. Reintges et al. [55] argue that models with a weakly (strongly) stable stratification tend to have a shorter (longer) re-occurrence time of deep convection and demonstrate the importance of sea-ice volume on the average length of both non-convective and convective periods.
Cheung et al. [3] investigate intrinsic variability both from the control and historical simulations, where they filter the forced signal in the latter simulations by using the multimodel ensemble mean. Their results show that many CMIP5 models clearly underestimate intrinsic variability, both in the North Atlantic and in the North Pacific, a point further stressed by Mann et al. [44] by looking at MV in the atmospheric surface temperature field.

Robust MV: collective interactions
In all CMIP5 model results discussed above, the spatial resolution of the ocean model component does not allow for ocean eddies to form and interact. Eddies arise usually in ocean models through a mixed barotropic/baroclinic instability mechanism, whose proper representation requires a horizontal resolution smaller than the local internal Rossby radius of deformation, which is about 25-50 km in the midlatitude oceans. In CMIP5 models, the ocean eddies are not captured but their net effect on the transport of heat and salt is captured in most cases through the GM [22] parameterization. However, the nature of diffusive parameterizations, such as GM, may act to damp MV in low-resolution models [46,67].

MV in eddying ocean models
Early simulations with adiabatic eddying models, such a multi-layer quasi-geostrophic models and shallow-water models have indicated that MV can occur through a collective interaction between ocean eddies and the time-mean flow. In the midlatitude wind-driven ocean circulation problem, such variability was studied and indicated by the 'Turbulent Oscillator' [1]. Collective interactions between eddies can lead to low-frequency variability in the surface ocean circulation relevant to, for example, the bi-modality of the Kuroshio. In Southern Ocean channel models, MV variability has also been found [27] and explained through the overall energy balances.
In Le Bars et al. [40], MV was found in a multi-century simulation with the Parallel Ocean Program (POP). It was named the Southern Ocean Mode, an example of SMV, and it is remarkably regular with a dominant 45 year period. Such variability, but with a slightly smaller period (about 25 year), was also found in a multi-century ocean eddying version of Community Earth System Model (CESM) as described in van Westen and Dijkstra [66]. The representation of mesoscale flows may be the key to capturing MV adequately in GCMs. Jüling et al. [30] show that both AMV and PMV are better represented (both temporally and spatially) when ocean eddies are explicitly captured in the CESM ocean model component.
Including ocean eddies also changes temporal scaling properties in the MV regime. In Fig. 7, results from multifractal temporally weighted detrended fluctuation analysis [32,49] on 250 model year time series of the AMV, PMV and SMV are shown for both the eddying and non-eddying version of CESM [30]. Here, the fluctuation function F 2 (s) is shown versus time scale s, and the scaling exponent (F 2 (s) ∼ s 1−H , where H is again the Hurst exponent [49]) are indicated by dashed lines. In particular, for the SMV in the multidecadal regime, the non-eddying model clearly shows a white noise scaling regime (H ∼ 0.5) whereas in the eddying model, pink noise is found (H ∼ 0).
To understand the mechanism of this emerging MV, one often resorts to energetic descriptions where the production and dissipation of and conversion between available potential energy and kinetic energy is studied over a specific ocean domain. A full energy analysis can be performed when a reference state of minimal potential energy can be computed [29]. Alternatively, an approximate method is the Lorenz Energy Cycle (LEC) approach [68]. An LEC analysis of the POP results in [40] was performed by Jüling et al. [31], supporting the basic elements of the SOM mechanism, already provided in Hogg and Blundell [27] and Le Bars et al. [40].

Conceptual model
The energy analysis provides a macroscopic view of MV where variability emerges through the interactions of eddies with the mean flow. We define K m , K e , P e and P m to be the volume integrated mean kinetic energy, eddy kinetic energy, mean available potential energy and eddy available potential energy. The equations of the LEC are then given by as sketched in Fig. 8a. Here, the generation of X is denoted by G(X ), conversion from X to Y by C(X, Y ) (and hence C(Y, X ) = −C(X, Y )) and dissipation of X by D(X ); the boundary terms (see [31]) are not considered here. Explicit expressions of the terms can be found in von Storch et al. [68] and Jüling et al. [31]. Following Sinha and Abernathey [60], the terms in (3) can be approximated to develop conceptual models of MV. In the following paragraphs, we introduce three models that describe the transfer of energy between the different reservoirs. These transfers can generate an oscillatory response due to phase differences between the different terms. Conditions for oscillatory behavior and the corresponding frequency can be derived. Typical values will be calculated for the Southern Ocean using the POP model output [40], noting the a priori limitations of the LEC analysis that were extensively discussed in Jüling et al. [31]. First of all, convection occurs which affects the energy budgets but is not well captured in the LEC. Second, the dissipation terms were not calculated explicitly in Jüling et al. [31], but adjusted to obtain equilibrium energy budgets. Despite these limitations, the LEC results can be used to estimate where conditions of multidecadal oscillatory behavior occur.
The first step in deriving the conceptual models is to approximate the full LEC in (3) as follows [60]: • The baroclinic instability process dominates eddy production, which occurs in the C(P e , P m ) terms, so that the Reynolds' stress related terms C(K e , K m ) is neglected; • The production of energy is dominated by the input to K m (such as the wind work), so that G(K m ) is the only significant generation term; • The dissipation of potential energy, D(P m ) and D(P e ) is neglected; and • It is assumed that dK m /dt ≈ 0.
These approximations can be well justified using the POP results [31]. With P ≡ P m + P e , Eq. (3) then reduce to

Model 1
The key term in (4) is the conversion term, C(P e , K e ) which represents the production of eddy kinetic energy through baroclinic instability. This term can, under conditions of adiabatic eddy stresses and Gent-McWilliams closure (see derivation in Sinha and Abernathey [60] with the supporting references), be approximated by C(P e , K e ) = α P, where α is a constant.
, and using E = K e we then obtain the reduced equations (model 1) Here, the 'useful' wind work W [60] and the dissipation of eddy kinetic energy D are given by (Fig. 8b) where τ is the stress and the subscripts b and s indicate bottom and surface, respectively. In addition, u is the horizontal velocity field and tilde and bar indicate eddy and mean flow, respectively. The dissipation term, D, is generated by bottom drag, which is often represented by a linear or quadratic drag law and hence D can be approximated to be a function of E (D ∼ −r E). W is also dependent on E, as an increase in E indicates an enhanced baroclinic instability, a more meandering mean flow and hence a lower useful wind work (when the wind stress τ s is fixed), so as a first approximation, one can take W (E) ∼ β E. Estimates from Sinha and Abernathey [60] provide α = K GM /L 2 ≈ 10 −9 s −1 where K GM ≈ 1000 m 2 s −1 is the Gent-McWilliams isopycnal diffusion coefficient and L ≈ 10 6 m a typical meridional length scale. Using the POP results [31,40] for the region SO30 (south of 30 • S), values of W , D, C, E and P normalized to the interval [−1, 1] are plotted in Fig. 9a. Linear regressions of W vs E, C vs P and D vs E (which have large error bars because POP data consist of only one MV cycle) give estimates of α = 7.7 × 10 −9 s −1 , r = 2.0 × 10 −7 s −1 and β = −1.9 × 10 −7 s −1 . The value of α appears large compared with the estimate of Sinha and Abernathey [60] and the value of β is probably too small. In this conceptual model, feedback processes drive an oscillatory flow. A low potential energy is associated with a smaller slope in isopycnals leading to less baroclinic instability to develop (phase A in Fig. 9b); hence, the minimum of E follows the minimum of P. The mean flow changes are determined by rectification and hence W already increases while E is still decreasing (which is measured by β). The increase in W leads to a stronger slope of the isopycnals and hence P increases also when E is still decreasing (phase B in Fig. 9b). Once isopycnal slopes are sufficiently strong, baroclinic instability occurs leading to an increase of E (phase C in Fig. 9b). Energy is drawn from the potential reservoir and the flow becomes less zonal such that W decreases before E becomes maximal. The isopycnal slopes reduce and P decreases, followed by a decrease in E (phase D in Fig. 9b).

Model 2
The results plotted in Fig. 9 show that C is not a sole function of P. The Gent-McWilliams isopycnal diffusion coefficient, K GM , may depend upon both P and E through eddy feedbacks. An extension of model 1 is to consider (as suggested by [60]) the closure C(P e , K e ) = k P E 1 2 , in which case the model equations of model 2 become Again linearizing around a mean state (P,Ē), we obtain a Jacobian matrix where the coefficients are given by a = −kĒ Complex conjugate eigenvalues, with period T 2 , are found when Using the POP results [31,40], again for the region SO30 (south of 30 • S), the relation between C/P vs E 1 2 provides a rough estimate of k = 2.0 × 10 −18 m −1 . With estimates P = 2.4 × 10 20 J andĒ = 10 18 J from the POP results, the period T 2 is plotted versus β as the red curve in Fig. 10. For this model, the value of β 2c = −2.28 × 10 −6 s −1 shifts to less negative values, but the range of MV in β is still very limited. However, the result is sensitive to the choice of the parameters used and can shift either way (to smaller or larger values of β).

Model 3
A third model is motivated by the apparent delay of a few years between the useful wind work and the eddy kinetic energy (Fig. 9). When we model the complex processes setting this delay by a single constant delay time δ, hence W (t) = β E(t − δ), then Eq. (5) of model 1 can be used to derive model 3: where the time argument is now made explicit. The characteristic equation of this model [61], is given by where the complex numbers σ are the eigenvalues; for δ = 0, this reduces to model 1. The effect of the delay (of 1 and 2 years) on the period of oscillation versus β is plotted in Fig. 10 (magenta and black curves). Compared to the earlier models 1 (Eq. 5) and 2 (Eq. 9), the range of β for which MV occurs is much broader and longer periods occur at smaller values of β, in particular for δ = 2 years.

Summary, discussion and outlook
The focus of this paper was on the understanding of one of the fundamental and open issues in climate science: multidecadal variability (MV) in sea surface temperature (SST), in particular in the North Atlantic, North Pacific and Southern Ocean. Although the instrumental record cannot fully characterize the SST patterns, together with proxy data there is ample evidence that MV exists, certainly in the North Atlantic and North Pacific. A first issue is whether MV can be attributed to specific large-scale processes, leading to quasi-periodic behavior (such as is the case for ENSO) or whether it is part of the background that emerges through the interaction of small-scale processes. In the quasi-periodic case, we have seen that studies with a hierarchy of models have suggested a number of possible mechanisms for MV, going from atmosphere-forced, via internal ocean variability to coupled ocean-atmosphere variability. Specific examples of mechanisms are the thermal Rossby wave for the AMV [63] or the noise-affected Rossby basin modes for the PMV [69]. One would hope that such mechanisms could be falsified by observations, but that is not the case yet. In an alternative paradigm, MV is characterized by scaling behavior in the climate system [42] or by specific stochastic processes [56]. It is interesting that multidecadal variability is just on the boundary of the macro-weather regime, with a break in the power-law spectral behavior and actually seems to follow pink noise behavior [49].
A possible mechanism for generating MV is the action of small-scale processes, i.e., eddy mean-flow interaction. This mechanism has been presented with specific focus on the SMV. Here, on the macroscopic scale, large-scale oscillations occur because the small-scale processes interact with the larger scales through instabilities and rectification. When these processes are out of phase, an oscillation appears which can be described by changes in the integral energy budgets of the potential and kinetic energy reservoirs and the fluxes between them.
The conceptual models introduced here, using crude approximations to the energy equations, illustrate that multi-decadal regimes of variability exist once there are eddy-induced changes in useful wind-energy input due to changes in the surface flow. The main process here is rectification, i.e., the effect of the mean flow due to the nonlinear interactions of eddies. It is interesting that a 1-2 years delay can lead to a long-term variability in the Southern Ocean, and this is likely connected to memory processes associated with eddy-mean flow interactions [46]. The differential delay model is an interesting prototype to describe the effects of these interaction processes and, while it is an ad-hoc model at the moment, its justification can probably be done using weakly nonlinear theory [65].
A current problem of concern is that many CMIP5 models underestimate MV with possible serious consequences for projections of future climate under greenhouse forcing [44]. CMIP5 models do not capture the processes at the Rossby deformation radius which are responsible for the MV in eddying ocean models [31,40]. One would then expect that MV is better represented in climate models with eddying ocean components, as is indeed found for CESM [30]. As we have shown here, the scaling regime of MV changes with CESM ocean model resolution, tending to pink noise for the eddying model. Such a change in scaling can possibly be explained by studying stochastic differential delay models.
In this respect, only an incremental improvement can be expected regarding MV in GCMs in the CMIP6 project. A new generation of eddying climate models, with at least 0.1 • degree horizontal ocean model resolution and maybe also with better parameterizations of submesoscale processes, is needed to capture MV much better than the CMIP5/6 models. It is when these models are available that the new leap in understanding of MV is expected to occur. As a final remark, considering the great importance of the availability of highresolution climate model simulations, it is remarkable that computing facilities for running such models are not more easily available to the wider community. Ideally, CMIP6 should not have been a minor extension of CMIP5, but should have induced a major leap in modeling the climate system. PC: Principle component PMV: Pacific multidecadal variability POP: Parallel ocean program SST: Sea surface temperature SMV: Southern ocean multidecadal variability