Improved estimates of ocean heat content from 1960 to 2015

See allHide authors and affiliations

Science Advances  10 Mar 2017:
Vol. 3, no. 3, e1601545
DOI: 10.1126/sciadv.1601545


Earth’s energy imbalance (EEI) drives the ongoing global warming and can best be assessed across the historical record (that is, since 1960) from ocean heat content (OHC) changes. An accurate assessment of OHC is a challenge, mainly because of insufficient and irregular data coverage. We provide updated OHC estimates with the goal of minimizing associated sampling error. We performed a subsample test, in which subsets of data during the data-rich Argo era are colocated with locations of earlier ocean observations, to quantify this error. Our results provide a new OHC estimate with an unbiased mean sampling error and with variability on decadal and multidecadal time scales (signal) that can be reliably distinguished from sampling error (noise) with signal-to-noise ratios higher than 3. The inferred integrated EEI is greater than that reported in previous assessments and is consistent with a reconstruction of the radiative imbalance at the top of atmosphere starting in 1985. We found that changes in OHC are relatively small before about 1980; since then, OHC has increased fairly steadily and, since 1990, has increasingly involved deeper layers of the ocean. In addition, OHC changes in six major oceans are reliable on decadal time scales. All ocean basins examined have experienced significant warming since 1998, with the greatest warming in the southern oceans, the tropical/subtropical Pacific Ocean, and the tropical/subtropical Atlantic Ocean. This new look at OHC and EEI changes over time provides greater confidence than previously possible, and the data sets produced are a valuable resource for further study.

  • Earth’s energy imbalance
  • ocean heat content
  • sampling
  • mapping
  • subsample test
  • global warming
  • decadal variability
  • hiatus


Global warming is driven by Earth’s energy imbalance (EEI). The EEI is likely forced to first order by a combination of greenhouse gas and aerosol forcing, which shapes the timing and magnitude of global warming (1). It is also linked to the internal variations of the climate system and episodic volcanic eruptions; the latter may provide episodic strong radiative forcing to the Earth system (2, 3). By definition, radiative forcing is the change in the net radiative flux due to a change in an external driver of climate change, such as greenhouse gas concentrations. More than 90% of EEI is stored in the ocean, increasing ocean heat content (OHC), while the residual heat is manifest in melting of both land and sea ice, and in warming of the atmosphere and land surface (1, 4, 5). It is therefore essential to provide estimates of OHC changes over time with high confidence to improve our knowledge of EEI and its variability (4).

How much has Earth really warmed in recent decades? The magnitude and location of the ocean warming have become an area of active research, because of the large historical uncertainty in estimated OHC changes (3, 610). For instance, tracking Earth’s heat and ocean heat is one of the key topics of the so-called “global warming hiatus” research surge (1120). There are several sources of uncertainty in OHC estimates. A detailed discussion of uncertainty can be found in section S1, supplementing the brief overview here. A key source of uncertainty is instrumental in nature, relating to the measurement process; another originates from methodological choices, for example, in dealing with gaps in sampling (6, 8, 2127). The community has made progress in detecting the systematic errors in expendable bathythermograph (XBT) data and has provided recommendations to correct the associated errors (25). These recommendations have markedly reduced the impact of XBT biases on multidecadal OHC estimation (26). Another major uncertainty arises from insufficient data coverage, mainly during the pre-Argo era (before 2005), that has led to spatial sampling errors in global and regional OHC estimation (2830).

There have been many efforts to improve estimates of historical OHC changes including reduction of sampling errors using various methodologies (3, 6, 8, 28, 29, 3137). One group is referred to as mapping methods, which statistically infill data gaps by using available observations, a prior guess, and an adjustment term based on nearby observations and covariance characteristics, because grid points within the ocean are often mutually correlated (38). For example, a widely referenced method used by the National Centers for Environmental Information (NCEI) (34, 35, 39) objectively analyzed gridded temperature anomalies on the basis of the distances between analysis grid boxes and nearby observations, with zero anomalies used as a first guess. Domingues et al. (36) presented the main ocean spatial variability by means of an empirical orthogonal function analysis. Recently, Cheng and Zhu (40) used an ensemble optimal interpolation method (CZ16) combined with covariances from Coupled Model Intercomparison Project phase 5 (CMIP5) multimodel simulations to provide an improved prior guess (section S2). The other group of gap-filling techniques is referred to as dynamic gap filling, which constrains an ocean/climate model by observations via a data assimilation method, for instance, Ocean Reanalysis System 4 (ORAS4) (3). In this case, the prior guess is produced by bringing forward past information using a model. The accuracy of this method can be limited by model error, which varies from model to model (9), although this limitation may be partly mitigated by model bias corrections (for example, based on assessed bias during the Argo era). This study will focus on the first group—that is, mapping methods—because these have been the main method used for historical OHC estimation.

