[1] The elusive nature of the post‐2004 upper ocean warming has exposed uncertainties in the ocean's role in the Earth's energy budget and transient climate sensitivity. Here we present the time evolution of the global ocean heat content for 1958 through 2009 from a new observation‐based reanalysis of the ocean. Volcanic eruptions and El Niño events are identified as sharp cooling events punctuating a long‐term ocean warming trend, while heating continues during the recent upper‐ocean‐warming hiatus, but the heat is absorbed in the deeper ocean. In the last decade, about 30% of the warming has occurred below 700 m, contributing significantly to an acceleration of the warming trend. The warming below 700 m remains even when the Argo observing system is withdrawn although the trends are reduced. Sensitivity experiments illustrate that surface wind variability is largely responsible for the changing ocean heat vertical distribution.

1 Introduction [2] Increasing greenhouse gases have apparently led to an increasing radiative imbalance at the top of the atmosphere (TOA) of order 0.5–1 W m–2 in the past two decades, based on observational [Trenberth, 2009; Trenberth et al., 2009; Murphy et al., 2009] and model [Hansen et al., 2011] estimates. Over the past 50 years, the oceans have absorbed about 90% of the total heat added to the climate system [Bindoff et al., 2007], while the rest goes to melting sea and land ice, and warming the land surface and atmosphere. Although the Mount Pinatubo eruption in 1991 caused a short‐term reduction in TOA radiation [Trenberth and Dai, 2007], increasing greenhouse gases should have led to increasing warming. However, sea surface temperature (SST) increases stalled in the 2000s and this is also reflected in upper ocean heat content (OHC) for the top 700 m in several analyses [Levitus et al., 2009, 2012; Lyman et al., 2010]. Although the energy imbalance from 1993 to 2003 could be accounted for, it was not possible to explain the energy imbalance from 2004–2008. This led to the concept of “missing energy” [Trenberth and Fasullo, 2010]. [3] Several recent modeling studies [Easterling and Wehner, 2009; Palmer et al., 2011; Katsman and van Oldenborgh, 2011; Meehl et al., 2011] advocate for the role of the deep ocean in the heat uptake. However, there are challenges in sampling the deeper reaches of the ocean although previous estimates [Purkey and Johnson, 2010; Kouketsu et al., 2011] found the deep ocean (below 1000 m) gaining heat at rates of less than 0.1 Wm−2 (95% confidence interval; global average). The most recent estimate of the ocean warming in 0–2000 m depth range of 0.39 W m–2 (per unit area of the World Ocean) from 1955 to 2010 emphasized the dominant role of the 700–2000 m depth range in the heat uptake (about one‐third of the total warming) [Levitus et al., 2012; L12 hereafter] . [4] A key question is how recent and robust is the role of deep ocean in the heat uptake, because the advent of the Argo ocean observing system is known to have a profound impact in ocean state estimations [Balmaseda et al., 2007]. Another question is how much disruption there is of the warming trends by natural variability and from the volcanic eruptions, which contributed to the TOA imbalance but have not been factored into most analyses of OHC variations. Several decadal hiatus periods of the upper ocean warming associated with natural decadal and interannual variability, such as La Niña events, have been found in a model [Meehl et al., 2011]. Recently, the 2001–2010 interannual variations of TOA radiation and OHC have been associated with El Niño‐Southern Oscillation (ENSO) [Loeb et al., 2012]. [5] Here we investigate the time evolution of the global OHC at different depth ranges for the period 1958–2009 using the latest European Centre for Medium‐Range Weather Forecasts ocean reanalysis system 4 (ORAS4) [Balmaseda et al., 2013, BMW13 in what follows], which provides a continuous record of the history of the global ocean by combining a wealth of observational information (section 2). Section 3 presents the time evolution of ORAS4 OHC at different depth ranges, the geographical distribution of the warming trends, and several sensitivity experiments, including the impact of Argo. The main findings are summarized in section 4.

