1 Motivation

A central issue in climate research in recent years has been the apparent paradox that global surface temperature has not been increasing in tandem with increasing emissions of human-induced greenhouse gases. This paradox has generated a lot of attention among researchers (e.g., Easterling and Wehner 2009; Foster and Rahmstorf 2011; Katsman and van Oldenborgh 2011; Santer et al. 2011, 2014; Trenberth and Fasullo 2013; Huber and Knutti 2014; Maher et al. 2014; Risbey et al. 2014; Watanabe et al. 2014) and the general public (Tollefson 2014). It has also been used to cast doubts about the reliability of climate research in general and climate models in particular (Showstack 2014) since the ensemble mean of the models do not reproduce the so-called surface temperature hiatus.

A large suite of possible factors influencing the recent hiatus in the global surface temperature have been identified, including: the possibility of too-high model sensitivity to increased greenhouse gas (GHG) forcing (Flato et al. 2013); decadal-scale ocean uptake and storage of heat in the Pacific Ocean (e.g., Meehl et al. 2011; Kosaka and Xie 2013; Meehl et al. 2013; Trenberth and Fasullo 2013; England et al. 2014; Hu et al. 2015), in the Atlantic Ocean (Chen and Tung 2014; Drijfhout et al. 2014), in the Indian Ocean (e.g., Lee et al. 2015), and in the Southern Ocean (Drijfhout et al. 2014); decreased solar activity (Schmidt et al. 2014); uncertainty in the natural and human-induced input of aerosols and particles (e.g., Fyfe and Gillett 2014; Santer et al. 2014); and reduced concentration of stratospheric water vapour (Solomon et al. 2010).

In addition, Cowtan and Way (2014) argues that available global surface temperature records underestimate the recent global warming signal due to sparse data coverage in the Arctic region. This was also supported by Saffioti et al. (2015), where a dynamical adjustment method was used to show that the wintertime cooling reported during recent years has been overestimated, partly due to missing data and partly due to anomalous atmospheric circulation. In a recent reevaluation of observed surface temperature records, Karl et al. (2015) conclude that the rate of global warming during the last 15 years has been as fast as that seen during the latter half of the last century. The surface temperature hiatus might therefore be of a smaller magnitude than commonly reported or even statistically absent.

In the pioneering days of global climate modelling, Charney et al. (1979) noted that surface warming from GHG emissions might not be uniform, due to absorption and storage of heat in the ocean. The notion of decadal-scale periods with non-increasing surface temperature under global warming is therefore far from new. Based on observations of the global surface air temperature record and CMIP3 model projections, Easterling and Wehner (2009) found that decadal-scale temperature hiatus periods are not unprecedented in the observational record and will likely also occur in the future. Roberts et al. (2015) even suggest that the resent hiatus period could extend for another few years. Maher et al. (2014) continued the analysis from Easterling and Wehner (2009) to include model projections from phase five of the Climate Model Intercomparison Project (CMIP5; Taylor et al. 2012). They found that the occurrence of hiatus periods in the future depends on the rate of anthropogenic forcing and on the strength and occurrence of volcanic eruptions.

Inspired by the decadal analysis of Easterling and Wehner (2009), we use the latest suite of state-of-the-art climate models to address the likelihood of short- and medium-term (up to 30 years) temperature variations, particularly non-warming or cooling periods during longer-term (multi-decadal to century-scale) warming. As an addition to the Maher et al. (2014) analysis, we also include two additional forcing scenarios from the CMIP5 database, the RCP 2.6 and RCP 6.0 (Sect. 2.2). Central to our analysis is the location, magnitude and duration of decadal-scale heat anomalies in the coupled atmosphere-ocean climate system, and whether the simulated variability resembles observation-based variability. To shed light on where the heat may go in hiatus-like periods we address variations in global surface and upper and deep ocean temperature with particular focus on decadal time scales.

2 Data and models

2.1 Observations

Our analysis make use of three publically available observation-based reconstructions of global surface temperature. For global land and ocean surface temperatures: NASA Goddard Institute for Space Studies (GISS) (Hansen et al. 2010); and comparable compilations by the Met Office Hadley Centre and the Climatic Research Unit at the University of East Anglia, (HadCRUT3; Brohan et al. 2006, and HadCRUT4; Morice et al. 2012). For land surface temperatures, we use HadCRUT4.

2.2 Models

Seventeen global climate models participating in CMIP5 have been used (see Table 1). We use the four Representative Concentration Pathway (RCP) scenarios RCP2.6, RCP4.5, RCP6.0 and RCP8.5, which start in 2006 (Meinshausen et al. 2011). These scenarios include annual prescriptions of the solar forcing, atmospheric concentration of the major greenhouse gases, tropospheric ozone, stratospheric ozone, aerosol abundances, and land-use patterns (Taylor et al. 2012).

Table 1 Overview of the model systems, by institution, acronym and experiment, used in this study
Table 2 The correlation between decadal trends of the global surface temperature and the listed variables after removal of long-term (>31 years) variations in global surface temperature and the listed variables for each class of model simulations
Table 3 Summary of the minimum ocean heat values (\(10^{21} \hbox { J}\)) from Fig. 5 (at zero lag), presented as an average for 0–700 m, depths below 700 m, and the full ocean depth

Two additional model experiments are used. One is the so-called control integration (hereafter piCtrl), where the models are run with fixed, pre-industrial (i.e., year 1850) composition of atmospheric greenhouse gases and particles, and with no volcanic eruptions and a constant solar forcing. In addition, the models are run with observed composition of atmospheric greenhouse gases and particles, including the radiative effect of volcanic eruptions and solar variability for the period 1850–2005 (hereafter Hist).

For the analysis of hiatus periods, we have randomly picked one realization from each model for each of the six model experiments to the extent these are available (Table 1). Surface temperature, including land surface and sea surface temperatures, TOA incoming solar and TOA incoming and outgoing long-wave radiation, and global ocean temperature over the full depth range are used. Land and sea surface temperatures, the annual mean net TOA radiation, and the heat storage at different ocean depth intervals (upper 100, 300, 500 and 700, and the deep ocean) are derived from these data.

Table 4 The number of periods based on surface temperature forming Warm, Cold, Pre-Warm and Pre-Cold across all models and from the six types of model runs
Table 5 Globally integrated, ensemble mean ocean heat anomalies (\(10^{21} \hbox { J}\)) from piCtrl for the decadal-long Warm and Cold phases, and the 3–7 year long Pre-Warm and Pre-Cold phases
Table 6 The magnitude of Warm and Cold heat anomalies (\(10^{21} \hbox { J}\)) from piCtrl for 0–100, 0–700 m and the full ocean depth

To account for any spurious drift in the simulated globally averaged time series, a linear trend equal to that in piCtrl of the respective model has been removed from the Hist and RCP-scenarios, using identical periods from the piCtrl and the model runs in question. To remove short and long-term variations in correlations when we investigate relationships between the global temperatures and heat storage in the ocean, the model variables are typically bandpass filter using a third order Butterworth filter with cutoff at 5 and 31 years prior to the correlation analysis.

3 Results

3.1 Global temperature variations

