Observation-Based Estimates of Global and Basin Ocean Meridional Heat Transport Time Series

Ocean meridional heat transports (MHTs) are deduced as a residual using energy budgets to produce latitude versus time series for the globe, Indo-Pacific, and Atlantic. The top-of-atmosphere (TOA) radiation is combined with the vertically integrated atmospheric energy divergence from atmospheric reanalyses to produce the net surface energy fluxes everywhere. The latter is then combined with estimates of the vertically integrated ocean heat content (OHC) tendency to produce estimates of the ocean heat divergence. Because seasonal sea ice and land runoff effects are not fully considered, the mean annual cycle is incomplete, but those effects are small for interannual variability. However, there is a mismatch between 12-month inferred surface flux and the corresponding OHC changes globally, requiring adjustments to account for the Earth's global energy imbalance. Estimates are greatly improved by building in the constraint that MHT must go to zero at the northern and southern extents of the ocean basin at all times, enabling biases between the TOA and OHC data to be reconciled. Zonal mean global, Indo-Pacific, and Atlantic basin ocean MHTs are computed and presented as 12-month running means and for the mean annual cycle for 2000–16. For the Indo-Pacific, the tropical and subtropical MHTs feature a strong relationship with El Niño–Southern Oscillation (ENSO), and in the Atlantic, MHT interannual variability is significantly affected by and likely influences the North Atlantic Oscillation (NAO). However, Atlantic and Pacific changes are linked, suggesting that the northern annular mode (as opposed to NAO) is predominant. There is also evidence of decadal variability or trends.

Keywords: Atmosphere-ocean interaction; ENSO; Ocean circulation; Climate variability; Energy budget/balance; North Atlantic Oscillation

For many years, the observed average oceanic meridional heat transport (MHT) has been best estimated as a residual of the top-of-atmosphere (TOA) radiation and the atmospheric energy budget ([44], [45]; [47]; [50]; [25]). [8] made an early attempt that subsequently has been substantially refined. In [47], the estimates were based upon the Earth Radiation Budget TOA measurements of the mid-1980s, which have been supplanted by higher-quality Clouds and the Earth's Radiant Energy System (CERES) observations since March of 2000 ([23]). At that time, atmospheric energy transports and their divergence were estimated from two different atmospheric reanalyses that have now been much improved with later generations; in particular, the European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERA-Interim) has proven to be the most stable and best quality overall, at least until now (2018) ([55]; [51], [52]). Further improvements in the methodology, in particular adjustments for the inevitable spurious mass imbalance and allowance for the enthalpy associated with precipitation ([53]), have changed the results in small ways. As well as global estimates, the net surface flux of energy can be estimated locally and combined with estimates of changes in ocean heat content (OHC) to allow for an estimate of the movement of ocean heat around. A schematic of the residual methodology is given in Fig. 1.

Graph: Fig. 1. Schematic of the main flows of energy through the climate system for the (top) atmosphere and (bottom) ocean, including radiation at the top and surface R T and R s , surface sensible heat flux H s , and surface latent heat flux L E . The last three terms combine to give the surface flux F s (also a small precipitation enthalpy term is included). Latent heat is realized in the atmosphere as precipitation LP. The vector transports of total vertically integrated energy in the atmosphere F A and ocean F O are indicated.

The challenge is to produce these energy flows annually and even monthly in order to fully understand and appreciate the interannual and longer-term variability and trends of energy on Earth. How the atmospheric and oceanic energy disposition and transports vary and change is an important question because it influences where the surplus heat resides and potentially influences the subsequent atmospheric state and climate ([33]).

Tremendous progress has been achieved with the global ocean observing system ([1]). Argo is a system of autonomous profiling floats that drift with the currents at 1-km depth and provide profiles of temperature and salinity every 10 days or so for the top 2 km ([32]). Argo became global after about 2005 with over 3500 routinely deployed. Accordingly, the knowledge of OHC and its changes in space and time has increased enormously and enabled better reconstructions of the past ([10]). Nevertheless, observations of currents and heat transports remain spartan, although a major advance occurred in April 2004 with the establishment of RAPID-MOCHA (Rapid Climate Change–Meridional Overturning Circulation and Heatflux Array) across about 26.5°N in the Atlantic (henceforth simply the RAPID array) ([5]; [12]; [18]; [17]; [29], [30]; [41]). The variability in the Atlantic MHT arises mainly from changes in the AMOC volume transports ([17]) and is surprisingly large ([41]). Several other smaller arrays have also been established (e.g., [31]), but the expense and difficult logistics limit the capabilities for direct ocean observations.

Accordingly, complementary approaches that enable global coverage and extensions of the records back in time can help build understanding of variability and trends. Indeed, the zonal mean oceanic MHT can be readily computed from the TOA budget plus vertically integrated atmospheric energy transports, although some adjustments are necessary to deal with Earth's energy imbalance (EEI) and the changes in storage of OHC ([52]).

[57] examined six OHC analyses plus one dynamic ocean reanalysis and found that the implied change from month to month was not physically possible—energy was not conserved—as it was far too large compared with the changes inferred from CERES. The CERES 12-month running mean EEI variations had a standard deviation of 0.4 W m−2 (global), but all OHC analyses had values over 3.6 W m−2 with the sole exception being ORAP5 (Ocean Re-Analysis Pilot v5) from ECMWF ([62]), which was 1.4 W m−2 ([57]). In several products, the monthly anomaly standard deviation was analyzed to be >10 W m−2 versus the CERES estimate of 0.64 W m−2. Hence OHC analyses remain a substantial source of uncertainty. [52] combined estimated surface fluxes with OHC changes to deduce time series of vertically integrated ocean MHT throughout the Atlantic that could be verified by direct ocean observations at 26.5°N from RAPID. [33] used gridded analyses of OHC along with surface flux estimates to explore the main drivers of OHC variability, and also compared with the RAPID results.