2 The Ocean Reanalysis [6] ORAS4 has been produced by combining, every 10 days, the output of an ocean model forced by atmospheric reanalysis fluxes and quality controlled ocean observations. These consist of temperature and salinity (T/S) profiles from the Hadley Centre's EN3 data collection [Ingleby and Huddleston, 2007], which include expendable bathythermographs (T only, with depth corrections from Table 1 of Wijffels et al. [2008]), conductivity‐temperature‐depth sensors (T/S), TAO/TRITON/PIRATA/RAMA moorings (T/S), Argo profilers (T/S), and autonomous pinniped bathythermograph (or elephant seals, T/S). Altimeter‐derived along track sea level anomalies from AVISO are also assimilated. Gridded maps of SST from NOAA are used to adjust the heat fluxes via strong relaxation, and altimeter global mean sea‐levels are used to constrain the global average of the fresh‐water flux. The ocean model horizontal resolution is approximately 1°, refined meridionally down to 1/3° at the equator. There are 42 vertical levels with separations varying smoothly from 10 m at the surface to 300 m at the bottom, with partial cell topography. [7] A model bias correction [BMW13] is used to reduce potential spurious variability resulting from changes in the observing system. The bias correction first guess—a seasonal cycle of 3‐D model error—is estimated from the data‐rich Argo period, and applied to ORAS4 from the beginning of the record. This is updated as the analysis progresses via an adaptive scheme (see BMW13 for details; see also Figure S3 of the auxiliary material). The five ensemble members of ORAS4 sample plausible uncertainties in the wind forcing, observation coverage, and the deep ocean. The uncertainty is probably underestimated in ORAS4, because the uncertainty in observations and their quality control [Lyman et al., 2010] is not sampled. Quality improvements in ORAS4 relative to earlier ocean reanalyses stem from the use of improved atmospheric surface fluxes, improved data assimilation, and more comprehensive quality‐control of the observation data set, with important corrections to the ocean observations. [8] The methods section S01 in the auxiliary material provides more specific information on the model, surface forcing, observation data sets, bias correction and ensemble generation. A detailed description and evaluation of ORAS4 is given in BMW13, and a discussion of the sensitivity of the reanalysis to several aspects not included in the ensemble generation.

