Abstract

This study investigates the multidecadal warming and interannual-to-decadal heat content changes in the upper ocean (0–700 m), focusing on vertical and horizontal patterns of variability. These results support a nearly monotonic warming over much of the World Ocean, with a shift toward Southern Hemisphere warming during the well-observed past decade. This is based on objectively analyzed gridded observational datasets and on a modeled state estimate. Besides the surface warming, a warming climate also has a subsurface effect manifesting as a strong deepening of the midthermocline isopycnals, which can be diagnosed directly from hydrographic data. This deepening appears to be a result of heat entering via subduction and spreading laterally from the high-latitude ventilation regions of subtropical mode waters. The basin-average multidecadal warming mainly expands the subtropical mode water volume, with weak changes in the temperature–salinity (θS) relationship (known as “spice” variability). However, the spice contribution to the heat content can be locally large, for example in Southern Hemisphere. Multidecadal isopycnal sinking has been strongest over the southern basins and weaker elsewhere with the exception of the Gulf Stream/North Atlantic Current/subtropical recirculation gyre. At interannual to decadal time scales, wind-driven sinking and shoaling of density surfaces still dominate ocean heat content changes, while the contribution from temperature changes along density surfaces tends to decrease as time scales shorten.

1. Introduction

The global ocean has been warming with an almost linear trend since the 1950s (Levitus et al. 2012). While the overall observational record is sparse during the 1950s and 1960s, some regions, notably the North Atlantic, have relative dense early observations, and the nearly ubiquitous presence of warming lends confidence to the overall result. Most attention has focused on the upper 2 km of the water column, yet warming is present also in the Southern-origin deep and bottom water masses, connected with deepening of the isopycnal surfaces (Purkey and Johnson 2012; Gille 2008). The upper ocean heat content trend has not been distributed evenly around the global ocean, with the North Atlantic Ocean experiencing an especially concentrated warming along the Gulf Stream and North Atlantic Current (Levitus et al. 2012). North Atlantic heat content analysis from several ocean state estimates (Häkkinen et al. 2015, hereafter HRW15) shows that the increase in ocean heat content has been associated with increasing volume of Subtropical Mode Waters (STMW; potential density 26.0 < σ0 < 27.0), while the northeastern Atlantic Subpolar Mode Waters (SPMW; 27.0 < σ0 < 27.6) have lost volume. The subpolar Atlantic as a whole has a strong decadal variability mode superimposed on a weak but significant warming trend (HRW15).

The best observational data coverage is limited mostly to the upper 700 m over the last 50+ years; this depth range covers subtropical mode waters adequately, but subpolar waters less optimally (Lyman and Johnson 2008, 2014). After 1990 the global observing system has greatly improved during the WOCE (1990–98) and ARGO (2006–present) programs. In the Southern Hemisphere (SH) the 0–700-m depth range is occupied by southern subtropical waters, Subantarctic Mode Waters (SAMW) (26.0 < σ0 < 27.1), while the subpolar mode water range includes the Antarctic Intermediate Waters (AAIW) (27.1 < σ0 < 27.6). The lighter SAMW densities (<26.6) are also called subtropical mode waters. Both SAMW and AAIW are formed from upwelled upper Circumpolar Deep Waters within the Antarctic Circumpolar Current (ACC). Exposed to surface heating, northward Ekman transport and mixing with fresh and cold Antarctic surface waters, these upwelled waters are then subducted north of the ACC.

Here our goal is to evaluate the historical patterns of heat content and temperature changes since the 1950s. We will characterize climate warming effects entering the thermocline from the subduction regions and advecting laterally, which is a separate mechanism from vertical diffusion of surface heating or deep convection. One signal of subduction of warmed upper-ocean waters is the change of isopycnal depth (“heaving” of isopycnals; Bindoff and McDougall 1994, hereafter BM). Hence besides focusing on heat content variability, we will analyze subsurface isopycnal depth changes. Warm subduction is a likely explanation for results by Wang et al. (2015), who showed that simulated future global warming manifests in a deepening of the 26.8 potential density surface in CMIP5 results. We also earlier described multidecadal deepening of isopycnals in the North Atlantic Ocean with a maximum at densities 27.0–27.3 (HRW15).

Our main focus is on global patterns of heat content variability with limited discussion of regional variations. These can be affected by very low-frequency atmospheric and ocean variability (the prime example is the North Atlantic) and their full discussion is beyond the scope of this study. We utilize widely available gridded datasets for our heat content analysis and invoke a decomposition of temperature variability into changes of isopycnal depth (heaving) and of the water mass properties (“spice” as in BM). (This decomposition implicitly involves also salinity fields through changes on neutral density.) We will show that on multidecadal scales, heaving of isopycnals dominates 50-yr trends of the heat content variability in all basins. While we are interested in multidecadal variability, we also explore interannual-to-decadal heat content variability, by making use of satellite altimeter data. The role of atmospheric forcing in the multidecadal deepening of subtropical mode water isopycnals is reviewed at the end.

2. Datasets and methods

For our analysis, we apply widely used, objectively mapped gridded hydrographic datasets. From the NOAA–NODC observational database we utilize annual in situ temperature T and salinity S anomalies interpolated from pentadal averages (1957–present) with 1° spatial resolution (Levitus et al. 2012). (NOAA and NODC also have an annual temperature dataset 1955–present, which was not used here). The interpolated data will cause weak biases in computing potential and neutral density but we focus on the top 700 m, where density anomalies are mostly determined by temperature anomalies outside the high latitudes. As described by Levitus et al. (2012), missing data regions, especially widespread in the 1950s to 1970s are filled by assuming zero T and S anomalies relative to a climatology estimated for 1955–2006.

The UK Met Office observational potential temperature θ and salinity S EN4.0.2 dataset (monthly; 1° resolution) is analyzed from 1955 to the present (Good et al. 2013). All monthly data are converted to annual data before computations. In EN4, missing data points are assigned a climatological value causing considerable small-scale variability, which is difficult to remove without extensive smoothing. Its SH data are biased to the more recent ARGO era with a high data coverage, but analyzed trends for the period ending in 2003 are shown in our supplemental material and confirm that the improved coverage by ARGO observations does not alter the conclusions.