The estimates obtained on the basis of all of these methods can reveal substantial differences in OHC changes on monthly, interannual, and interdecadal time scales (9, 10, 21, 41), although all of them still indicate a robust ocean warming in recent decades. Therefore, designing the best way to infill the data gaps in both space and time while minimizing sampling error is a key task in global OHC assessment that demands a comprehensive analysis. Quantification of the reliability of obtained OHC estimates across temporal and spatial scales is also essential.

The success of a mapping method can be judged by how accurately it reconstructs the full ocean temperature domain. When the global ocean is divided into a monthly 1°-by-1° grid, the monthly data coverage is <10% before 1960, <20% from 1960 to 2003, and <30% from 2004 to 2015 (see Materials and Methods for data information and Fig. 1). Coverage is still <30% during the Argo period for a 1°-by-1° grid because the original design specification of the Argo network was to achieve 3°-by-3° near-global coverage (42). However, we typically use temporal persistence from 1 month to the next as well as spatial covariance (as carried out using all the mapping methods cited above). Because the early observations were mainly carried out on commercial and scientific research vessels, their locations are limited to the regions around developed countries and along shipping routes (6, 28, 43). From 1990 to 1997, the World Ocean Circulation Experiment (44, 45) extended the in situ data to a mostly one-time global network. Since late 1992, altimetry from space has provided essential complementary information on sea surface height. This century (mainly since 2005), the deployment of the Argo floats into a global ocean network markedly increased data coverage (46, 47), with the exception of coastal, marginal sea, and ice-covered regions (4, 48). Now, the Argo community is extending the Argo float deployment under sea ice and into marginal seas (49).

Fig. 1 Fractional coverage of the mapping method used in this study.

(A) Averaged fraction (0 to 700 m) of the global ocean considered sampled for temperature data in each month for this study. (B) Same as (A) but for 700 to 2000 m. The mean fractional coverage of observations is shown in blue if the global ocean is divided into 1°-by-1° grids, and fractional coverage for NCEI yearly and the 5-year method is shown in green and orange, respectively. This figure starts from 1940. (C) Fractional coverage of the global ocean for layers within 0 to 700 m. (D) Same as (C) but for 700 to 2000 m.

In this paper, we extend and improve a recently proposed mapping strategy (CZ16) to provide a complete gridded temperature field for 0- to 2000-m depths from 1960 to 2015. We carefully evaluate the estimate by comparing the fully analyzed temperature signals with those based on a “subsample test,” in which subsamples of data in the Argo era are colocated with historical ocean observations and used to assess sampling error. The time scales for which reliable estimates are obtained are quantified. Using this improved reconstruction, we derived an updated historical (1960–2015) ocean energy budget and contributions to Earth’s total energy budget.


The efficiency of a mapping method

We first address the fractional coverage of monthly ocean temperature data achieved with reconstructed signals, which is defined as the fraction of total ocean area obtained by the mapping method (50). The large fractional coverage helps to ensure that a near-complete global reconstruction can be reached. A mapping method uses the data only near the analyzed grid within an assumed area to perform the reconstruction at individual grid cells, with the size of the area defined by the influencing radius or spatial correlation length scales. The choice of influencing radius mainly governs the extent of fractional coverage.

When the present method (CZ16; see Materials and Methods and section S2) is applied, the analyzed fractional coverage is larger than 90% from the late 1950s to 2015 from the sea surface to 2000 m (Fig. 1). For comparison, the NCEI method has a fractional coverage of 60 to 90% for its pentadal estimate (compositing 5 years of data) before 1990 and 90% in the upper ocean after 1990 (0 to 700 m). Below 700 m, the coverage for the NCEI method increases from 50 to 80% from 1960 to 1990 for its pentadal estimate and sharply increases from 30 to 90% during the 2000–2005 period for its yearly estimate. The present method has larger fractional coverage than NCEI and most of the other published methods (21) mainly owing to the difference in the assumed influencing radius. For example, the NCEI method (39) uses an 888-km (~8°) radius, and many other methods adopt a similar radius (3133, 37). These methods take account of the spectrum of ocean spatial variability and are effective owing to the large spatial correlations within 8° (>0.6; Fig. 2). Instead, CZ16 uses 20° for an influencing radius within 0 to 700 m (25° for 700 to 2000 m) by using the relatively weak correlations (0.0 to 0.6; Fig. 2) for longer spatial distances (see section S3.1 for the choice of influencing radius in the CZ16 method). The application of a large influencing radius helps to ensure a high fractional coverage and also filters out the smaller-scale signals (that is, smaller than 20°), resulting in a smooth spatial pattern (an example is shown in fig. S2). However, the ocean shows spatial variability at multiple length scales (33, 38). To improve methods for assessing this variability, we adopted an iterative strategy similar to that adopted by NCEI (39). Three iterative scans were performed successively using 20° (25° below 700 m), 8°, and 4° to effectively include ocean variability across those spatial scales. Quantitative assessment of the iterative strategy (section S3) shows a significant reduction of analysis error. However, because this fractional coverage does not give information regarding how well a mapping method can reconstruct temperature signals, quantificative assessment of the accuracy of the reconstruction is also needed.