As well as short-term weather-related variations, interannual variations in EEI associated with El Niño–Southern Oscillation (ENSO) are substantial, and typically are on the order of ±0.5 W m−2. They are associated with fluctuations in global mean surface temperature (GMST) ([54], [56]; [26]) as heat is stored in the oceans before being redistributed and some is released back into the atmosphere with an El Niño event ([11]).

The energy imbalance clearly varies over time ([56], [57]). The only way to determine the EEI is to perform a detailed inventory of the changes in energy in various forms in the climate system, the dominant component of which is changes in OHC ([59]). Because of changes in atmospheric temperatures, sea and land ice, land temperatures, and snow and water on land, there is no requirement for a match between TOA radiation and OHC tendencies on a monthly basis, but globally the mismatch is strongly constrained to be less than about 0.2 W m−2 [see [56] for an estimate both from observations and a climate model]. But the TOA radiation once calibrated plus the atmospheric reanalyses and OHC then provide the best information on the variability of the EEI and its regional manifestations. As shown here, it is necessary to adjust the OHC tendencies to match the EEI values, and by doing so, the results are "cleaned up" enormously.

In this paper we compute monthly MHT estimates, first globally, and then for two ocean basins. The first is the joint Indo-Pacific region, which is combined through the Indonesian Throughflow (ITF). Independent estimates of ITF would allow the results to be further broken out into the contributions from the Indian and Pacific Oceans, but this will be done elsewhere. The second is the combined Arctic and Atlantic Oceans. In both cases, as discussed below, the flow through the Bering Strait is sufficiently small that we can ignore it. We briefly consider the Indo-Pacific results in the context of ENSO, and we update the Atlantic results and compare with MHT from the RAPID array at 26°N and further consider the influence of the North Atlantic Oscillation (NAO).

We use monthly TOA CERES Energy Balanced and Filled (EBAF) Ed. 4.0 radiation on 1° × 1° grids ([23], [24]), from Langley Atmospheric Science Data Center (http://ceres.larc.nasa.gov/order%5fdata.php). Observations from CERES begin in March 2000, and have been extended back in time ([2]) using model results and other constraints. Nevertheless, as questions remain about the earlier reconstruction, those values are used here only for January–February 2000. The global mean net TOA radiation RT is too small to measure directly from satellite, and raw CERES analyzed values of global EEI of 6.5 W m−2 were much larger than the estimated 0.85 W m−2 owing to what are thought to be primarily systematic errors that were adjusted in a somewhat ad hoc manner based on overall OHC changes ([23]). However, instruments are far more stable than they are absolutely accurate, with calibration stability <0.3 W m−2 per decade (95% confidence) ([23]), and hence there is considerable confidence in the changes from year to year. Therefore, the long-term average absolute value of global mean RT is established from an inventory of the energy and, in particular, estimates of mean ocean heat uptake ([24]; [46]; [57]).

The atmospheric computations here all utilize only the ERA-Interim (ERA-I) ([13]) (http://data-portal.ecmwf.int/data/d/interim%5fdaily/), as they are superior in several assessments and much improved over earlier reanalyses (e.g., [55]; [51]). They have been comprehensively evaluated for conservation properties by [6] and for air temperatures and humidity ([35], [36]) and the water and energy cycles ([55]; [51]). ERA-I did not include comprehensive TOA forcings and volcanic aerosols, such as those from the eruption of Mount Pinatubo in 1991, and the TOA radiation is biased ([51]). Accordingly, we have here confined diagnostics to after 2000. The budget computations for the atmosphere are performed on 60 model levels every 6 h at T255 (about 79 km) resolution, and results are mass corrected to ensure that the atmospheric mass is conserved (see [48],[49]; [53]). All datasets were averaged or interpolated to a 1° resolution for the computations presented.

Here we mainly make use of ocean reanalyses from ECMWF, using ORAS5 ([63]), which was developed from ORAP5 ([52]), which was the best ocean reanalysis in the Arctic ([58]). ORAS5 uses the same ocean and sea ice model as for ORAP5 and was produced using the V3.4.1 of the NEMO ocean model at a resolution of 0.25° in the horizontal and 75 nonuniformly spaced levels in the vertical. ORAS5 uses four-dimensional data assimilation of multivariate fields, including sea surface height from altimetry, with a full global ocean circulation model used to carry information from past observations forward in time ([3],[4]; [63]; [28]). ORAS5 is an eddy-permitting ocean reanalysis with a prognostic thermodynamic–dynamic sea ice model with assimilation of sea ice concentration data and surface forcing from ERA-I. [63] note advances and improvements in a number of areas for ORAS5 versus ORAP5. It has been run in reanalysis mode through 2014 with consistent data streams, quality control, and forcing fluxes and has been extended into an operational setting in January 2015. For ORAS5, five ensemble members are generated by perturbing both observations and forcing fields, to reflect the main uncertainties in both, and we have analyzed results using all five, as well as the ensemble mean.

A concern with ORAP5, noted by [52], was a problem area in the midlatitudes of the North Atlantic, evidently associated with Mediterranean exchanges, that developed after about 2007. The problem was only evident below 1000-m depths with a vertical dipole structure (hence cancelling in part in the OHC) and is still present in ORAS5 but has been reduced. ORAS5 experiments ([63]) reveal maximum analyzed RMS temperature error at 1000 m associated with spurious convection between 1000 and 2000 m due to warm and salty Mediterranean outflow. A new "capping" procedure for salinity helps mitigate this effect. We computed all results using OHC down to both 1000 and 2000 m and, as the problem appears to have little effect on our results, we use the latter as it is more inclusive over the rest of the domain.

While other ocean analyses have been analyzed, we do not report the results here because the OHC is very noisy and the noise is spurious ([57]). Because the computation of OHC tendency inflates the noise (it is effectively a high-pass filter), we first take centered differences (thereby keeping the tendency centered on the time of interest), and then smooth the result with a 1/12 [1–3–4–3–1] filter that removes fluctuations shorter than 4 months.

For the global ocean, transports of energy (or any quantity) necessarily go to zero at the boundary. This applies universally at all times. Hence, an immediate check and constraint on the results of ocean heat transports is that they must go to zero at the North Pole and the edge of Antarctica. In general, we compute MHT results by integrating southward, so that values start at zero in the north, and the accumulated inferred value at Antarctica is the error. The errors primarily stem from the OHC mismatch with the CERES and Fs values, and we adjust the discrepancy on an annual basis to ensure that this constraint is satisfied. This is done uniformly over the ocean, so that a discrepancy at the Antarctic coast of 0.3 PW converts to an adjustment of about 0.8 W m−2 over the global ocean. This actually cleans up results from different OHC datasets by effectively bias correcting them. The documentation of this correction is given below (see Fig. 7). In reality, a small portion of the error may arise from CERES, and errors in energy divergence in the atmosphere ([52]), and there may also be small contributions from changes in sea ice, land ice, and E − P over land, as documented for a model by [57].

It is worth emphasizing that it is essential to carefully first balance the atmospheric mass budget, because small imbalances are magnified in any energy budgets. The procedures have been revised and improved by recognizing that the mass imbalance is identified with precipitation and moisture content of the atmosphere, so that rather than the previous barotropic adjustment, the profile is now weighted by the moisture profiles ([53]). The overall framework is illustrated schematically in Fig. 1.

The computation of the divergence of the vertically integrated atmospheric energy ([53]) includes a new formulation of the energetics and an allowance for enthalpy associated with precipitation. When combined with RT, the net surface energy fluxes are computed as a residual; see [20][20], [21]) and [52], [53]). The energy and water budget analyses can be assessed by how well closure is achieved and are ideally suited to analysis of large-scale variability, because errors from the method decrease with spatial averaging. In general, over the oceans, the errors in Fs, the net surface flux upward, are less than 10 W m−2 over 1000-km scales. Values are somewhat higher over land owing to complications with orography and the way surface moisture is restored in the reanalyses (providing spurious sources or sinks of water and thus energy). [20][20], [21] chose a different approach to adjust the land values but their method can lead to inconsistencies in the matching of the ocean plus atmospheric energy transports to the TOA radiation implied values.

