Internal variability of Earth's energy budget simulated by CMIP5 climate models

Developments in Earth system observing capability and the unexpectedly modest rise in global surface temperature over the last decade or so have prompted a number of research efforts to better understand Earth's energy budget, and its relationship to changes in near surface temperature and upper ocean heat content (e.g. Palmer et al 2011, Katsman and van Oldenborgh 2011, Meehl et al 2011, Loeb et al 2012). Satellite observations from CERES (Clouds and Earth's Radiant Energy System, Wielicki et al 1996) allow us to monitor changes in the net radiation at top-of-atmosphere (TOA) from 2001 onwards. In addition, the Argo array of ocean profiling floats provide estimates of near-global ocean heat content (OHC) for the upper 2 km of the open ocean since about 2005 (Roemmich et al 2009). While inter-annual variability in TOA (from CERES) and OHC changes (from Argo) show some agreement (Loeb et al 2012), it is difficult to assess from observations how important deep (>2 km) OHC changes are, both in terms of annual-to-decadal variability and the long-term response to anthropogenic forcing of the climate system (Purkey and Johnson 2010).

Palmer et al (2011, hereafter referred to as 'P11') used pre-industrial control simulations from three generations of climate model from the Met Office Hadley Centre to show that unforced decadal trends in sea surface temperature are only weakly indicative of TOA over the same period. This result arises because the ocean is capable of vertically re-distributing substantial quantities of heat on decadal timescales. In contrast, P11 found that trends in total OHC strongly constrain TOA, since the ocean represents Earth's primary heat store on decadal and longer timescales (Bindoff et al 2007). Meehl et al (2011) showed that decades of near-zero surface temperature rise in simulations of future climate were associated with enhanced heat uptake below 300 m and a transition to La Niña conditions in the tropical Pacific. Recent simulations from an ocean data assimilation model suggest a recent acceleration in ocean heat uptake below 300 m that seems to be partly driven by surface wind forcings and associated changes in upper ocean circulation (Balmaseda et al 2013). In addition, hindcasts from a decadal prediction system were able to reproduce the onset of the recent pause in observed surface warming (Easterling and Wehner 2009, Knight et al 2009, Foster and Rahmstorf 2011), which was accounted for by an increase in ocean heat uptake over the upper 700 m in the Pacific and Atlantic oceans (Guemas et al 2013). The same study found no role for external forcings in simulating this event, implying that internal climate variability is the primary explanation.

In this work, we apply a similar analysis to that conducted by P11 to a large number of pre-industrial control simulations from the CMIP5 model archive. In addition to consideration of a much larger number of models than the previous study, we address the following questions related to Earth's energy budget: (i) on what timescale does TOA come into balance with changes in total OHC?; (ii) what is the potential role of internal variability in TOA for the recent pause in surface warming?; (iii) what is the potential role of ocean heat redistribution for the recent pause in surface warming? In section 2 we describe the model simulations and the pre-processing methods. The statistical methods and main results are presented in section 3. Finally, in section 4, we discuss the wider implications of our results and summarize our findings.

The model data analysed here are multi-century pre-industrial control simulations ('piControl') from the CMIP5 archive (Taylor et al 2012). The piControl simulations fix atmospheric constituents at pre-industrial levels, circa 1850, and here their primary use is to assess the role of unforced internal variability on Earth's energy budget on inter-annual to multi-decadal timescales. The models include all of the major terms in the energy budget (Bindoff et al 2007), including the cryosphere, ocean, atmosphere and land heat storage, but generally have only a crude representation of ice sheets (e.g. Collins et al 2011). Mountain glaciers are not represented in the CMIP5 models and continental heat storage may be underestimated (Stevens et al 2007). However, observation-based estimates suggest that these two components account for only 1% and 2% of the change in Earth's heat storage over the period 1972–2008, respectively (Church et al 2011).