The construction of the above NODC and EN4 datasets involves estimates for filling missing data regions, particularly widespread in the SH during the 1950s and 1960s. Abraham et al. (2013) have extensive discussion of the temperature observations and treatment of missing data, and corrections made to mechanical and expendable bathythermographs (MBTs and XBTs), which constitute a large portion of the pre-ARGO observations. Several choices of climatologies or estimated anomalies used to filling in missing data are discussed by Lyman and Johnson (2008, 2014).

To support the results from the observational datasets we analyze the Simple Ocean Data Assimilation (SODA) reanalysis for 1955–2011 (v. 2.2.8, monthly, data available from J. Carton; ½° resolution; here data have been converted to annual averages) (Carton and Giese 2008). SODA fills in the missing data by optimizing the model physics and forcing and available ocean observations. Overall the benefit of using SODA is to uncover robust large-scale patterns based on sparse observations, which is particularly important in the early period in the Southern Hemisphere. Also SODA has a smooth transition to the ARGO era compared to the observational datasets, which is especially noticeable in the Southern Ocean, where the global trend is concentrated in the ARGO era (Roemmich et al. 2015).

Additionally for recent interannual heat content investigation we utilize AVISO altimetric sea surface height data (1993–present; daily with 0.25° resolution). These data are converted to monthly and/or annual data for the analysis. We also provide limited analysis of atmospheric wind stress and derived variables based on the NCEP–NCAR reanalysis.

In the analysis we consider ocean temperature change at a fixed depth due to heaving and spice variability: heaving is an Eulerian measure of the temperature change at fixed depths, due to vertical migration of isopycnals, either adiabatically (as with changes in wind forcing of a gyre circulation) or through diabatic heat flux divergence. Spice variability is a shift of θS profiles, hence a change in θ and S upon a fixed neutral-density surface (BM). In the BM decomposition of potential temperature θ (similarly for the salinity), the changes at depth z (/dt|z)

 
formula

are divided into a change along the neutral density surface, spice (the first term on right), and, a change due to vertical movement of the neutral density surface, heave (dz/dt) (the second term on right). An example sketch of pure warming leading to negative spice and positive heave is shown in Fig. 1a on the θ–S plane. Both θ and S have positive gradients (i.e., warm/salty over cold/fresh), which is the typical stratification in the subtropical gyres. Water is warmed at a fixed depth, its θ moving from point 1 to point 2; this is the sum of vector 1 to 3, which is heaving along the θS profile, then corrected by vector 3 to 2, which is spice change at constant potential density σ (light solid lines). Conditions at the original potential density σ1 become colder and fresher (point 4). By contrast, subpolar regions with fresher water above more saline deep waters tend to have warm heaving accompanying a warm spice change: this is readily seen by altering the slope of the θS profile in Fig. 1 from positive to negative. The sign of spice is determined by the slope of the temporal changes at fixed pressure and depth, δθ/δS, compared to the background slope of the θS profile.

Fig. 1.

(a) In a θ–S diagram, pure warming, increasing θ from 1 to 2 at a fixed depth, is depicted as the vector sum of movement of fixed depth points along the original θ–S profile without altering it (“heave,” 1 to 3) plus coordinated change in θ and S without change in potential density (“spice” change, 1 to 4). Now parcel 2 appears cooler (and fresher) than the parcel with the same density (point 3) on the original thermocline (solid blue). (b) The same warming effect in a z plane. Warming (1 to 2) moves the density from 3 to 4. Parcel 5 with the same density as before warming (parcel 3) is much cooler than the original parcel (at point 1): Parcel 5 has a corresponding temperature at point 6. Connecting the lines from 1 and 6 to the salinity curve would show that the parcel 5 is also fresher than the parcel 3.

Fig. 1.

(a) In a θ–S diagram, pure warming, increasing θ from 1 to 2 at a fixed depth, is depicted as the vector sum of movement of fixed depth points along the original θ–S profile without altering it (“heave,” 1 to 3) plus coordinated change in θ and S without change in potential density (“spice” change, 1 to 4). Now parcel 2 appears cooler (and fresher) than the parcel with the same density (point 3) on the original thermocline (solid blue). (b) The same warming effect in a z plane. Warming (1 to 2) moves the density from 3 to 4. Parcel 5 with the same density as before warming (parcel 3) is much cooler than the original parcel (at point 1): Parcel 5 has a corresponding temperature at point 6. Connecting the lines from 1 and 6 to the salinity curve would show that the parcel 5 is also fresher than the parcel 3.

Another view of the pure warming process at constant pressure in the case of warm/salty over cold/fresh stratification is displayed in Fig. 1b. Warming at fixed depth (1 to 2) moves the density from 3 to 4. A parcel at 5 with the same density as before warming (parcel at 3) appears now much cooler than the original parcel (at 3) with temperature at point 6. Connecting the lines from 1 and 6 to the salinity curve would show that the parcel at 5 is also fresher than the parcel at 3. The density change at fixed depth is small but the destabilizing effect of the salinity gradient creates a large negative θ change on isopycnal. After the warming, the parcel at 5 with the original density is found at much greater depth than the parcel at 3. The depth change, or “warm heave,” represents global warming effects in the subsurface waters.

In theory the decomposition equation is exact, but in practice the irregular vertical sampling/gridding and treatment of surface changes creates a residual (e.g., it is linearized with respect to vertical displacements and requires an approximation of /dz from the gridded data, which are widely spaced in the deeper ocean). To compute the neutral density values for each dataset, we used software available online at http://www.teos-10.org/preteos10_software/neutral_density.html (Jackett and McDougall 1997; BM). Total heat content change and its heave and spice components are computed by multiplying potential temperatures by the gridpoint volume, density, and heat capacity of seawater, and integrated over 700 m (and 2000 m).

3. Multidecadal ocean heat content evolution and its decomposition