It is impossible to accurately project future human-induced emissions of \(\hbox {CO}_2\) and other greenhouse gases and particles. As such, one cannot confine key constraints on global climate in the twentyfirst century and beyond. We therefore include four main emission scenarios for the twentyfirst century. These range from a high, business-as-usual, scenario (RCP8.5), two intermediate scenarios with maximum emissions occurring around 2080 (RCP6.0) and 2040 (RCP4.5), and a low-emission scenario (RCP2.6). RCP2.6 is constructed to be in line with the so-called two-degree target (Rijsberman and Swart 1990; Jaeger and Jaeger 2011). The observed \(\hbox {CO}_2\)-emissions indicate that we are currently tracking the RCP 8.5 scenario (Quéré et al. 2014). If continued, this will lead to almost a threefold increase of present day \(\hbox {CO}_2\)-emissions by the end of the century.

In addition to the long-term warming trend, the observed global surface temperature records show substantial short-term variations. These variations are most pronounced over land (Fig. 1a; Morice et al. 2012), but are also present in the combined land and ocean records (Fig. 1b; Jones et al. 2012b). For both types of time series, several 5–15 year periods without warming occur, despite a clear, long-term warming trend.

Fig. 1
figure 1

Annual global mean land surface temperature (a) and land and ocean surface temperature (b) relative to the period 1850–1900. Observed temperature for 1850–2014 in magenta from Morice et al. (2012) (land) and Jones et al. (2012b) (land and ocean). The remaining curves and shadings are obtained from the models listed in Table 1. Blue colours show Hist, red colours RCP8.5, and green colours RCP2.6. Solid/dashed blue, red and green lines show the ensemble mean value, the corresponding colour shadings show the full range of the models, and the grey shadings show the model spread expressed as 1 std. For completeness, the ensemble mean for of RCP6.0 and RCP4.5 for the last 30 years are shown as thin, black lines. The black, vertical line at 2006 denotes the transition from Hist to the RCP integrations

The future global surface temperature based on the RCP emission scenarios are also displayed in Fig. 1. Despite reduced greenhouse gas emissions in the second half of the twentyfirst century, both of the intermediate scenarios, RCP4.5 and RCP6.0, show surface warming with time. Only RCP2.6 gives a slight cooling in the second half of the twentyfirst century, with a temperature maximum occurring in the 2040’s. For the different RCP-scenarios, variations in the model specific climate sensitivities (e.g., Andrews et al. 2012) lead to long-term differences in the simulated surface warming, whereas internal variations lead to short-term divergences in surface temperature, including periods with non-increasing temperature.

3.2 Likelihood of hiatus periods

To examine the likelihood of hiatus periods in surface temperature, we compute probability density functions (PDFs) for non-positive (i.e., zero or negative) linear trends for periods of up to 30 years duration (see Fig. 2a). The probabilities are generated by calculating the running 2–30 year linear trends for the individual time series, divided by the total number of periods of the given length, similar to Easterling and Wehner (2009). Here we use the three datasets of the observed surface temperature (hereafter Obs) plus the simulated land and ocean surface temperature fields from the 17 CMIP5-models and all 6 scenarios. For Hist, the PDF’s are calculated for the period 1911–2005 in order to have the same number of years as the RCP-simulations spanning the period 2006–2100. piCtrl is computed for the period equivalent to 1911–2100.

Fig. 2
figure 2

Ensemble mean probability for having a negative trend of the length indicated on the abscissa (2–30 years) for a global mean surface temperature, b net TOA incoming radiation, and c 0–700 m ocean temperature. The probabilities are calculated according to the method of Easterling and Wehner (2009). The coloured lines show the model mean for piCtrl (black dashed), Hist (black), RCP2.6 (green), RCP4.5 (orange), RCP6.0 (blue) and RCP8.5 (red), with the shading showing the total spread among the models. The magenta line and shading in (a) shows the mean PDF of the observed (GISS, HadCRUT3 and HadCRUT4) global mean surface temperature. Note different scaling in (a) and (b). d The probability of having a 10 year negative temperature trend over land as a function of latitude for the first (thin lines) and second (thick lines) 45 year periods of the scenarios. Open circles show piCtrl, dashed lines Hist, and the fully drawn lines the RCP scenarios with colouring as above. Hist starts in 1911 to ensure consistency in length with the RCP runs

For piCtrl, negative and positive surface temperature trends are equally likely for all periods up to 30 years (see Fig. 2a). The magnitude of the trends are also equal whether positive or negative (not shown). This is the same as Easterling and Wehner (2009) found for the 10 year period. For Hist, the likelihood of negative temperature trends drops with increasing duration. In this case, negative surface temperature trends have probabilities ranging from 0.36 for a 10 year hiatus period to about 0.23 for a hiatus period of 25–30 years.

The ensemble mean probability of having a hiatus period in Hist is similar to that in Obs for short and long periods. However, for the intermediate length periods, around 10 years, the probabilities in Obs are below all of the simulated ones. This indicates either an overestimation of modelled global-scale decadal variations or an underestimation of the twentieth century global warming trend compared to the decadal scale variability.

Figure 3 shows the temporal evolution of the observed and simulated global temperature trends between 5 and 30 years used in Fig. 2a. By comparing Fig. 3a, b, we see that the majority of the models tend to simulate larger decadal variations than in Obs. This overestimation of the modeled decadal temperature variations relative to Obs could be caused by the prescribed concentrations of (short-lived) greenhouse gases, aerosol loadings or the solar cycle, impacts of interactions between these forcings and naturally occurring variability modes in the models (e.g., Otterå et al. 2010; Chylek et al. 2014), model deficiencies like the troposphere-stratosphere coupling (e.g., Manzini et al. 2012), or the ocean uptake and release of heat.

Fig. 3
figure 3

Global temperature trends (in \(^\circ \hbox {C}\) per decade) of 5–30 years duration, from a three observational time series, and from all models for the experiments, b Hist, c RCP2.6, d RCP4.5, e RCP6.0, and f RCP8.5. The colouring shows the actual occurrence of the trends, denoted by the first year of the period. For increased clarity, GISS (Hansen et al. 2010) is shifted −0.3 years, HadCRUT3 (Brohan et al. 2006) is plotted at the actual time, and HadCRUT4 (Morice et al. 2012) is shifted +0.3 years in (a) and the models have been evenly distributed between −0.3 years, the actual year, and +0.3 years. Note the different range in (a) compared to the other panels

For the intermediate and high RCP-scenarios (Figs. 2a, 3d–f), there is a gradual decrease in the occurrence of periods with negative trends, reflecting the increasing rate of global warming (Fig. 1). For RCP6.0 and RCP8.5, the trends with strongest warming are found towards the end of the integrations. Negative trends with a duration of 10 years are found in all scenarios, but tend to occur early in the time series when the global warming signal is still relatively weak. Even for the aggressive RCP8.5 scenario, there is a tiny—one percent—chance of having a negative temperature trend of 10 years, despite the rapid warming (Figs. 2a, 3f). The RCP2.6 scenario (Figs. 2a, 3c) diverges from the other scenarios in that all temperature trends show increased occurrence of negative trends in the second half of the twentyfirst century, in accordance with Fig. 1.

