The Southern Hemisphere oceans absorb most of the excess heat stored in the climate system due to anthropogenic warming. By analyzing future climate projections from a large ensemble of the CMIP5 models under a high emission scenario (RCP8.5), we investigate how the atmospheric forcing and ocean circulation determine heat uptake and redistribution in the Southern Hemisphere oceans. About two-thirds of the net surface heat gain in the high-latitude Southern Ocean is redistributed northward, leading to enhanced and deep-reaching warming at middle latitudes near the boundary between the subtropical gyres and the Antarctic Circumpolar Current. The projected magnitudes of the ocean warming are closely related to the magnitudes of the wind and gyre boundary poleward shifts across the models. For those models with the simulated gyre boundary biased equatorward, the latitude where the projected ocean warming peaks is also located farther equatorward and a larger poleward shift of the gyre boundary is projected. In a theoretical framework, the subsurface ocean changes are explored using three distinctive processes on the temperature–salinity diagram: pure heave, pure warming, and pure freshening. The enhanced middle-latitude warming and the deepening of isopycnals are attributed to the pure heave and pure warming processes, likely related to the wind-driven heat convergence and the accumulation of extra surface heat uptake by the background ocean circulation, respectively. The equatorward and downward subductions of the surface heat and freshwater input at high latitudes (i.e., pure warming and pure freshening processes) result in cooling and freshening spiciness changes on density surfaces within the Subantarctic Mode Water and Antarctic Intermediate Water.
More than 90% of the excess heat stored in the climate system from anthropogenic greenhouse gas emissions is stored in the ocean, leading to global ocean thermal expansion and sea level rise (Church et al. 2011; Rhein et al. 2013). The Southern Hemisphere (SH) oceans are one of the key regions in absorbing and storing the anthropogenic heat. For example, the oceans south of 30°S, while occupying only 30% of the global ocean surface area, account for 75% of ocean surface heat uptake over the historical period of 1861–2005 in the global climate models of phase 5 of the Coupled Model Intercomparison Project (CMIP5) (Frölicher et al. 2015). Recent Argo observations demonstrated that most of the global ocean warming since 2006 is found in the extratropical SH oceans, with the strongest ocean warming centered at middle latitudes (Roemmich et al. 2015). Although historical estimates of the SH ocean warming might be biased low due to sparse coverage of observational samples (Durack et al. 2014), the fast ocean warming at the SH middle latitudes observed since the middle twentieth century is well reproduced by the climate model historical simulations (Cai et al. 2010; Durack et al. 2014) and has been primarily attributed to the anthropogenic greenhouse gas forcing (Swart et al. 2018; Bilbao et al. 2019; Tokarska et al. 2019). The enhanced SH middle-latitude ocean warming is also a robust feature in the twenty-first-century projections (Sen Gupta et al. 2009; Collins et al. 2013; Durack et al. 2014).
It is of great interest to understand the drivers and mechanisms that are involved in the SH ocean heat uptake and storage resulting from the anthropogenic climate change. One major feature of SH climate change is a poleward shift and strengthening of the westerly winds in both observations and climate model simulations (Swart and Fyfe 2012; Collins et al. 2013). Consistent with the changing winds, a strengthening and poleward expansion of the SH wind-driven subtropical ocean gyres has been observed (Roemmich et al. 2016) and is also projected to continue in the future (Sen Gupta et al. 2009; Meijers et al. 2012; Zhang et al. 2014), although a shift of the Antarctic Circumpolar Current (ACC) remains uncertain (e.g., Meijers et al. 2019). Observed and simulated zonal mean structures of the SH ocean temperature change can be largely reproduced through a simple poleward shift of the ocean temperature climatological field (Gille 2008; Cai et al. 2010), implying a linkage between the ocean warming and a poleward shift of the ocean circulation. Diagnosis of climate model outputs indicated that heat is accumulated at SH middle latitudes through lateral advection of warm waters and downwelling, which can be explained by the wind-driven circulation change and resulting heat convergence (Exarchou et al. 2015; Kuhlbrodt et al. 2015). Over the recent Argo period 2005–15, the wind-driven thickening of Subantarctic Mode Water (SAMW) is a key process in the SH ocean warming (Gao et al. 2018). These results suggest that the wind changes could play an important role in driving the SH ocean warming patterns across different time scales.
Cai et al. (2010) pointed out that while the wind forcing drives an adiabatic redistribution of heat, extra surface heat flux input is required to explain the SH ocean warming, suggesting that both the wind and surface heat flux must be at work. Numerical model experiments using different forcing fields and model configurations have been conducted to address the relative importance of different forcing. Fyfe et al. (2007) examined ocean responses to changes in the westerly winds and direct radiative effect of CO2 concentration increases separately in a climate model of intermediate complexity. Their model results confirmed the importance of both forcing and suggested the wind changes are critical to determine the latitude of enhanced warming through increased downwelling. Similarly, recent work by Liu et al. (2018) separated the wind-induced ocean changes and direct CO2 effects in a coupled climate model. They found that most of the ocean warming could be attributed to the direct CO2 effects when the winds are kept unchanged in the model. In contrast, when the CO2 is fixed at its preindustrial level, the wind change itself contributed only 20% of the SH middle-latitude warming. Based on passive-tracer simulations, Armour et al. (2016) suggested that the northward advection of anomalous warming by climatological ocean currents plays a more dominant role than changes in ocean circulation.
The Flux-Anomaly-Forced Model Intercomparison Project (FAFMIP; Gregory et al. 2016) was planned to address these issues by examining the responses of coupled climate models to individual surface flux perturbations derived from multimodel ensemble projections under doubled CO2 concentration. Preliminary results from the FAFMIP generally support the above findings that the uptake of surface heat input as a passive tracer could explain most of the ocean warming under anthropogenic climate change, with smaller contributions from the wind and freshwater flux forcing (Gregory et al. 2016).
While conducting numerical model experiments is an effective way to examine ocean responses to different forcing, the model results may depend on the choice of the model (coupled or ocean only), model parameters (Kuhlbrodt and Gregory 2012; Huber and Zanna 2017), and the way surface forcing fields are modified. Attributing the ocean changes to the change in each individual surface flux may not be achieved easily in these model experiments, as other flux fields may also change in response to the applied change. This may occur through bulk formulas in an ocean-only model (Garuba and Klinger 2018) or in a more complicated way in coupled climate models (Gregory et al. 2016).
While the ocean circulation and different atmospheric forcing should all be at work in the SH ocean heat uptake and redistribution, a better understanding of their individual roles in determining patterns and magnitudes of the ocean warming would help to reduce future projection uncertainties. To highlight the key role of ocean redistribution in shaping the ocean warming patterns, we first quantify the SH ocean surface heat uptake and ocean heat storage in a large ensemble of the CMIP5 models (section 3a). We then link the SH ocean warming structures and magnitudes to the background wind-driven ocean gyre circulation and its poleward shifts (section 3b). To better identify signatures of ocean responses to changes in air–sea fluxes (i.e., heat, freshwater, and momentum), we separate and examine temperature changes along isopycnals and those due to the movement of isopycnals (section 3c). Finally, we attempt to identify connections between subsurface ocean changes and individual surface flux changes using the theoretical framework of Bindoff and McDougall (1994) (section 3d).
2. Data processing of the CMIP5 model outputs
We analyzed historical and future simulations under the representative concentration pathway (RCP) 8.5 from 29 CMIP5 models (Table 1) that also provide the corresponding preindustrial simulations. Only the first realization of each model ensemble is used. The variables used in this study are the three-dimensional ocean potential temperature, salinity, zonal velocity, and two-dimensional fields of the surface wind stress, net heat flux, net freshwater flux, and ocean barotropic streamfunction. The full-depth ocean heat content (OHC) at each ocean grid point was calculated as the vertical integration of ocean temperature multiplied by constant values of density (1035 kg m−3) and specific heat capacity (3985 J K−1 kg−1). Different choices of these constants within a reasonable range results in an uncertainty of only 3% (Hobbs et al. 2016). The ocean barotropic streamfunction, if not provided as model direct output, was calculated by cumulatively summing the vertically integrated zonal velocities from Antarctica northward (Griffies et al. 2016). The calculated streamfunctions were adjusted to ensure that the zero contour passes through the northern limit of the Drake Passage (Sen Gupta et al. 2009; Meijers et al. 2012).
The projected changes under RCP8.5 are defined as differences between the current and future climatologies, calculated as averages over the last two decades of the historical (1986–2005) and the RCP8.5 simulations (2081–2100), respectively. The model drifts have been corrected by subtracting corresponding changes estimated from the parallel preindustrial control runs (Zhang et al. 2014). For a given variable X, at each grid point the de-drifted future change is thus
Compared to the polynomial fitting, this is a relatively convenient method to remove model drifts, especially when dealing with three-dimensional ocean data. However, for the purpose of calculating ocean surface heat uptake, the full time series of de-drifted ocean surface heat flux anomalies over the whole period is required. These were derived by removing a cubic polynomial fit to the corresponding preindustrial simulations. Using a linear fit of the preindustrial simulations to de-drift gives virtually identical results. Note that with this procedure a constant time average from the preindustrial simulations is also removed so the surface heat flux anomalies are relative to the preindustrial climatology. The surface heat flux anomalies were then integrated over time to calculate cumulative ocean heat uptake over 1861–2100. Our calculations of the surface heat uptake over the SH oceans are very close to those reported in recent studies (Frölicher et al. 2015; Shi et al. 2018) despite the use of different groups of CMIP5 models. Results from each individual model were remapped to the same 1° × 1° horizontal grid and depth level for multimodel ensemble average calculations.
a. Identifying the main features of SH ocean heat uptake and storage
We first look at global zonal mean sections of ocean temperature changes from available observations and CMIP5 models to show where heat is stored in the SH oceans (Fig. 1). Both recent Argo observations (2006–17) and an estimate of historical ocean changes based on observations over 1950–2008 (Durack and Wijffels 2010) show enhanced and deep-reaching ocean warming centered at ~45°S (Figs. 1a,b). The CMIP5 multimodel ensemble historical simulations over 1950–2005 show similar change patterns and magnitudes compared to observations over the same period, including slight cooling to the north of the enhanced warming at middle latitudes (Fig. 1c). The consistency between historical observations and climate model simulations indicates that the observed SH ocean warming over the historical period is mainly driven by the external forcing from greenhouse gases and other emissions (Swart et al. 2018; Bilbao et al. 2019; Tokarska et al. 2019). The rate of warming observed over the recent Argo period is considerably larger than changes estimated and simulated over the longer historical period. In addition to possible contributions from the acceleration of anthropogenic warming, the internal climate variability from decadal anomalies in the wind forcing and ocean circulation could temporally reinforce (or offset) the SH ocean warming over such a short period (Roemmich et al. 2016; Gao et al. 2018). Based on future projections from the CMIP5 multimodel ensemble, the enhanced ocean warming centered between 40° and 45°S is projected to continue under RCP8.5, with warming larger than 1°C reaching down to around 800 m by the end of the twenty-first century (Fig. 1d).
Hereafter our analyses will be focused on future projections under RCP8.5, given the larger signal-to-noise ratio between climate change signal and internal climate variability than in historical simulations and climate change mitigation scenarios (Lyu et al. 2015). Corresponding to the enhanced ocean warming, the depth-integrated OHC changes also show the maximum increase at middle latitudes (Figs. 2a,b). While the enhanced ocean warming and the largest OHC increase are located at middle latitudes, the ocean surface heat uptake (i.e., positive and downward surface heat flux into the ocean) mainly occurs in the SH high latitudes with the global zonal mean maximum at ~58°S (Figs. 2c,d), where the upwelled cold waters come into surface and meet with the warmer atmosphere. Spatially integrated from the Antarctic coastline northward to 50°S, the implied OHC increase by the net surface heat flux into the ocean is 8.53 ± 1.27 × 1023 J (multimodel mean and one intermodel standard deviation). However, the projected OHC increase in this area (i.e., the heat stored locally) is only 2.88 ± 0.73 × 1023 J. That is, the relatively slow warming of the high-latitude Southern Ocean results from the residual between surface heat uptake and northward advection across 50°S of about two thirds of this heat uptake, as pointed out by Saenko et al. (2015) and Armour et al. (2016).
The ocean surface heat uptake in the middle-latitude band predominantly occurs in the Indian Ocean sector (Fig. 2c). The net heat loss to the atmosphere (i.e., negative surface heat flux change) is located to the southeast of Australia and in the Atlantic basin sector at about 40°S where a larger sea surface temperature (SST) warming is projected (not shown). The OHC increase poleward of 30°S (7.49 ± 1.08 × 1023 J) accounts for much of the surface heat uptake (11.59 ± 1.79 × 1023 J), with the remaining heat (~36%, 4.10 ± 1.22 × 1023 J) transported northward across 30°S. In the CMIP5 historical simulations, a larger portion (48%) of the surface heat uptake by the SH oceans south of 30°S is transported northward across 30°S (Frölicher et al. 2015). As a result of this northward ocean heat transport, there was little interhemispheric difference in historical ocean heat storage (Irving et al. 2019). In summary, the ocean heat redistribution in the SH oceans determines the regional and global ocean warming patterns.
In addition to the surface heat flux changes, the ocean also responds to changes in the wind field and the surface freshwater flux. The poleward intensifying westerly winds lead to a meridional dipole structure of wind stress curl changes over the SH middle-to-high latitudes (Figs. 3a,b). The surface freshwater flux changes, with the sea ice contribution included, show net freshwater input into the ocean poleward of 45°S and net freshwater loss to atmosphere farther to the north (Figs. 3c,d). In the following, we will examine how ocean circulation and these various atmospheric forcing contribute to the ocean heat redistribution and attempt to identify distinct ocean responses to each individual forcing.
b. Linking enhanced middle-latitude warming to ocean gyre circulation and wind forcing
The depth-integrated OHC change patterns show zonal variations in the latitudes of maximum ocean warming in the SH middle-latitude band (Fig. 2a). Given that the SH ocean warming has been shown to be accompanied by a poleward shift of the subtropical gyres (Cai et al. 2010), here we examine how the SH ocean warming patterns are related to the background ocean gyre circulation and its change. In particular, we identified the boundary between the subtropical gyres and the northern edge of the ACC as the zero contour of the barotropic streamfunction based on the time-averaged streamfunction map over 1986–2005 (Fig. 2a). The largest full-depth OHC increases at middle latitudes are generally located around (e.g., the Indian Ocean) or slightly north of (e.g., the Atlantic Ocean) the climatological gyre boundary. To the southeast of Australia, there is a second warming center farther north that is linked to the increased poleward transport by the poleward shift and strengthening of East Australian Current (EAC) extension (Oliver and Holbrook 2014; Zhang et al. 2014).
To determine the depth-dependent ocean warming structures near this gyre boundary, we first identified the boundary location at each longitude and then averaged the ocean temperature data as a function of the meridional distance from this gyre boundary (Fig. 4). Averaging across all longitudes (or in each ocean basin sector) reveals a narrow band of deep-reaching warming located close to the gyre boundary. For example, the enhanced warming of over 1.8°C extends deeper near the gyre boundary (400–700 m in different basins) than inside the subtropical gyres (only the upper 200 m). Consistent with the depth-integrated OHC change patterns (Fig. 2a), the strongest upper-ocean warming is centered at the boundary in the Indian Ocean (Fig. 4b) but is mainly located north of the boundary in the Pacific (Fig. 4c) and Atlantic (Fig. 4d) Oceans.
The strong spatial correspondence between SH ocean warming patterns and background ocean gyre circulation can also be seen from an intermodel relationship. For those models in which the climatological gyre boundary is located farther north, the largest OHC increase is also projected to occur farther north (Fig. 5). When averaged across all basins, the latitudes of the climatological gyre boundary are correlated (r = 0.57) with the latitudes of projected largest OHC increase in individual models (Fig. 5a). This intermodel correlation is slightly weaker in the Pacific Ocean (r = 0.58) and Atlantic Ocean (r = 0.52) basins than in the Indian Ocean basin (r = 0.70), where the largest OHC increases occur closer to the gyre boundary (Fig. 2a). In short, the varying locations of maximum ocean warming in the models are largely related to the background ocean circulation simulated in the individual models.
We also derived the gyre boundary locations based on ocean velocities averaged over 1986–2005 from the ECMWF Ocean Reanalysis System 4 (ORAS4; Balmaseda et al. 2013). When averaged across all basins, most of the models (23 out of 29) simulate the gyre boundary more equatorward than the ORAS4 reanalysis (Fig. 5a). Using ORAS4 ocean reanalysis as a reference, the CMIP5 multimodel mean location of the gyre boundary displays the smallest bias in the Indian Ocean sector (Fig. 5b), but the largest bias in the Pacific sector among three basin sectors (Fig. 5c). As the ocean gyre circulation is mainly set up by the large-scale wind fields, the equatorward bias of the ocean gyre boundary may be related to the similar equatorward bias that has been found in the location of the SH westerly wind jet (Swart and Fyfe 2012; Bracegirdle et al. 2013). We identified the zero wind stress curl line in the Southern Ocean under current climate 1986–2005 from the individual CMIP5 models to represent the simulated locations of the westerly winds. As expected, we found good intermodel correlations between the climatological latitudes of the zero wind stress curl line and the ocean gyre boundary (Fig. 6). For the same period (1986–2005), most of the CMIP5 models show an equatorward bias of the climatological location of the zero wind stress curl line compared to the Japanese 55-Year Reanalysis (JRA-55; Kobayashi et al. 2015) (Fig. 6), especially in the Pacific basin (Fig. 6c). Therefore, the climatological locations of both the westerly winds and the ocean gyre circulation in the SH are biased equatorward in the climate models relative to current reanalyses. It is likely that the locations where the projected SH ocean warming peaks tend to follow the ocean gyre boundary (Fig. 5) and might also be biased equatorward.
A comparison between the gyre boundary locations under current climate (1986–2005) and future climate (2081–2100) indicates that the boundary between the subtropical gyres and the ACC is projected to be displaced poleward (Fig. 2a). It has been shown that the SH ocean warming patterns are consistent with a poleward shift of the subtropical gyres (Cai et al. 2010). Here we found that the projected magnitudes of ocean warming and gyre boundary poleward shift are well correlated across the models. A larger poleward shift of the gyre boundary is generally accompanied with a larger OHC increase zonally integrated along the historical position of the boundary (Fig. 7). Zhang et al. (2014) has shown that the poleward expansion of the South Pacific subtropical gyre is featured by enhanced ocean warming and large sea level rise through thermal expansion around its poleward edge. Therefore, this intermodel relationship essentially indicates the dominant contribution from thermosteric effect (i.e., ocean temperature change) to the gyre circulation change. In summary, while the locations of SH maximum ocean warming are largely determined by the climatological gyre circulation (Fig. 5), the magnitudes of warming are closely related to the gyre circulation changes (Fig. 7).
The ocean gyre circulation change could be driven by the projected changes in the westerly winds through a Sverdrup-type ocean response. Consistent with wind-driven dynamics, the models with a larger poleward shift of the zero wind stress curl line also tend to have a larger poleward shift of the ocean gyre boundary (Fig. 8). Similarly, Bouttes et al. (2012) also found that the wind could account for most of the intermodel spread in magnitudes of projected sea level change in the Southern Ocean. The enhanced middle-latitude warming could also be explained by the corresponding wind changes (Fig. 3a). The westerly wind anomalies over the SH middle-to-high latitudes imply enhanced Ekman transport to move heat northward across the ACC, where the positive wind stress curl anomalies help to transfer heat downward into the ocean interior through enhanced Ekman convergence and the resultant downwelling. Therefore, the magnitudes of the poleward shift of the zero wind stress curl line and of the ocean warming are also well correlated across the models (Fig. 9). In short, a larger change in the westerly winds tends to induce a larger poleward shift of the ocean gyre boundary and a stronger ocean warming, further supporting the wind-driven mechanisms.
The projected magnitude of the poleward shift in the westerly winds has been found to depend on the current climatological westerly winds position (Kidston and Gerber 2010; Bracegirdle et al. 2013). Specifically, the models having a larger equatorward bias for the climatological westerly winds would project a larger poleward shift of the westerly winds. For the ocean counterpart, here we also found a similar intermodel relationship that a more equatorward location of the climatological ocean gyre boundary is associated with a larger poleward shift of the gyre boundary in future projections (Fig. 10). This relationship between current climatology and future change of the ocean gyre boundary is stronger in the Pacific sector (r = 0.65) than in the Indian (r = 0.44) and Atlantic (r = 0.39) sectors. Similarly, Bracegirdle et al. (2013) also found that the intermodel relationship between climatological locations and future shifts of the westerly winds is strongest over the Pacific. Given that the climatological gyre boundaries simulated in most of the CMIP5 models are biased equatorward compared to the current reanalysis, particularly in the Pacific Ocean (Fig. 10c), the projected poleward shift of the gyre boundary might be overestimated. In the Pacific Ocean, the equatorward bias of the mean fields in the CMIP5 models is also associated with stronger ocean warming in future projections (Fig. 11c). However, this relationship does not hold in the Indian and Atlantic Oceans (Figs. 11b,d). It is worth noting that the model mean state bias is only one of the factors that contribute to the uncertainties in future projections. For example, Chen et al. (2019) suggested that a larger decline of the Atlantic meridional overturning circulation (AMOC) could remotely induce a stronger poleward shift of the SH westerly winds.
c. Decomposing ocean temperature changes into spiciness and heave components
Wind is not the only forcing that drives the ocean and determines the magnitude of regional ocean warming. The wind-driven mechanisms highlighted in the previous section correspond to the redistribution of the ocean’s existing heat through changing winds and ocean circulation. Meanwhile, changes in the surface buoyancy forcing, including the heat and freshwater fluxes, could also force the ocean heat redistribution and dynamical adjustments of the ocean circulation. Although the heat enters the ocean at high latitudes (Fig. 2c), the absorbed heat at the surface is preferentially transported to the north and accumulated in the subduction regions at middle latitudes, where the climatological meridional surface transport is convergent (Morrison et al. 2016). That is to say, even without changes in wind, it is still possible to see the enhanced ocean warming at middle latitudes due to the convergence of heat by the background ocean circulation.
We conducted similar intermodel analyses as in the previous section but in terms of the relations between surface heat flux change and ocean warming. The regional ocean warming magnitudes across the models are significantly correlated with the local surface heat flux changes both poleward of 50°S (r = 0.72) and between 50° and 30°S (r = 0.54). In contrast, the magnitudes of surface heat uptake poleward of 50°S are not significantly correlated with either the magnitudes of the middle-latitude ocean warming between 50° and 30°S (r = 0.30) or the magnitudes of the gyre boundary shift (r = 0.27) across the models. Therefore, although most of the ocean warming is originated from the high-latitude surface heat uptake in the individual models (Fig. 2), the varying magnitudes of the middle-latitude ocean warming across the models are mainly determined by the different magnitudes of local surface heat fluxes and ocean heat transport.
To better reveal the subsurface ocean responses to the wind and buoyancy forcing, here we take a different approach by decomposing the temperature changes at depth levels into two components (Fig. 12a): spiciness change along isopycnals and change due to the movement of isopycnals known as heave (Bindoff and McDougall 1994). The decomposition of subsurface ocean changes into spiciness and heave has been widely applied to historical ocean observations or reanalysis datasets of different time lengths (e.g., Durack and Wijffels 2010; Köhl 2014; Häkkinen et al. 2016; Desbruyères et al. 2017). The spiciness reflects change in water mass properties due to the subduction of surface temperature and salinity anomalies as well as modifications by the interior mixing. We first calculated neutral densities (Jackett and McDougall 1997) based on the temperature and salinity data averaged over 1986–2005 and 2081–2100. Then the spiciness changes were derived as temperature and salinity differences on the same density surface between these two periods. The zonally averaged spiciness changes on density surfaces show enhanced warming toward the upper ocean at isopycnal outcrops (Fig. 13a). To the north of 40°S, cooling spiciness anomalies occur on the density surfaces between 25.4 and 27.4 kg m−3, which generally correspond to the multimodel mean density ranges for the SAMW and the Antarctic Intermediate Water (AAIW) (Downes et al. 2010; Sallée et al. 2013). Although the density structures and ranges for water masses vary across the models, the large-scale features of the spiciness changes seen in the multimodel ensemble mean are robust in the individual models (not shown). The maximum cooling of up to 1°C appears near the density surface of 26.1 kg m−3, which is close to the average density of the SAMW in the CMIP5 model ensemble (Sallée et al. 2013).
When the spiciness changes are remapped onto depth levels using the mean depth of each density surface, the major cooling occurs roughly between 200 and 800 m (Fig. 13b). On the density surface of 26.8 kg m−3 (i.e., roughly the CMIP5 ensemble average density of the AAIW; Sallée et al. 2013), there exist warm anomalies in the outcrop regions between 40° and 50°S and widespread cooling anomalies to the north in all basins (Fig. 13c). The density-compensated nature along the same density surface implies that a warm (cold) spiciness anomaly is accompanied by a salty (fresh) signal (Fig. 13d). Similar cooling and freshening spiciness changes have been observed from historical hydrographic measurements on density surfaces at the SH middle latitudes (e.g., Bindoff and Church 1992; Banks and Bindoff 2003; Aoki et al. 2005; Durack and Wijffels 2010; Helm et al. 2010). The observed cooling and freshening spiciness changes are found at denser surfaces and deeper depths than model simulations (Fig. 13e), as the CMIP5 models simulate the Southern Ocean water masses at lighter densities (Sallée et al. 2013).
It has been suggested that under anthropogenic warming, the subduction of warmed surface water would lead to a cooling and freshening of the subsurface thermocline water along a given density surface, although the waters still become warmer at constant depth (Church et al. 1991). Bindoff and McDougall (1994) further explained that the sign of spiciness change depends on the slope of the temperature–salinity (θ–S) relation relative to the isopycnals on the θ–S diagram, which is equivalent to a stability ratio defined by the vertical gradients of mean temperature and salinity [see Eq. (A3) in the appendix]. Above the salinity minimum of the AAIW, the θ–S relation has a positive slope; that is, the waters in the upper layer are warmer and saltier than waters below. In such a case, the new water properties after either warming or freshening at constant depth level (from point 1 to point 2 on Figs. 12b and 12c) would appear cooler and fresher compared to the original water properties along the same density surface (point 2 compared to point 3 on Figs. 12b and 12c). Therefore, the cooling and freshening spiciness changes north of 40°S are consistent with warming and freshening at depth levels, which could be explained as signatures due to the equatorward subductions of warm and/or fresh waters originated from the surface net heat and freshwater input in the high-latitude Southern Ocean.
The heave component, derived as the residual between the total change at depth level (Fig. 1d) and the spiciness component (Fig. 13b), is dominated by cooling at high latitudes poleward of 60°S and warming to the north (Fig. 14a). The enhanced warming at middle latitudes can be primarily attributed to the heave component related to isopycnal movements rather than spiciness. Zonally averaged around the globe, the isopycnals in the SH oceans are projected to displace both downward and poleward. Due to the slope of the isopycnals, either a poleward shift or a downward shift of the isopycnals could induce a similar feature of isopycnal changes. The temperature changes due to the isopycnal movements can be defined following Köhl (2014) as
in which γ′|z is the density change at depth levels and and are gradient vectors of the mean temperature and density, respectively. Because the ocean’s horizonal scale is much larger than its vertical scale, in Eq. (3) the horizonal gradients are several orders smaller than the vertical gradients and thus the related terms could be omitted. By doing so the temperature change due to the vertical motion of isopycnals (Fig. 14b) can be directly calculated as the density change γ′|z (Fig. 14c) multiplied by the vertical gradient of the mean temperature against density (Fig. 14d) (Durack and Wijffels 2010),
with the assumption that all the density changes are related to isopycnal vertical motions. This assumption holds well, as through Eq. (2) the isopycnal vertical motions can generate similar patterns and magnitudes as the heave-related temperature change (Figs. 14a,b). The isopycnal deepening generally induces warming except cooling at high latitudes where the near-surface waters are cooler than the waters below (see contours in Fig. 1 for the climatological temperature distribution). Therefore, a poleward shift of the isopycnals should be understood as a consequence of isopycnal deepening of the tilted isopycnals.
The wind-driven heat convergence and isopycnal deepening at middle latitudes can explain the poleward shift of the isopycnals and the subtropical gyres. However, the isopycnals would be expected to deepen at middle latitudes but shoal at high latitudes as a result of the meridional dipole structure of wind stress curl anomalies if the heave process is controlled by the wind forcing alone. Rather, the large-scale deepening of isopycnals (Fig. 14) is caused by the large-scale ocean warming due to the added heat at the surface being subducted into the ocean interior and redistributed (Wang et al. 2015; Liu et al. 2018). The vertical heat mixing is another important process to transfer the heat downward and deepen the isopycnals (Wang et al. 2015).
d. Roles of atmospheric forcing in the SH ocean heat redistribution: Insights from a theoretical framework
While the decomposition of temperature change into spiciness and heave components represents a transformation from depth to density coordinates, neither of these two components corresponds to a single process or a specific atmospheric forcing. A theoretical framework proposed by Bindoff and McDougall (1994) associates the subsurface temperature and salinity changes to three processes on the θ–S diagram, namely, pure warming, pure freshening, and pure heave. Here we attempt to attribute the subsurface temperature and salinity changes to these three processes and explore their possible linkages to changes in different types of air–sea fluxes.
The pure warming process corresponds to temperature change without salinity change at a given depth level (point 1 to point 2 on Fig. 12b). In a global ocean model, Lago et al. (2016) showed that applying a surface warming induces much weaker subsurface salinity changes at depth levels than applying a surface salinity change, suggesting that the warming at constant salinity (i.e., pure warming process) is a good approximation of the subsurface ocean response to the surface heat uptake. The pure freshening process is similar to the pure warming process but allows only salinity change at constant depth (point 1 to point 2 on Fig. 12c). The pure freshening process might be related to the surface freshwater flux change, given that the subsurface salinity changes at depth levels are primarily induced by the surface salinity changes (Lago et al. 2016) due to an amplification of the global hydrological cycle (Durack et al. 2012). The pure heave process corresponds to the adiabatic heave of density surfaces without any spiciness change, which could be due to changes in the rate of subduction between two density surfaces, likely driven by the winds, eddies, wave propagation, or the large-scale gyre circulation change (Bindoff and McDougall 1994, 2000).
The mathematical definitions of each process and their relations are shown in the appendix. These mathematical constraints can be put together into the following matrix form (Bindoff and McDougall 1994):
where are the unknown coefficients for three processes and the other quantities can be calculated from the ocean temperature and salinity data. On the left-hand side of Eq. (5), is the mean density, γ′|z is the density change at depth levels, and Rρ is a stability ratio [Eq. (A3) in the appendix]. On the right-hand side of Eq. (5), the six observables (i.e., temperature and salinity changes at depth levels and their spiciness and heave components) are not independent and only provide two pieces of independent information. That is, there is not enough observational information available to constrain this linear system and uniquely solve for the three processes, which means this linear system is underdetermined and there are many possible solutions with different combinations of the three processes. Following Bindoff and McDougall (1994), we solve this linear system using the singular value decomposition (SVD) method, which provides the least squares solution with the minimum length. The SVD is applied at each grid point and depth level in each model and then the multimodel ensemble mean is derived (Fig. 15). The large-scale features seen in the multimodel mean are also common in the individual models (not shown).
The temperature changes at depth levels (Fig. 1d) are attributed to the pure warming (Fig. 15a) and pure heave (Fig. 15b) processes. By definition the pure freshening process does not contribute to temperature change at constant depth. The pure warming process could be understood as the redistribution of added heat into the ocean from the surface heat flux input. The main feature that the heat tends to accumulate at middle latitudes in the pure warming process is consistent with the redistribution of added heat tracers in the climate model simulations, which is largely determined by the background ocean circulation (Armour et al. 2016; Gregory et al. 2016). For the pure heave process, the enhanced warming and isopycnal deepening at middle latitudes seem to be consistent with the local positive wind stress curl changes there (Fig. 3b) and thus enhanced Ekman downwelling. However, the simulated large-scale deepening of the isopycnals is inconsistent with the expected ocean responses to the wind forcing that include the rising of isopycnals at high latitudes with the negative local wind stress curl changes (Fig. 3b). Also, there should be no net heat gain or loss for wind-driven heat redistribution, while the pure heave process shows warming in most of the global ocean. Such disagreement arises because this linear system is solved for each location individually and the constraint of zero global integral for the pure heave process is not explicitly included.
In addition to the pure heave process, both the pure warming and pure freshening processes can also induce the heave of isopycnals as a result of the redistribution of added heat or freshwater at the surface which changes the internal density field. The pure freshening process, despite causing no temperature change at depth levels, still induces the heave of isopycnals via salinity change and contributes to the heave-related temperature change (Fig. 15d). The heave component induced by the pure warming process, featured by isopycnal deepening and large warming at middle latitudes (Fig. 15c), resembles the observed total heave component (Fig. 14a). In contrast, the pure freshening process shows different patterns of isopycnal movements, with the induced warming centered farther northward at ~30°S (Fig. 15d), which implies that the redistribution of added freshwater is not the dominant process to explain the observed isopycnal movements and related temperature changes (Fig. 14a). It thus supports our understanding that both the redistribution of surface heat input and the wind-driven redistribution (i.e., the pure warming and pure heave processes) contribute to the isopycnal deepening and enhanced warming at middle latitudes. However, the exact contributions from these two processes cannot be partitioned in a definite way here because of the underdetermined nature of this linear system [Eq. (5)].
The spiciness changes (Fig. 13b) are separated into contributions from the pure warming (Fig. 15e) and pure freshening (Fig. 15f) processes. The cooling spiciness changes northward of 40°S under the pure warming process (Fig. 15e) can be explained by the equatorward subduction of surface heat uptake and a positive slope of the θ–S relation in this region (Fig. 12b; Church et al. 1991; Bindoff and McDougall 1994). The spiciness changes due to the pure freshening process (Fig. 15f) largely resemble the observed historical salinity changes on density surfaces (Fig. 13e; Durack and Wijffels 2010). The cooling and freshening spiciness signals under the pure freshening process correspond to the subduction of net freshwater input at high latitudes (Fig. 3c). The near-surface warm and salty spiciness anomalies within the subtropical gyres north of 40°S due to the pure freshening process correspond to the net freshwater loss to the atmosphere and the salinification of surface waters (Fig. 3c). Therefore, the spiciness changes associated with the pure warming and pure freshening processes are reasonably consistent with changes in the surface heat and freshwater fluxes, respectively. However, we would like to highlight again that this is an underdetermined problem. The decomposition results shown here only illustrate the contributions from different processes qualitatively and should not be viewed in a quantitative way. Potentially the underdetermined problem might be addressed by adding additional constraints to the system or predefining one process based on the model perturbation experiment.
It should be noted that this approach is kinematic rather than dynamic. In our case it is applied to the mean states of ocean temperature and salinity under current and future climate and cannot replicate the dynamical adjustments in the model. Also, no dynamic constraint has been considered to ensure correspondence between each idealized pure process and atmospheric forcing (Bindoff and McDougall 1994, 2000). For example, the surface heat flux or freshwater flux changes should drive both subsurface temperature and salinity anomalies, although the pure warming or freshening might be still good approximation. Rintoul and England (2002) suggested that the wind can drive the interannual variability of the SAMW (i.e., spiciness signal) through the anomalous advection of temperature and salinity in the Ekman layer. Durack and Wijffels (2010) pointed out that a wind-driven lateral isopycnal migration may lead to surface waters of different properties being subducted along a given isopycnal and thus subsurface spiciness anomalies. Therefore, the pure heave without spiciness signal might be an idealized process. The wind forcing may indirectly generate subsurface spiciness signals rather than a pure heave scenario.
4. Summary and discussion
In this study, we analyzed future climate projections from a large ensemble of the CMIP5 models to investigate features and mechanisms for the heat uptake, redistribution, and storage in the SH oceans under anthropogenic warming. The CMIP5 models project an enhanced warming in the SH middle latitudes roughly between 35° and 50°S during the twenty-first century. The SH ocean warming centered at middle latitudes is not only a robust feature of the CMIP5 historical simulations and future projections but is also a feature observed in the recent decade-long Argo observations and the longer historical observations over several decades (Fig. 1). The surface heat input into the SH oceans mainly occurs in the high-latitude Southern Ocean south of 50°S, with the ocean transporting much of the absorbed heat northward and accumulating the heat at middle latitudes (Fig. 2).
We found that the most deep-reaching warming in the SH oceans is located around the boundary between the subtropical gyres and the ACC (Fig. 4). The linkages between ocean warming and ocean gyre circulation were explored based on their intermodel relationships across the CMIP5 models. The maximum ocean warming is located farther equatorward in models that have the climatological gyre boundary at a more equatorward location. The projected magnitudes of ocean warming are larger in models in which there is a larger poleward shift of the gyre boundary. Both the current climatology and future change of the ocean gyre circulation could be linked to the wind forcing. There is a good correspondence between the westerly winds and ocean gyre boundary in terms of their climatological locations. For future projections, the varying magnitudes of both the ocean warming and gyre boundary poleward shift among models are closely related to the projected magnitudes of the poleward shift of the westerly winds. These findings extend earlier studies on the importance of wind-driven gyre circulation change in the SH ocean warming (Gille 2008; Cai et al. 2010; Zhang et al. 2014) by revealing how this mechanism helps to explain the projection uncertainties across the models.
While the wind forcing is important to drive the ocean gyre circulation change and explain the heat convergence and enhanced warming at middle latitudes, the large-scale ocean warming requires additional heat input into the ocean (Cai et al. 2010). The heat and freshwater added at ocean surface is redistributed, which further modifies the ocean density fields, ocean circulation, and water mass properties. To better reveal ocean responses to the wind and buoyancy forcing, we decomposed the subsurface ocean temperature changes at depth levels into the spiciness and heave components (i.e., changes along isopycnals and changes due to movements of isopycnals, respectively). A major feature of the simulated (and observed) spiciness changes is the cooling and freshening along density surfaces within the SAMW and AAIW, a signature that is consistent with the warming and freshening at depth levels. The enhanced warming at middle latitudes is mainly related to the heave component through the deepening of isopycnals, which also explains the poleward shift of isopycnals as the isopycnals are titled in the SH middle-to-high latitudes.
A theoretical framework by Bindoff and McDougall (1994) was used to interpret the projected ocean temperature changes and their spiciness and heave components. Based on three “pure” processes on the θ–S diagram, this relatively simple framework, as an alternative way from the numerical model solution, provides new insight to understand ocean responses to the different types of air–sea fluxes. The pure warming and pure heave processes, which could be roughly understood as the redistribution of surface heat uptake by the background ocean circulation and the wind-driven redistribution, are the two dominant processes to explain the isopycal deepening and enhanced warming at middle latitudes. The spiciness changes are separated into two parts due to the pure warming and pure freshening processes, driven by the subductions of surface heat and freshwater input in the high-latitude Southern Ocean, respectively. It will be interesting to check in future work the similarities and differences between these idealized three processes and the results of carefully designed model experiments in which only one air–sea flux is changed at a time and the other two fluxes are kept unchanged.
We also found that, in addition to the uncertainties in air–sea fluxes, the model mean state biases are also responsible for part of the future projection uncertainties. For example, the magnitudes of ocean gyre boundary shift in future projections could depend on the current position of the simulated gyre boundary. The models in which the climatological gyre boundary is located more equatorward tend to project a larger poleward shift of the gyre boundary. Given that the climatological gyre boundary simulated in most models is equatorward biased compared to the ocean reanalysis, an important implication would be that the climate model projections might overestimate the poleward shifts of the SH subtropical gyres. In the Pacific Ocean the equatorward bias of the gyre boundary is also accompanied by stronger ocean warming. Therefore, another important outcome from this study is that the future projection uncertainties across the models would be potentially reduced by correcting model mean biases relative to current observations. It would be interesting to examine if the model biases (e.g., the equatorward bias of the westerly winds) will be reduced or not in the next generation of climate models from phase 6 of the Coupled Model Intercomparison Project (CMIP6; Eyring et al. 2016) or in the high-resolution coupled models.
The SH oceans will continue to be a key region in absorbing and storing heat from anthropogenic warming. This study reveals how the ocean circulation and various atmospheric forcing work individually, and together, to determine the structures, magnitudes, and intermodel uncertainties of the projected SH ocean warming. Further understanding and quantifying ocean responses to changes in different air–sea fluxes may help to reduce uncertainties in future projections of the ocean climate change, including global sea level rise due to thermal expansion and regional sea level change patterns (Gregory et al. 2016). While the coarse-resolution CMIP5 models analyzed in this study show partial eddy compensation in the Southern Ocean (Downes and Hogg 2013), it is worth investigating what role eddies play in the ocean heat uptake and redistribution by analyzing high-resolution models that resolve mesoscale eddies (Griffies et al. 2015; Morrison et al. 2016).
This study was supported by the Centre for Southern Hemisphere Oceans Research (CSHOR), jointly funded by the Qingdao National Laboratory for Marine Science and Technology (QNLM, China) and the Commonwealth Scientific and Industrial Research Organisation (CSIRO, Australia), and the Australian Research Council’s Discovery Project funding scheme (project DP190101173). Q.W. was supported by the National Natural Science Foundation of China (91958203). We are grateful to Damien Irving, Annie Foppert, and three anonymous reviewers for their constructive comments. We thank Trevor McDougall for helpful discussions on applying and interpreting the theoretical framework. We acknowledge the climate modeling groups (listed in the Table 1), the World Climate Research Programme’s (WCRP) Working Group on Coupled Modelling (WGCM), and the Global Organization for Earth System Science Portals (GO-ESSP) for producing, coordinating, and making available the CMIP5 model output.
Pure Warming, Pure Freshening, and Pure Heave
In this appendix, we provide some equations that could help to understand the definitions of these three processes. More details about this theoretical framework can be found in Bindoff and McDougall (1994, 2000).
The observed potential temperature and salinity changes on height surfaces (θ′|z and S′|z) can be decomposed into spiciness and heave components,
in which θ′|n and S′|n are spiciness changes along neutral density surfaces, N′ is the neutral density surface height change (up is positive), and and are vertical gradients of mean temperature and salinity. The neutral density definition requires (Jackett and McDougall 1997)
where α and β are thermal expansion and saline contraction coefficients. The stability ratio Rρ is defined by the vertical gradients of mean temperature and salinity,
Therefore, the pure warming process has
Similarly, from the geometry on θ–S space of Fig. 12c, the pure freshening process has
The pure heave process requires no spiciness change on density surfaces and thus