Earth's energy imbalance since 1960 in observations and CMIP5 models.

Observational analyses of running 5 year ocean heat content trends (Ht) and net downward top of atmosphere radiation (N) are significantly correlated (r ∼ 0.6) from 1960 to 1999, but a spike in Ht in the early 2000s is likely spurious since it is inconsistent with estimates of N from both satellite observations and climate model simulations. Variations in N between 1960 and 2000 were dominated by volcanic eruptions and are well simulated by the ensemble mean of coupled models from the Fifth Coupled Model Intercomparison Project (CMIP5). We find an observation-based reduction in N of - 0.31 ± 0.21 W m-2 between 1999 and 2005 that potentially contributed to the recent warming slowdown, but the relative roles of external forcing and internal variability remain unclear. While present-day anomalies of N in the CMIP5 ensemble mean and observations agree, this may be due to a cancelation of errors in outgoing longwave and absorbed solar radiation. KEY POINTS: Observed maximum in ocean heat content trend in early 2000s is likely spuriousNet incoming radiation (N) reduced by 0.31 ± 0.21 W m-2 during the warming pausePresent-day estimates of N may contain opposing errors in radiative components.

1. Introduction

The recent slowdown in the rate of global surface temperature rise [Easterling and Wehner, 2009, Knight et al., 2009] has highlighted a need for a better quantification of Earth's energy imbalance [Hansen et al., 2011; Katsman and van Oldenborgh, 2011; Balmaseda et al., 2013; Trenberth et al., 2014] and improved understanding of the flow of energy in the climate system [Trenberth, 2009; Meehl et al., 2011; England et al., 2014]. Earth's energy imbalance is defined here as the net downward radiative energy flux (N, expressed in W m−2 averaged over the entire surface of the Earth) through the top of the atmosphere (TOA), representing the residual between incoming solar radiation and outgoing reflected (mainly short-wavelength visible) and emitted (mainly long-wavelength infrared) radiation. N varies naturally on all timescales due to internal variability of the climate system, volcanic eruptions, and variations in incoming solar radiation. However, increased levels of greenhouse gases initially increase N by reducing outgoing longwave radiation, driving climate change until radiative equilibrium is restored and N returns to zero on average. Humans also influence N by changing outgoing shortwave radiation through emissions of aerosol particles and changes to land usage which affect surface albedo. Furthermore, subsequent changes in N in response to initial perturbations may be amplified or reduced by feedbacks within the climate system [Boucher et al., 2013; Myhre et al., 2013]. Quantifying the time mean and variability of N is therefore fundamental for understanding climate variability and change [Trenberth, 2009; Hansen et al., 2011; Stephens et al., 2012; Trenberth et al., 2014].

Allan et al. [2014] produced a homogenized satellite data set of N covering the period since 1985. Satellites provide reasonably accurate estimates of interannual variations in N [Loeb et al., 2012], but absolute values cannot be determined accurately [Stephens et al., 2012]. However, on multiyear timescales, N is expected to be largely balanced by changes in global ocean heat content [Palmer et al., 2011; Palmer and McNeall, 2014], since changes in energy of the atmosphere, land, and cryosphere are much smaller [Levitus et al., 2001; Hansen et al., 2011; Trenberth et al., 2014]. Observations of trends in ocean heat content (hereafter Ht) therefore potentially enable the satellite observations to be anchored to obtain absolute values [Loeb et al., 2012; Allan et al., 2014], as well as providing independent estimates of variations in N. However, the subsurface ocean was poorly observed historically, with less than 20% of 1° latitude by 1° longitude regions sampled yearly in the upper 700 m before 1990 [Abraham et al., 2013], and estimates of Ht from different ocean analyses appear to show little consistency, at least on interannual timescales [Trenberth et al., 2014]. Observational coverage increased to 30–40% of 1° regions after 1990, but observation-based estimates of Ht for the period 1998 to 2008 vary by a factor of 2, from 0.31 (full depth ocean [Church et al., 2011]) to 0.64 W m−2 (upper 700 m only, [Lyman et al., 2010]). Even in the period since 2005 in which more than 70% of 1° regions are sampled by Argo floats [Roemmich et al., 2009], estimates of Ht in the upper 700 m range from 0.16 to 0.39 W m−2 over the period 2005–2012 [Abraham et al., 2013]. Furthermore, annual values of Ht appear to be inconsistent with satellite observations of N in some years [Trenberth et al., 2014], suggesting that Earth's energy budget cannot be closed with observations even over the relatively well-observed period since 2005 [Trenberth and Fasullo, 2010].