For all scenarios, the computed probability of negative trends is likely underestimated since these simulations use a repeated solar cycle (cycle number 23) and do not take into consideration variations in volcanic eruptions. These two factors contribute to externaly-driven, short-term fluctuations in global surface temperatures (Lean and Rind 2008, 2009; Feulner and Rahmstorf 2010; Foster and Rahmstorf 2011; Jones et al. 2012a; Rypdal 2012). As demonstrated by Maher et al. (2014), adding volcanic forcing to RCP4.5 dramatically increases the likelihood of a decadal hiatus period by the end of the 21st century. In their analysis, the RCP8.5 greenhouse gas forcing even masks the effect of the cooling from volcanoes in the likelihood of hiatus decades by the end of the century.

Figure 2b shows how the probability of having negative trends in the net radiative balance at the top of the atmosphere (net TOA incoming radiation) decreases with increasing emissions. This is in line with increased trapping of heat in an atmosphere with elevated levels of greenhouse gases (Hansen et al. 2010; Trenberth and Fasullo 2012). The only exception is scenario RCP2.6, which shows an increased likelihood of a trend in net loss of heat at TOA.

When trends longer than 31 years are removed from the time series, no significant short-term relationship between the net TOA incoming radiation and the global surface temperatures is found in piCtrl or in any of the scenario runs. This implies that a heat loss at the TOA cannot explain decadal-scale periods that show no warming. This is similar to the findings of Palmer and McNeall (2014). However, Hist shows a positive but low correlation between TOA and surface temperature, indicating that solar variations and volcanic eruptions (and possibly aerosol emissions) influence this model experiment.

Since the 1980’s there has likely been a net gain of heat at TOA (Hansen et al. 2005; Trenberth et al. 2009; Trenberth and Fasullo 2010; Loeb et al. 2012). This is also the case for all of the scenarios, except for RCP2.6 (Fig. 2b). With the assumption of no loss of heat at TOA or, alternatively, a net gain of heat at TOA, conservation of heat implies that periods with non-positive trends in surface temperature can only occur in tandem with redistribution of heat within the atmosphere–cryosphere–land–ocean system. Of these components, the ocean is by far the largest, dynamically active reservoir of heat (e.g., Charney et al. 1979; Levitus et al. 2005; Church et al. 2011; Katsman and van Oldenborgh 2011).

To check the above assumptions, Fig. 2c shows the probability of having a negative temperature trend for the upper 700 m of the global ocean. There is a small probability of a temperature hiatus in the upper ocean in Hist and, to some extent, in RCP2.6. For the other experiments the surface ocean accumulates heat at all time scales between 2 and 30 years.

3.3 Global decadal-scale temperature variations

As discussed in the previous section, the ocean is the prime candidate for temporary uptake and storage of heat (e.g., Levitus et al. 2005; Church et al. 2011; Katsman and van Oldenborgh 2011; Meehl et al. 2011, 2013; Trenberth and Fasullo 2013; Kosaka and Xie 2013; England et al. 2014). Table 2 shows that decadal trends in the upper ocean heat content (0–700m) coincide with decadal trends in the global surface temperature in all of the scenarios investigated, where the correlations between the two variables range between 0.31 and 0.56 for the different scenarios.

Sun and Trenberth (1998) found that the upper ocean heat content imbalance on inter-annual time scales probably occurs in relation to El Niño-La Niña or El Niño-Southern Oscillation (ENSO) variations. The ensemble mean Niño 3.4 index shows a correlation of about 0.60 with the global surface temperature across all scenarios (Table 2). A similarly strong relationship holds for the Niño 3.4 index and the upper ocean heat content. One can therefore expect ENSO, or variability modes with related surface temperature characteristics, to be a central component in decadal-scale surface temperature variations.

The surface ocean alone cannot sequester heat for a very long time without releasing it back to the atmosphere. We therefore investigate the role of the sub-surface ocean. By comparing the global upper (<700 m) versus deep (>700 m) ocean heat content trends, a robust, inverse relationship is found across all simulations, with a correlation ranging between \(-0.53\) and \(-0.66\) for piCtrl and the RCP scenarios. For Hist, a slightly weaker relationship is found (\(r=-0.44\)). This indicates that external forcings, like aerosols or volcanic eruptions, may interact with the phasing or strength of the naturally occurring variations in the models (e.g., Pacific Decadal Oscillation (PDO) or the Atlantic Multidecadal Oscillation/Atlantic Meridional Overturning Circulation; Otterå et al. 2010; Chylek et al. 2014).

In the following sections, the location, magnitude and duration of upper and deep ocean heat anomalies are examined. The analysis is focussed on vertically integrated heat contents rather than the sea surface temperature since the former is less influenced by short-term fluctuations, yielding more robust patterns.

3.4 Latitudinal variations in the ocean heat content

In this section, heat anomalies per unit area are first calculated based on the temperature anomaly from each model, a constant heat capacity (\(4168 \hbox { J kg}^{-1} \hbox { K}^{-1}\)) and water density (\(1027 \hbox { kg m}^{-3}\)), and the bathymetry from the respective models. The anomalies are then mapped onto a standard 1-by-1 degree latitude−longitude grid and multiplied by the area of the remapped grid cell yielding the total heat content of each grid cell.

The latitude−time distribution of the inter-model spread, represented by 1 standard deviation, of the 5–31 year bandpass filtered heat anomalies above and below 700 m is provided in Fig. 4. The global ocean, upper heat anomalies are most pronounced in a band just north of equator with center of action between \(10^\circ \hbox {N}\) and \(20^\circ \hbox {N}\), followed by latitudinal variations with weaker magnitude just south of equator and at \(\pm 40^\circ\) latitude (panel a). At depth (panel b), the Southern Ocean shows by far the largest anomalies, followed by enhanced values at subtropical latitudes in the northern hemisphere.

Fig. 4
figure 4

Latitude−time (Hovmöller) plot of 1 std of the model response in zonally integrated heat anomalies (\(10^{20} \hbox { J}\)). Left panels heat anomalies integrated above 700 m and right for below 700 m. Shown are the global ocean (a, b), Atlantic Ocean between 20 and \(70^\circ \hbox {N}\) (c, d), Pacific Ocean between 10 and \(60^\circ \hbox {N}\) (e, f), and Southern Ocean between 70 and \(20^\circ \hbox {S}\) (g, h) from piCtrl. All models are re-gridded to 1-by-1 degree resolution before the sum for each latitude was made for each model. The time series have been bandpass filtered (5–31 years), then averaged over all models. To guide visual inspection, all of the ocean basins are plotted with a latitudinal range of 50°

The corresponding basin specific ocean heat anomalies (Fig. 4c–h) show that in the upper portion of the water column (left panels), largest variations are found in the Pacific Ocean basin between 10 and \(20^\circ \hbox {N}\), followed by the Southern Ocean at about \(40^\circ \hbox {S}\). At depth (right panels), variations are clearly largest in the Southern Ocean with the largest amplitudes near \(40^\circ \hbox {S}\).