3 The ORAS4 Ocean Heat Content 3.1 Time Evolution and Trends [9] Figure 1 shows the time evolution of the global OHC at different depth ranges (upper 300 m, upper 700 m and total column) as represented by the five ensemble members of ORAS4. The 300 and 700 m depth ranges are chosen to be consistent with previous OHC reconstructions. In this context, depths below 700 m are sometimes referred as “deep ocean”, not to be confused with “deep waters” commonly used in oceanography. A large uncertainty (more than 5 × 1022 J ) in the first two decades in the total OHC arises from the sparse observations that do not constrain the values well. The values for the upper 300 and 700 m, where there are more observations, converge faster (within the first 5 years of the reanalysis). The evolution of the OHC is dominated by a clear warming trend starting around 1975, with two distinct cooling episodes in the mid‐1980s and in the mid‐1990s. These follow the volcanic eruptions of El Chichón in March–April 1982 and Pinatubo in June 1991, respectively [Trenberth and Dai, 2007]. The sharp cooling in 1963—attributed to the eruption of Mount Agung in the same year—is followed by a weak cooling trend until around 1975, although uncertainty in this period is large. Aside from the volcanic cooling episodes, there is an additional cooling episode following the huge 1997–1998 El Niño event after 1998, which mainly affects the upper 700 m. The event led to a global warming of the atmosphere and made 1998 the warmest year on record to that point as heat came out of the ocean, largely through evaporative cooling [Trenberth et al., 2002]. After 1998, there was a rapid exchange of heat between the regions above and below 700 m (Figure S1 in the auxiliary material). The heat exchange between the layers above and below 700 m during 1998 is consistent with a recent study based on Argo data for more recent events [Roemmich and Gilson, 2011]. Then after 1999 the warming starts again dramatically, this time also involving all depth ranges. This signals the beginning of the most sustained warming trend in this record of OHC. Indeed, recent warming rates of the waters below 700 m appear to be unprecedented. Figure 1 Open in figure viewer PowerPoint OHC integrated from 0 to 300 m (grey), 700 m (blue), and total depth (violet) from ORAS4, as represented by its 5 ensemble members. The time series show monthly anomalies smoothed with a 12 month running mean, with respect to the 1958–1965 base period. Hatching extends over the range of the ensemble members and hence the spread gives a measure of the uncertainty as represented by ORAS4 (which does not cover all sources of uncertainty). The vertical colored bars indicate a 2 year interval following the volcanic eruptions with a 6 month lead (owing to the 12 month running mean), and the 1997–1998 El Niño event again with 6 months on either side. On lower right, the linear slope for a set of global heating rates (W m–2) is given. [10] The overall magnitude of the ocean warming in ORAS4 in the upper 700 m is comparable with other observational estimates [BMW13 and Figure S2 (left panel)]. The total depth warming is also comparable with that in the global OHC from L12 (Figure S2, right panel). Compared with the 5 year L12 product, ORAS4 OHC shows more distinctive interannual variability, making it easier to identify the impact of volcanoes and ENSO. Other differences arise from the analysis procedures and data sets used, as discussed in the methods section S01. Table 1 provides a summary of OHC trends at different depth ranges and for different periods, in W m–2 per unit area of the World Ocean. The trends and corresponding standard error are estimated using a linear trend fitting of the ORAS4 ensemble mean. The 1975–2009 trend is 0.39 ± 0.03 W m–2 and 0.47 ± 0.03 W m–2 for the upper 700 m and total column, respectively. For 1960–2009, the corresponding values are 0.22 and 0.29 W m–2, slightly lower than those reported by L12 for 1955–2010 (0. 27 and 0.39 W m–2 for the upper 700 and 2000 m). As in L12, the contribution to warming of depths below 700 m is stronger after the 1990s. Table 1. OHC Linear Trends (W m–2) 1975–2009 1980s 1990s 2000s NoArgo 2000s Total 0.47 ± 0.03 0.58 ± 0.15 −0.26 ± 0.13 1.19 ± 0.11 0.82 ± 0.10 Upper 300 m 0.26 ± 0.02 0.46 ± 0.11 −0.04 ± 0.12 0.45 ± 0.08 0.35 ± 0.08 Upper 700 m 0.38 ± 0.03 0.59 ± 0.13 −0.22 ± 0.10 0.90 ± 0.13 0.67 ± 0.09 Below 700 m 0.10 ± 0.01 −0.01 ± 0.02 −0.04 ± 0.05 0.30 ± 0.02 0.15 ± 0.01 3.2 Sensitivities and Diagnostics [11] The deployment of the Argo floats from 2000–2004 is a revolution in the ocean observing capabilities and it is only after 2003 that regular and spatially homogenous temperature soundings of the upper 2000 m are available, giving new confidence in the OHC assessment in the recent years, but rendering comparison of the OHC in the pre‐Argo and post‐Argo difficult. To address the impact of Argo in the ORAS4 OHC estimate, an additional experiment (NoArgo), equivalent to ORAS4 but withdrawing Argo data, was carried out for 2001–2009. NoArgo will also exclude the elephant seal data (treated as Argo in ORAS4), but will include the changes experienced in expendable bathythermographs observations. Withdrawing Argo influences the total OHC, but the gap of about 5 × 1022 J between the OHC in the top 700 m and total column remains (Figure 2, top panel). Figure 2 Open in figure viewer PowerPoint (top) OHC evolution of the NoArgo experiment (thick single lines) compared with the OHC in ORAS4 (multiple lines hatched). The layers are 0–300 m, 0–700 m and total column. (bottom) Total (thick lines) and upper 300 m (thin lines) OHC anomalies (respect to the 1990–1993 period) in ocean only simulations: the control experiment (CNTL black) is forced by daily heat and momentum fluxes from atmospheric reanalyses; ClimWind (blue) is equivalent to CNTL but using a daily wind‐stress climatology. (Atmospheric reanalyses include a climate change signal, both in atmospheric composition and observed state variables, and in the SST forcing). [12] In addition to the NoArgo experiments, other experiments have been carried out to establish the sensitivity of results to different aspects. Experiments separately withdrawing the altimeter and mooring data did not show any significant variation from the heat budget above. An experiment similar to ORAS4 but without bias correction shows a similar OHC evolution, but with larger decadal variability and, in particular, a much stronger contribution from the depths below 700 m (Figure S3). Ocean‐only simulations forced by atmospheric reanalysis fluxes and constrained by SST, without any data assimilation, exhibit stronger warming trends (and stronger than observational estimates), probably because the ocean model diffuses too much heat vertically, but it shows similar cooling episodes to ORAS4, except for the El Chichón event. The post‐1999 hiatus in the upper 300 m OHC (Figure 1) and the two volcanic cooling episodes are also present in these ocean‐only simulations (albeit with smaller magnitude), suggesting that the signal is not a spurious artifact of the ocean observing system. Another ocean‐only simulation where the time‐varying forcing was replaced by climatological daily winds shows a much weaker heat uptake by the ocean, since most of the warming is confined within the upper 300 m (bottom panel in Figure 2), although both experiments exhibit an accelerated warming around 2000–2006. In the experiment where the wind interannual variability is suppressed, the heat is mostly confined to the upper 300 m, the ENSO signal disappears, and the warming stalls (in fact it reverses) after 2006. This result suggests that changes in the atmospheric circulation are instrumental for the penetration of the warming into the ocean, although the mechanisms at work are still to be established. One possibility suggested by Lee and McPhaden [2008] is related to the modified subduction pathways in response to changes in the subtropical gyres resulting from changes of the trade winds in the tropics (Figure S4), but whether as low frequency variability or as a longer term trend remains an open question. The 2000–2006 warming trend is likely associated with the weakening of the Atlantic Meridional Overturning Circulation (MOC) in both experiments (see BMW13). [13] The role of the different terms contributing to variations of OHC in ORAS4 (surface fluxes, assimilation increment, SST relaxation and bias correction) has been quantified. The surface fluxes from the reanalyses show large spurious changes over time especially with the switch from ERA‐40 to ERA‐I in 1989, but these are compensated for by the bias correction and the relaxation to observed SSTs (see auxiliary Figure S5 and section S02). The analysis shows that information on the El Chichón cooling comes mainly from subsurface observations, while the Mt Pinatubo and Mt Agung coolings are largely provided by the SST relaxation. The El Niño signature is present in the ERA‐I fluxes, and this explains why it is not so apparent in other subsurface observation‐based ocean analyses. The subsequent warming is both in the ERA‐I fluxes and subsurface observations (Figure S5). 3.3 Geographical Distribution [14] The OHC has also been examined to find where the main variations occur (Figure 3) and it is clear that while the long term trend is present in all the regions, the dominant contributions to the low frequency variability arise from the Tropics (±30°) (especially the Pacific Ocean). Note especially the extreme amount of cooling of the tropical Pacific Ocean following the 1997–1998 El Niño event that released heat from the ocean [Trenberth et al., 2002] to such an extent that a major recovery period was required. There is also a net poleward heat transport during the discharge phase of ENSO as can be seen by the exchange of heat between tropics and extratropics, which is likely favored by the intensification of the trades after 1998 (Figure S4). The drop in the OHC around 2007, probably associated with an El Niño, was also noted by Roemmich and Gilson [2011] using Argo‐only data. The spatial patterns of OHC trends in ORAS4 (Figure S6) differ over time. Trends above 700 m clearly show the influence of ENSO and PDO, in the Pacific and Indian oceans, and consistent warming in the Atlantic Ocean. There is a clear change after the year 2000 in the structure of the upper 700 m Pacific OHC trends. There are also changes in the trend patterns of the OHC below 700 m. The tropical Pacific and Indian oceans show intensified positive trends below 700 m for the last decade. The trend patterns below 700 m in the North Atlantic are suggestive of changes in the ocean circulation. After the year 2000, there is a pronounced warming along the path of waters associated with the Mediterranean outflow. This warming is reduced in the NoArgo experiment (not shown). Figure 3 Open in figure viewer PowerPoint Contributions to the total‐column OHC anomalies by regions. Northern (NXTRP) and Southern extratropics (SXTRP), and Tropics (TROP) where the divide occurs at ±30°. For the tropics, the regions are further subdivided into the Atlantic (TRATL), Pacific (TRPAC) and Indian (TRIND) ocean contributions to the total. Units 1022 Joules. For clarity only one ensemble member is shown, but the dominant features (interannual variability and trend) are robust across ensemble members.