The following monthly mean model fields are used in our investigation, with the relevant CMIP5 variable name included in parentheses: (i) top-of-atmosphere incoming shortwave radiation ('rsdt'); (ii) top-of-atmosphere outgoing shortwave radiation ('rsut'); (iii) top-of-atmosphere outgoing longwave radiation ('rlut'); (iv) near surface air temperature ('tas'); and (v) ocean potential temperature on all model levels ('thetao'). Variables (i)–(iii) are globally integrated and combined to give TOA and then time integrated to give time series of the total energy in the Earth system (TE). Note that the trend in TE is essentially synonymous for the time-averaged TOA over a given period and we will use the terms interchangeably. Variable (iv) is spatially averaged to provide time series of global surface temperature (GST). Variable (v) is converted to ocean heat content time series for each model vertical level by spatially integrating and multiplying by reference values for sea water density and heat capacity of 1025 kg m³ and 3985 J kg⁻¹ K⁻¹. In order to carry out spatial integration or averaging grid cell areas for the atmosphere ('areacella') and the grid cell volumes for the ocean ('volcello') were also retrieved from the CMIP5 archive. Where volcello files were not available we made use of the ocean grid cell areas ('areacello') and used the nominal layer thicknesses to estimate the cell volumes. See supplementary materials (available at stacks.iop.org/ERL/9/034016/mmedia) for further details.

As discussed by Sen Gupta et al (2013), climate model simulations can be subject to 'drift': spurious long-term changes that are unrelated to either changes in external forcing or internal low-frequency variability. Since it can take many centuries for the deep ocean to reach equilibrium, it is important to remove any associated model drift from our analyses. Following the approach of P11, we use a Butterworth filter with a cut-off frequency of 100 years to high-pass the global time series of GST, TE and OHC (figure S1 available at stacks.iop.org/ERL/9/034016/mmedia). This process inevitably lead to a trade-off between elimination of data due to filter edge effects versus elimination of genuine low-frequency variability from the models. Since we focus mainly on decadal timescales, our results are not sensitive to using a longer cut-off, such as 150 or 200 years. The resulting high-passed time series are effectively anomaly time series relative to a 'time local' reference state. In all subsequent analysis, we remove 50 years from each end of the time series to avoid filter edge effects. Finally, we remove the mean annual cycle from each time series to prevent any aliasing of trends in the data at sub-annual timescales. Before any further analysis of the data is undertaken, we check the consistency between the high-passed time series of TE and total OHC (figure S2 available at stacks.iop.org/ERL/9/034016/mmedia).

3.1. Surface temperature, ocean heat content and Earth's energy budget

To investigate the relationships between TE, GST and total OHC on decadal timescales we compute linear trends for overlapping 10-year periods and plot the regressions (figures 1 and 2). As previously demonstrated by P11, we find that in general decadal trends in GST generally provide only a weak constraint on the trend in TE for the same period (figure 1). The correlation coefficients range from 0.08 (GISS-E2-R) to 0.63 (HadGEM2-CC). Those models with a larger range of trends in GST tend to show a stronger correlation between the two variables. This relationship is significant at the 5% level with an average increase in correlation of 0.25 for each 0.1 K increase in GST trend range. In contrast, decadal trends in total OHC generally provide a strong constraint on corresponding trends in TE (figure 2); consistent with the ocean being Earth's primary heat store. In this case the correlation coefficients range from 0.81 (IPSL-CM5A-LR) to 0.99 (GISS-E2-R). The lower correlations may arise partly from inaccuracies in ocean grid cell volumes (see supplementary materials available atstacks.iop.org/ERL/9/034016/mmedia) or inadequate removal of model drift. Two-thirds of the models show a correlation between TE and total OHC of 0.95 or greater on decadal timescales. It is implicit that for those models with the lowest correlation coefficients, such as the MPI models, that other terms in the planetary heat budget, such as the cryosphere, atmosphere and land surface, play a more important role on decadal timescales. Given the strong relationship between TE and total OHC on decadal timescales, a relevant question for climate monitoring is: on what timescale does the ocean become the dominant term in Earth's energy budget?

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Following a similar approach to the analysis above, we find all overlapping trends of length 2–36 months and fit a linear model for the relationship between trends in TE and total OHC. The resulting correlations between trends in TE and OHC for a range of periods provides an estimate of the timescale on which the ocean becomes the dominant term in the planetary energy budget (figure 3). For most models, the correlation between trends in TE and OHC increases to near saturation over the first 12 months or so. The MPI models take the longest for trends in TE and OHC to become highly correlated and also show some of the weakest relationships between the two variables on decadal timescales. This seems a logical result if for those models other energy storage terms are more important, as discussed earlier. In any case, the indicative timescale of approximately 12 months for most models supports the attempts to cross-validate annual averages of TOA with annual changes in ocean heat storage (e.g. Loeb et al 2012).