Fig. 2 Observational correlations with distance.

(A and C) Zonal-mean (Z) and (B and D) meridional-mean (M) correlation as a function of distance calculated using the ORAS4 reanalysis data. A linear trend and the mean seasonal cycle have been removed when calculating the correlation. The ORAS4 1° by 1° and monthly data were used.

Quantification of the global mean sampling error

The first step in quantifying the global mean sampling error is to investigate its contribution to both global temperature and OHC changes. The sampling error is assessed on the basis of a subsample test, which is defined as the difference between the reconstructed fields and “truth fields.” The truth field is taken to be a set of the gridded averaged temperature anomalies during the Argo era, where the anomalies are relative to 1997–2005 climatology (see Materials and Methods). Each truth field is subsampled according to the locations of historical observations and mapped to obtain the reconstructed fields. The sampling error is estimated from the difference between the reconstructed and truth fields (see Materials and Methods).

Figure 3A presents the global averaged sampling errors since the late 1950s for different choices of truth fields, with the uncertainty estimated as 2 SDs (2σ). The global mean error due to the historical sampling is around 0°C from the late 1950s to 2014, which is significantly different from the temperature signals of the truth field (Fig. 3A, stars). This indicates that no significant sampling bias exists despite the marked evolution of the observation system over the past 56 years.

Fig. 3 Global and basin-averaged sampling error compared with reconstructed temperature change.

(A to G) The sampling errors for different truth fields subsampled by historical observation locations are shown as green dots, with the green solid line and the error bar for the mean and ±2σ, respectively. Blue stars denote the 16 different truth fields. The gray line is the reconstructed monthly temperature anomaly time series from 1960 to 2015 with ±2σ interval in gray shading, and the dark line is the time series after a 7-year low-pass filter (see Materials and Methods). The orange curve is the NCEI result, along with a 1-SE bar (dashed orange). (A) Global, (B) tropical/subtropical Pacific, (C) North Pacific, (D) Indian, (E) tropical/subtropical Atlantic, (F) North Atlantic, and (G) southern oceans. (H and I) S/N ratio of the temporal variability of our reconstruction on (H) decadal scales (solid lines) and (I) interannual scales (triangle lines) compared to the sampling error.

Comparison between the sampling errors and ocean temporal variability on different time scales provides an indication of whether the variability from sampling error is detectable. The temperature time series (Fig. 3A, gray lines) for decadal and multidecadal scales exhibit variability that is much larger than the uncertainties induced by sampling error (Fig. 3A, green dots and 2σ error bars). The decadal/multidecadal variability, which is defined as 2σ of the time series after applying a 7-year low-pass filter (see Materials and Methods), is 0.076°C, at least two times larger than the 2σ sampling error, with a signal-to-noise (S/N) ratio larger than 3 (Fig. 3H). This indicates a robust reconstruction of ocean decadal/multidecadal variability.

The temporal variability of ocean temperature change on interannual scales is 0.008°C after 1990, comparable with the sampling error (S/N ratio between 0.5:1 and 3:1; Fig. 3I). Errors in the signals in both the truth and reconstructed fields are unavoidable because of the imperfect removal of sampling error, imperfect quality control processes, the impact of mesoscale signals, and the residuals of instrumental biases. The year 1990 is chosen as a boundary here to better represent interannual variability, because associated errors in the reconstructed field are reduced (but not removed) in the most recent periods owing to improvements in data quality and density. On the other hand, the estimated sampling error also reveals uncertainty from other error sources in recent decades. Therefore, the comparable S/N ratio between sampling error and interannual variability implies that the reconstructions of year-to-year global temperature changes remain uncertain and can easily be affected by observational error. A recent study (10) based on a comparison of month-to-month variations of global OHC for the upper 2000 m with radiative observations at the top of atmosphere (TOA) noted that short-term fluctuations in OHC contain large noise.

An analysis of these temperature reconstructions across different layers (section S3.2) addresses uncertainty in interannual variability estimates for all ocean layers. In the upper ocean (that is, 20 m in fig. S7), the interannual variability is significant and larger than the sampling error (S/N ratio >2) because the upper ocean experiences strong interannual variability mostly dominated by El Niño–Southern Oscillation (ENSO) (3). However, the S/N ratio decreases markedly below 300 m (shown in figs. S8 to S10), and interannual variability is of comparable magnitude to sampling error in the deep ocean, simply because there is much weaker variability (signal) than in the upper 300 m.

Furthermore, the S/N ratio of both decadal and interannual variability increases with time, starting in the 1960s, and the interannual variability becomes slightly stronger than the sampling error during the Argo era (S/N ratio between 2:1 and 5:1), suggesting the need to maintain and extend the integrated ocean observation system. The latter includes the current Argo network in the open ocean, a ship-based observation system, glider and coastal moorings, the coastal regions, the moored array in the tropics, instrumented pinnipeds in polar regions, and modified Argo floats and ice-tethered profilers in ice-covered regions.