4 Summary and Conclusions [15] The time evolution of the global OHC for the period 1958–2009, as estimated by the ORAS4 ocean reanalysis, is dominated by a warming trend and pronounced cooling episodes, and shows an increasing warming trend at depths below 700 m. The cooling episodes correspond to cooling seen in SSTs in response to the El Chichón and Mt Pinatubo eruptions, and the radiative imbalance associated with the latter [Trenberth and Dai, 2007] is consistent with the cooling found here. More surprising is the extra cooling following 1998, a likely consequence of the ocean heat discharge associated with the massive 1997–1998 El Niño event [Trenberth et al., 2002]. Meehl et al. [2011] have demonstrated in a model study how La Niña events and negative PDO events could cause a hiatus in warming of the top 300 m while sequestering heat at deeper layers. This mechanism can also explain the increasing role of the depths below 700 m after 1999 in the ORAS4 OHC, consistent with La Niña‐like conditions and a negative phase of the PDO which has dominated the last decade. The deep ocean warming, which mostly involves the depth range 700–2000 m, may also be related to the weakening of the MOC after 1995, which is present in ORAS4 [BMW13]. Possibly changes in MOC and PDO are connected through changes in the atmospheric circulation patterns. [16] The deep ocean has continued to warm, while the upper 300 m OHC appears to have stabilized. The differences in recent trends among the different ocean layers are profound. The small warming in the upper 300 m is belied by the continuing warming for the ocean as a whole, with considerable warming occurring below 700 m. However, this raises the question of whether this result is simply because of the new Argo observing system? The results shown here suggest otherwise, although Argo clearly is vitally important quantitatively. Instead changes in surface winds play a major role, and although the exact nature of the wind influence still needs to be understood, the changes are consistent with the intensification of the trades in subtropical gyres. Another supporting factor is the uniqueness of the radiative forcing associated with global warming. [17] The magnitude of the warming trend is consistent with observational estimates, being equivalent to an average 0.47 ± 0.03 W m–2 for the period 1975–2009. There is large decadal variability in the heat uptake, the latest decade being significantly higher (1.19 ± 0.11 W m–2) than the preceding record. Globally this corresponds to 0.84 W m–2, consistent with earlier estimates [Trenberth et al., 2009]. In an observing system experiment where Argo is withdrawn, the ocean heating for the last decade is reduced (0.82 ± 0.10 W m–2), but is still significantly higher than in previous decades. The estimation shows depths below 700 m becoming much more strongly involved in the heat uptake after 1998, and subsequently accounting for about 30% of the ocean warming. [18] The analysis of ORAS4 OHC shows some interesting signals. In particular, the prolonged and intense cooling events during the 1980s and 1990s are not as distinct in other observation‐only analyses [BMW13], and the rapid involvement of the deep ocean starting around the 1998–1999 La Niña needs further investigation. Sensitivity experiments indicate that these features are robust, and suggest that changes in the atmospheric circulation play an important role in the heat uptake. Detecting, understanding and modeling the processes that lead to the vertical distribution of heat within the ocean is a key for the correct initialization of decadal predictions, because the trends in forecasts of the SST will likely depend on whether the ocean is in a recharge (low stratification) or discharge (high stratification) mode.

Acknowledgments [19] The work of Trenberth was supported by NASA award NNX09AH89G. The National Center for Atmospheric Research is supported by the National Science Foundation. The production of ORAS4 has been partially funded by the EU FP7 COMBINE project.

Supporting Information Filename Description grl50382-sup-0001-2013GL055587_AuxiliaryMaterial_14March2013.pdfPDF document, 1.3 MB Supporting information grl50382-sup-0002-55587Rreadme.docxPDF document, 27.6 KB Supporting information Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.