Zoom In Zoom Out Reset image size

3.2. Multi-annual to decadal variability in GST, TE and OHC

In order to address the potential role of internal climate variability for the observed pause in global surface temperature rise over the last decade or so, we quantify the range of trends of a number of key variables for trend periods between 2 and 20 years. For all variables we determine the 2.5th and 97.5th percentile trends empirically and also compute the maximum overall trend for any model for each period. These are approximately equivalent to a 1 in 50-year event and roughly commensurate with one warming pause event in about 50 years that isn't associated with major radiative forcing, such as a volcanic eruption (Morice et al 2012). In the following discussions all energy fluxes expressed in W m⁻² are relative to the surface area of the Earth.

Similar to previous studies (Easterling and Wehner 2009, Knight et al 2009), we find considerable variability in GST for the CMIP5 models at decadal timescales (figure 4(a)). As one might expect, the GST trend magnitudes initially decay rapidly as the period increases from 2 to 5 years (the typical timescale of ENSO variability) but there remains substantial variability out to 20 years. For 10-year and 15-year periods the mean 2.5th (97.5th) percentile trend is −0.29 (+0.28) and −0.18 (+0.17) K dec⁻¹, respectively. The largest negative (positive) trends we see for any CMIP5 model at these periods is −0.60 (+0.79) and −0.37 (+0.34) K dec⁻¹, respectively.

Zoom In Zoom Out Reset image size

Compared to GST, trends in TE decay less rapidly as a function of time period and there seems to be greater spread among the models (figure 4(b)). Even at periods of 20 years, GFDL-CM3 shows trend magnitudes greater than 0.15 W m⁻². This model shows quasi-oscillatory behaviour in TE with a period of approximately 100 years (see supplementary materials available at stacks.iop.org/ERL/9/034016/mmedia). For 10-year and 15-year periods the mean 2.5th (97.5th) percentile trend is −0.14 (+0.13) and −0.11 (+0.10) W m⁻², respectively. The largest negative (positive) trends we see for any model in the archive at these periods is −0.28 (0.28) and −0.23 (0.22) W m⁻², respectively.

In addition to trends in GST and TE, we investigate the variability in heat fluxes across two isobaths: 700 and 1800 m (figures 4(c) and (d)). The fluxes are computed as the trends in OHC from 700 m-to-bottom and 1800 m-to-bottom. The depths are informed by the historical in situ observations of ocean temperature, which are dominated by expendable bathythermographs (XBTs) from the 1970s to mid-2000s and the Argo array of profiling floats since the mid-2000s (Abraham et al 2013). The depth sampling associated with these instruments means that most estimates of multi-decadal OHC change are for the 0–700 m layer (e.g. Lyman et al 2010) and estimates since the early-to-mid-2000s are generally limited to the 0–1800 m layer, or thereabouts (e.g. von Schuckmann et al 2009). Positive trends represent an export of heat to the deeper ocean and are therefore the primary focus in the context of the recent pause in surface temperature rise.