We first focus on the total 0–700-m ocean heat content which we call simply OHC, computed as annual values for latitude range 60°S–65°N from SODA, EN4, and NODC (Fig. 2). Note the latitude limits: the calculation does not include the entire World Ocean. Overall, the OHC in all basins shows a well-known long-term increase, with the variance explained by a linear trend for each dataset listed in Table 1. In the NODC data 67% or more (over 88% in the Atlantic sector) can be explained by a statistically significant linear trend that is consistent with the Levitus et al. (2012) values. The variances explained by a linear trend in EN4 are smaller, 58%–84% for the Atlantic and Indian Oceans, while the Pacific sector variance explained is less than 16% (barely significant at 95%). In SODA the linear trend explains 67%–91% of variance. In the 2000s the North Atlantic and North Pacific OHC increases have weakened according to NODC and EN4 data. Because of the pentadal derived temperature anomalies, the NODC evolution is very smooth and removes the hiatus signal visible in the EN4 OHC. The OHC evolution shows the largest basin changes in the North Atlantic peaking in 2004, but the Indian Ocean has continued a rapid increase after 2003, with its OHC change exceeding all other basins. The recent period with the most complete global data shows a dramatic shift of warming to the Southern Hemisphere [Roemmich et al. (2015) based on ARGO data (2006–13)]. EN4 is similar to NODC in the Atlantic and Indian Oceans but displays large decadal fluctuations in the Pacific basins. The Pacific basins and Indian Ocean in SODA reanalysis show the largest trends compared to NODC and EN4.

Fig. 2.

(top three rows) Ocean heat content (1022 J) integrated over 0–700 m in NH and SH basins and their heave and spice components in NODC, EN4, and SODA, respectively. Heat contents are referenced to 1955. (bottom) The global heat content (1022 J) in 0–700-m and 0–2000-m columns, with and without the Indian Ocean contribution; heat contents are referenced to 1955, except NODC to 1957. Values are smoothed by one binomial filter for SODA and NODC, twice in EN4.

Fig. 2.

(top three rows) Ocean heat content (1022 J) integrated over 0–700 m in NH and SH basins and their heave and spice components in NODC, EN4, and SODA, respectively. Heat contents are referenced to 1955. (bottom) The global heat content (1022 J) in 0–700-m and 0–2000-m columns, with and without the Indian Ocean contribution; heat contents are referenced to 1955, except NODC to 1957. Values are smoothed by one binomial filter for SODA and NODC, twice in EN4.

Table 1.

Variance explained by a linear trend in the basin average OHC and its heave and spice components in each basin.

Variance explained by a linear trend in the basin average OHC and its heave and spice components in each basin.
Variance explained by a linear trend in the basin average OHC and its heave and spice components in each basin.

To address the changes in the ARGO era since 2003 and discussion of a possible hiatus in global warming (Easterling and Wehner 2009; Balmaseda et al. 2013; Kosaka and Xie 2013), we combine OHC from all basins for a global trend (bottom row of Fig. 2). NODC and EN4 continue to show positive global trends, albeit weaker, since the early 2000s. A previous global hiatus (in NODC data or a decline in EN4 data) occurred in the 1980s, thus hiatus events are not unusual. Both NODC and EN4 suggest that warming in the Pacific and Atlantic sector started to recover after 2010. The difference between the global heat content with and without the Indian Ocean increases significantly after 2003. This implies that the weak global OHC trend since 2003 originates mainly from the Indian Ocean warming, while warming in the Pacific and Atlantic sectors has stalled or declined.

Using the BM decomposition of potential temperatures to divide OHC changes into heaving and spice components, we see that in all basins, in all datasets, there is a positive trend in the heaving component similar to the trend of the total OHC (Fig. 2). The positive heaving denotes deepening of the isopycnal depths. In all basins except in the SODA North Pacific the 50-yr trend is largely determined by heaving. The maximum heaving occurs in the Southern Hemisphere, either in the Indian Ocean in SODA and NODC, or in the South Pacific in EN4. The vertical structure of heaving and spice contributions to warming is discussed in section 6.

The variance explained by a linear trend in the OHC heave and spice components is shown in Table 1. Only NODC and SODA have a large fraction of the OHC heave component explained by a linear trend. The linear trend fraction in the spice components is more variable with only the NODC and EN4 Pacific basins reaching highly significant values (and the SODA North Pacific, although there is a question of the sign of the trend).

In basin averages the spice component is of secondary importance compared to the heave in all basins except in the SODA North Pacific (Fig. 2). A cooling/freshening spice trend also opposes the heaving in all southern basins. NODC data show robust changes in all basins compared to EN4 where the amplitude of the spice components is small and noisy (Fig. 2). Unlike the other basins in NODC and EN4 data, the North Atlantic spice displays quasi-decadal fluctuations, with a warmer/saltier spice trend prevalent since the 1990s. In SODA, the North Pacific spice component is an outlier as it has a positive trend (i.e., waters are warming on the midthermocline isopycnals, and hence adding to the positive heaving trend).

The heat content for the 0–2000-m column contains the same trends as OHC for 0–700 m but with enhanced amplitudes, shown for EN4 and SODA (see Fig. S1 in the supplemental material). Their residual OHCs, also in Fig. S1, obtain their amplitude (1 × 1022 J) mostly from the thin surface layer. In NODC the global heat content change 1957 to 2014 in the latitudes 60°S–65°N is 18.0 × 1022 J for 0–700 m and 27.6 × 1022 J for 0–2000 m (i.e., about 35% of global heat content change in the 0–2000-m column is occurring in the 700–2000-m depth range). The global average heat flux required for the 50+ years of warming in NODC is 0.3 and 0.45 W m−2 for the 700- and 2000-m columns, respectively. The corresponding values of the 700- and 2000-m columns for EN4 are 0.2 and 0.3 W m−2 and for SODA are 0.41 and 0.55 W m−2. These heating rates are comparable with estimates of top-of-atmosphere radiation imbalance of ~0.5 W m−2 (Stephens et al. 2012).