The atmospheric energy transport ${\mathbf{F}}_A=1/g{\displaystyle {\int }_0^{p_s}\left(h+k\right)}\mathbf{v}dp$ , where h is the moist static energy h = s + Lq and k is the kinetic energy; s is the dry static energy computed relative to a reference temperature T0, s = cp(T − T0) + Φ; T0 is set to 0°C; and Φ is the geopotential. Here q is the specific humidity and L is the latent heat of vaporization. Setting the atmospheric energy AE = cp(T − T0) + Lq + k + Φs, where the subscript s refers to the surface, then

(1) ${\displaystyle \frac{\partial }{\partial t}}{\displaystyle \frac{1}{g}}{\displaystyle {\int }_0^{p_s}{A}_{E}dp+\nabla \cdot {\mathbf{F}}_A=R_T+F_s},$

where RT is the TOA radiation downward, and Fs is the total surface flux upward (Fig. 1). This consists of the surface sensible heat flux, the surface net radiation, and the surface evaporative cooling (latent heat flux), plus a small term from the enthalpy associated with transport of moisture and associated precipitation ([53]). Sensitivity to the reference temperature was analyzed in detail in [53] where they show benefits of using 0°C through a reduction in potential error sources (also [27]).

The atmospheric energy budget is used to compute the divergence of the total transport of atmospheric energy FA, which is balanced by the vertically integrated diabatic heating ${\tilde{Q}}_1$ , and the atmospheric moisture budget is used to compute the column latent heating ${\tilde{Q}}_2$ , ([48],[49]; [52]). Subtracting these two removes the dominating effects of precipitation and replaces it with the surface moistening (evaporation). The frictional heating ${\tilde{Q}}_f$ is very small and, as it is included in ${\tilde{Q}}_1$ , ${\tilde{Q}}_1\mathrm{-}{\tilde{Q}}_f$ is the nonfrictional heating. Although we include the frictional heating and the atmospheric tendencies, for simplicity these are ignored in the following equations.

(2) $\nabla \cdot {\mathbf{F}}_A={\tilde{Q}}_1-{\tilde{Q}}_2=R_T+F_s.$

Within the ocean, the net surface heat flux Fs is balanced either by changes in OHC or a divergence of vertically integrated ocean heat transport:

(3) $F_s=–d\text{OHC}/\mathit{dt}–\nabla \cdot {\mathbf{F}}_O,$