The mean 2.5th (97.5th) percentile trends in heat flux across the 700 m isobath −0.09 (+0.09) and −0.07 (+0.07) W m⁻², for a 10-year and 15-year period, respectively. The largest negative (positive) trends we see for any model at these periods is −0.25 (+0.26) and −0.19 (+0.19) W m⁻², respectively. For the heat flux across the 1800 m isobath the 10-year and 15-year period the mean 2.5th (97.5th) percentile trends are −0.06 (+0.05) and −0.05 (+0.04) W m⁻², respectively. The largest negative (positive) trends we see for any model at these periods is −0.22 (+0.16) and −0.22 (+0.11) W m⁻², respectively. We note that there is a strong asymmetry between the maximum and minimum heat fluxes across 1800 m, which indicates that the dominant processes favour heat release from the deep ocean rather than sequestration on decadal timescales. An example could be enhanced (suppressed) convection in the North Atlantic subpolar gyre, which tends to cool (warm) the deeper ocean (Katsman and van Oldenborgh 2011).

3.3. OHC as a predictor for TE on decadal timescales

Following P11, we carry out similar regressions to that illustrated in figure 2 for decadal trends in OHC integrated over successively deeper model levels against corresponding decadal trends in TE. For each model we treat OHC trends as a predictor of the total energy trend over the same decade and express this relationship in terms of uncertainty (90% prediction interval) against OHC integration depth (figure 5). We use non-overlapping decades, in order to minimize the effect of autocorrelation on the estimate of the prediction intervals. Visual inspection of the residuals of a large sample of the regressions confirms little or no structure in the residuals. Part of the uncertainty in the predictions arises from ocean heat rearrangements, which determines the vertical structure of uncertainty as a function of integration depth. There is also a depth independent contribution to uncertainty from the other components of the energy budget—such as the cryosphere, land surface and atmosphere—that remains even when integrating OHC over the full column (indicated by the values intersecting the x-axis on figure 5).

Zoom In Zoom Out Reset image size

We see a wide variety of model behaviour for the predictive uncertainty as a function of depth (figure 5). All models show a constant or near-constant uncertainty in the uppermost model levels. This result is consistent with our expectations of an ocean mixed layer, where trends should be well correlated with depth (and therefore fairly constant for the prediction interval), despite the fact that we are probably integrating over large local variations in the simulated ocean mixed layer. The depth of this 'effective' mixed layer is strongly model dependent. For example, the strong reduction in uncertainty deeper than 100 m seen for the GFDL models suggests that the effective mixed layer depth in these models is approximately 100 m. Conversely, the IPSL models suggest that the effective mixed layer depth is more like 400 m.

As we integrate OHC with depth many of the models show intervals for which prediction interval increases with depth. This can be explained by many models exhibiting a number of correlated and anti-correlated layers in their OHC variability (figures S3–5 available at stacks.iop.org/ERL/9/034016/mmedia). Thus, as one integrates over just part of a vertical dipole (or multipole) in OHC anomaly, the uncertainty can actually increase with depth. The wide range of behaviours seen in both figures 5 and S3–5 (available at stacks.iop.org/ERL/9/034016/mmedia) is a striking result and illustrates how differently coupled model systems can behave despite being built on the same fundamental physics. Even those models from the same modelling centre can show substantial differences in the profile and magnitude of OHC variability.

Despite the range of model behaviours, all models support the notion that more complete sampling over the full ocean depth can dramatically reduce the uncertainty in predicting changes in TOA—consistent with the findings of P11. On a 10-year timescale, the 90% prediction interval for the trend in TE given a trend in OHC integrated over the 0–700 and 0–1800 m layers ranges from 0.07–0.16 to 0.05−0.12 W m⁻², respectively. As an example, for a 90% prediction interval of 0.1 W m⁻² the true value of trend in TE should lie within approximately ± 0.05 W m⁻² for 9 out of 10 decades.

In section 3 we have presented a number of results that describe the internal variability in Earth's energy budget as simulated by CMIP5 pre-industrial control simulations. One of the primary motivations for this work is to assess the potential role of unforced internal variability as an explanatory factor during periods of reduced surface temperature rise. Similar to previous studies (Easterling and Wehner 2009, Knight et al 2009), we find that typical 1-in-50-year trends in GST arising from internal variability (∼0.3 K dec⁻¹) are large compared to the 'background' warming rate of approximately 0.2 K dec⁻¹.