Quantification of the regional sampling error

It is also worthwhile examining the sampling error of basin-scale averages because different mechanisms may be responsible for OHC changes in different ocean basins. Moreover, ocean heat transport within and between the basins is of great interest in understanding Earth’s key energy flows. Here, we divide the oceans into the North Pacific Ocean (30° to 65°N), the North Atlantic Ocean (30° to 65°N), the tropical and subtropical Pacific Ocean (30°S to 30°N), the tropical and subtropical Atlantic Ocean (30°S to 30°N), the Indian Ocean (30°S to 20°N), and the southern oceans (south of 30°S). There are continual observations for the North Atlantic and North Pacific oceans in middle latitudes from mechanical bathythermographs (1940s to 1960s) and XBTs (late 1960s to 2000), mainly distributed along cruise ship lines. In the tropics, the Tropical Atmosphere Ocean moorings started in 1979, although they were only fully deployed by about 1992. The least sampled area is south of 30°S, especially in the winter half-year. Even in the most recent decade, there are minimal data in the ice-covered regions.

In the six ocean basins, the mean sampling error is approximately zero (Fig. 3), indicating a robust reconstruction of temperature signals. Five ocean basins (not including the Indian Ocean) have experienced significant decadal-scale temperature changes, which are all larger than the 2σ sampling error (S/N ratio is larger than 2), confirming that long-term basin-averaged temperature changes are robust and insignificantly affected by sampling error. In the Indian Ocean (where there is evidently low signal), decadal variability is not significant before 1970, with an S/N ratio between 1.6:1 and 2:1. In particular, the southern oceans and the Atlantic Ocean have experienced the most marked decadal variations, with larger S/N ratios than the other basins. The magnitude of interannual variability in all six basins is comparable to the 2σ sampling error (S/N ratio between 0.3:1 and 5:1; Fig. 3I), as was found for the global ocean.

In addition, a geographical analysis of sampling error is presented in section S3.3. The sampling error for each 1°-by-1° grid has mean errors centered mostly about zero (fig. S11), suggesting that no significant regional bias in the reconstructed field exists over most ocean regions (see also fig. S12 for the 20- and 1600-m layers). However, in boundary current systems, Antarctic Circumpolar Current regions, and the Southern Hemisphere midlatitudes, there are larger 2σ sampling errors than in other regions (fig. S11B).

Global OHC change for the upper 2000 m

On the basis of reconstructed temperature fields and associated error bars, monthly OHC changes within 0 to 700 and 700 to 2000 m (Fig. 4A) show significant warming in the past 56 years. A stronger ocean warming trend has existed since the late 1980s for both 0- to 700- and 700- to 2000-m depths compared with the 1960s to the 1980s. The linear trend of OHC at 0- to 700-m depth is 0.15 ± 0.08 × 1022 J/year (0.09 ± 0.05 W/m2) during 1960–1991 and 0.61 ± 0.04 × 1022 J/year (0.38 ± 0.03 W/m2) during 1992–2015, a warming trend four times stronger than the 1960–1991 period. The linear trend of OHC at 700- to 2000-m depth is 0.04 ± 0.08 × 1022 J/year (0.02 ± 0.05 W/m2) in the period 1960–1991 and 0.37 ± 0.02 × 1022 J/year (0.23 ± 0.02 W/m2) from 1992 to 2015 (nine times stronger than that from 1960–1991) (table S1). This indicates an accelerating heat input into both the 0- to 700-m and 700- to 2000-m layers. The acceleration is most probably linked to the increasing EEI with time. In particular, a new study (51) shows that the recent acceleration is partly due to recovery from the volcanic eruption in 1992 (Mt. Pinatubo), which led to strong ocean cooling.

Fig. 4 Global OHC change time series.

(A) OHC from 0 to 700 m (blue), 700 to 2000 m (red), and 0 to 2000 m (dark gray) from 1955 to 2015 as obtained by this study, with the uncertainty of the ±2σ interval shown in shading. All time series of the new analysis are smoothed by a 12-month running mean filter, relative to the 1997–2005 base period. (B) The new estimate is compared with an independent estimate from NCEI with its SE as dashed lines. Both OHC 0 to 700 m and OHC 700 to 2000 m are shown from 1957 to 2014. The baseline of the time series from NCEI is adjusted to the values of the current analysis within 2005–2014.