As an example of the meridional structure of the OHC 50-yr trend and its decomposition, we construct 15-yr averages for 1997–2011 and 1957–71 to suppress interannual variability, and hence enhance multidecadal signal. We show results for NODC in the Atlantic, Pacific, and Indian sectors and their global sum for (1997–2011) − (1957–71) (Fig. 3a). (EN4, very similar to NODC, and SODA decompositions are shown in Fig. S2.) The largest changes in OHC are confined to subtropical latitudes in both hemispheres as already shown in Levitus et al. (2012). These are here seen to be associated with heaving of isopycnals, yet with regionally significant spice contribution. The trends in the Indian Ocean contribute the most to the SH zonal average at latitudes 35°–50°S; however, that may be due to inclusion of ARGO data. In the zonal average frame, the Southern Hemisphere and North Pacific negative spice anomaly shows a strong damping effect on the heaving contribution, consistent with the subtropical water column illustrated in Fig. 1. A notable exception to this isopycnal cooling/freshening is the subtropical North Atlantic, where from 20° to 40°N the spice is positive; at 20°N it is even the dominant contributor to OHC. This is consistent with the large evaporative salinity increase directing the thermohaline trend vector along isopycnals (Curry et al. 2003). The meridional structure for the ARGO-era OHC change (2003–08 to 2009–14) from EN4 (NODC update not available yet) is shown in Fig. 3b, where the Pacific and Atlantic sectors display large changes across latitudes averaging to (almost) zero net change. The Indian Ocean shows a positive net gain over its north–south extent. Because of the short averaging span the OHC differences may contain heaving arising from low-frequency Rossby waves.

Fig. 3.

The zonal OHC (0–700 m) difference (1020 J per one degree of latitude) vs latitude and its decomposition are shown for the Atlantic, Pacific, and Indian sectors and globally: (a) (1997–2011) minus (1957–71) from NODC and (b) (2009–14) minus (2003–08) from EN4. Total OHC is marked black; heave and spice components are red and green respectively.

Fig. 3.

The zonal OHC (0–700 m) difference (1020 J per one degree of latitude) vs latitude and its decomposition are shown for the Atlantic, Pacific, and Indian sectors and globally: (a) (1997–2011) minus (1957–71) from NODC and (b) (2009–14) minus (2003–08) from EN4. Total OHC is marked black; heave and spice components are red and green respectively.

4. Spatial variability of heat content from decadal to multidecadal scales

a. OHC on multidecadal time scales

The spatial variability of OHC over the 50+ years as a difference from the early (1957–71) to the late (1997–2011) years is depicted in Fig. 4 for NODC, EN4, and SODA (for an enlarged image of the NODC panels, see Fig. S3a). Their common pattern shows decreased OHC in the tropical Pacific and increased OHC in the subtropical gyres in both hemispheres. We have added the long-term mean of the zero-wind stress curl lines to the spice and heave panels (the zero-curl lines show minimal movement; Fig. S4a). Most of the OHC increase is derived from increased downward heaving of the isopycnals as shown in Fig. 2 and in the second column of Fig. 4. The spice component tends to oppose heaving with the most negative values located in the poleward side of the subtropical gyres, downstream of the subduction/ventilation regions (in vicinity of the zero wind stress curl), as would be expected based on warm subduction (or warm heave) (Church et al. 1991; BM). In NODC and EN4 datasets the zero-curl lines enclose the enhanced heave and negative spice in the subtropical gyres in the Pacific and Indian Oceans. The Atlantic Ocean shows increased heave almost everywhere, and the North Atlantic subtropics shows mainly positive spice in contrast to the South Atlantic and other basins. The SODA heave component is very similar to the observational datasets, but its spice component differs from observations especially in the Pacific sector. Examples of θS diagrams in locations representative of subtropical regions with strong spice and heave are shown in Fig. S5. They show a θS shift toward colder and less salty at the subtropical mode waters densities, with an exception of the central North Atlantic, where the θS shift is toward warmer and more saline. The magnitude of heave over the Gulf Stream appears to be an outlier.

Fig. 4.

Spatial pattern of heat content change [109 J m−2 (OHC divided by the grid area)] from 1957–71 to 1997–2011 in NODC, EN4, and SODA. The black curves in the heave and spice components represent the long-term mean zero line of the wind stress curl. Stippling denotes 95% significance for the difference with 28 degrees of freedom.

Fig. 4.

Spatial pattern of heat content change [109 J m−2 (OHC divided by the grid area)] from 1957–71 to 1997–2011 in NODC, EN4, and SODA. The black curves in the heave and spice components represent the long-term mean zero line of the wind stress curl. Stippling denotes 95% significance for the difference with 28 degrees of freedom.

On multidecadal scales heaving could originate from low-frequency ocean dynamics, such as strengthening of the gyres, and the meridional overturning circulation (MOC) in the North Atlantic; however, we have no means of assessing these effects directly from hydrographic observations. Second, these processes compete with the main cause for heaving, which is horizontal advection and mixing of warmed surface layers and their subduction into subtropical gyres as originally proposed by Church et al. (1991). In section 7 we will discuss the role of surface forcing in light of the recent modeling results by Wang et al. (2015) and Marshall et al. (2014), which suggest that surface heating effects of global warming on ocean are of primary importance with other forcing components, such as winds and precipitation minus evaporation (PE), being of secondary importance.

b. OHC on interannual to decadal time scales; comparison with altimetry

To explore the more recent changes we first investigate whether satellite altimetry will support the spatial patterns seen in OHC derived from observations. The recent decade is also of interest because of the apparent hiatus in global surface air temperature warming and its relationship with OHC (Balmaseda et al. 2013; Kosaka and Xie 2013; England et al. 2014). Figure 5 displays the sea surface height (SSH) and EN4 OHC trend from 1993 to 2014. Most notably, the western tropical Pacific has experienced a large increase in both OHC and SSH compared to the eastern tropics, creating an east–west dipole. The trends are significant at 95% level, particularly in the western tropical Pacific. This large-scale east–west dipole in the Pacific OHC and SSH resembles the negative phase of the Pacific decadal oscillation (PDO) in sea surface temperature (SST) anomalies. This mode is a contribution to the temporary OHC hiatus in the Pacific as described by Kosaka and Xie (2013). We also show the trend of OHC heave and spice components, where the 22-yr trend in the Pacific OHC dipole arises from heaving of isopycnals, that is, from the acceleration of the easterly tropical winds over this period as described in England et al. (2014). This dipole is an example of the OHC heave component driven by winds, which are the dominant factor at interannual to decadal scales. Away from the tropics, the OHC spice component, while opposing the heave component and SSH trend, shows surprising strength on decadal time scales as shown by the significance of its trends in Fig. 5. Notably the spice component shows strength in the same areas as the multidecadal change in Fig. 4, in the poleward sides of the subtropical gyres (except the North Atlantic).

Fig. 5.

AVISO sea surface heights trend (computed from monthly data; cm decade−1) and EN4 OHC 0–700-m trend and their heave and spice component trends (computed from annual data; 1018 J decade−2). Stippling denotes 95% significance.