Here we investigate Earth's energy imbalance over the period since 1960 in observations, atmospheric models forced by observed sea surface temperatures and radiative changes (Atmosphere Model Intercomparison Project, AMIP), and fully coupled model simulations (Coupled Model Intercomparison Project, CMIP5). We build on Allan et al. [2014] by including estimates of Ht and developing an observation-based record of N extending back to 1960. We show that two reanalyses of Ht in the early 2000s are inconsistent with satellite observations and AMIP simulations of N and are therefore likely unreliable, helping to explain our inability to close the energy budget. However, prior to 2000 variations in Ht in ocean reanalyses are in reasonable agreement with observation-based estimates of N and with CMIP5 models, as Earth's energy imbalance was dominated by volcanic eruptions over this period. We also highlight an observation-based reduction in N that likely contributed to the recent slowdown in surface warming.

2. Data and Methods

Observational estimates of N since 1985 are based on the Clouds and the Earth's Radiant Energy System (CERES) instruments since 2000 and the Earth Radiation Budget Satellite (ERBS) wide field of view nonscanning instrument from 1985 to 1999 [Allan et al., 2014]. These two data sets were homogenized, and gaps filled, using the ERA Interim reanalysis [Dee et al., 2011] and high-resolution atmosphere model simulations driven by observed sea surface temperatures and sea ice concentrations (see Figure S1 in the supporting information and Allan et al. [2014] for full details). This satellite-based N is adjusted to be 0.58 W m−2 over the period July 2005 to June 2010 [Loeb et al., 2012] based on ocean heat content changes (Ht) of 0.47 W m−2 in the upper 1800 m obtained from a linear fit through Argo data [following Lyman and Johnson, 2008] and 0.07 W m−2 below 1800 m [Purkey and Johnson, 2010], plus 0.04 W m−2 estimated for heating and melting of ice and heating of the land and atmosphere.

We also use N from AMIP model simulations form (i) the UPSCALE project [Mizielinski et al., 2014] over the period 1986 to 2011 and (ii) multimodel simulations of the climate of the twentieth century project (C20C) [Scaife et al., 2009] for the period 1950 to 2002. Three of the 16 C20C models were eliminated from our study due to unrealistic anomalies of N. AMIP simulations cannot provide absolute values of N because the sea surface temperature (SST) forcing potentially provides unrealistic sources or sinks of energy. The AMIP fluxes were therefore anchored to the Allan et al. [2014] analysis with a constant adjustment to make their global means equal during overlapping periods.

AMIP simulations are in good agreement with satellite observations (r > 0.6) for monthly deseasonalized variability in N over the period 1986 to 2011 (Allan et al. [2014] and Figure S1), showing that SST observations provide useful additional constraints on N. We therefore create an extended observation-based estimate of N (hereafter No) by averaging the satellite observations with the ensemble means of the C20C and UPSCALE simulations (with each data set given an equal weight where available). We use No even when satellite data are available since errors are likely to be reduced by including independent information from the SST observations (via the AMIP simulations). We compute uncertainties from differences between UPSCALE, C20C, and satellite data (0.12 W m−2) as well as the average spread of the C20C simulations (0.09 W m−2) added in quadrature to the 1 sigma anchoring uncertainty (0.24 W m−2) estimated in Allan et al. [2014] to obtain a total uncertainty on No of 0.28 W m−2.