The new estimates for OHC 700 to 2000 m show a stronger long-term trend before 2005 than NCEI estimates, with 10 to 20% differences, although the two estimates are not distinguishable within the margin of error (Fig. 4B). For both the upper 700-m and 700- to 2000-m layers, our current OHC estimates show larger long-term changes in the period 2005–2010 than NCEI estimates. To see where the “extra heat” is coming from, we include the NCEI estimates in Fig. 3 (A to G, orange lines) for comparison. The major differences appear in the southern oceans and the tropical/subtropical Pacific Ocean: The current OHC estimate shows quicker warming since 1960 than do NCEI estimates, albeit with sparse data coverage before the Argo era. This is probably due to the prior assumption of zero anomalies used in NCEI, which biases the reconstructed field toward zero in data-sparse regions. Rapid warming in the southern oceans has been observed in all Argo-based data sets since 2004 (10) and in the altimetry sea level (ASL) change since 1993 (fig. S18).

Observational estimates (both NCEI and this study) show a very weak (insignificant) deep (700 to 2000 m) ocean warming before 1980, in apparent contradiction to the strong radiative forcing since the 1960s and many model simulations (1). There are several hypotheses for this moderate deep ocean warming that require a more careful analysis in the future: (i) The warming signal is small during the 1960s and can be influenced by sampling error (note that the S/N ratios during the 1960s are small on global and regional scales in Fig. 3); (ii) there is a robust ocean warming in the upper 700 m, but the warming does not penetrate to the deep ocean below 700 m. By contrast, climate models may be too diffusive and thus may prematurely warm the deep ocean (5153).

Regional OHC change in recent decades

Regional OHC changes are of great interest to the broad climate community because they provide a basis for understanding the mechanisms of ocean energy flow, especially since 1998 (the so-called “hiatus” period), as previous studies have sometimes shown contradictory results based on different data sets (14, 16, 17, 54, 55). Our studies show that there has been no slowdown in global OHC change since 1998 compared with the previous decade (Fig. 5): There is an acceleration of ocean warming (for both OHC 0 to 700 m and OHC 0 to 2000 m), consistent with a previous study based on model-based data assimilation (3). Moreover, the percentage of 700- to 2000-m ocean heating relative to 0- to 2000-m OHC is 32% (32%) for 1960–1998 (1960–2005), but 38% (49%) for 1998–2015 (2005–2015), indicating that the deep ocean has played an increasingly important role in the ocean energy budget since 1998. The total OHC increase from 1998 to 2015 is 15.2 × 1022 J in the upper 2000 m, with 17% stored in the Pacific Ocean, 24% in the Indian Ocean (30°S northward), 31% in the Atlantic Ocean, and 28% in the southern oceans (south of 30°S). Again, the total OHC change calculated here is not well characterized by a linear trend because of the relatively short time period considered and the presence of strong decadal variability. It is evident that all six ocean basins have experienced significant warming since 1998 but that heat was mainly stored in the southern oceans, the tropical/subtropical Pacific Ocean, and the tropical/subtropical Atlantic Ocean from 1960 to 1998 (Fig. 5). Understanding how this heat has been transported or redistributed in the ocean continues to be an important research topic.

Fig. 5 OHC changes from 1960 to 2015 for different ocean basins.

(A) For 0 to 2000 m, (B) 0 to 700 m, and (C) 700 to 2000 m. All the time series are relative to the 1997–1999 base period and smoothed by a 12-month running filter. The curves are additive, and the OHC changes in different ocean basins are shaded in different colors.

Atlantic Ocean warming shows the largest OHC 0 to 2000 m increase (about 3.5 times larger than the Pacific Ocean), despite having an area only 47% as large as that of the Pacific (Figs. 3 and 5). The OHC change in the North Atlantic Ocean shows strong decadal variability and is likely linked to the strengthening of Atlantic meridional overturning circulation up to the middle 1990s and a subsequent weakening during the 2000s (Fig. 3) (56, 57), although the reasons for these changes are still debated. Other regions with significant decadal variability are the North Pacific Ocean and the tropical/subtropical Pacific Ocean. However, these changes are mainly limited to the upper 700 m (with very weak changes below 700 m; Fig. 5C and figs. S7 to S10), showing the importance of the shallow subtropical circulation (11, 58). The southern oceans and the tropical/subtropical Atlantic Ocean have experienced continuous and monotonic long-term warming since the 1960s, revealing a robust footprint of global warming (59). The Atlantic Ocean and the southern oceans are the major new heat reservoirs (59%) even though their total area is just 48% that of the global ocean. This implies that the meridional overturning circulation is a key mechanism involved in the sequestering of heat in the deeper (below 700 m) ocean.

For the upper ocean (0 to 700 m), the OHC increase from 1998 to 2015 is 10.11 × 1022 J, with 21% occurring in the Pacific Ocean, 21% in the Indian Ocean, 28% in the Atlantic Ocean, and 30% in the southern oceans, generally consistent with the OHC 0 to 2000 m. Our data reveal significant heat input into each basin in the upper ocean in the recent decade (Figs. 3 and 5), including the Indian Ocean (10).

Implications for the planetary energy imbalance