where FO is the transport of oceanic energy (Fig. 1). In areas involving sea ice, this should also include terms related to the transport of sea ice (a latent heat of fusion term effectively) as well as sea ice tendency. Over land, there may be transports of water and ice in rivers that could play a role, but these are relatively small. Seasonal snowfall can also play a significant role at high latitudes ([27]). Integrating (3) the zonal integral [⋅] of ocean variables from the North Pole southward gives the zonal mean MHT at latitude ϕ

(4) $\text{MHT}\left(\varphi \right)={\displaystyle {\int }_{\varphi }^{90}\left[F_s+d\text{OHC}/\mathit{dt}\right]}ad\varphi ,$

where a is Earth's radius. The zonal integral includes an a cosϕ factor.

Given estimates of RT from satellite measurements and using computed values for terms on the left-hand side of (1), Fs can be estimated as a residual. Over the oceans, Eq. (3) then allows estimates to be made of the vertically integrated divergent ocean heat transport component, given the changes in OHC. Equation (4) provides the meridional heat transport in the ocean.

In practice, we integrate from the north, so that the Arctic and Atlantic are included together. The net northward heat flux by the ocean through the Bering Strait varies substantially from year to year reaching a maximum in 2007 of 15–20 TW ([60], [61]), but even then the values are a factor of 20 smaller than the zonal mean transport and they are ignored here. As a result, for the Pacific, we integrate southward from the Bering Strait. However, the Indonesian Throughflow means that the Pacific and Indian Oceans are inextricably linked and here we have combined them. The ITF appears to be stronger during La Niña, influencing the heat budget in both Pacific and Indian Oceans on interannual time scales ([39], [40];?. [22][22]; [26], [28]) but owing to complications in leads and lags of ITF versus Niño indices, this aspect is dealt with elsewhere ([11]). We performed all calculations over the ocean at 1° resolution and took care to include all partially filled grid squares, but found overall that this refinement made very little difference.

It is desirable but very difficult to perform an error analysis. [23], [24]) extensively discussed errors in TOA radiation, and why their product is somewhat arbitrarily calibrated overall with a value that is too low ([57]; [10]) and that should vary in time. CERES assigns a single global value of EEI, while the procedures we use insist on constraints at all times, so that any bias is removed. We choose to use the best reanalysis available, rather than assign errors based upon comparisons with other reanalyses that have demonstrated problems (e.g., [7]). The derived surface flux has errors over the ocean estimated to be ~10 W m−2 on 1000-km scales (106 km2) ([53]). Similarly, rather than using multiple OHC analyses to assign the spread, because that result would be unduly large, we choose to use what we estimate to be the best product available based on evaluations ([57]; see also [58]. In addition, the smoothing we have adopted has knocked down the noise somewhat in the tendency term. The mismatch between the EEI inferred from OHC tendencies versus the CERES values is determined and removed, as documented below. The magnitude of this adjustment is likely a useful indicator of uncertainty. [52] made estimates of the uncertainty for the MHT at 26°N and assessed the contributions to uncertainty from both structural and temporal sampling of the TOA radiation, the atmospheric divergence, surface fluxes, and OHC tendency (see the online supplemental material), and similar values apply here. In the current analysis we use improved OHC analyses (ORAS5 vs ORAS4 or ORAP5) and we apply the physical constraint over the global ocean to further reduce the uncertainty.

The details of the TOA radiation and atmospheric diabatic heating and energy transports are given in [52], [53]) and the latter presents results for 2000 to 2016. The net annual mean result deduced as a residual of the TOA net radiation and vertically integrated atmospheric energy divergence for the net surface flux of energy Fs (Fig. 2) reveals where there is a net flux of energy into the ocean (in blue areas) versus the losses of ocean heat (in red areas) so that in the absence of changes in storage of heat within the ocean, the ocean must transport heat via the ocean currents from the blue to the red areas. This is not easy to measure because here we are dealing with only the divergent transport, while in reality there is also a substantial rotational transport whereby the heat just moves around. This is especially the case for the Antarctic Circumpolar Current, for instance. The tropical Pacific is a large heat source region for the ocean while the largest sinks lie off the coasts of Asia in the Pacific and North America in the Atlantic in the vicinity of the Kuroshio and the Gulf Stream. These fluxes are largest in winter.

Graph: Fig. 2. Annual average net surface energy flux F s for 2000–16 in W m−2. Positive values are upward out of the ocean.

The fluxes of energy into the atmosphere over the Arctic and the Atlantic cannot be met by the flux of energy into the Atlantic tropical ocean and, as a result, there has to be a transport from the Pacific and Indian Oceans into the South Atlantic and a northward MHT throughout the Atlantic to balance the energy flows. In the Indian Ocean, there is a strong source region of heat in the equatorial zones that exceed transports out of the southern Indian Ocean and, combined with the net westward flow through the ITF, contributes to a net southward flow of heat across the southern boundary of the Indian Ocean.

The long-term 2000–16 annual average zonal mean MHTs (Fig. 3) show the overall energy transports that result from the integrated divergences. On the left-hand side, the CERES TOA radiation provides the total transport once the EEI is removed. The atmospheric reanalyses in turn provide the atmospheric transports [see [48],[49]) for details]. Because they involve the divergence, they automatically integrate to zero globally. The ocean MHTs are computed as a residual. Contributions from land are assumed to be small: the only meridional transports come from north–south flowing rivers, as diffusion within land is extremely small. The right-hand side (Fig. 3) then breaks the ocean down into the contributions from the basins, where use is now made of the OHC tendencies to match the TOA radiation and remove the energy imbalance or, equivalently, spurious OHC storage. The results are artificial for the Indian and Pacific between about 7° and 40°S because the ITF is ignored. The results are almost identical using all 5 ORAS5 ensemble members. Here the Arctic is combined with the Atlantic, and south of 35°S the ocean basin contributions become meaningless as the zonal transports become large (south of Africa), and indeed continuous around the globe near 60°S through the Drake Passage between South America and the Antarctic Peninsula.