Fig. 5.

AVISO sea surface heights trend (computed from monthly data; cm decade−1) and EN4 OHC 0–700-m trend and their heave and spice component trends (computed from annual data; 1018 J decade−2). Stippling denotes 95% significance.

Strong Southern Hemisphere warming is evident in EN4 and the altimetric SSH rise along the rim of the ACC and neighboring subtropical gyres. However, correlations between SSH and EN4 OHC are strongest outside the Southern Ocean where the geostrophic dynamical contribution to SSH is strongest (Fig. S6). The global ARGO profiling float array was fully populated by about 2006, and warming trends during 2006–13 (Roemmich et al. 2015) show the concentration of ocean warming in the Southern Hemisphere to be even more extreme.

5. Multidecadal vertical migration of isopycnals

Deepening of isopycnals dominates the heat content evolution in all basins on multidecadal scales and in HRW15 the maximum deepening was found to be at the subtropical mode water densities. To find if this is valid globally, we first show the vertical zonal average sinking of potential density surfaces in Fig. 6 (negative values denote shoaling) as a trend over 1957–2011 from all three datasets. [The trends are not affected if the ARGO era is excluded (Fig. S7a) or if early data, such as the trend for 1970–2011, are neglected (Fig. S7b).] NODC trends below 27.2 are less reliable due to the use of 5-yr mean salinities. In all basins the largest deepening, over 50 m per 50 years, occurs in the 26.5 to 27.3 isopycnal range extending from the poleward edge of the density range (i.e., from the ventilation region) deep into the subtropics. This density range consists of the subtropical mode waters (STMW and SAMW) and the lightest subpolar mode waters (SPMW and AAIW). However, it is the sinking of the STMW (<27.0) and of SAMW (<27.1) layers that is pushing down the SPMW and AAIW isopycnals. The hemispheres are separated by a shoaling trend in the equatorial regions with stronger shoaling south of the equator than north of it in all datasets. The basin average sinking is 20–30 m from the early period (1955–69) to the late period (1997–2011) peaking at densities σ0 = 27.0–27.2 (Fig. S8), as also found for the North Atlantic in HRW15.

Fig. 6.

Zonal average sinking/shoaling trends [m (50 yr)−1; positive downward] of potential density surfaces for the (left to right) Atlantic, Pacific, and Indian Oceans for 1957–2011 shown for (top to bottom) NODC, EN4, and SODA. Trend values significant at 95% level are stippled, the black contour is 50 m per 50 yr.

Fig. 6.

Zonal average sinking/shoaling trends [m (50 yr)−1; positive downward] of potential density surfaces for the (left to right) Atlantic, Pacific, and Indian Oceans for 1957–2011 shown for (top to bottom) NODC, EN4, and SODA. Trend values significant at 95% level are stippled, the black contour is 50 m per 50 yr.

The spatial variability of the isopycnal heaving trend for 1957–2011 is shown for the depth of the 26.0, 27.0, and 27.3 potential density surfaces in Fig. 7. [The trends before the ARGO era (Fig. S9a) or for 1970–2011 (Fig. S9b) are very similar.] The EN4, NODC, and SODA results in Figs. 6 and 7 (see also Figs. S7 and S9) are highly consistent, giving mutual credibility to both the ocean state estimate and the observational datasets, which were constructed from very limited data in the early part of the analysis period. While SODA relies on the same sparse observations, the model physics were adequately constrained to recreate the ocean state and evolution.

Fig. 7.

Sinking/shoaling trends (positive downward) of potential density surfaces (left to right) 26.0, 27.0, and 27.3 for 1957–2011 shown for (top to bottom) NODC, EN4 and SODA [m (50 yr)−1]. Trend values significant at 95% level are stippled.

Fig. 7.

Sinking/shoaling trends (positive downward) of potential density surfaces (left to right) 26.0, 27.0, and 27.3 for 1957–2011 shown for (top to bottom) NODC, EN4 and SODA [m (50 yr)−1]. Trend values significant at 95% level are stippled.

At density σ0 = 26.0 surface, NODC, EN4 and SODA subtropical gyres show almost uniform sinking of isopycnals by 20–25 m per 50-yr period (Fig. 7). The equatorial regions display a shoaling trend in all basins as seen also in the meridional section (Fig. 6). The STMW–SPMW (SAMW–AAIW) interface isopycnal of 27.0 shows a strong deepening in all basins. In the SH, deepening rates can reach 100 m in 50 years, compared to about half that rate elsewhere (Fig. 7). The strong deepening over the North Atlantic Current over the past 50 years (Fig. 7) and associated heat content changes have been attributed to the positive shift in the North Atlantic Oscillation (NAO) index from the negative values in the 1950s and 1960s to mostly positive values since the 1990s (Leadbetter et al. 2007; Lozier et al. 2008; Williams et al. 2014). This shift is associated with increased Ekman pumping leading to downward heaving of isopycnals (see the references above), and large-scale gyre and MOC circulation changes (Zhai and Sheldon 2012). However, deepening isopycnals elsewhere do not necessarily require change in Ekman pumping, as will be discussed in section 7.

The middle isopycnal 27.3 of the SPMW–AAIW layer also shows a widespread sinking trend in the same SH and North Atlantic locations, but weaker than for 27.0, because it is pushed down by the layer above. As a new feature the tropical Atlantic (and also the EN4 tropical Pacific) has increased sinking. Multidecadal changes in isopycnal surfaces in EN4 and NODC display robust changes in large regions, which are captured by the SODA results, in terms of both regional patterns and magnitude. In all datasets, sinking of σ0 = 27.0 and 27.3 surfaces appears to be maximum in the Southern Hemisphere downstream from their ventilation region. The only notable sinking activity in the Northern Hemisphere is along the Gulf Stream and North Atlantic Current.

6. Temperature changes at midthermocline isopycnals

Here we explore the vertical profiles of basin averaged temperatures and their heave and spice components from the early period (1957–71) to the late period (1997–2011) in NODC, EN4, and SODA (Fig. 8). The total warming and its heave and spice component amplitudes weaken but maintain the same structure if the data after 2003 (ARGO era) are excluded (Fig. S10). In basin averages, the heave contribution is at maximum in the top 200–800 m with larger basin amplitudes in SODA compared to the observational datasets. This depth range defines subtropical mode waters, SAMW and STMW, in every basin. The basin average spice components below 200 m are cooling in all basins, with two exceptions; the North Atlantic has near zero spice change and the SODA has a positive spice change in the North Pacific.