To infer the magnitude and duration of the upper and deep ocean heat anomalies, integrated over the ocean basins, 10 year non-overlapping periods from all models are binned by anchoring the time of maximum heat anomaly at year zero. For the global analysis, this is done for the entire ocean domain, whereas the basin analyses are restricted to the actual basin domains. The resulting negative heat anomalies are presented in Fig. 5 and Table 3 (the positive anomalies have similar shapes and absolute values, and are not shown). The integrated negative heat anomalies in the upper 700 m of the Southern and Atlantic Oceans are 75 and 63 % of that in the upper Pacific Ocean, respectively. At depths greater than 700 m, the anomaly in the Southern Ocean is about twice as large as those in the Atlantic and Pacific Oceans.

To put these anomalies in an observation-based perspective, the linear rate of change of the 0–700 and 700–2000 m ocean heat contents since the mid 1950’s are \(3.0\times 10^{21}\) and \(1.3\times 10^{21} \hbox { J yr}^{-1}\), respectively (Levitus et al. 2012). The global full depth and the upper Pacific Ocean values from Table 3 are large, roughly 2.9 and 2.3 times greater than the long-term, 0–2000 m yearly trend in ocean warming reported by Levitus et al. (2012). The ensemble mean point-by-point anomalies are, however, rather short lived; accumulation and release of heat occurs on a pentadal time scale (mainly from −2 to +2 years in Fig. 5).

Fig. 5
figure 5

Heat anomaly composites (\(10^{21} \hbox { J}\)) of all non-overlapping periods where the minimum heat anomaly for the individual model, averaged over 10 years, is less than \(-1\) std. The upper (<700 m) and deep (>700 m) heat anomalies are identified at each latitude band and then summed over all latitudes. The heat anomalies are binned together by anchoring the time of minimum heat anomaly at year zero, and are shown from 3 years before to 3 years after the minimum anomaly. All of the ocean basins span a latitude band of 50°; the Atlantic Ocean covers the area between 20 and \(70^\circ \hbox {N}\), the Pacific Ocean between 10 and \(60^\circ \hbox {N}\), and the Southern Ocean between 80 and \(30^\circ \hbox {S}\). The time series for each latitude and model are bandpass filtered (5–31 years) prior to making the composites

3.5 Spatial relationship between upper and deep ocean heat content

To further examine the relationship of the heat contents between the upper and the deep ocean and regional distributions thereof, local (point-by-point) and global correlation analyses are presented below. The local analysis reveals the vertical transfer of heat without, or with only weak, horizontal transport processes, for instance related to deep winter-time mixing at mid to high latitudes or the net effect of seasonal ventilation of the upper water masses. The global analysis will, on the other hand, incorporate all variations in the ocean heat content, irrespective of the physical process involved.

3.5.1 Local and global correlations

Figure 6a, c show the ensemble mean point-by-point correlations of the integrated heat content above and below 300 and 700 m, respectively. The upper and deep waters in the North Atlantic and the North Pacific subpolar gyres, and in the Southern Ocean/southern subpolar gyres, show prominent co-variability (\(r\,\gtrsim\,0.6\)). This pattern is in line with the large variability signals at \(\pm 40^\circ\) in Fig. 4. Negative correlations are mainly found for the 700 m threshold depth, where weaker (\(r\,\lesssim\, 0.3\)), negative correlations are mainly confined to the region west of Australia, in the tropical East Pacific Ocean and partially in the subtropical gyres in the Atlantic Ocean.

Fig. 6
figure 6

Local (point) correlation between the 5 and 31 year bandpass filtered mean heat content above and below 300 m (a) and above and below 700 m (c) (in colour). The model spread of the correlations is represented by 1 std (black contours), computed for each model and then averaged across all model experiments. b, d Correspond to a, c, but show the correlations between the bandpass filtered (5–31 years) time series of the global mean upper ocean heat content above and the grid-point ocean heat content below 300 and 700 m, respectively

The overwhelmingly positive correlations in Fig. 6a, c are in stark contrast to the general negative correlation between the upper and deep ocean heat content shown in Table 2. Local correlations, and consequently pure vertical water mass exchange processes, are therefore not suited to account for the decadal-scale variations in the ocean heat content.

In contrast to the local correlations in Fig. 6a, c, the ensemble global mean upper ocean heat content in the upper 300 and 700 m show large regions of anti-correlation with the local, sub-surface ocean heat content at depths greater than 300 and 700 m, respectively (Fig. 6b, d). Except for the negative correlation in the Amundsen and Ross Seas of the Southern Ocean, the single most prominent, large-scale region with negative correlation covers most of the eastern Pacific Ocean, particularly at low and mid latitudes. This region is the focus for the remainder of Sect. 3. It also follows that the global correlation patterns (Fig. 6a, c) are similar, indicating that the obtained signal is a robust feature in the uppermost 700 m of the water column.

3.6 Decadal-scale temperature composites

To investigate sub-decadal preconditioning phases, if any, and fully developed decadal-scale warm and cold periods, two sets of composites are constructed and analysed.

The first set, hereafter named the Warm and the Cold phases, or composites, denote non-overlapping 10 year periods with global mean surface temperature above or below a given threshold value \(\overline{T}\) relative to a smoothed version of the unfiltered time series. Since the global mean surface temperature time series based on Hist for 1910–2005 and the four RCPs for 2006–2100 exhibit quite different long-term evolutions (see Fig. 1), the regression method Locally Weighted Scatterplot Smoothing (LOWESS; Cleveland 1979, 1981) is used for the temporal filtering. The fraction of time series points used to compute each fitted value is set to 0.2, corresponding to a smoother window of about 38 points (years) in our case. LOWESS is viewed as a robust, local smoother (e.g., Foster and Brown 2015).

Figure 7a shows unfiltered and LOWESS smoothed time series of the global surface temperature anomalies based on a merge of Hist and RCP4.5 for the 15 models with ocean data (Table 1). By subtracting the filtered from the unfiltered time series, residual time series are obtained. The residual time series might have a very weak trend, so these are furthermore linearly detrended. The resulting time series are displayed in Fig. 7b.

Fig. 7
figure 7

a Unfiltered (black) and LOWESS filtered (gray) time series of the global mean temperature anomalies (K) from model 1–15 (see Table 1) based on a combination of Hist and RCP4.5. In (b), the linearly detrended residual time series based on (a) is shown in black. Warm, Cold, Pre-Warm and Pre-Cold are shown in deep red, deep blue, magenta and cyan colours, respectively. See text for details. The ordinate shows the temperature anomaly for every second model

Based on the individual residual time series \(T_n\) (n denotes the year) in Fig. 7b, all 10 year periods with an average temperature above \(\overline{T}\) for Warm (below \(-\overline{T}\) for Cold), with \(0.05\,\mathrm{K}\le |\overline{T}|\le 0.09 \hbox { K}\) in intervals of 0.01 K, are extracted. Warm (Cold) phases are then chosen from all of the warmest (coldest), non-overlapping periods.