This updated OHC estimate offers the opportunity to revisit the planetary energy imbalance over the historical period, since more than 90% of EEI is stored in the ocean. The revised ocean energy budget (Fig. 6) based on OHC in the upper 0 to 2000 m from this study is combined with OHC below 2000 m adapted from Purkey and Johnson (60) (0.11 ± 0.10 × 1022 J/year), which has been updated by Desbruyères et al. (61). It is extended to 1960 by assuming a zero heating rate before 1990 below 2000 m, because OHC changes are found to be very small for 700 to 2000 m before 1990 (with large uncertainty), consistent with the assumption adopted in previous sea level analyses (62). The OHC 700 to 2000 m in this study (NCEI) indicates a one-ninth (one-third) smaller warming rate for the 1960–1991 versus the 1992–2015 period. Therefore, the maximum temperature change below 2000 m before 1991 should be one-third of its change after 1991 (that is, 0.04 × 1022 J/year). To be conservative, the uncertainty of the OHC change below 2000 m is set to ±0.04 × 1022 J/year.

Fig. 6 Estimate of the ocean energy budget.

The three major volcanic eruptions are marked. The energy budgets are relative to the 1958–1962 base period. The integrated net radiative imbalance from Allan et al. (65) estimated from the TOA is included in yellow and is multiplied by 0.93 to be comparable with the ocean energy budget. The TOA radiation is adjusted to the value of OHC within 2013–2014. The dashed gray lines encompass the 95% confidence interval.

The new result (Fig. 6) suggests a total full-depth ocean warming of 33.5 ± 7.0 × 1022 J (equal to a net heating of 0.37 ± 0.08 W/m2 over the global surface and over the 56-year period) from 1960 to 2015, with 36.5, 20.4, 30.3, and 12.8% contributions from the 0- to 300-m, 300- to 700-m, 700- to 2000-m, and below 2000-m layers, respectively. Here, we prefer to use the total energy change in the budget analyses rather than the linear trend because the change is not linear owing to the general increase in radiative forcing with time. The new reconstruction confirms the previous finding that the upper ocean experiences the most statistically significant warming, while the 0- to 2000-m layer contributes to the vast majority of the ocean warming since 1960 (1). By comparison, the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (1) reported an ocean heat gain of 25.5 ± 6.1 × 1022 J from 1971 to 2010, whereas the results of the current study are somewhat greater at 28.8 ± 4.4 × 1022 J for the same period, at the upper end of their uncertainty range (table S1). Furthermore, the full-depth OHC for 1970–2005 (26.5 ± 4.8 × 1022 J) agrees with a recent analysis (63) based on adjusted observation-based estimates, a simple transparent gap-filling method, and ORAS4 reanalyses (data assimilation method) (28.3 ± 1.8 × 1022 J); it also agrees with the ensemble mean of CMIP5 simulations (26.6 ± 4.4 × 1022 J) (table S1) (63). The consistency among these independent studies indicates the strength of the assessment of the ocean energy budget since 1960.

In addition to warming the oceans, the energy imbalance of the Earth system has also melted ice and warmed the atmosphere and land (4, 5, 64). These events account for ~7% of the EEI (1, 5, 64) but may be slightly smaller in the earlier years (64). Therefore, the total planetary energy imbalance incorporating all appreciable thermal reservoirs since 1960 is 36.0 ± 7.5×1022 J (equal to a net heating of 0.40 ± 0.09 W/m2 over the global surface during the 56-year period).

Although the EEI can best be estimated by OHC changes, the radiation measurements from space provide an important complement (4, 5). Tracking changes in the energy inventory of the planet gives an estimate of the energy imbalance. Figure 6 incorporates the accumulated net downward radiative flux imbalance at the TOA reconstructed by Allan et al. (65) since 1985. Because 7% of EEI is absorbed into other energy reservoirs, we multiply the TOA net radiation by 0.93 to estimate the ocean component in Fig. 6. The EEI based on OHC and TOA radiation shows consistent long-term changes since the late 1980s, further confirming the robustness of this new assessment of Earth’s energy inventory (tables S1 and S2). The inconsistency of the interannual variations of EEI based on TOA and OHC reveals the impacts of sampling errors in OHC, deeper OHC change below 2000 m, interannual variations in atmosphere/land/ice heat change, and errors in TOA radiation observations (66). Although there are errors in the various observation systems, comparison and reconciliation among independent observations (that is, OHC, TOA radiation, sea level, and also surface flux data) and models (that is, CMIP) will promote a more complete understanding of climate variability change and provide an opportunity for independent validation (1). Improving the existing integrated observation system, including the ocean and TOA observations, is therefore deemed essential for future progress.


Documentation of a detailed and reliable structure of ocean temperature changes is highly valuable in understanding the mechanisms of ocean heat variability driven by both internal processes and external forcing. This study provides a new reconstruction of the ocean temperature field (0 to 2000 m) and a quantification of its uncertainty on global and regional scales and interannual/decadal time scales based on subsample tests. These tests suggest a small sampling error related to temperature signals on decadal time scales. This new reconstruction provides a basis for accurately detecting and quantifying both regional and global variability across recent decades while minimizing error due to irregular ocean sampling by in situ observations.