Fig. 8.

Total potential temperature change (°C) 1957–71 to 1997–2011 divided into heave and spice for (left to right) NODC, EN4, and SODA for the SH and NH basins. (Residuals dominating the top 50 m not shown.)

Fig. 8.

Total potential temperature change (°C) 1957–71 to 1997–2011 divided into heave and spice for (left to right) NODC, EN4, and SODA for the SH and NH basins. (Residuals dominating the top 50 m not shown.)

In the examples of BM a negative spice trend (i.e., cooling on isopycnals) is expected when parcels (depicted in Fig. 1) with a warm or fresh anomaly subducts into subtropical regions stratified stably by temperature and unstably by salinity (climatological θS relation with positive slope). In general, isopycnal changes in θ and S depend on the “thermohaline vector,” which on the θS plane expresses the combined changes in θ and S at a fixed pressure. The direction of this vector on the θS plane compared with the background θS profile determines whether a water parcel will gain or lose spice as seen in Fig. 1.

The spatial variability of the isopycnal cooling/freshening trend close to the maximum heaving, at potential density surface σ0 = 27.0, is shown in Fig. 9 for 1957–2011. (The trend with ARGO era excluded is shown in Fig. S11). The large blue regions of cooling/freshening on this potential density surface are notable. The strongest cooling (a consequence of warm heave) is located in the southern basins north of the ACC and gives an appearance in all datasets that cooling is spreading northward away from the ACC. SODA fields are similar to EN4 and NODC except in the North Pacific, where SODA displays an opposite sign. In the dataset constructed by Durack and Wijffels (2010) the salinity trend (1950–2008) at σ0 = 26.75 (to be compared with Fig. S12 for our analyzed temperature trend at σ0 = 26.75) has a similar pattern of negative trends as in Fig. 9, with a strong freshening, associated with cooling, dominating the Southern Hemisphere.

Fig. 9.

Potential temperature trend [°C (50 yr)−1] on potential density surface = 27.0 for NODC, EN4, and SODA for period 1957–2011. Trend values significant at 95% level are stippled.

Fig. 9.

Potential temperature trend [°C (50 yr)−1] on potential density surface = 27.0 for NODC, EN4, and SODA for period 1957–2011. Trend values significant at 95% level are stippled.

7. Atmospheric forcing and the multidecadal deepening of subtropical isopycnals

The roles of wind changes and of global warming effects on the ocean on multidecadal time scales have been assessed in the recent model simulations and analysis by Wang et al. (2015) by comparing CMIP5 historical and future CO2 warming simulations. These simulations use coarse-resolution ocean components, which are usually too diffusive and do not resolve western boundary processes well. Interestingly, as one of the global warming signals, the CMIP5 ensemble of simulations and offline ocean model with CMIP5 derived forcing predicted a deepening of the 26.8 potential density surface. Wang et al.’s Figs. 2b and 6a have a great similarity with our Fig. 7 for the 27.0 surface with largest deepening north of the ACC and dominance of sinking in the Atlantic. An uncoupled ocean model forced by globally uniform SST anomaly alone reproduces subsurface warming similar to the spatially varying SST of the full CMIP5 experiment. They suggest a weak role for the wind anomaly compared to that of globally uniform SST anomaly.

In a somewhat similar approach, Marshall et al. (2014) used estimates of surface heating from CMIP5 doubled CO2 warming simulations. Their ocean model was forced by a globally uniform heating of 4 W m−2 with climatological wind and PE forcing. This uniform heating anomaly reproduced the surface SST structure of the CO2 warming in the ensemble CMIP5. They show that warming reaches almost to 500–700 m in the southern and northern subtropics. Despite being a coarse-resolution model, it shows that the background ocean circulation, including MOC, can create large heat content anomalies in the presence of a uniform surface heating anomaly. They also note, like Wang et al., that wind and PE anomalies are of secondary influence in their simulation of effects from CO2 warming. The conclusions from these two studies only apply to multidecadal time scales and not to interannual to decadal time scales, when winds are expected to be important as seen in section 4b. As earlier noted, a multidecadal shift in the NAO index and possible changes in MOC likely impacts our North Atlantic results, but detailed discussion is beyond the scope of this study.

In contrast to the Wang et al. and Marshall et al. results, some studies have found that multidecadal wind changes over the ACC can be important globally. The southern annular mode (SAM) index describes broadly the Southern Hemisphere atmospheric variability, and the index intensified since the 1950s and particularly from the 1970s onward (Marshall 2003). Ocean modeling studies, both coarse and eddy resolving, and analytical studies (Klinger and Cruz 2009; Wolfe and Cessi 2010; Radko and Kamenkovich 2011) show that adjustment to increased wind stress over ACC will result in a downward adjustment of isopycnal surfaces globally (i.e., changes in the meridional overturning circulation). However, the intensification of SAM does not translate to a continuous intensification of the westerlies [weakened until 1978 and increased after 1980, as shown by Fan et al. (2014)], and hence these results are not readily applicable for the changes seen in the last 50+ years.

The Wang et al. and Marshall et al. results imply that changes in wind stress curl (i.e., in Ekman transport and pumping) are secondary in the oceanic response to multidecadal global warming. Taken together with our result regarding deepening of midthermocline isopycnals, similar to Wang et al. any trend in Ekman pumping is unlikely to be an important factor. Indeed, the long-term trend in Ekman pumping (Fig. S4b) does not appear to have a strong one-to-one correspondence with maximum sinking trend. Figure 10 shows the deepening of the 27.0 surface from Fig. 7 for NODC, now overlaid by the outcropping of densities 26.0 and 27.0 from the early and late period and by the mean zero-line of wind stress curl (the zero-curl line shows very little movement; Fig. S4a). The almost stagnant outcrop lines define the ventilation regions, which overlap with the largest deepening (except in the Gulf Stream region) and with mean anticyclonic wind stress curl (Fig. S4c). Stagnancy of outcrop lines is consistent with Gille (2014), who finds no significant movement of ACC and associated frontal system in the altimetry era. Another common feature of the ventilation regions is that they and the adjacent poleward areas are all located in westerly wind regimes and have equatorward Ekman transport (Fig. S4a). Thus at high latitudes, warmed-up surface parcels are moving with equatorward Ekman transport toward the subduction region of the subtropical mode waters. This subduction of the warmed parcels would occur even without any change in winds since the ventilation region is quite stagnant, as is the collocated zero line of the wind stress curl. The anomalous heat enters subsurface via subduction and spreads laterally into the thermocline as envisioned by Church et al. (1991). One signal of this process is the warm heave, the deepening of isopycnals, within the subducted water masses as seen in Fig. 1b.