Graph: Fig. 3. Zonal-mean annual long-term average meridional transports in PW for 2000–16. (left) Those inferred from TOA radiation (black), within the atmosphere from ERA-I (red), and the residual for the oceans (blue). (right) The breakdown for the oceans for the Atlantic (purple), Pacific (red), and Indian (green) and combined for the Southern Ocean south of 35°S (blue) in PW. (bottom) The domains used and the standard errors are 1 σ. The ITF transport is not included here.

The long-term mean annual cycle of these values (Fig. 4) based upon the adjusted transports (see discussion related to Fig. 7) does not include the contributions from sea ice melt or freeze, nor does it include the effects from the strongly seasonal runoff from land. [51] estimated Fs for the large continents and showed that the annual cycle of land surface fluxes can have amplitudes of 20 W m−2, peaking in northern summer for Eurasia and North America and with opposite sign in winter. However, these are mostly accounted for, because we use Fs values only over the ocean. [50] made the first attempt to estimate the annual cycle of ocean MHT using indirect methods.

Graph: Fig. 4. The inferred mean annual cycle of ensemble mean meridional heat transport (MHT) by the oceans in PW for the (top) global, (middle) Indo-Pacific, and (bottom) Arctic–Atlantic ocean domains in PW. These do not allow for seasonal influences of runoff, ice thaw, and freeze.

In the Atlantic at 26°N for the RAPID array, [17] report on the mean annual cycle of MHT based only on the 3.5-yr time series, and the error bars are quite large. Using slightly revised and updated RAPID estimates of MHT (courtesy of W. Johns), Fig. 5 gives updated RAPID results. For comparison, our results include those for the five ensemble members of ORAS5 as well as the ensemble mean. Here we used the 1σ error bars from [17] and a rough estimate of the error bars on our results ([52]). We also suspect that there is an offset of 0.15 PW for the RAPID values, as discussed below.

Graph: Fig. 5. Estimates of the mean annual cycle of MHT at 26°N, from RAPID based upon April 2004–February 2017 updated from Johns et al. (2011) (red); our initial result (blue) including the ensemble mean and all 5 members for ORAS5 (thin lines); an approximate estimate of the effect from melting and freezing of sea ice from PIOMAS (gold) as the long-term average plus the values for each year from 2000 to 2016 (gray); and an adjusted result including the latter (black). The error bars are 1 σ.