This study aims to provide reconstructions on larger spatial and longer time scales (than mesoscale signals, eddies, and internal tides). Small-scale signals have been taken into account by smoothing in two ways: (i) collecting data within a time window for estimates in a given month and (ii) collecting data within the influencing radius, as carried out by most mapping methods (that is, the current method and NCEI), which acts as a spatial smoother. However, the small-scale signals could become errors in the reconstruction if they are spuriously mapped onto the large scale. This is where quality control becomes crucial. It is inherently subjective and based on knowledge of the ocean circulation and temperature/salinity patterns, but it remains the primary means to removing aliased mesoscale signals. Mesoscale eddies can be important in moving heat and fresh water around but do not show up in mean fields. The subsample test used in this study partly takes account of this error because gridded averaged data are used as truth without applying a spatial smoothing. One deficiency of this method is that the historical sampling error can be underestimated. We have tested for this effect and have shown that it does not affect the key conclusions of this study (section S4). It is difficult to fully address this issue because the Argo network does not fully capture mesoscale signals and lacks observations in coastal areas. High-resolution data, such as output from eddy-resolving ocean models, sea surface temperature, ocean color, and high-resolution ASL data (but only after 1993), will be helpful in further investigations of the impact of mesoscale signals, although the instrumental error of in situ profile observations may not be fully represented by these data.

Note that the new OHC analysis indicates only a weak ocean cooling after the Mt. Pinatubo eruption event in 1992 (Figs. 3 and 4), in contrast to the imbalance suggested by satellite data. Our subsample test suggests that interannual variations during the pre-Argo era are of comparable magnitude to sampling errors, but the current Argo system begins to allow for a resolution of interannual OHC variability. Therefore, if this volcanic eruption event were to occur with today’s observing network, we expect to be able to resolve its impact on OHC with greater confidence. Further reduction in the error in OHC estimates on interannual scales likely requires international coordination, including extending the current integrated observation system and understanding the influence of quality control processes, instrumental bias corrections, and ocean sampling. In addition, the improvement of methods that synthesize all available observations, such as reanalyses, will continue to be essential.



In situ temperature data from 1960 to 2015 were from the World Ocean Database (WOD) from the NCEI (67) and were prepared following CZ16. The quality flags provided by WOD were used to remove erroneous profiles and measurements. There are 41 standard depths in the upper 2000 m: 1 m, 5 m, 10 to 100 m in 10-m intervals, 120 to 200 m in 20-m intervals, 250 to 900 m in 50-m intervals, and 1000 to 2000 m in 100-m intervals. In the upper 700 m, we used data collected over a 3-month window to calculate monthly means, which helped to increase the data coverage without losing interannual signals, such as ENSO variations. To achieve improved spatial coverage before 2005 below 700 m, the data within 4-, 8-, 12-, and 18-month bins around the selected month were averaged together for the 750- to 1000-m, 1100- to 1500-m, 1600- to 1900-m, and 2000-m ocean layers, respectively, using moving averages. This choice was based on the amount of data available before 2005 in these layers. For models, we assembled 40 outputs of historical runs from 1940 to 2005 and 31 outputs using the Representative Concentration Pathways 4.5 scenario from 2006 to 2015 (68). Here, we examined the temperature reconstruction for 70°S to 65°N. For both observations and models, we constructed the climatology by using the data for the 1997–2005 period (see section S1 for a discussion of climatology).


The mapping method adopted in this study was proposed by Cheng and Zhu (40) and is called the ensemble optimal interpolation method with dynamic ensemble (EnOI-DE) provided by CMIP5 historical simulations. We modified this method to better represent the ocean variability on various spatial scales by using iterative strategies. Three iterative runs were carried out, where the analysis fields of the previous scan provided the prior estimate for the next scan (a thorough description can be found in section S2). Different influencing radii were set in each iterative run, that is, 20° in scan 1 (25° for 700- to 2000-m layers), 8° in scan 2, and 4° in scan 3. Furthermore, the present study extended the method to ocean depths below 700 m, where the observations were sparser than the upper 700 m. The temporal variations of deep ocean temperatures were much smaller relative to the upper ocean (37, 39). Therefore, to obtain a reliable assessment of the decadal and multidecadal variations of ocean temperature below 700 m, it is desirable to use the data within adjacent months to increase the spatial data coverage.

Subsample test