The second set, hereafter named the Pre-Warm and the Pre-Cold composites, is constructed to identify characteristic features when Warm (Cold) follow after relatively cold (warm) years. Pre-Warm is defined as the 3–7 year time period prior to Warm, satisfying at least \((T_{m-1}+T_{m-2})\,<\,\overline{T}\) and \(T_{m-3}\,<\,\overline{T}/2\) and, if present, also \(T_{m-i}\,<\,\overline{T}/3\) (\(i=4, \ldots ,7\)), where m denotes the first year of Warm. Corresponding conditions define Pre-Cold.

Finally, if the first or the first and second year of Warm are relatively cold, these years (\(T_m\), or \(T_{m}\) and \(T_{m+1}\)) are also included in Pre-Warm. Similarly for Pre-Cold. In any case the length of the precondition phases does not exceed 7 years.

Based on the above definitions, Pre-Warm can be expected to share similarities with Cold, and Pre-Cold with Warm. In situations when Warm is followed by Cold or vice versa, the preconditioning phases will partly overlap with the previous Cold and Warm periods.

The resulting Warm, Cold, Pre-Warm and Pre-Cold periods are all marked in Fig. 7a, b. The selection criteria are clearly arbitrary, and other definitions may work equally well. We find that a threshold of \(|\overline{T}|=0.06 \hbox { K}\) yields a reasonable number of non-overlapping Warm and Cold decades. The remaining analyses, including the time series shown in Fig. 7, are based on this threshold value. The other values of \(\overline{T}\) yield similar results.

The two sets of composites are computed for the depth intervals 0–100, 0–300, 0–500, 0–700 m and over the full water depth (Sect. 4.1). The total number of analysed model years is 9810, and the individual Warm (Pre-Warm) and Cold (Pre-Cold) cases occur 17 % (6 %) and 14 % (6 %) of the total integration time, respectively (see Table 4).

Without filtering, long-term temperature trends and variations, whether governed by internal variability or caused by model drift due to numerical formulations, model resolution or choices of physical parameterisations, will necessarily influence the decadal-scale Warm/Cold composites. Therefore, long-term local (grid-point by grid-point) trends are removed from all model runs and the considered depth intervals, prior to forming these composites. For piCtrl a linear trend has been removed. For the other scenarios, a LOWESS filter has been applied and removed as explained above. A comparison between the original and detrended Warm and Cold composites from piCtrl shows negligible differences in the upper 300–500 m, ventilated portion of the water column. Detrending is, however, required for depths greater than about 500 m, since most models show a slow warming or cooling trend in the abyss.

3.6.1 Ocean state during Warm and Cold phases

The ensemble mean SST anomalies for Warm and Cold are depicted in Fig. 8a, b, whereas Fig. 8c–f, shows the corresponding 0–100 and 0–300 m ocean heat content anomalies, respectively. The heat anomalies are first calculated based on the temperature anomaly from each model mapped onto a standard 1-by-1 degree latitude−longitude grid. For simplicity, a constant heat capacity of \(4168 \hbox { J kg}^{-1} \hbox { K}^{-1}\), a water density of \(1027 \hbox { kg m}^{-3}\) and a common ETOPO5 bathymetry are used. These anomalies are then averaged over all models. By not using the model’s own bathymetry, the values given in the figures will only give a rough estimate of the real heat content in the models. However, when calculating the heat content for individual basins or for the total global ocean, the model’s own bathymetry is used.

Fig. 8
figure 8

Ensemble mean SST anomaly (K, upper panels), and 0–100 and 0–300 m vertically integrated ocean heat content anomaly (×\(10^8 \hbox { J m}^{-2}\), mid and bottom panels) for Warm (left column) and Cold (right column) composites based on all model runs. Regions with large, small and no dots denote regions where <70, 70–90 and >90 % of the models agree with the sign of the ensemble mean

The spatial signature of the SST and the ocean heat content anomalies for the Warm and Cold phases are inversely related, but with a tendency of larger amplitudes for the Warm phase. Large regions in the Atlantic, the Central and East Pacific and the Indian Oceans have anomalies of the same sign, with clearly the largest geographical extent and largest amplitudes at low to mid latitudes in the Central and East Pacific.

Despite differences in the modelled SST anomalies and, notably, in the sub-surface heat anomalies (next section), the sign of the majority (90 %) of the models agree throughout the Central and East Pacific Ocean, at the poleward rims of the Pacific basin, in the western Indian Ocean and, in general, at low latitudes in the Atlantic Ocean. These patterns can therefore be viewed as robust features across the model ensemble. The obtained patterns resemble, as expected, those of Fig. 6b, d. The spatial patterns of Pre-Warm and Pre-Cold resemble those of Cold and Warm, respectively.

For comparison, the globally integrated ocean heat content associated with the two sets of composites are provided in Table 5. The magnitudes of Warm and Cold show small changes below 300 m, implying that most of the decadal-scale heat anomalies, when globally averaged, are confined to the upper 300 m of the water column.

The magnitude of the upper 700 m ocean heat anomalies in Warm and Cold are about \(7.5\times 10^{21} \hbox { J}\), a value 2.5 times larger than the observed annual increase in the 0–700 m global ocean heat content of \(3.0\times 10^{21} \hbox { J}\) between 1955 and 2010 (Levitus et al. 2012). Between 1961 and 2008, the observed increase in ocean heat content has roughly been an order of magnitude larger than the increased heat content of the atmosphere (Church et al. 2011). This indicates that the ensemble mean ocean heat content anomalies in Warm and Cold, when compensating for the steady increase in the observed ocean heat content since the 1960s, have the duration (by construction) and magnitude (Table 5; 0–700 m) roughly in line with decadal-scale periods with non-increasing global surface temperature.

Furthermore, Table 5 shows that decadal-scale periods with anomalously low surface temperatures (Cold) are associated with an anomalously cold upper and an anomalously warm deep ocean, and vice versa for Warm. Integrated over the full ocean depth, the heat anomalies have, however, the sign of the 0–100 m anomaly.

The magnitude of the full depth heat anomalies for Pre-Warm and Pre-Cold are a factor 2–4 smaller than those of Warm/Cold. For both Pre-Warm and Pre-Cold, the largest heat anomalies are located between 100 and 300 m depth. The ocean heat anomalies in the depth range 100–300 m have the sign of the following Warm and Cold phases, respectively. The sign of the 0–100 m anomalies are, however, opposite to the following Warm and Cold phases. So despite the spatial similarities between Warm and Cold, and Pre-Cold and Pre-Warm, the vertical distribution differs (Sect. 3.6.2).

Due to the drift in the deep ocean in most of the models (e.g., Sen Gupta et al. 2013), the column integrated heat anomalies have been linearly detrended before making the total global and regional heat content anomalies for piCtrl in Table 5. Since the drift affects the integrated heat content differently if the ocean is split into two depth intervals, removing the trend will result in a slight inconsistency when adding upper and deep heat anomalies compared to full depth calculations (Tables 56). Regional differences and suitable choices for the depth horizon separating the upper and deep ocean are elaborated on in Sect. 4.

3.6.2 Associated thermocline variations

Composites of vertical 0–300 m distribution of temperature anomalies along the equatorial Pacific for the four phases are presented in Fig. 9. Also shown in these panels is the mean depth of the thermocline, sloping upwards from about 150 m in the west to less than 50 m in the east.