To roughly estimate the effects of freezing and melting sea ice, we use the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) estimates of volume changes from [34] (http://psc.apl.uw.edu/research/projects/arctic-sea-ice-volume-anomaly/). Although uncertainties exist, as discussed by Schweiger et al. these estimates are credible and seem to be among the best available. The annual cycle of Arctic sea ice volume has a range from about 27 to 11 × 103 km3 or 16 × 1012 m3. Assuming a latent heat of fusion of 3.34 × 105 J kg−1 and ice density of 900 kg m−3 results in a change in MHT that occurs almost entirely in the Atlantic. The maximum sea ice volume occurs in April and the minimum in September, implying a smaller northward midlatitude MHT from October to March and a larger value from May to September (Fig. 5). The values vary from year to year and we have taken the monthly mean values from 2000 to 2016 to estimate the monthly changes in volume converted to an implied MHT associated with the melting and freezing of sea ice via the latent heat of fusion component. This assumes that the OHC accounts fully for the changes in ocean temperature in the vicinity of the sea ice, which is unlikely to be the case, given poor observations among sea ice. All 16 years are included in Fig. 5.

The influence of the change in phase with sea ice forming or melting is quite large (Fig. 5), and ranges over 1 PW in MHT with a maximum MHT northward of 0.6 PW in July as sea ice melts (relative to the annual mean), and a broader minimum of −0.4 PW peaking in December. Our result based on the residual technique has a maximum in March and a minimum in July, but with an average of 1.06 PW, somewhat smaller than observed by RAPID of 1.21 PW. If we add our result to that from PIOMAS, the peak moves to June and July, and the annual cycle is more in accord with that from RAPID although with the largest discrepancy in October. At that time, sea ice is beginning to freeze and latent heat of fusion effects reduce MHT, and the main discrepancy could be OHC tendencies that are relatively too negative. We consider it more likely that OHC errors occur in the season of extensive sea ice when observations are few and altimetry is not available (i.e., October departures are influenced by the annual mean).

The RAPID values are slightly higher from June to November and are on the order of 1.3 to 1.4 PW and lowest in March near 1.0 ± 0.2 PW. In Fig. 4 for the Atlantic, our values are positive year-round but slightly higher for November through June at 26°N and with values of 0.9 (July, September, October) to about 1.3 PW in March–April. Hence it is worth noting that the mean for all months from both estimates lies between 0.9 and 1.4 PW, well within the 2σ error bars. The inclusion of the sea ice contribution improves agreement considerably.

Accordingly, as a first approximation, we adopt the mean annual cycle effects from changes in sea ice in Fig. 5 and apply them to Fig. 4 for the Atlantic and global domains everywhere south of 70°N ([34], to give Fig. 6.

Graph: Fig. 6. The inferred mean annual cycle of ensemble mean MHT by the oceans in PW for the (top) global, (middle) Indo-Pacific, and (bottom) Arctic–Atlantic domains in PW adjusted to include changes in sea ice in the Atlantic south of 70°N.

Returning to Fig. 4, we note the large annual cycle of MHT across the equator into the winter hemisphere but entirely arising from the Indo-Pacific oceans. The new results refine those of [50] although the broad picture given by the latter remains correct. MHT in the Atlantic is northward year-round north of 10°N once sea ice effects are accounted for (Fig. 6), but weak southward transport is evident from September to December in the southern Atlantic.

The results of the raw global MHT for the 12-month running mean total and anomalies using the ensemble mean ORAS5 OHC (Fig. 7) reveal the problems over the southern oceans, referred to above. The cleaned-up version is given in Fig. 8 and the problem is revealed in Fig. 9. The latter shows the implied values in Fig. 7 at the Antarctic coast near 78°S (which includes the sea ice zone) for the ensemble mean OHC and the individual members. In Fig. 9 we see the results of the energy imbalance where the OHC tendency has not matched the Fs and implied RT at TOA. [14] note that ERA-I surface fluxes are biased by an amount on the order of 5 W m−2 globally and the ability of an analysis to correct that bias will depend on the ocean observations. The mismatch is mostly positive from 2000 to 2005, negative for 2006 through 2014, and wildly negative in 2015. Also shown in Fig. 9 are some key factors about why the values may have changed. Before about 2006, the ocean observations were much more limited in the pre-Argo era, and in January 2015 ORAS5 switched to an operational setting, with clear consequences for the overall EEI implied by the changing OHC. For the Argo era, implied imbalances produce a spurious MHT ranging from about 0 to −0.5 PW; 0.3 PW globally is equivalent to about 0.8 W m−2 imbalance over the global ocean—and on average may be attributable to the way CERES was adjusted. This spurious imbalance, which changes over time, is what has been applied as a uniform correction for each year over the global ocean to produce Fig. 8. The comparison of Figs. 7 and 8 then reveals the effects of the adjustments, which is mostly imperceptible north of about 40°S, but large over the southern oceans.

Graph: Fig. 7. Latitude–time series of zonal mean MHT for the global ocean as the (top) total and (bottom) anomalies (departure from annual cycle of Fig. 4) in PW. Values are 12-month running means computed using 13 terms with the first and last point weighted by 0.5. This shows the result for the ensemble mean for the raw calculation.

Graph: Fig. 8. Latitude–time series of zonal and ensemble mean MHT for the global ocean as the (top) total and (bottom) anomalies (departure from annual cycle of Fig. 4) in PW. Values are 12-month running means adjusted to remove the bias at 78°S.

Graph: Fig. 9. Zonal mean raw MHT at 78.5°S in Fig. 7 for the global ocean in PW for the ensemble mean (red) and 5 members (gray). Values are 12-month running means. The Argo period spun up after about 2005, and the ORAS5 analyses changed to an operational setting in January 2015.

Also shown in Fig. 9 are the results for the individual ensemble members from ORAS5, and the spread is often about 0.4 PW, and at times higher still. The cleanup procedure described above is very beneficial in narrowing the spread among the ensemble member results.

The MHT time series in Fig. 8 may seem unduly spotty in the deep tropics, but most of this comes from the Indo-Pacific (Fig. 10) and is real. Indeed, the variability is associated with ENSO events. In the Indo-Pacific there is consistent northward MHT by the ocean north of about 5°N which is the location of the intertropical convergence zone (ITCZ). The largest value and farthest south extension of the northward MHT is in 2015–16 in association with the major El Niño event. The anomaly field also reveals a penchant for values from the Southern Hemisphere to progress northward over about 6 months to the equatorial region, especially from 2008 on.

Graph: Fig. 10. Latitude–time series of zonal and ensemble mean MHT for the Indo-Pacific oceans as the (top) total and (middle) anomalies in PW. Values are adjusted 12-month running means. (bottom) ONI index smoothed with 12-month running mean.

The Indian Ocean MHT (not shown) is southward everywhere, because of the continental boundaries. It appears that the variability in the Indian Ocean is mostly negatively correlated with the Niño-3.4 sea surface temperature (SST) index, also called the ONI (Ocean Niño Index; Fig. 10), while the Pacific Ocean MHT is pronouncedly positively correlated. However, the ITF adds ambiguity to this and instead we focus on the combined result (Fig. 11). The lead and lag correlations between ONI and MHT for the global ocean, Indo-Pacific, and Atlantic (Fig. 11) show the dominance of the Indo-Pacific. On the equator, the strongest correlations of about 0.65 occur about 3 months after the ONI. However, much stronger correlations are evident throughout the southern tropics with a maximum near 20°S of 0.76, and with the MHT leading ONI by some 6 months or so. Values greater than 0.46 are significant at the 5% significance level.

Graph: Fig. 11. Correlations between the oceanic Niño index (ONI) and zonal and ensemble mean MHT for (top) the global ocean, (middle) the Indo-Pacific, and (bottom) the Atlantic as a function of lead and lag. Negative values mean that MHT leads ONI, while for positive values ONI leads. The 5% statistical significance level is 0.46 and values greater than 0.55 are outlined in white.

The details of movement of heat around will mostly be covered elsewhere, including in [11], but we have examined correlation and regression maps for these times (6 months before and 3 months after the ONI) (not shown). For 6 months prior to the El Niño peak there is a huge dipole structure to the ocean heat divergence with divergence west of the date line and convergence to the east, so that a great deal of this cancels in the MHT. However, there is also a region of divergence of heat from 15° to 30°S, east of the date line, linked to the convergence farther north. By 3 months after the ONI, the dipole has largely abated and the region east of the date line from 10°S to 5°N is a region of divergence of ocean heat with strong convergence from about 5° to 20°N, extending up the coast of the Americas to Alaska ([54]).

In El Niño events, the deep warm water stored in the warm pool in the tropical western Pacific shoals as warm waters surge eastward in association with a downwelling Kelvin wave. The latter deepens the thermocline in the central and then eastern Pacific and results in warmer waters being upwelled, greatly changing the SSTs. Subsequently, warm waters spread poleward along the coast of the Americas. In the Northern Hemisphere warm waters can spread all the way to the Gulf of Alaska. The regional details are further documented by [11], who also perform a detailed energy budget analysis and show that the tropical Pacific Ocean loses heat during the later stages of an El Niño event to other parts of the tropics, and then heat moves out of the tropics. The biggest losses of ocean heat are through evaporation and manifested as increases in rainfall over the central and eastern Pacific.

Clearly, much more can be done with these results, but here the focus is on MHT.

New results for the Atlantic (Fig. 12) update those from [52], [53]) with the ORAS5 dataset. There are a number of modest changes as a result. The main picture is much the same, however, with coherent fluctuations from the equator to about 45°N. MHT tends to be somewhat higher in the first part of the record until 2009 when values drop precipitously, and with a smaller repeat in 2012. These fluctuations can be seen more clearly in the anomalies in MHT (Fig. 12), which are quite a lot weaker than for the Pacific.

Graph: Fig. 12. Latitude–time series of zonal and ensemble mean MHT for the Atlantic Ocean as the (top) total and (bottom) anomalies in PW. Values are adjusted 12-month running means.

In the RAPID array, MHT is computed from several components from the Florida Straits, western boundary array, interior Atlantic mean, interior eddy gyre, and Ekman transports. Changes from the methodology used in [17], described in [30], involve use of ERA-I surface winds to compute the Ekman transports, assumed for the upper 50 m of the column. The midocean eddy heat flux is derived from an objective analysis of interior Argo profiles and temperature/salinity profiles from the RAPID moorings. Meridional velocity anomalies across the section are derived from this analysis using a geostrophic approximation relative to 1000 m. The time-varying interior zonal average temperature transport is derived from Argo and mooring data, and merged into a seasonal temperature climatology below 2000 m ([17]).

At 26°N (Fig. 13), our MHT is slightly higher than in [52], [53]) using the OHC down to 2000 m, and the offset for the RAPID values is 0.15 PW. [14], [15]) and [42] also find values 0.2 PW lower at 26.5°N, with much lower values in 2004, as we find, and perhaps this is associated with sampling issues arising from 2004 hurricane activity damage in western boundary current regions. Questions have been raised about whether the RAPID array properly resolves the transport processes near the western boundary and whether the assumption of geostrophy falls short in the presence of ocean recirculation ([42]). [37] explore uncertainties in RAPID estimates and note the dependence on the reference level used.