Ocean heat content trends (Ht) are computed from the slope of a least squares regression using annual time series of heat content and expressed as a flux relative to the total surface area of the Earth [following Lyman and Johnson, 2008]. We compute Ht over 5 year periods and compare with rolling 5 year mean N, thereby smoothing out interannual variability related to El Niño. Absolute values of N cannot be determined accurately from satellite observations [Stephens et al., 2012]. We therefore also assess rates of change of Ht and N (Htt = dHt/dt and Nt = dN/dt, again computed with least squares regression over 5 year periods), since Nt does not depend on absolute values.

Ocean data are taken primarily from the Met Office statistical ocean reanalysis (hereafter MOSORA, updated from Smith and Murphy [2007]; see supporting information). Several published analyses of Ht for the upper 700 m appear to capture signals associated with major volcanic eruptions, but the timing is not always consistent, and there is little agreement at other times [Trenberth et al., 2014]. However, the ORA-S4 analysis, in which temperature and salinity observations are assimilated into a dynamical ocean model [Balmaseda et al., 2013], was found to be in reasonable agreement with model simulations of volcanic impacts and also simulated plausible changes in Ht associated with El Niño [Trenberth et al., 2014]. For these reasons, and because full-depth values are available, we also include ORA-S4 in our study.

Time-varying biases in expendable bathythermographs (XBTs) and mechanical bathythermographs were corrected in MOSORA based on an average of adjustments proposed by Wijffels et al. [2008], Levitus et al. [2009], and Ishii and Kimoto [2009]. MOSORA uses global covariances to interpolate between observations (Table S1). This approach was tested, and uncertainties estimated, by performing data withholding experiments for the data-rich period 2007–2013, which was not used for the covariances. In these experiments, analyses created using all available observations were taken as the reference and compared with analyses created using subsampled observations at locations typical of historical periods. A total of 12 tests were performed using data locations from 1950 to 1956, 1955 to 1961, 1960 to 1966,…, 2005 to 2011. Patterns of Ht reconstructed with historical observations agree well with the reference (Figure1). For example, the spatial correlation is 0.88 and 0.76 using the 1990s and 1960s observations, respectively (compare Figure1b with Figures1d and 1f). Uncertainties (which we quote as 1 sigma throughout) in global heat content range from around 28 ZJ (ZJ = 1021 J) using 1950–1956 observations to 8 ZJ using 2005–2011 observations, corresponding to uncertainties in Ht of 0.7 to 0.2 W m−2 (Table S1). We stress that these estimates only assess uncertainties arising from the observational coverage relative to the recent Argo era and do not include unknown instrument errors or the potentially large contributions from regions not sampled adequately by Argo floats, such as the shelf seas, ice-covered regions, and ocean below 2000 m [e.g., von Schuckmann et al., 2014]. Furthermore, our subsampling experiments suggest that the magnitude of Ht is not fully captured using the historical data, especially in the Southern Ocean (Figure1f) so that Ht might be underestimated in the historical period prior to Argo, as suggested by Durack et al. [2014].

Figure 1

Assessment of Met Office ocean reanalysis. (left) Total number of temperature observations at 1000 m depth during the period 2008 to 2012 using (a) all data and subsampled data typical of (c) the 1990s and (e) the 1960s. (right) Patterns of ocean heat content trends (Ht, W m−2) for the period 2008 to 2012 computed using (b) all observations and subsampled observations typical of (d) the 1990s and (f) the 1960s.

We compare observational estimates of Ht and N with 21 CMIP5 coupled model simulations (Taylor et al. [2012] and Table S2) driven by observed variations in greenhouse gases, aerosols, solar radiation, and land use until 2005 and following representative concentration pathway (RCP4.5) thereafter. There are long-term drifts in the CMIP5 models because the preindustrial control simulations were not necessarily in radiative equilibrium [Sen Gupta et al., 2013]. We therefore compute anomalies in the historical and RCP4.5 simulations relative to preindustrial control simulations (using at least 100 years of the control simulations).