Fig. 9
figure 9

Ensemble mean temperature for a Warm, b Cold, c Pre-Warm and d Pre-Cold composites from piCtrl averaged between \(5 \, ^\circ \hbox {S}\) and \(5 \, ^\circ \hbox {N}\) in the Pacific Ocean. The ensemble mean thermocline, here indicated by the model mean \(20^\circ \hbox {C}\) isotherm, is shown in gray, and the 1 std model spread is given by black contours. Local (column-by-column) linear trends are removed from each run prior to generating the ensemble mean. Note different colour scales in the upper and lower figures

East of about \(140^\circ \hbox {W}\), all of the phases show anomalies of the same sign extending from surface to a depth of about 300 m. In the East and Central Pacific, a pronounced sub-surface wedge of opposite sign is found. The latter is positioned around the depth of the mean thermocline east of \(150^\circ \hbox {W}\).

The temperature signals in the Pre-Warm and Pre-Cold phases (Fig. 9c, d) are stronger and, consistent with Table 5, penetrate deeper than in Warm and Cold (Fig. 9a, b). The actual magnitude could partly be explained by the 3–7 year duration of the pre-cases, compared to the 10 year duration of the Warm and Cold phases. The signal in the pre-phases is, however, prominent. Furthermore, Pre-Warm and Cold, and Pre-Cold and Warm, show similar features as expected based on the construction of the pre-phases (Sect. 3). Thus, a rapid transition from one ocean state to another is found in the models, partly linked to the inter-annual ENSO variability.

The surface and sub-surface signals seen in the four phases, and then particularly in the Pre-Warm and Pre-Cold phases, are related to the ventilation of subsurface waters along the equator (e.g., England et al. 2014). In Pre-Cold (and Warm), the thermocline moves upward in the west and downward in the east, leading to a flattening of the main thermocline, with a resulting surface cooling in the west and warming in the east. In the Pre-Warm (and Cold) phase, the thermocline steepens, yielding an inverse surface temperature response.

Although the model spread in Fig. 9 is large and comparable in magnitude to the ensemble mean signal, a comparison between variations in the ensemble-mean and observation-based ocean heat content and related fields (Sect. 4.3) gives credibility to the presented ensemble-mean results.

3.6.3 Surface wind stress anomalies

A prime candidate for changes in the zonal slope of the mean thermocline in the equatorial region is the changes in the surface wind stress (England et al. 2014; Watanabe et al. 2014). Figure 10 displays changes in the mean sea level pressure (slp) and surface winds in the two sets of composites. It follows that Warm and Cold, and Pre-Warm and Pre-Cold, are inversely related. Furthermore, the Warm/Pre-Cold and the Cold/Pre-Warm phases are very similar except for the Atlantic and the Arctic regions. The largest amplitudes are, however, clearly found during the Pre-Warm and Pre-Cold phases.

Fig. 10
figure 10

Ensemble mean sea level pressure (hPa, contoured) and surface wind (vectors) anomalies from piCtrl for a Warm, b Cold, c Pre-Warm and d Pre-Cold. The contouring interval is 0.1 hPa with warm colours showing positive slp anomalies. Reference vectors are shown over South America

In the Pacific Ocean, Warm and Pre-Cold are characterised by weakened Walker and Hadley cells, similar to the findings of Kosaka and Xie (2013) and England et al. (2014). The tendency for reduced trade winds and a reduced zonal pressure gradient in the tropics imply reduced Ekman-induced upwelling of colder sub-surface waters in the eastern part of the basin and reduced transport of surface water towards the west. The net result is surface cooling in the west and warming in the east. The opposite response, with accelerated Walker and Hadley cells, seen in Cold and Pre-Warm.

The zonal component of the surface wind anomalies in the four phases in Fig. 10 change sign at, or west of, the date line at low latitudes. This shift is associated with the Interdecadal Pacific Oscillation (IPO, England et al. 2014), with positive IPO-anomalies related to weakened trade winds and vice versa for negative IPO-anomalies (Meehl et al. 2011, 2013; Maher et al. 2014; Trenberth and Fasullo 2013). Positive (negative) IPO-anomalies therefore coincide with our Warm/Pre-Cold (Cold/Pre-Warm) phases.

The IPO-mode/variations in the trade winds, related to low-frequency ENSO variations, is in line with the high correlation between the decadal trends of the ensemble mean Niño 3.4 index and both the global surface temperature and the upper ocean heat content (Sect. 3.1). High-resolution, observation-based evidence from analysis of corals in the West Pacific, extending back to the late nineteenth century, indicates that decadal-scale variations in the Pacific trade winds influence decadal-scale variations in the global surface temperature (Thompson et al. 2014), illustrating the key role of ENSO-related variations on global temperature.

The most pronounced slp anomalies in Fig. 10 are related to the Aleutian Low in the central North Pacific and the Amundsen-Bellingshausen Seas low off the Antarctic continent, both with total changes in amplitude of 0.6 hPa (1.2 hPa) between Warm and Cold (Pre-Warm and Pre-Cold). The intensified Aleutian low in the Warm/Pre-Cold phases imply cooling of the upper ocean in the western North Pacific caused by (cold) northerly winds and warming of the upper ocean towards the North American coast (left panels in Fig. 8), and an anomalous equatorward Ekman transport south of the Aleutian low (Ceballos et al. 2009). The weakened Aleutian low in the Cold/Pre-Warm phases lead to a reversed forcing of the upper ocean (right panels in Fig. 8).

3.7 Simulated SST pattern linked to changes in the thermocline temperature

To further explore the large-scale surface imprint of changes in the thermocline in the Pacific Ocean, the region with the largest decadal-scale variations in the thermocline temperature is roughly bounded by 125–\(95^\circ \hbox {W},\, 5^\circ \hbox {S}\)\(5^\circ \hbox {N}\) and 50–130 m depth (Fig. 9). Temperature variations in this region, which can be viewed as a sub-surface counterpart to the Cold Tongue Index region (\(180^\circ\)\(90^\circ \hbox {W},\, 6^\circ \hbox {S}\)\(6^\circ \hbox {N}\); Mantua et al. 1997), have a surface temperature imprint as shown in Fig. 11. The pattern resembles the Pacific Decadal Oscillation (e.g., Mantua et al. 1997; Mantua and Hare 2002) and variations in slp associated with the Aleutian low (Liu 2012, see also Sect. 3.6.3). The obtained pattern also resembles the ensemble mean SST and for 0–100 m heat content anomalies (Fig. 8a, c). The similarity is not limited to the Pacific Ocean, but holds for most of the global ocean (with a notable exception for the northern North Atlantic). Variations in the thermocline depth and thermocline temperature in the eastern Pacific, governed by changes in the surface wind stress, are therefore a proxy for near global, decadal-scale surface temperature anomalies.

Fig. 11
figure 11

Ensemble mean correlation between time series of the East Pacific thermocline temperature and local SST across all piCtrl runs. The model spread of the correlations is represented by 1 std (black contours). The thermocline temperature is defined as the average temperature in a box bounded by the depth interval 50–130 m and confined by 125–\(95^\circ \hbox {W}\) and \(5^\circ \hbox {S}\)\(5^\circ \hbox {N}\)