Graph: Fig. 13. The 12-month running mean northward heat transports across 26°N (blue) from Fig. 12 compared with results from the RAPID array (red) in PW, with an offset of 0.15 PW (right scale). The error bars are ±1 σ , for RAPID in pink and for current results in light blue. Also in black is the NAO index (right scale).

The large variability in the MHT zonal mean time series is primarily associated with the atmospheric circulation and, in particular, with modes of natural variability such as the NAO (Fig. 14) and the closely related northern annular mode (NAM) ([16]). The latter is the dominant mode of wintertime variability for sea level pressure for the Northern Hemisphere as a whole. It has a strong signature into the stratosphere and connects the Pacific to the Atlantic. In Fig. 14, positive correlations in the Atlantic of up to 0.6 between the equator and about 45°N are verified with the NAO leading by about 2 months. At higher latitudes between 50° and 60°N the correlation reverses and is up to about −0.4, marginally significant.

Graph: Fig. 14. Correlations between the NAO index and zonal and ensemble mean MHT for (top) global ocean, (middle) Indo-Pacific, and (bottom) Atlantic as a function of lead and lag. Negative values mean that MHT leads NAO, while for positive values NAO leads. The 5% statistical significance level is 0.46 and values greater than 0.55 are outlined in white.

Globally, however, the relationships are dominated by the Pacific, reinforced by the Atlantic, with highest correlations of 0.6 between 25° and 30°N, with the NAO leading by 2 to 3 months. The predominance of positive correlations between MHT and NAO at multiple lags and latitudes suggests decadal or lower frequency variability. In Fig. 13 the NAO has large variability and perhaps a slight upward trend, but the MHT at 26°N is distinctly lower, with a small downward trend ([38]), especially relative to the NAO.

Given the obvious NAO interannual effects on MHT, we have performed a regression of MHT on the zero-lag NAO index and removed that component from the anomalies in Fig. 12 (see Fig. 15). The main reduction in variance is in midlatitudes. For instance, at 30°N the original standard deviation of 0.08 PW is reduced to 0.06 PW after removal of the linear NAO contribution of up to 0.05 PW standard deviation. However, exploring the northern annular mode (as opposed to the NAO) and leads and lags further may be worthwhile.