3. Observation-Based Estimates

Global ocean heat content, Ht and Htt from MOSORA and ORA-S4 are in good agreement (r = 0.94, 0.80, and 0.82 respectively for running 5 year averages) over their overlapping period (compare red and dashed magenta curves in Figure2), despite the lack of historical ocean observations and additional uncertainties involved in computing derivatives for Ht and Htt. Furthermore, over the period 1960–1999, Ht and Htt are in broad agreement with N (r = 0.58 and 0.56 between MOSORA and ORA-S4 and No) and Nt, (r = 0.63 and 0.69). These correlations are significantly greater than zero with 99% confidence based on bootstrapping in 5 year blocks to take serial correlation into account [e.g., Goddard et al., 2012].

Figure 2

Observed energy budget. (first panel) Annual mean total ocean heat content anomalies (relative to 1961–2010, ZJ = 1021 J) from MOSORA (red) [Smith and Murphy, 2007] and ORA-S4 (dashed magenta) [Balmaseda et al., 2013]. Black asterisks show MOSORA omitting Argo data. (second panel) Five year running mean ocean heat content trends (Ht, red and magenta) and net incoming radiation (N) at the top of the atmosphere from satellites [Allan et al., 2014] (blue curves) and atmosphere models driven by observed surface temperatures for the climate of the twentieth century (C20C) [Scaife et al., 2009] (cyan) and UPSCALE [Mizielinski et al., 2014] (blue shading) projects. (third panel) Trends in Ht (Htt, red and magenta) and N (blue and cyan). Red shading shows uncertainties (1 sigma) in MOSORA arising from observational coverage (see text and Table S1). Grey shading highlights potentially spurious ocean reanalyses (see text), with significant differences (90% level) between MOSORA and N shown by black squares. In Figure2 (second to third panels) the year represents the center of the time period used to compute trends.

After 2000, Ht and Htt show pronounced variability that is not present in N and Nt (Figure2, second and third panels). This is largely caused by the rapid increase in ocean heat content between 2001 and 2004 resulting in a peak Ht in 2002 followed by a trough in 2006–2007 in both MOSORA and ORA-S4 (Figure2, second and third panels). A peak in Ht around 2002 is also evident in other ocean analyses [Trenberth et al., 2014]. However, such an increase in heat content is not supported by thermosteric sea level observations, because total sea level increased at a fairly constant rate [Church et al., 2013] while freshwater input from glaciers and ice sheets likely increased [Marzeion et al., 2012; Shepherd et al., 2012; Church et al., 2013], implying a reduction in the rate of thermosteric sea level rise (and hence ocean heat content trend) over the period 2001 to 2004. Although interpretation of sea level changes is complicated by the nonuniformity of thermal expansion coefficients, warming is not consistent with a reduction in thermosteric sea level trends. Given this, and the agreement between independent estimates of N (and Nt) from satellite observations and AMIP simulations (Figure2), and from forced ocean model experiments (Figure S2), we suggest that there are likely to be errors in the observational reanalysis of ocean heat content between 2001 and 2007 (grey-shaded region in Figure2). Cheng and Zhu [2014] drew similar conclusions from a different approach, suggesting that these errors are caused by the transition from XBT to Argo data. However, this does not appear to be the sole cause in our study since MOSORA shows a similar peak in Ht in 2002 even when Argo observations are omitted (Figure2, black asterisks), although the subsequent trough is reduced. Further work is therefore needed to understand the cause of these errors in ocean reanalyses. Possible explanations include inadequate vertical or horizontal sampling and/or unresolved biases in the instruments themselves.

4. Comparison With CMIP5 Model Simulations