To illustrate the relevance of the present study to the question of Earth's energy budget, we consider the following example. For the period 2000–2010, Stephens et al (2012) have used satellite observations to estimate an energy imbalance at top-of-atmosphere of 0.6 ± 0.4 W m⁻², constrained by estimates of global ocean heat uptake based on the results of Lyman et al (2010). Over this same period it has been estimated that changes in stratospheric water vapour account for − 0.1 W m⁻² to the rate of change of TE in the Earth System (Solomon et al 2011). Kaufmann et al (2011) estimated a − 0.06 W m⁻² contribution to the planetary heating associated with increased sulfate emissions and a further − 0.18 W m⁻² associated with the prolonged solar minimum, for the period 2002–2007 (we note that the solar cycle transitioned into the increasing phase of solar output before 2010). Thus the order-of-magnitude variations in some of the key external radiative forcings during the first decade of the 21st century are 0.1 W m⁻². Our analysis of CMIP5 models shows that a representative 2.5th and 97.5th percentile decadal trend in TE has a magnitude of approximately 0.1 W m⁻², and the largest trends in any model having a magnitude of up to 0.3 W m². Thus, taking the CMIP5 results at face value, we suggest that internal variability in TOA fluxes on a 10-year timescale may be equally important as some of the radiative forcing agents that have been invoked as contributory factors to the recent pause in global surface warming.

Estimates of OHC change for the 0–700 m or 0–1800 m layers provide an important constraint on the rate of accumulation of energy in the Earth system (e.g. Murphy et al 2009, Stephens et al 2012) because the ocean represents the Earth's primary heat store (Bindoff et al 2007). Ocean temperature observations are very limited below 700 m prior to the mid-2000s, and since then observations still remain sparse deeper than 1800 m or so (Palmer et al 2010, Abraham et al 2013). Observational estimates of ocean heat uptake for the ocean below these depths are based on sparse time-and-space-sampled hydrographic cruise data (Purkey and Johnson 2010). The CMIP5 simulations analysed here suggest that fluxes into the deeper ocean across the 700 and 1800 m isobaths arising from internal variability alone can be of order 0.1 W m⁻² on decadal timescales. The largest fluxes we find for any model across these isobaths exceed 0.2 W m⁻², which is about one third of the top-of-atmosphere radiation imbalance estimated by Stephens et al (2012) for the period 2000–2010.

The analysis presented here highlights the potential importance of internal variability in Earth's energy budget on multi-annual to decadal timescales—both in terms of changes in the top-of-atmosphere energy budget and ocean redistribution of heat. We find that the ocean becomes the dominant energy store in the Earth System on a timescale of about one year and changes in global OHC provide a more reliable indication of TOA variability than GST. However, the spatial patterns and ultimately the processes controlling this internal variability warrant further attention. In particular, understanding of the spatial patterns and associated mechanisms from model simulations would promote greater insight into processes behind the recent pause in surface warming. Consideration of these spatial patterns may also help to account for deficiencies in model simulations by the scaling up or down of model signals to match observed amplitudes. Ultimately, we would like to understand the mechanisms for ocean heat rearrangement in the models (e.g. Meehl et al 2013) and relate this back to quantitative statements about the real world. This is a challenging problem to tackle, and must be done with a good understanding of the limitations of current climate models to simulate key ocean processes, such as Antarctic Bottom Water formation (Heuze et al 2013).

We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modelling groups (listed in Table S1 in supplementary materials available at stacks.iop.org/ERL/9/034016/mmedia) for producing and making available their model output. For CMIP the US Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. Ian Edmond and Jamie Kettleborough provided scripts to help with downloading and archiving of CMIP5 climate model data. We are grateful to Rowan Sutton, Jonathan Gregory and John Church for useful discussions on this work. Two anonymous reviewers provided helpful comments on an earlier draft of the manuscript. This work was supported by the Joint UK DECC/Defra Met Office Hadley Centre Climate Programme (GA01101) and represents a Met Office contribution to the Natural Environment Research Council DEEP-C project NE/K005480/1.