Graph: Fig. 15. As in Fig. 12 for the anomaly, but with the linear NAO effect at zero lag removed.

While it is evident that an NAO contribution remains, the modest downward trend, especially north of about 30°N, is more evident in Fig. 15. Although the trend in our MHT product is not as large as for RAPID in Fig. 13, the residual after removing the linear NAO influence does indeed suggest decadal variability in AMOC. Nevertheless, any trend is also complicated by trends in Arctic sea ice, but this effect seems to be quite small. For 2000 to 2016, the net trend in latent heat of fusion from melting Arctic sea ice (about 4000 km3 decade−1 from PIOMAS) is about 0.004 PW on a global annual basis, much smaller than the total MHT at 26.5°N. This is only 1% to 2% of the total increase in OHC. Therefore, the impact of sea ice melt is negligible in our calculation for long-term change.

In the South Atlantic, our results are quite compatible with those of [31] based on ocean reanalyses, but there is some evidence of a downward trend, although with large uncertainty (Fig. 12).

In the previous section, the results from this work were all for the product using ORAS5 to provide the OHC tendencies, including all five ensemble members. We have performed similar calculations using other OHC products. However, the results are not as good or credible. As shown by [57], other OHC datasets violate energy constraints in terms of the tendencies much more than for ORAS when compared with TOA radiation changes, and the variance is therefore a lot greater. Much of it is demonstrably spurious. Nonetheless, many aspects of the analyses in the different basins are confirmed. There is better agreement among products in the Pacific, including many regional aspects of the ocean heat divergence with ONI and thus the general patterns of change with ENSO ([11]). In part this is likely due to the larger signal-to-noise ratio. At the same time, there is clearly scope to do better. Also, there are a number of other OHC datasets that could be worth exploring ([14]; [33]; [19]; [43]; [9]).

In this paper we have focused on MHT. A key point is that by putting together the latest and best datasets for TOA radiation, atmospheric reanalyses, and OHC, we have been able to construct time series of MHT for the global, the Indo-Pacific, and the Atlantic–Arctic ocean domains as a residual. In [47] we were able to construct time average annual mean MHT for the oceans, and in [50] we constructed an approximate annual cycle, and only now can we progress to MHT as 12-month running means, and refine the mean annual cycle. The latter is recognized as somewhat incomplete owing to lack of consideration of sea ice melt and freeze, other than rather crudely, and also effects from river runoff. Examining the departures from the mean annual cycle largely removes those influences, and the estimates have been greatly improved by building in the constraint that MHT must go to zero at the northern and southern extents for the ocean basin at all times. This enables the biases between the CERES and OHC data to be reconciled.

For the Pacific, we have presented new results for MHT in association with El Niño, and for the Atlantic we can compare with the results from direct ocean estimates at the RAPID array. A major advantage over in situ observations is the global mapping and extensions in time. By estimating the TOA radiation further back in time ([2]) estimates can be made of surface fluxes ([20]) and the potential is to be able to further extend these results. However, limitations with OHC analyses and sea ice contributions remain.

In the Atlantic, our results suggest that the mean MHT is about 0.15 PW less than given by RAPID and there are good reasons why this may be correct. In addition, there are some homogeneity concerns with RAPID related to changes in instrumentation and loss of data, especially in the first 2 years when hurricanes were especially active. But there are also remaining concerns with all datasets used in our product. The atmospheric reanalyses are being redone (e.g., ERA-5), with some notable upgrades. How CERES data are adjusted to get the correct EEI remains an issue, and this then relates to the monthly reconciliation of the global mean energy budget. This can only get better as these data are utilized together. However, even in the current form, the information on AMOC is invaluable.

In the Pacific, regional and global aspects of the energy budget with El Niño are given in [11]. Our results can be improved by splitting out the Indian from the Pacific data using information on the ITF. Moreover, the dominant east–west movement of heat can be explored using our datasets. New information on ENSO evolution is thereby emerging.

The results we have generated can clearly be exploited much further to explore regional influences and questions such as what really is the ocean heat transport across the equator. There are challenges in using the methods outlined here. The biggest relate to the need to carefully balance mass budgets in the atmosphere, making adjustments for EEI in TOA radiation, closing the energy and water budgets over land, including effects of snow and ice, and accounting for precipitation enthalpy better along with liquid and frozen water in the atmosphere. In addition, the magnitude of the latent heat of fusion effects for changing phase of seawater to and from sea ice is substantial yet quite uncertain, and the OHC in the vicinity and under sea ice is likely poorly known. Use can also be made of the multiple surface flux fields, which have their own strengths regionally (e.g., [33]; [58]). Ocean reanalyses provide estimates of currents and heat transports, and deserve to be evaluated. Finally, climate models can be used to produce similar plots to those shown here and they can be evaluated in new ways, leading to future improvements.

The EEI perspective has implications for the energy accumulation and regional differentials, which thus affect the future evolution of the system and risk of extremes, and there is a lot of information in the coupled climate system not being utilized in many analyses. The water and energy constraints should be built into or at least checked for ocean analyses of OHC, such as has been shown in proof of concept by [43].

NCAR is sponsored by the National Science Foundation. Data from RAPID-MOCHA are funded by the U.S. National Science Foundation and U.K. Natural Environment Research Council and were downloaded from https://mocha.rsmas.miami.edu/results/index.html and we thank Bill Johns for the updates. We thank the reviewers, especially reviewer 1 for thoughtful comments. The values will be made available at NCAR at https://doi.org/10.5065/9v3y-fn61.