CMIP5 coupled models simulate variations in N and Ht that are in good agreement with observations (r = 0.82 between 5 year running mean CMIP5 ensemble mean and No, and 0.68 with MOSORA excluding the period after 2000, Figure3a). Since internal variability is largely removed in the ensemble mean of the CMIP5 models, this agreement suggests that much of the multiyear variability in Earth's energy imbalance from 1960 to 2000 was externally forced by the volcanic eruptions of Agung (1963), El Chichon (1982), and Pinatubo (1991), as noted by Trenberth et al. [2014]. Furthermore, the CMIP5 models do not simulate a peak in Ht in 2002, and the observed peak is significantly different from the CMIP5 models, consistent with our earlier conclusion that estimates of Ht based on ocean analyses between 2001 and 2007 (grey shading in Figure2) are likely unreliable.

Figure 3

Earth's energy imbalance. (a) Time series of 5 year running mean N and Ht (as Figure2, second panel) for 21 CMIP5 coupled model simulations (N in green, Ht in orange, ensemble mean in thick lines) compared with Ht from MOSORA (red) and No (blue, see text). Black squares (diamonds) show where differences between MOSORA and No (CMIP5) are significant with 90% confidence. (b) N averaged over different periods in No (blue, with 1 sigma uncertainties) compared to the CMIP5 models (green, box showing the mean ±1 sigma and whiskers showing the range) and estimates from the IPCC fifth assessment (red) [Rhein et al., 2013, Box 3.1]. Numerical values are given in Table S3.

Averaged over different periods, No is in broad agreement with the Intergovernmental Panel on Climate Change (IPCC) estimates [Rhein et al., 2013, Box 3.1] and the CMIP5 models (Figure3b), but it is difficult to draw firm conclusions because uncertainties in the absolute values of No are large. These uncertainties are dominated by the anchoring estimate of Ht and could potentially be reduced by averaging Ht over longer periods. We do not attempt this because ocean observations in the early 2000s are potentially unreliable (as discussed above), historical analyses do not capture the full magnitude of Ht, especially in the Southern Ocean (Figure1) [Durack et al., 2014], and uncertainties are evident even in the most recent 5 year means which show a decrease in No but an increase in Ht (Figure3a). We therefore investigate changes in N and anomalies relative to the period 1960–2011 since these are less affected by anchoring uncertainties.

In general, there is good agreement between the CMIP5 ensemble mean and observation-based anomalies of N (r = 0.82 for 5 year running means) and its components, absorbed shortwave (ASR, r = 0.87) and outgoing longwave (OLR, r = 0.80) radiation (Figure4). The CMIP5 ensemble mean therefore appears to simulate a realistic magnitude of variability associated with volcanic eruptions. The largest discrepancy between models and observations occurs during the 1970s (Figure4). The mid-1970s shift in the Pacific Decadal Oscillation [Mantua et al., 1997; Power et al., 1999] may have contributed to No. However, an increase in N is also apparent in the CMIP5 ensemble mean in the late 1970s (especially in ASR; red curves in Figure4), suggesting that external forcing may have played a role. Further work is therefore needed to understand these apparent differences between models and observations.

Figure 4

Components of radiative fluxes. Time series of 5 year running mean anomalies (relative to 1960 to 2011) of TOA-absorbed shortwave radiation (ASR, red), outgoing longwave radiation (OLR, green), and net radiation (N = ASR-OLR, blue) in No (dashed) and the CMIP5 models (solid, with thick line showing the ensemble mean).

Models and observations both show a rapid increase in N during the recovery from the eruption of Mount Pinatubo, reaching a peak in the late 1990s. However, there is a subsequent decline of −0.43 ± 0.21 W m−2 in No between 1999 and 2005 in the early part of the global warming slowdown which is larger in magnitude than the model ensemble mean (−0.12 ± 0.06 W m−2) and at the lower end of individual models (+0.38 to −0.32 W m−2). Both the satellite observations and the AMIP simulations show a decline over this period (Figures2, middle, 2 bottom, and S3), showing that the observed signal is not an artifact of our averaging procedure. It is also not caused by the transition from ERBS to CERES satellites since all six UPSCALE AMIP simulations show reductions, ranging from 0.30 to 0.52 W m−2. Furthermore, a similar decline in surface flux between 1999 and 2005 is simulated in forced ocean model experiments (Figure S2). However, the maximum decline in No could be subject to noise. We therefore compute the difference averaged over 2004–2006 minus 1998–2000, giving a reduction of −0.31 ± 0.21 W m−2 between the late 1990s and mid-2000s.