Fig. 10.

The heaving trend [m (50 yr)−1] of the NODC 27.0 isopycnal from Fig. 7 with 26.0 and 27.0 outcrops; solid black line is the winter average 1957–71, dashed is 1997–2011. Green lines are the zero contour of the long-term mean wind stress curl.

Fig. 10.

The heaving trend [m (50 yr)−1] of the NODC 27.0 isopycnal from Fig. 7 with 26.0 and 27.0 outcrops; solid black line is the winter average 1957–71, dashed is 1997–2011. Green lines are the zero contour of the long-term mean wind stress curl.

8. Summary

Our analysis of historical ocean heat content reveals how the global warming effect during the last 50+ years is distributed in the upper ocean above the permanent pycnocline, where the largest heat gain has occurred. The vertical and lateral structure of θ variability, decomposed into heave and spice after Bindoff and McDougall (1994), provide a framework to analyze both multidecadal ocean state estimates and ocean observations. In separating trends of the θ–S relation from the vertical displacement of isopycnals, the decomposition relates closely to water mass evolution and to air–sea forcing of the thermohaline shifts, as well as trends in general circulation and winds.

In basin averages the dominant contribution to the fixed depth heat content comes from deepening of isopycnals: near-surface deepening is caused by downward diffusion of heat, but maximum deepening occurs at depth at densities of the subtropical mode waters. This subsurface maximum is a result of heat entering thermocline from the high-latitude subduction regions and advecting laterally. This deepening reaches nearly 100 m in 50 years on the poleward side of the subtropical gyres. Confidence in the basin-average and zonal-average results arises from the monotonic nature of the OHC increase and its heave component in most of the World Ocean.

This Bindoff and McDougall approach applies also to high-resolution ARGO data in order to analyze interannual-to-decadal scale variability, which can aid interpretation of altimetric data. We show that during the last 20 years the Pacific OHC contains large changes from isopycnal sinking and shoaling (heaving) episodes due to PDO variability. On decadal and longer time scales the spice component shows strength away from equatorial regions in all basins.

At the same time as the multidecadal OHC evolution is dominated by the heave component, the OHC spice component shows local importance by opposing the heave in the subtropical regions. An exception is the strongly evaporative central North Atlantic subtropics with shoaling isopycnals leading to positive OHC spice contribution. Cooling at midthermocline isopycnals is consistent with Church et al. (1991) and Bindoff and McDougall (1994). They suggest that in a typical subtropical stratification (warm/salty over cold/fresh), the mixed layer under a warming climate will subduct at a lighter isopycnal, making the new subsurface waters appear cooler (and fresher) on their original isopycnal. Bindoff and McDougall also predicted the sinking of isopycnals due to warming (and freshening), which we here have shown to be active globally and reaching a maximum at around σ0 = 27.0–27.2 isopycnals. Since the deepening of isopycnals is not constant in the vertical, the volume of waters above maximum deepening (i.e., subtropical mode waters) increases and volume of waters below (i.e., subpolar mode waters) decreases.

Our analysis cannot immediately discriminate between dynamically induced adiabatic vertical heaving, changes in lateral advection, and thermodynamic forcing such as diabatic change in water-mass renewal rates due to either warmer temperatures or weaker winter winds. On multidecadal scales water mass renewal is expected to be the primary contributor to heaving, but additional dynamic impacts may arise from very low-frequency changes in ocean (e.g., MOC and gyre circulation). However, heaving is dominating the multidecadal basin-average OHC in every basin with heaving (deepening) at maximum at subtropical mode water densities. This global linkage supports a single source of forcing (i.e., global warming), which can efficiently heat thermocline waters through subduction at high latitudes. These conclusions are consistent with Wang et al. (2015), who show that the future CO2 warming leads to the deepening of midthermocline isopycnals. While Wang et al. focus on future warming, our results from the last half-century show that the global warming effect on the upper ocean structure is already visible as deepening of the subtropical mode water isopycnals.

Acknowledgments

We acknowledge the support of NASA Headquarters by the Physical Oceanography Program (SH, DLW) and its Ocean Surface Topography Mission (SH, PBR, and DLW). We also thank the reviewers for useful and constructive comments.

REFERENCES