The subsample test was designed as follows. The 1°-by-1° gridded averaged temperature anomalies during the Argo period were used for each selected month (January and August from 2007 to 2014) as truth values (for a total of 16 truth fields). Each truth field was then subsampled according to the history of the observational locations (every 5 years from 1955 to 2014 for January and 1957 to 2014 for August). Thus, 24 subsampled fields were generated and then mapped using the CZ16 method and compared with the truth fields. The data during the Argo era not only observed ocean variability in the real world in each ocean layer but also contained instrumental errors in observations and insufficient coverage in coastal areas. However, because one could never repeat ocean observations made before 2000, the data during the Argo era were the best for examining the performance of the reconstruction for the historical period. One study (50) has subsampled ASL data to examine the sampling error in the OHC estimate. Using ASL has benefits for its complete global ocean coverage, and provides an important insight into the sampling error. However, ASL is a combination of full-depth thermal expansion and ocean mass change and is not directly representative of OHC change. Further, ASL is a two-dimensional metric without depth information, but the ocean sampling changes markedly with depth.

Statistical analysis

The S/N ratio for the ocean temporal variability compared to the sampling error is defined asEmbedded Imagewhere the signal is defined as twice the SD of the filtered temperature time series (Thistorical), which represents the magnitude of the OHC temporal variation. The noise is also twice the SD of the sampling errors. When S/N > 2, it indicates a significant signal given the sampling error; this is similar to a t test estimate of the >95% confidence level. When the ocean decadal/multidecadal variation was examined, a 7-year low-pass filter was applied to remove the ocean variability with frequencies higher than 1/7 per year. The method filtered a time series via the Lanczos method (cosine filter) (69). When the interannual variations were examined, the 7-year high-pass filter was applied to remove the ocean variability with frequencies lower than 1/7 per year.

Error bars of the OHC estimates (in Figs. 3, 4, and 6) are given as the ±2σ range, two SDs of the ensemble members obtained by the mapping method. ±2σ is equal to an approximate 95% confidence interval when assuming a Gaussian distribution, which means that there is only a 5% chance that the OHC falls outside the given interval. The distribution of ensemble anomalies at 20 and 1200 m is shown in fig. S19; although not strictly Gaussian, ±2σ captures more than 95% of this distribution.


Supplementary material for this article is available at

Supplementary Materials and Methods

fig. S1. Illustration of the iterative EnOI-DE/CMIP5 method used in the current study.

fig. S2. An example of the reconstructed fields after three iterative scans.

fig. S3. Reconstruction of temperature field at 1200 m in August 2011 for the historical sampling (in February).

fig. S4. Reconstruction of temperature field at 1200 m in August 2011 for the historical sampling (in September).

fig. S5. Mean temperature error as a function of different choices of the influencing radii between the reconstructed and truth fields.

fig. S6. Six major ocean basins defined in this study.

fig. S7. Global and basin-averaged sampling error compared with reconstructed temperature change at 20 m.

fig. S8. Global and basin-averaged sampling error compared with reconstructed temperature change at 300 m.

fig. S9. Global and basin-averaged sampling error compared with reconstructed temperature change at 800 m.

fig. S10. Global and basin-averaged sampling error compared with reconstructed temperature change at 1200 m.

fig. S11. Geographical distribution of mean and 2σ sampling error for 0- to 2000-m average.

fig. S12. Geographical distribution of mean and 2σ sampling error in 1°-by-1° grid at 20 and 1600 m.

fig. S13. 2σ sampling error for different scans at different depths.

fig. S14. Global OHC time series for the reconstructions after scan 1, scan 2, and scan 3.

fig. S15. Frequency distribution of temperature anomalies for the years 1986 and 2015.

fig. S16. Sampling error as calculated by two subsample methods.

fig. S17. S/N ratio analysis for two methods of subsample test.

fig. S18. Comparison between OHC and sea level change since 1993.

fig. S19. Distribution of the ensemble anomalies.

table S1. OHC trends obtained in this study for the 1960–1991 and 1992–2015 periods.

table S2. Net OHC and EEI changes obtained in the current study compared with some independent estimates.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank the three anonymous reviewers for their insightful comments on the previous version of this manuscript. Funding: L.C. and J.Z. were supported by the Chinese Academy Sciences’ Project “Western Pacific Ocean System: Structure, Dynamics and Consequences” (XDA11010405), the Natural Science Foundation of China (41506029 and 41476016), and 2016YFC1401800. K.E.T. and J.F. were partially sponsored by the U.S. Department of Energy (grant DE-SC0012711). J.F. was sponsored additionally by NSF Award ID 1243125 and NASA Award Number NNH11ZDA001N. The National Center for Atmospheric Research is sponsored by the U.S. NSF. T.B. was funded by the National Oceanographic and Atmospheric Administration Climate Program Office, OAR-CPO-NES-NCEI SLA 17-1-001, from 2016 to 2018. Author contributions: L.C. and J.Z. initiated this study, designed the experiments, and prepared figures and the first manuscript. K.E.T., J.F., T.B., and J.A. contributed to the design of the experiments and analysis of the results and refined the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The WOD13 data set is available at the National Oceanic and Atmospheric Administration/NCEI website The reconstructed data set is available at Argo data were collected and made freely available by the International Argo Program and the national programs that contribute to it (, The Argo Program is part of the Global Ocean Observing System. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article