4 Discussion

4.1 Penetration depth of the heat anomalies

The ensemble mean heat content anomalies for Warm and Cold for 0–500 and 0–700 m are very similar to the 0–300 m anomalies shown in Fig. 8e. This, and the quantification given in Table 5, indicates that a global-scale threshold depth in the range 300–700 m is appropriate for a simplified, two-layer analysis as presented here. An implication of this is that the presented correlation analyses with a threshold depth of 700 m (Tables 23) is also representative for threshold depths in the range 300–700 m.

The main difference between the heat anomalies in the upper 100 m and at greater depths is that the predominantly uniform positive (negative) surface signal in Warm (Cold), see Fig. 8c (Fig. 8d), is split into an east-west dipole at depth (Fig. 8e, f). This dipole-structure is characterised by negative (positive) heat anomalies in Warm (Cold) in the tropical West Pacific.

Outside of the Pacific sector, the heat anomalies are generally weak (Fig. 8), with the exception of the North Atlantic sub-polar region and the western half of the Indian Ocean, both with anomalies of similar sign as the East Pacific anomaly.

Most of the remaining discussion is confined to piCtrl since this is the only model run where climate variations are governed by internal, stochastic processes only, though the results are similar for all of the scenarios when the LOWESS fit has been removed.

4.2 Quantification of the ocean basin heat anomalies

The magnitudes of the decadal-scale heat anomalies in Warm and Cold for the main ocean basins are given in Table 6. In the upper 0–100 m, about 47/17/21 % (48/25/14 %) of the Warm (Cold) heat anomalies are found in the Pacific/Atlantic/Southern Ocean. For 0–700 m, the corresponding distributions are 19/24/33 % (34/25/16 %). The contribution to the global, decadal-scale heat anomalies from the Pacific Ocean is thus reduced with depth, going from being the dominant contributor in the upper 100–300 m of the water column to become on par with the other ocean basins integrated over greater depths.

4.3 Observation-based evidence

Trenberth and Fasullo (2013) examined variations in the ocean heat content of the recent global hiatus (1999–2012) minus global warming (1976–1998) period, using the ORAS-4 ocean reanalysis from ECMWF (Balmaseda et al. 2013). When compensating for the observed, century-scale global surface warming, by removing the long-term surface temperature trend, the hiatus period can be viewed as a decadal-scale cooling episode. The Trenberth and Fasullo analysis, i.e., hiatus minus warming, therefore shares similarities to our Cold minus Pre-Cold phases. Based on this, the difference between Cold and Pre-Cold is shown in Fig. 12 for the depth intervals 0–100, 0–700 m and the full ocean depth.

Fig. 12
figure 12

The simulated, vertically integrated ocean heat content difference (\(\times 10^8 \hbox { J m}^{-2}\)) between Cold and Pre-Cold based on all model runs. From top: 0–100, 0–700 m and over the full ocean depth. Regions with large, small and no dots denote regions where <70, 70–90 and >90 % of the models agree with the sign of the ensemble mean. Note that the contouring range is half of those in Fig. 12 in Trenberth and Fasullo (2013)

Visual inspection of the ORAS-4 ocean reanalysis (Fig. 12 in Trenberth and Fasullo 2013) and our Fig. 12 uncovers both similarities and differences. In the upper 0–100 m, the main differences are that the ensemble mean heat content shows more extensive negative anomalies in the Central and East Pacific Ocean, and weaker positive anomalies west of the date line, and that the Atlantic Ocean has an overall negative, rather than a positive, heat content anomaly (see below).

The upper 700 m and the full depth analyses in the Pacific sector show, however, a perhaps surprisingly close overall agreement between the ensemble mean fields both with respect to the location of the most prominent positive and negative anomalies, and their spatial extensions. Both analyses show a tongue of negative heat anomalies extending westward along the equator to about the date line, and bands of positive heat anomalies extending northeastwards and southeastwards from the equatorial West Pacific poleward of the Central/East Pacific cold wedge. In the South Pacific, the positive anomaly approaches the southern tip of South America, as is the case in the reanalysis. It is also worth noting that most of the CMIP-models agree on the sign where ORAS-4 shows the largest amplitudes.

The ensemble mean amplitudes are typically a factor of 2–3 smaller than in ORAS-4. Obviously, shorter-term variability modes like the actual El Niño-La Niña states during the 1976–1998 and 1999–2012 periods are not present in the ensemble mean, so a close match between the two analyses cannot be expected.

Furthermore, the slp-field from the Cold and Pre-Cold composites displayed in Fig. 10b, d, and the related Cold–Pre-Cold slp (not shown), resemble the reanalysed slp difference between the recent global hiatus (1999–2012) minus global warming (1976–1998) period (Fig. 11b in Trenberth and Fasullo 2013). This is particularly the case for the Pacific Ocean, the Indian Ocean and the Indo-Pacific portion of the Southern Ocean.

The slp-fields in Fig. 10 (this paper) also share similarities with the 1992–2011 trend analysis presented by England et al. (2014). Together with the global imprint of heat anomalies in the tropical East Pacific that are in line with the analysis of Kosaka and Xie (2013), we conclude that the simulated ensemble mean decadal-scale temperature anomalies in CMIP5 shares resemblance to several large-scale, observation-based fields.

The ensemble mean ocean heat anomalies in the Atlantic Ocean are, however, largely opposite to those in ORAS-4. It is beyond the scope of the present analysis to elaborate on this difference. It is, however, well-known that decadal and multi-decadal variations in the Pacific and Atlantic Oceans are driven by different physical mechanisms and that they exhibit variations on different time scales. It is therefore possible that the actual ocean state during the analysis period of Trenberth and Fasullo (2013), including the Pinatubo eruption in 1991 and the unprecedented El Niño-event in 1997–1998, in combination with global warming, contributed to the actual ocean state during the last decades.

4.4 Regional temperature variations and trends

The presented analysis has addressed basin to global-scale temperature variations, and particularly the likelihood of non-increasing decadal-scale surface temperatures under global warming. On a more local scale, and presumably closer to an individual’s perception of variations in temperature, decadal-scale temperature variations will generally differ from those being reported on a global scale, and certainly from an ensemble of model simulations.

The probability of 10 year periods without warming on land, as a function of latitude for the analysed model experiments, is displayed in Fig. 2d. It follows that the smallest probability of non-warming periods occur between 20° and \(40^\circ \hbox {N}\) and that decadal-scale variations increase towards both poles. It also follows from Fig. 2d that one can expect reduced probability of non-warming periods in the second half of the century for scenarios with strong greenhouse-gas emissions (RCP6.0 and RCP8.5), with an opposite response for scenarios with greatly reduced emissions (RCP2.6 and RCP4.5).

The spatial representation of the likelihood of a negative temperature trend with 10 years duration in the four RCP experiments are displayed in Fig. 13. For piCtrl (not shown), the likelihood is around 0.5 everywhere, as expected. For the RCP-scenarios as well as for Hist (not shown), the probability for non-warming decades is lowest in the region extending from the equatorial Atlantic towards, and including, Indonesia.

Fig. 13
figure 13