5. Summary and Conclusions

We find significant correlation (r ∼ 0.6, p < 0.01) between rolling 5 year mean ocean heat content trends (Ht) in two ocean analyses and net incoming radiation at the top of the atmosphere (N) from satellite observations and AMIP model simulations over the period 1960 to 1999. Ocean reanalyses of this historical period are therefore potentially useful for further studies of the physical processes of ocean heat uptake. However, ocean reanalyses of Ht in the early 2000s are inconsistent with observation- and model-based constraints on N. This helps to explain our inability to close Earth's energy budget [Trenberth and Fasullo, 2010; Trenberth et al., 2014] and suggests that observed estimates of Ht covering this period [e.g., Lyman et al., 2010, Levitus et al., 2012, Abraham et al., 2013; Rhein et al., 2013] should be treated with caution. It also calls into question whether ocean reanalyses can be used to robustly attribute increased Ht below the ocean mixed layer in the early 2000s as an explanation for the onset of the recent slowdown in surface warming [Guemas et al., 2013].

We create an observation-based estimate of N (No) extending back to 1960, using satellite data and AMIP model simulations, and compare with 21 CMIP5 coupled climate model simulations. Variations in N over this period are dominated by the volcanic eruptions of Agung in 1963, El Chichon in 1982, and Pinatubo in 1991 [Trenberth et al., 2014], with good agreement (r = 0.82) between No and the CMIP5 ensemble mean. Absolute values of No are in broad agreement with the models, but uncertainties are large and improved estimates of Ht are needed to provide a more accurate anchor for satellite observations [Loeb et al., 2012].

We find a reduction of 0.31 ± 0.21 W m−2 in No between the late 1990s and the mid-2000s which may have contributed to the recent slowdown in global surface warming. This is consistent with minor volcanoes [Solomon et al., 2011; Fyfe et al., 2013; Haywood et al., 2014; Santer et al., 2014], an extended and deeper solar minimum [Lean, 2009; Kaufmann et al., 2011], and possible nitrate and indirect aerosol effects [Shindell et al., 2013; Bellouin et al., 2011] that were not included in the CMIP5 models [Flato et al., 2013, Box 9.2; Schmidt et al., 2014; Huber and Knutti, 2014]. However, the reduction in No is caused by an increase in OLR rather than a reduction in ASR which would be expected from these factors (compare green and red dashed curves in Figures4 and S3), though ASR may have been reduced relative to the models. Internal variability may also reduce N during cooling decades [Katsman and van Oldenborgh, 2011; Palmer and McNeall, 2014; Brown et al., 2014], but the CMIP5 ensemble mean does simulate a weak reduction (−0.12 ± 0.06 W m−2, mainly in ASR), suggesting a potential role for external forcing. Further work is therefore needed to unravel the relative roles of internal variability and external factors on Earth's energy imbalance during the recent warming pause.

While present-day anomalies in N in the CMIP5 ensemble mean are in good agreement with observations, this is potentially achieved through a cancelation of errors (in models and/or observations) since both ASR and OLR are larger in almost all of the models than the satellite observations (Figures4 and S3). Both ASR and OLR are projected to increase further over the coming decades, with anomalies in ASR about twice as large as those in OLR on average (in agreement with CMIP3 models [Trenberth and Fasullo, 2009]). Further analysis of the roles of forcing and feedback on these signals is needed, along with continued satellite and in situ observations of Earth's energy imbalance to monitor future climate variability and change.