REFERENCES
Abraham
,
J. P.
, and Coauthors
,
2013
:
A review of global ocean temperature observations: Implications for ocean heat content estimates and climate change
.
Rev. Geophys.
,
51
,
450
483
, doi:.
Balmaseda
,
M. A.
,
K. E.
Trenberth
, and
E.
Källen
,
2013
:
Distinctive climate signals in reanalysis of global ocean heat content
.
Geophys. Res. Lett.
,
40
,
1754
1759
, doi:.
Bindoff
,
N. L.
, and
T. J.
McDougall
,
1994
:
Diagnosing climate change and ocean ventilation using hydrographic data
.
J. Phys. Oceanogr.
,
24
,
1137
1152
, doi:.
Carton
,
J. A.
, and
B. S.
Giese
,
2008
:
A reanalysis of ocean climate using Simple Ocean Data Assimilation (SODA)
.
Mon. Wea. Rev.
,
136
,
2999
3017
, doi:.
Church
,
J. A.
,
J. S.
Godfrey
,
D. R.
Jackett
, and
T. J.
McDougall
,
1991
:
A model of sea level rise caused by ocean thermal expansion
.
J. Climate
,
4
,
438
456
, doi:.
Curry
,
R.
,
R.
Dickson
, and
I.
Yashayev
,
2003
:
A change in the freshwater balance of the Atlantic Ocean over the past four decades
.
Nature
,
426
,
826
829
, doi:.
Durack
,
P. J.
, and
S. E.
Wijffels
,
2010
:
Fifty-year trends in global ocean salinities and their relationship to broad-scale warming
.
J. Climate
,
23
,
4342
4362
, doi:.
Easterling
,
D. R.
, and
M. F.
Wehner
,
2009
:
Is the climate warming or cooling?
Geophys. Res. Lett.
,
36
,
L08706
, doi:.
England
,
M. H.
, and Coauthors
,
2014
:
Recent intensification of wind-driven circulation in the Pacific and the ongoing warming hiatus
.
Nat. Climate Change
,
4
,
222
227
, doi:.
Fan
,
T.
,
C.
Deser
, and
D. P.
Schneider
,
2014
:
Recent Antarctic sea ice trends in the context of Southern Ocean surface climate variations since 1950
.
Geophys. Res. Lett.
,
41
,
2419
2426
, doi:.
Gille
,
S. T.
,
2008
:
Decadal-scale temperature trends in the Southern Hemisphere ocean
.
J. Climate
,
21
,
4749
4765
, doi:.
Gille
,
S. T.
,
2014
:
Meridional displacement of the Antarctic Circumpolar Current
.
Philos. Trans. Roy. Soc. London
,
372A
, 20130273, doi:.
Good
,
S. A.
,
M. J.
Martin
, and
N. A.
Rayner
,
2013
:
EN4: Quality controlled ocean temperature and salinity profiles and monthly objective analyses with uncertainty estimates
.
J. Geophys. Res. Oceans
,
118
,
6704
6716
, doi:.
Häkkinen
,
S.
,
P. B.
Rhines
, and
D. L.
Worthen
,
2015
:
Heat content variability in the North Atlantic Ocean in ocean reanalyses
.
Geophys. Res. Lett.
,
42
,
2901
2909
, doi:.
Jackett
,
D. R.
, and
T.
McDougall
,
1997
:
A neutral density variable for the world’s oceans
.
J. Phys. Oceanogr.
,
27
,
237
263
, doi:.
Klinger
,
B. A.
, and
C.
Cruz
,
2009
:
Decadal response of global circulation to Southern Ocean zonal wind stress perturbation
.
J. Phys. Oceanogr.
,
39
,
1888
1904
, doi:.
Kosaka
,
Y.
, and
S.-P.
Xie
,
2013
:
Recent global-warming hiatus tied to equatorial Pacific surface cooling
.
Nature
,
501
,
403
408
, doi:.
Leadbetter
,
S. J.
,
R. G.
Williams
,
E. L.
McDonagh
, and
B. A.
King
,
2007
:
A twenty year reversal in water mass trends in the subtropical North Atlantic
.
Geophys. Res. Lett.
,
34
,
L12608
, doi:.
Levitus
,
S.
, and Coauthors
,
2012
:
World ocean heat content and thermosteric sea level change (0–2000 m), 1955–2010
.
Geophys. Res. Lett.
,
39
,
L10603
, doi:.
Lozier
,
M. S.
,
S.
Leadbetter
,
R. G.
Williams
,
V.
Roussenov
,
M. S. C.
Reed
, and
N. J.
Moore
,
2008
:
The spatial pattern and mechanisms of heat content change in the North Atlantic
.
Science
,
319
,
800
803
, doi:.
Lyman
,
J. M.
, and
G. C.
Johnson
,
2008
:
Estimating annual global upper-ocean heat content anomalies despite irregular in situ ocean sampling
.
J. Climate
,
21
,
5629
5641
, doi:.
Lyman
,
J. M.
, and
G. C.
Johnson
,
2014
:
Estimating global ocean heat content changes in the upper 1800 m since 1950 and the influence of climatology choice
.
J. Climate
,
27
,
1945
1957
, doi:.
Marshall
,
G. J.
,
2003
:
Trends in the southern annular mode from observations and reanalyses
.
J. Climate
,
16
,
4134
4143
, doi:.
Marshall
,
J.
,
K. C.
Armour
,
J. R.
Scott
,
Y.
Kostov
,
U.
Hausmann
,
D.
Ferreira
,
T. G.
Shepherd
, and
C. M.
Bitz
,
2014
:
The ocean’s role in polar climate change: Asymmetric Arctic and Antarctic responses to greenhouse gas and ozone forcing
.
Philos. Trans. Roy. Soc. London
,
372A
, 20130040, doi:.
Purkey
,
S. G.
, and
G. C.
Johnson
,
2012
:
Global contraction of Antarctic Bottom Water between the 1980s and 2000s
.
J. Climate
,
25
,
5830
5844
, doi:.
Radko
,
T.
, and
I.
Kamenkovich
,
2011
:
Semi-adiabatic model of the deep stratification and meridional overturning
.
J. Phys. Oceanogr.
,
41
,
757
780
, doi:.
Roemmich
,
D.
,
J.
Church
,
J.
Gilson
,
D.
Montelesan
,
P.
Sutton
, and
S.
Wijffels
,
2015
:
Unabated planetary warming and its ocean structure since 2006
.
Nat. Climate Change
,
5
,
240
245
, doi:.
Stephens
,
G. L.
, and Coauthors
,
2012
:
An update on Earth’s energy balance in light of the latest global observations
.
Nat. Geosci.
,
5
,
691
696
, doi:.
Wang
,
G.
,
S.-P.
Xie
,
R. X.
Huang
, and
C.
Chen
,
2015
:
Robust warming pattern of global subtropical oceans and its mechanism
.
J. Climate
,
28
,
8574
8584
, doi:.
Williams
,
R.
,
V.
Roussenov
,
D.
Smith
, and
M. S.
Lozier
,
2014
:
Decadal evolution of ocean thermal anomalies in the North Atlantic: The effects of Ekman, overturning, and horizontal transport
.
J. Climate
,
27
,
698
719
, doi:.
Wolfe
,
C. L.
, and
P.
Cessi
,
2010
:
What sets the strength of the mid-depth stratification and overturning circulation in eddying ocean models?
J. Phys. Oceanogr.
,
40
,
1520
1538
, doi:.
Zhai
,
X.
, and
L.
Sheldon
,
2012
:
On the North Atlantic Ocean heat content change between 1955–70 and 1980–95
.
J. Climate
,
25
,
3619
3628
, doi:.

Footnotes

Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-15-0607.s1.

Supplemental Material