Ensemble mean probability for local (grid point) 10 year negative surface temperature trends for the period 2006–2100 for a RCP2.6, b RCP4.5, c RCP6.0 and d RCP8.5

The regions with highest probability for non-warming decades are the northern North Atlantic and the Southern Ocean, with a probability of around 0.40–0.50 depending on the actual scenario. In RCP4.5 (RCP8.5), the lowest probability of a local, non-warming decades is in the range 0.20–0.25 (0.05–0.10). The regions of lowest probability are consistent with the finding that the largest signal-to-noise ratio in surface temperature (e.g., Räisänen 2001), and the first detection of the global-warming signal (e.g., Kattsov and Sporyshev 2006), occurs at low latitudes, with the exception of the equatorial East Pacific where the pronounced ENSO-signal dominates the signal.

5 Summary and concluding remarks

We have analysed nearly 10,000 model years of data with particular focus on non-increasing decadal-scale trends in surface temperature. The data comes from a control integration (piCtrl), an integration with prescribed composition of atmospheric greenhouse gases and particles for the period 1850–2005 (Hist) and four RCP scenarios for the period 2006–2100, carried out with 17 state-of-the-art climate models.

For piCtrl, negative and positive surface temperature trends are equally likely, as expected. The probability of non-warming periods between 2 and 30 years in Hist is within the uncertainty of the observed surface temperature analyses with the exception for periods near 10 years. For this time window, the observed probability is slightly less than that found in the models. When the models are run with anthropogenic forcings, there is a low (a few percent) likelihood of non-warming periods lasting more than about 10, 15 and 30 years for RCP8.5, RCP6.0 and RCP2.6, respectively. Consequently, non-warming surface temperature is found across all model experiments up to about 10 years, even for the aggressive greenhouse-gas emissions in the business-as-usual scenario RCP8.5.

Variations in decadal-scale surface temperatures, and in upper and deep ocean temperatures, can be analysed by removing longer-term trends. It follows that decadal-scale, non-warming surface temperature periods under long-term global warming cannot be explained by enhanced loss of heat at TOA, but by anomalous uptake and storage of heat in the ocean. Comparison of global mean heat anomalies above and below 700 m show an inverse relationship; when the upper ocean warms the deep ocean cools, and vice versa. Furthermore, a non-increasing surface temperature coincides with anomalously low heat content in the upper ocean and warming in the abyss.

In apparent contradiction to the above, local point-by-point correlations of the heat content above and below thresholds between 300 and 700 m depth show strong co-variability throughout most of the ocean, particularly at mid to high latitudes, and in the Southern Ocean. The duration of the covariabilities are, however, rather short-lived, occurring on a pentadal time-scale. Local, vertical exchange of heat at mid to high latitudes can therefore contribute to, but cannot explain, decadal-scale periods with non-increasing surface temperature trends under global warming.

If, instead, the global mean heat content in the upper ocean is correlated with the local point-by-point heat anomalies in the sub-surface ocean, a prominent, large-scale region with negative correlations occur throughout most of the eastern Pacific Ocean, particularly at low and mid latitudes. To examine this pattern, 10 year periods with anomalously warm and cold global mean surface temperatures (Warm and Cold composites), as well as 3–7 year periods prior to Warm with relatively cold anomalies and vice versa for Cold (Pre-Warm and Pre-Cold composites), are examined.

When the surface temperature is anomalously high (Warm composite), the upper ocean is warmer than normal whereas the deep ocean has a negative heat anomaly. Summed up over the entire ocean domain, the ocean heat anomaly is, however, positive. The opposite situation is the case when the surface temperature is anomalously low.

The decadal-scale anomalies associated with Warm and Cold, integrated over the upper 700 m, have a magnitude of about \(7.5\times 10^{21} \hbox { J}\). Based on the observed 1960–2010 global warming, this amount of heat has the magnitude needed to maintain a 10 year period without increasing surface temperature under global warming.

The decadal-scale anomalies for Warm, Cold, Pre-Warm and Pre-Cold are related to variations in the slope and the temperature of the main equatorial thermocline in the Pacific Ocean. Anomalously high surface temperatures coincide with weakened trade winds, a flattening of the zonal slope of the main thermocline and reduced upwelling of cold sub-surface waters in the east, resulting in a positive (negative) upper ocean temperature anomaly in the east (west).

For the upper 0–100 m in Warm and Cold, about 47 % (23 %) of the total ocean heat anomalies are located in the Pacific (Atlantic) Ocean. For the 0–700 m heat anomalies, the relative importance of the Pacific Ocean is reduced with increasing contributions from the Atlantic and Southern Oceans. The contribution from the Indian Ocean is generally small.

The above findings are based on an idealised two-layer (upper versus deep ocean) analysis. More elaborate choices for the framework of the analysis can be made, for example based on density criteria or constrained by the actual ventilated water masses. For Warm and Cold, threshold depths between 300 and 700 m are suitable for a rough, two-layer analysis. For Pre-Warm and Pre-Cold, the most profound anomalies are present closer to the ocean surface, typically in the upper 100–300 m.

One key difference between the analysed simulations and reality is that nature shows one out of an infinite number of possible combinations of naturally occurring variability modes, whereas the simulations used in this analysis do not have the information about the actual state of, say, ENSO, the Southern Annular Mode or the Atlantic Meridional Overturning Circulation. One-to-one comparisons between observed and ensemble mean anomalies can, therefore, not be easily made.

Nevertheless, our analysis, averaged over 17 models, show some resemblance to observed, decadal-scale periods with weak or without increasing surface temperature under global warming. This is partly because of pattern similarities between the ensemble mean, decadal-scale ocean heat anomalies and sea level pressure anomaly distributions, and similar results from reanalysis products for the recent global hiatus (1999–2012) minus global warming (1976–1998) period.

Although this analysis indicates that the World Ocean can absorb and store sufficient amounts of heat to compensate for 10 year periods without surface warming under long-term global warming, it can only be taken as one of several contributing factors explaining the recent slow-down in global surface warming. A follow-up on this is that the surface temperature is certainly not the first choice for a reliable measure of the total heat content in the climate system. Here the ocean heat content is a much more robust quantity.

It is also worth mentioning that for all RCP scenarios, the computed probability of weak positive or even negative decadal-scale surface temperature trends is likely underestimated since these simulations do not take into consideration volcanic eruptions and variations between solar cycles, two factors contributing to externally-driven, short-term fluctuations in global surface temperature.

The low latitude region encompassing the tropical Atlantic, the Indian Ocean and the West Pacific is the first region to show a decreasing likelihood of decadal-scale non-warming periods with global warming. The North Atlantic and the Southern Ocean are the regions with large variability, and thus largest likelihood of non-warming decades in a warming world.

Finally, the generally large model spread in the tropical, upper ocean heat content is a challenge and should, to the extent possible, be understood and reduced. This is a well known problem (e.g., Michael et al. 2013; Taschetto et al. 2014), but is worth highlighting because of the leading role of the low latitudes in regulating the global heat transfer. Close interaction with the ongoing decadal predictability efforts (e.g., Goddard et al. 2013) might be a valuable avenue in this respect.