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.

Keywords

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.

( 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.

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 ).

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.

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 , 31 – 37 ). 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.

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 , 6 – 10 ). For instance, tracking Earth’s heat and ocean heat is one of the key topics of the so-called “global warming hiatus” research surge ( 11 – 20 ). 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 , 21 – 27 ). 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 ( 28 – 30 ).

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 ).

RESULTS

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 (31–33, 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 (51–53).

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).