Surface warming–induced global acceleration of upper ocean currents (original) (raw)

Abstract

How the ocean circulation changes in a warming climate is an important but poorly understood problem. Using a global ocean model, we decompose the problem into distinct responses to changes in sea surface temperature, salinity, and wind. Our results show that the surface warming effect, a robust feature of anthropogenic climate change, dominates and accelerates the upper ocean currents in 77% of the global ocean. Specifically, the increased vertical stratification intensifies the upper subtropical gyres and equatorial currents by shoaling these systems, while the differential warming between the Southern Ocean upwelling zone and the region to the north accelerates surface zonal currents in the Southern Ocean. In comparison, the wind stress and surface salinity changes affect regional current systems. Our study points a way forward for investigating ocean circulation change and evaluating the uncertainty.

Anthropogenic surface warming dominates and drives a global acceleration of the upper ocean currents in a warmer climate.

INTRODUCTION

The ocean circulation transports heat and nutrients and is important for climate change and the marine ecosystem response in the face of greenhouse warming (1). A recent study suggested that ocean currents have accelerated since the early 1990s (2), but uncertainties remain large because of changes in the sampling and type of measurements and to internal variability. A follow-up study using energy-conserving state estimates (3) confirmed a weak increasing trend in the surface mean kinetic energy (MKE) but not the acceleration in the deep ocean as originally suggested. Intensified winds were invoked to explain the global increase in total kinetic energy (2), but a global-scale wind increase was neither observed (4) nor simulated by climate models.

The upper ocean circulation is largely wind driven. Naturally, most studies have so far focused on the effect of wind change, but the projected change in the atmospheric circulation suffers from large uncertainties (fig. S1) (5, 6). Surface buoyancy (heat and freshwater) forcing could become important in regulating ocean circulation under global warming (7–9). For example, decreased sea surface salinity (SSS) in the subpolar North Atlantic is projected to slow the Atlantic Meridional Overturning Circulation (AMOC) (10–13). The impacts of increased sea surface temperature (SST) on ocean circulation, on the other hand, receive much less attention. The increased SST dominates the upper ocean heat content change (14, 15) and could intensify the surface subtropical gyre in a warmer climate (9, 16). Using an ocean general circulation model (OGCM), here we systematically evaluate the ocean current responses to the three components of surface forcing—changes in wind stress (∆Wind), SST (∆SST), and SSS (∆SSS). We show that the ∆SST effect dominates and causes a striking global acceleration of surface currents as reported in recent observations (2, 3). However, our results do not support the “deep-reaching” acceleration. Instead, the surface warming–induced upper layer acceleration in subtropical gyres and equatorial regions is accompanied by a deceleration in the lower thermocline.

RESULTS

We force the MIT general circulation model (MITgcm) with the Coupled Model Intercomparison Project Phase 6 (CMIP6) ensemble mean changes in wind, SST, and SSS at years 101–140 following an abrupt four-time CO2 increase (4xCO2) (fig. S1; see Materials and Methods). The ocean model under the full forcing (AllForc) broadly reproduces the CMIP6 results (Fig. 1A and fig. S2A), with a 14% increase in the globally integrated MKE in the upper layer (0 to 200 m) (Fig. 2A). As the response of major current systems is vertically coherent in the upper 400 m except for the equatorial current system with eastward change in the upper 200 m and westward change beneath (Fig. 2, C to H), here we focus on the top 200-m average currents to assess the global upper ocean current response to greenhouse warming. We have forced MITgcm with ∆Wind (∆WindForc), ∆SST (∆SSTForc), and ∆SSS (∆SSSForc) separately to evaluate their relative importance (see Materials and Methods). The results are largely linearly additive, and the sum of single component forcing runs well reproduces the all-forcing run (fig. S3). The validity of MITgcm decomposition is further verified by the Community Earth System Model (CESM) experiments (see details in Materials and Methods).

Fig. 1. OGCM simulated oceanic response to surface forcing.

Fig. 1.

Left panels: Annual mean changes in the ocean currents (vectors, m/s) and steric height (∆SH, color shading, m) (averaged in the top 200 m). Red (blue) arrows indicate that current anomalies strengthen (weaken) climatological current speed. Right panels: The zonal mean of the upper 200-m MKE changes per unit mass (∆MKE, J): (A and B) all forcings (AllForc), (C and D) ∆SST forcing (∆SSTForc), (E and F) ∆SSS forcing (∆SSSForc), and (G and H) ∆Wind forcing (∆WindForc).

Fig. 2. Relative importance and vertical structure of current responses to individual surface forcing.

Fig. 2.

(A) The changes of global upper 200-m total integral of MKE in each experiment. (B) Zonal mean changes in U currents (color shading, m/s) and ocean temperature (contours at intervals of 0.5°C; 2°C contour thickened; values smaller than 2°C were omitted) in ∆SSTForc. Vertical structure of current changes (cm/s) averaged along the section (C) Kuroshio (30°N, 130°E–134°E), (D) Brazil Current (35°S, 48°W–52°W), (E) ACC (zonal current averaged in 48°S–58°S, 180°W–180°E), (F) Gulf Stream (30°N, 76°W–80°W), (G) Agulhas Current (35°S, 26°E–30°E), and (H) equatorial Pacific (zonal current averaged in 2°S–2°N, 140°E–90°W). Black, red, blue, and gray curves denote the outputs from the AllForc, ∆SSTForc, ∆SSSForc, and ∆WindForc, respectively. All the changes are computed as the difference between each experiment and the control run (CTRL).

The resolution for the MITgcm is 1° in the zonal direction and 1/3° in the meridional direction at low and high latitudes, stretching to 1° at mid-latitudes. Unresolved eddy effects are parameterized (see Materials and Methods). Comparison of the eddy-resolving (~0.1°) CESM high resolution (CESM-HR) simulation with the standard resolution (1°) CESM1 Large Ensemble (CESM-LENS) shows a similar spatial pattern of the long-term MKE trend (1950–2050) (fig. S4), indicating that the standard resolution model with eddy parameterizations simulates the MKE response to anthropogenic climate change. In addition, we have carefully assessed the sensitivity of the eddy parameterizations and show that the results presented in this study are qualitatively insensitive to the choice of eddy-mixing coefficients (see details in Materials and Methods). The possible resolution impacts on the current changes, especially in the Southern Ocean, will be further evaluated by comparison with observations.

Surface warming–induced broad acceleration

We first examine the ocean current changes forced by the surface warming (fig. S1A). The model simulates a coherent acceleration of the upper ocean currents in 77% of the world ocean (Fig. 1, C to D). Our results show that the steric height (SH) and current changes in ∆SSTForc resemble those from the AllForc (Fig. 1, A and C) in a large portion of the world ocean, especially in the wind-driven current systems, including the subtropical gyres, the Antarctic Circumpolar Current (ACC), and the equatorial currents. Figure 2A shows that the globally integrated MKE in the upper layer (0 to 200 m) in ∆SSTForc increases by 21%, which is much larger than either the ∆SSS or ∆Wind forcing. Together, these results show that ∆SST plays a dominant role in modulating upper ocean currents globally under greenhouse warming, characterized by a notable intensification of the subtropical gyres, the equatorial currents, and the zonal currents in the Southern Ocean. Note that not all the currents accelerate in ∆SSTForc. For instance, the currents on the north flank of the ACC in the Indian sector, the Indonesia Throughflow (ITF), and the South Equatorial Current (SEC) in the Indian Ocean decelerate (Fig. 1C). The slowdown of the SEC in the upper 200 m is likely attributed to the uneven warming pattern in the Indo-Pacific Ocean: The strong (muted) warming in the northern (southern) tropical Indian Ocean gives rise to a large SH gradient across ~10°S (Fig. 3A) and, hence, an anomalous eastward geostrophic flow over the SEC region.

Fig. 3. The simulated current response to increasing vertical stratification.

Fig. 3.

(A) Upper (0 to 200 m) and (C) lower layer (500 to 1000 m) vertically averaged SH changes (m, color shading) and ocean current changes (m/s, vectors) from ∆SST4KForc; (B) and (D) same as (A) and (C) but for the outputs of first and second layers from the reduced gravity model experiment (STRAT_RG), respectively. (E) and (F) represent ocean temperature changes (°C, color shading; 3°C contour in black) meridionally averaged in 28°N–32°N from ∆SST4KForc and PassiveT, respectively.

Does the spatially varying ∆SST pattern affect the upper ocean current changes? Figure S5 shows that the current change in response to a uniform 4K surface warming (∆SST4KForc) is similar in magnitude and spatial pattern to the results of the ∆SSTForc, implying that the ocean current response is insensitive to the detailed spatial pattern of surface warming. This is consistent with a previous study that revealed the ocean heat uptake (OHU) is insensitive to the spatial pattern of ∆SST (17). To our knowledge, such coherent and worldwide acceleration of upper ocean currents due to surface warming has not been documented before.

To evaluate potential model dependency, we examine an ocean-only Flux-Anomaly-Forced Model Intercomparison Project (oFAFMIP) ensemble that includes four models (see details in Materials and Methods) (18). oFAFMIP aims to isolate the surface heat, freshwater, and momentum flux effects on the ocean response to greenhouse warming. The ocean current response to heat flux change from oFAFMIP resembles our ∆SSTForc results, including the intensification of the Kuroshio, the equatorial currents, and the Southern Ocean zonal flow (fig. S6). The North Atlantic response is notably different, where oFAFMIP heat flux forcing induces anomalous freshwater flux from the melt Arctic sea ice, resulting in an AMOC slowdown and thus a weakened Gulf Stream not seen in our ΔSSTForc run (the decreased SSS due to sea ice melt is included in the ΔSSS forcing in our decomposition; Fig. 1). In both oFAFMIP and ∆SSTForc, the zonal mean MKE in the top 200 m increases at most latitudes with sharp peaks at the equator and in the Southern Ocean (fig. S6, middle panels). The oFAFMIP response is smaller because its forcing is derived from CMIP5 1pctCO2 experiments. The broad agreement among oFAFMIP models and our MITgcm results suggests that the global surface current acceleration in ΔSSTForc is robust.

At nominally 1°, our ocean model does not resolve mesoscale eddies explicitly, so we only focus on the MKE response. In reality, the increased MKE would cause changes in eddy kinetic energy (EKE) as reported in the Southern Ocean (7, 19). Recent observations show a global increase in ocean eddy activity due to global warming, characterized by “eddy-rich regions become richer” (20). Further research is needed to document the detailed EKE response.

Subtropical gyres

Physical mechanisms for the current acceleration in ∆SSTForc differ in the subtropical gyres, the equatorial regions, and the Southern Ocean. The intensification of the subtropical gyres and western boundary currents (WBCs) is evident in all basins. Figure 2 (C, F, and G) shows the baroclinic nature of the WBC changes (except for the Brazil Current): The Kuroshio Current, Gulf Stream, and Agulhas Current accelerate in the upper ~500 m but decelerate beneath. The subtropical gyres are driven by large-scale wind stress curls between the tropical trades and mid-latitude westerlies. The vertically integrated volume transport follows the Sverdrup relation [V=curl(τ)βρ0]. Surface warming and the resultant increase in vertical stratification of the upper layer are robust changes in observations and models (figs. S1, A and B, and S2C). Consider an idealized three-layer ocean with the dense bottom layer at rest. In the first case, the top two layers share the same density, while in the second case, the first layer warms with a lower density. Assume no change in the wind. The density stratification between the top two layers allows the vertical shear to develop and intensifies the currents in the first layer. Since the vertical-integrated volume transport does not change, the flow in the second layer decelerates (see more details in the Supplementary Text). To demonstrate the effect of enhanced stratification, we have performed a sensitivity experiment with a 2.5-layer reduced gravity (RG) model by increasing the RG (_g_1′, _g_2′) from (0.02, 0.01) to (0.032, 0.012) m/s2 to represent strong (weak) warming in the upper (lower) layer (see details in Materials and Methods). The RG experiment (STRAT_RG-CTRL_RG) broadly reproduces the key regional patterns of the SH and baroclinic current changes in the upper/lower layers over the subtropical gyres (Fig. 3, A to D). The SH in the RG model is determined by both vertical stratification terms (_g_1′ and _g_2′) and layer depths (_h_1 and _h_2) (see Eqs. 9 and 10). The SH change pattern is such that it favors anticyclonic (cyclonic) current changes and accelerates (decelerates) the subtropical gyres in the upper (lower) layer in both hemispheres (Fig. 3, B and D). There are some discrepancies in the current change between MITgcm and RG experiment, especially in the lower layer of the North Atlantic Ocean (Fig. 3, C and D) associated with the slowdown of the AMOC (fig. S7, C and D).

The intensified surface subtropical gyre in STRAT_RG is accompanied by the shoaling of the subtropical gyre depth. To illustrate this dynamical constraint, we performed a passive tracer experiment (PassiveT) with the MITgcm (see Materials and Methods), which is the same as the ∆SST4KForc except that the “temperature-like” tracer does not alter the water density. The penetration depth of temperature increase is markedly shallower in ∆SST4KForc than in PassiveT (Fig. 3, E and F) because the enhanced density stratification limits the subtropical gyre to a shallower depth in a dynamic ocean. This dynamical constraint has received little attention in the literature but is important to limit the ocean warming (or heat content increase) in the upper layer.

Equatorial oceans

The equatorial current changes in the Pacific/Atlantic also feature a baroclinic vertical structure, with strong eastward anomalies in the 40 to 200 m but relatively weak westward flow anomalies above and beneath (Fig. 2, B and H). On the equator, the surface warming causes the equatorial undercurrent (EUC) to shoal in the absence of wind stress change, giving rise to strong eastward (westward) flow anomalies above (below) the mean EUC core (Fig. 4, A and B). The intensification of the upper part of the EUC could be linked to the enhanced vertical stratification (21) or the intensified boundary transports (22). Our experiments support the notion that the increased vertical stratification is crucial in regulating the equatorial current system in a warmer climate, consistent with (21). Dynamically, the enhanced surface stratification causes an upward displacement of the main pycnocline and contributes to an increase in the pressure gradient force in the upper ~150 m but a decrease underneath (fig. S8), accelerating (decelerating) the upper (lower) EUC. The equatorial current changes in the Atlantic are similar and not discussed in detail here.

Fig. 4. The equatorial Pacific Ocean and the Southern Ocean response to uniform 4K surface warming.

Fig. 4.

(A) Annual mean density changes (kg/m3, color shading) and EUC core depth (CTRL, dashed line; ∆SST4KForc, solid line) along the equator in terms of the contribution from uniform 4K warming (∆SST4KForc − CTRL). (B) same as (A) but for U changes (m/s, color shading) and U (0.1 m/s contour; CTRL in gray and ∆SST4KForc in black). The zonal mean of temperature changes (°C, color shading), geostrophic current changes (Ug, contours at 0.005 m/s intervals, positive in solid black line and negative in dashed gray line), as well as zonally integrated ocean heat uptake (OHU) changes (blue line) as well as full-depth ocean heat storage (OHS) changes per latitude (red line) obtained from (C) ∆SST4KForc and (D) PassiveT.

Southern Ocean

The Southern Ocean features a circumglobal acceleration of the zonal current velocity along the ACC between 45°S and 60°S that extends to ~2000 m depths (Figs. 1C and 2, B and E), with the upper 2000-m volume transport increasing at 1.2 Sv per decade in ∆SSTForc. In a warming climate, the upwelled cold water in the ACC flows northward in the Ekman layer and absorbs heat from the warming atmosphere. The warmed water subducts as mode and intermediate waters and carries heat into the ocean interior at ~40°S to 50°S, giving rise to muted (enhanced) warming in (on the north flank of) the ACC (Fig. 4C) (23, 24). This differential warming results in anomalous eastward geostrophic flow in 45°S to 60°S (Fig. 4C). The PassiveT experiment broadly reproduces the zonal mean warming pattern and the acceleration in the region, reaffirming the role of advection by the mean currents (Fig. 4D) (14, 25). Compared to PassiveT, the ocean warming in the dynamic model is confined to a shallow layer, and the resultant eastward “geostrophic current” is relatively weak (Fig. 4, C and D), possibly because mesoscale eddies (parameterized here) could partly flatten the steep isopycnals (17, 26).

The uneven warming pattern across the ACC is also evident in observations (Fig. 5E) (24, 27), accompanied by a large north-south sea level gradient between 45°S and 60°S (Fig. 5A) that drives anomalous eastward geostrophic currents and increased MKE in the Southern Ocean (Fig. 5, left panels, and fig. S9A). Here, we estimate that the upper 2000-m transport in 45°S to 60°S has increased by about (1.9 ± 0.6) Sv per decade for 1993–2017. The Southern Ocean zonal flow acceleration since the early 1990s is also apparent in (2). A recent study (28) revisited the Southern Ocean zonal velocity change based on historical hydrography and ocean reanalyses, Argo, and satellite altimetry. Because of marked internal variability, no reliable trends in the ACC were detectable before 2008 [Figure 2 of (28)], consistent with an influential study of Böning et al. (29). The extended observations since then now show a notable Southern Ocean zonal flow acceleration during 1993–2019 coupled with enhanced ocean warming on the northern flank of the ACC.

Fig. 5. Zonal mean sea level, currents, and ocean temperature trend from observations and simulations.

Fig. 5.

Left panels: Observed linear trend (1993–2017) in (A) zonal mean sea level (m/century, blue line, with the global mean trend subtracted) and geostrophic velocity (Ug, m s−1 century−1, red line) from satellite altimetry; (C) the zonal mean of the upper 200-m averaged MKE (10−3 J per century) and (E) zonally averaged U currents (color shading, m s−1 century−1) and ocean temperature (contours at 0.5°C intervals, positive in solid black line and negative in gray line, 0.5°C contour thickened). Right panels same as left panels but from CESM-LENS ensemble mean during 1993–2017. Shading in (B) is ±1 standard deviation (SD) of surface currents from all the ensemble members. The stippled areas denote the U current signals significant at 95% confidence level from the two-tailed t test. The observed ensemble mean current trend is computed from four reanalysis datasets: ECCO v4r4, ORAS4, GECCO3, and SODA3.4.2.

Modeling studies showed that the intensified westerly winds over the Southern Ocean energize mesoscale eddies but result in little changes in ACC transport, a mechanism known as eddy saturation (19, 29, 30). Vigorous eddy activities effectively compensate the momentum input by increased wind stress but cannot fully offset the uneven ocean uptake of anthropogenic heat. The eddy-resolving (0.1°) CESM simulation shows the characteristic ocean warming pattern and the baroclinic acceleration of Southern Ocean zonal currents (fig. S4, right panels) in support of observations. In summary, models and observations consistently show that the Southern Ocean zonal flow accelerates in a warmer climate, not because of wind intensification but because of meridional differential ocean warming. This result is robust among models (fig. S6) and qualitatively insensitive to model resolution (fig. S4) or eddy parameterization coefficients (fig. S10).

Effects of SSS and wind stress changes

Unlike the broad surface current acceleration forced by surface warming, ΔSSS and ΔWind are crucial for several regional current systems. The ∆SSS effect is small on the globally integrated upper ocean MKE, amounting to a ~1% decrease (Fig. 2A), but it causes a sizable deceleration of the Gulf Stream and an acceleration of the Brazil Current (Figs. 1E and 2, D and F). Figure S1C shows that the projected SSS changes are dominated by strong freshening in the North Atlantic-Arctic Ocean. To further isolate the impacts of this regional freshening, we run another experiment where only the SSS anomalies in the North Atlantic-Arctic Ocean are applied (∆SSSNAtlForc; see Materials and Methods). The regional ∆SSS forcing reproduces the deceleration of the Gulf Stream and the acceleration of the Brazil Current (fig. S11A). The strong freshening in the North Atlantic-Arctic Ocean reduces the formation of North Atlantic Deep Water (NADW) and slows down the AMOC (fig. S7, E and F) (11). Observations from the Rapid Climate Change (RAPID) array suggest a weakening trend in the AMOC (31) (~1.3 Sv per decade during 2004–2018), consistent with model results. The anomalous southward currents in the “upper limb” of the MOC (fig. S7F) weaken the Gulf Stream and strengthen the Brazil Current.

The ∆Wind effect is small at middle/high latitudes (Fig. 1, G and H) but important near the equator. Globally, it causes a ~3% decrease in upper ocean MKE (Fig. 2A). Near the equator where the Coriolis effect vanishes, the upper ocean currents are sensitive to wind stress changes (32, 33). In a warmer climate, the Walker circulation slows down (34, 35), with easterly and westerly wind anomalies over the equatorial Indian and Pacific, respectively (fig. S11B). The wind stress changes drive strong downwind surface current changes (Fig. 1G and fig. S11B). These current changes are largely confined in the upper 50 to 100 m and counter the background currents across the three basins (fig. S11B), leading to a net decrease in MKE in the upper 100 m near the equator. In the Southern Ocean, the marked westerly wind intensification between 40°S and 60°S has small impacts on the zonal currents (Fig. 1H) because of the eddy saturation.

DISCUSSION

Wind is the dominant driver of the interannual to decadal variations of upper ocean currents (36) and contributes to the variability of the global meridional overturning circulation (37). Our ocean model experiments show that globally, the surface warming (∆SST) effect accelerates and dominates the upper ocean circulation changes in a warmer climate. In comparison, wind stress and SSS changes are important in regional circulations, e.g., modulating the equatorial currents and decelerating the AMOC, respectively.

Surface warming and the resultant increase in upper ocean density stratification are the robust outcomes of greenhouse forcing (fig. S1). We show that the worldwide acceleration is insensitive to the spatial pattern of surface warming, suggesting that the warming-induced broad acceleration is robust. Unlike the focus of previous studies on wind change, the intensification of zonal currents in the Southern Ocean is associated with differential ocean warming, hence not fully compensated by eddies (fig. S4). Our results also show that the warming-induced acceleration of surface subtropical gyres and equatorial currents is associated with a slowdown at depth, a prediction that needs to be tested against future observations as climate warming amplifies.

How have ocean surface currents changed, and do the observed changes support model results? Figure 5 compares the linear trends in sea level, currents, ocean temperature, and upper 200-m MKE from observations and the CESM-LENS simulation during 1993–2017. Robust upper ocean warming with enhanced vertical stratification has emerged in the world ocean, with similar amplitude and vertical structure between observations and CESM-LENS (Fig. 5). Both observations and CESM-LENS show an intensification of zonal currents in the Southern Ocean and the equatorial Pacific Ocean (Fig. 5 and fig. S9), with greater magnitudes in observations. The spatial pattern of the upper ocean MKE trend (fig. S9A) is similar to that in (2) (their Figure 2A). Observations show that the globally integrated upper 200-m MKE increased by (24 ± 9)% per century for 1993–2017 (significant at 99% confidence level), lending support to the warming-induced acceleration of upper ocean currents. Not all the observed change is due to greenhouse warming. The anthropogenic aerosol forcing may have delayed the greenhouse warming–induced AMOC slowdown (38, 39), and the wind stress forcing is probably important for the surface current change in the Pacific during 1993–2017 when the Pacific Decadal Oscillation experienced a negative phase transition (40).

In the coupled climate system, changes in SST, SSS, and surface wind are not independent but are mutually interactive and coupled with the ocean circulation. As SST, SSS, and surface wind are being routinely measured from space and in situ, the OGCM-based decomposition we used is a useful framework for attribution studies to understand and interpret unfolding climate change. In this framework, the surface warming–induced global acceleration of ocean surface circulation as identified here is a simple baseline, against which one can evaluate other surface forcing effects. Insights gained into the dynamics of ocean circulation change help develop a fuller understanding of the coupled system and its response to radiative perturbations.

MATERIALS AND METHODS

Observations

The observed sea level trend during 1993–2017 is computed from the monthly satellite altimetry data. Ocean current trends (1993–2017) are computed from four reanalysis datasets: ECCO v4r4 (41), ORAS4 (42), GECCO3 (43), and SODA3.4.2 (44). All the current data are spatially smoothed by 2° × 2° to reduce the noise. Observed ocean temperature and salinity trends for 1993–2017 are computed from the above four reanalysis datasets and three independent objectively analyzed products: Ishii (45), EN4 (46), and the Institute of Atmospheric Physics (IAP) (47). Long-term SSS and SST trends during 1960–2019 are obtained from Ishii, EN4, and IAP together with other three SST datasets: Extended Reconstructed Sea Surface Temperature (ERSST) (48), Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST) (49), and Centennial Observation-Based Estimates of SST (COBE SST) (50). Wind stress trends are calculated using four reanalysis datasets: ERA5 (including the back extension) (1960–2019) (51), National Centers for Environmental Prediction (NCEP) reanalysis 1 (1960–2019) (52), 20CRv3 (1960–2015) (53), and ERA20C (1960–2010) (54). The observed AMOC stream function profiles and transports during 2004–2018 are obtained from the RAPID project (55).

CMIP6 simulations

The preindustrial control (piControl) and abruptly quadrupled CO2 (abrupt4xCO2) simulations of 30 models (table S1) from CMIP6 (56) are analyzed. The piControl is a control experiment, with constant radiative forcings that are representative of preindustrial. The abrupt4xCO2 is the standard CMIP6 simulation in which the CO2 concentration is instantaneously quadrupled from the preindustrial value and then held fixed, with all other gases set to the piControl conditions. The differences between abrupt4xCO2 (101–140) and piControl (101–140) are used to explore the ocean and atmosphere response (denoted by Δ) to increasing CO2. All variables are interpolated onto a regular 1° × 1° latitude-longitude grid for easy comparison. The first member (r1i1p1f1) of each model is analyzed.

CESM Large ensemble

We compute the simulated sea level, ocean current, and temperature trend for 1993–2017 from the 40-member CESM-LENS and compare them with observations. The CESM-LENS is a set of fully coupled simulations at nominal 1° standard resolution. Each member is integrated forward under historical (1920–2006) and future (RCP8.5 for 2006–2100) emissions scenarios (57) with different initial conditions for air temperature.

CESM high-resolution simulations

The CESM-HR simulation used in this study has a nominal 0.1° resolution in the ocean and 0.25° in the atmosphere, forced with historical forcings during 1950–2014 and high-emission scenario similar to RCP8.5 for 2016–2050. Mesoscale eddies in the ocean are resolved in this model. We compare the ocean current trend (1950–2050) of CESM-HR to that from the CESM-LENS simulation with a nominal 1° ocean resolution, where ocean eddies are parameterized (58). The comparison thus allows us to assess the possible impacts of ocean model resolution on our results.

Ocean model experiments

We use the MITgcm to quantify the relative roles of global warming–induced SST (ΔSST), SSS (ΔSSS), and wind (ΔWind) changes (fig. S1, left panels) in altering the ocean currents. The configuration is similar to the Estimating the Circulation and Climate of the Ocean Version 4 Release 4 (ECCO v4r4), which has been widely used to investigate the ocean dynamic process and current variations (3, 59). Specifically, there are 50 vertical levels, with layer thickness gradually increasing from 10 m near the surface to about 456 m in the deep ocean. The model resolution is 1° in the zonal direction and 1/3° in the meridional direction at low and high latitudes, stretching to 1° at mid-latitudes, which has limits in resolving mesoscale eddies, especially for the Southern Ocean. The model configuration uses an optimized and spatially varying specification of the Gent-McWilliams (GM)/Redi eddy parameterization scheme (58, 60). Our ocean model uses spatial-varying harmonic and bi-harmonic viscosity for lateral momentum exchange but harmonic and constant value for vertical momentum exchange. In the mixing layer, a so-called GGL90 mixing scheme (61) parametrizes the eddy-induced vertical mixing to replace the constant value mentioned above. We did not parametrize submesoscale in this standard resolution model. The surface forcing fields include 6-hourly winds, specific humidity, downward longwave and shortwave radiation, precipitation, and 2-m air temperature from the ECCO v4r4 (41, 62). The model restarts from an initial state obtained from ECCO v4r4 and is first integrated forward from 1 January 1992 to 31 December 2017. After the integration, the monthly air-sea fluxes including surface net heat flux, surface net short wave flux, and surface net freshwater flux together with SST and SSS are stored as diagnosed data, and the monthly present-day climatological forcing fields—including climatological wind field (WindClim), SST (SSTClim), and SSS (SSSClim)—are computed by using these diagnosed data. With these present-day climatological forcing fields, the model is further integrated for 100 years to reach a quasi-equilibrium state. Note that sea ice is not allowed to change during the integration even with perturbed forcing. The MITgcm is able to reproduce the main current systems in the world ocean such as the subtropical gyre, the zonal currents in the Southern Ocean, the EUC, and the AMOC. For instance, the simulated AMOC strength (heat transport) in our OGCM is 16.31 ± 1.24 Sv (1.10 ± 0.12 PW), which is comparable to the observed values 17.4 ± 1.65 Sv (1.25 ± 0.31 PW). Our OGCM also successfully captured the EUC in the Pacific and Atlantic Oceans, although the simulated EUC speed (~0.40 m/s, meridionally averaged in 2°S to 2°N) is weaker than observations (~0.55 m/s).

To isolate the effects of SST, SSS, and wind stress changes on the upper ocean currents, sensitive experiments were conducted (table S2). In these experiments, the surface forcing fields are perturbed by ΔSST, ΔSSS, and ΔWind, which are computed from the multimodel ensemble mean difference between abrupt4xCO2 (101–140) and piControl (101–140) (fig. S1, left panels). Each experiment is integrated forward for an additional 140 years; The average of the last 40 years (101–140) of each experiment is used in the analysis presented here. The control run (CTRL) is forced by present-day monthly climatological surface forcing, and the SST and SSS are strongly restored to SSTClim and SSSClim on a time scale of 10 days. In AllForc, the model is forced by WindClim + ΔWind, SSTClim + ΔSST, and SSSClim + ΔSSS. Its solution (AllForc- CTRL) contains the complete forcing impacts and is compared with CMIP6 results to evaluate the model performance. In the ∆SSTForc, the forcing fields are the same as CTRL, but SST is restored to the prescribed SSTClim + ΔSST. The difference, ∆SSTForc − CTRL, thus isolates the oceanic response to ΔSST. Likewise, the ∆SSSForc and ∆WindForc isolate the ocean current response due to ∆SSS (∆SSSForc − CTRL) and ∆Wind (∆WindForc − CTRL), respectively. The ocean response to buoyancy flux changes is computed by the sum of ΔSST and ΔSSS effects.

The standard resolution (~1°) model parameterizes mesoscale eddies, which could affect the ocean current response, especially in the Southern Ocean. To assess the possible eddy impacts, we compare our model simulation with the results from an eddy-resolving CESM model and observations. The broad agreement indicates that the eddy parameterization in our ocean model captures the mesoscale eddy effect (figs. S4 and S5). Further, we conducted two additional experiments, which are similar to the ∆SSTForc run but with a 25% increase (∆SSTForc_large) and 25% decrease (∆SSTForc_small) in the coefficients _K_Redi and _K_GM, to assess the sensitivity to eddy parameterizations. A ±25% is a reasonable range appropriate for global warming based on CMIP5 models (17, 63). The broad surface acceleration pattern in the world ocean is insensitive to the choice of eddy-mixing coefficients (fig. S10), although there are quantitative differences (8, 17). We recognize that the three components of the surface forcing are coupled in nature, not independent as in ocean-only models. For instance, surface wind changes affect ∆SST patterns, which could further affect ocean currents. As the broad acceleration pattern is insensitive to the ∆SST spatial pattern, the ocean-only model approach is useful in isolating the underlying physical mechanism and is a necessary step toward a fully coupled view.

To examine the sensitivity of ocean current response to the spatial pattern of ΔSST, we run the ∆SST4KForc. In the ∆SST4KForc, we use the same surface condition as the CTRL, except that a uniform 4K temperature anomaly (the average of the global SST warming during years 101–140 in the abrupt4xCO2 experiments) is added to the target temperature (SSTClim + 4K). The difference, ∆SST4KForc − CTRL, implies the ocean response to a uniform 4K warming. By comparing ∆SST4KForc − CTRL and ∆SSTForc − CTRL, we can assess the impacts of surface warming spatial patterns on the ocean currents. The ∆SSSNAtlForc is similar to the ∆SSSForc but is forced by SSS freshening only in the Northern Atlantic-Arctic Ocean (∆SSSNAtl), which is used to assess the impacts of ∆SSSNAtl. In addition, we introduce a passive tracer experiment (PassiveT) (64) to compare the passive advective and active dynamical effects. The passive tracer has the unit of temperature but does not affect ocean circulation in any way. This temperature-like tracer is initialized with the control temperature distribution (CTRL) and is subjected to the same surface conditions as in the ∆SST4KForc. Together with the ∆SST4KForc run, this tracer simulation helps reveal the dynamical constraints on the redistribution of temperature anomalies.

Ocean-only FAFMIP

We also analyzed the output of the oFAFMIP experiments from four different models: GFDL-MOM5, ACCESS-OM2, HadOM3, and NEMO3.4 (18). Unlike our experiments, the oFAFMIP directly prescribes surface heat, freshwater, and momentum flux perturbations for the ocean. The perturbations are monthly flux anomalies of the ensemble mean of the 61st to 80th years 1pctCO2 experiments from 13 CMIP5 models (65). The faf-heat, faf-stress, and faf-water experiments are carried out to assess the impacts of perturbed heat flux, wind stress, and freshwater flux, respectively. In this study, we compare the results of faf-heat with those from ∆SSTForc. As the experimental design and ocean models are distinct from our OGCM experiments, the comparison thus could test the robustness of our main conclusions. Note that oFAFMIP allows the sea ice to interact with the perturbed heat flux in the faf-heat experiment. Consequently, the strong freshening effects due to sea ice melting in the North Atlantic-Arctic Ocean are largely considered as the result of heat flux changes in oFAFMIP [Figure 3 of (18)]. However, in the MITgcm decomposition, the freshening impacts are primarily attributed to ∆SSS forcing. This difference in the experimental design could introduce inconsistent flow changes (especially for the Gulf Stream and AMOC) between oFAFMIP and MITgcm results.

Coupled general circulation model (CGCM) simulations

To verify the decomposition results from the above OGCM experiments, we use the CESM (66), version 1.0.5, from the National Center for Atmospheric Research (NCAR). The model comprises the Community Atmospheric Model version 5 (CAM5), the Community Land Model version 4 (CLM4), and the Parallel Ocean Program version 2 (POP2). The horizontal resolution of the atmosphere and land is 1.9° latitude × 2.5° longitude with 26 atmospheric layers in the vertical. The resolution of the POP2 is nominal 1° in the horizontal with 60 ocean layers in the vertical. We conducted four experiments (table S3) to quantify the role of buoyancy change and wind stress changes in modifying the ocean currents. The CTRL_C, acting as a reference simulation, is a fully coupled run with greenhouse gas (GHG) forcing corresponding to preindustrial values. The AllForc_C branches from CTRL_C, with CO2 instantaneously quadrupled. The difference, AllForc_C − CTRL_C, represents the overall impacts of quadrupled CO2, which can be used to evaluate the model performance. We then ran a partial coupling experiment (∆BuoyForc_C), in which the atmospheric CO2 level quadrupled but prescribing the surface wind from CTRL_C. We can isolate the ocean responses to buoyancy changes by subtracting CTRL_C from ∆BuoyForc_C (∆BuoyForc_C-CTRL_C). Similarly, in the ∆WindForc_C, we keep the atmospheric CO2 at the preindustrial level and prescribe the surface wind from AllForc_C. The ocean responses to wind change can be obtained by taking the difference between ∆WindForc_C and CTRL_C. All the experiments are integrated forward 90 years, and the last 50 years (41–90) of each experiment are analyzed in our study. More details about these CESM experiments can be found in (67).

Directly comparing the effects of buoyancy and wind changes from CESM and MITgcm experiments provides insight into confirming the validity of the OGCM experiments. The AllForc_C could reproduce the broad horizontal and vertical patterns of the CMIP6-projected current changes (figs. S2, S12, and S13), which resembles the outputs of AllForc (Fig. 1A). Figures S12 and S13 show that the impacts of buoyancy forcing (∆BuoyForc_C) and wind stress (∆WindForc_C) on ocean current from CESM experiments closely resemble those quantified from MITgcm (∆BuoyForc and ∆WindForc) horizontally and vertically. For instance, the buoyancy forcing in the CESM (∆BuoyForc_C) strengthens the Kuroshio, the Brazil Current, the upper equatorial Pacific current, and the ACC, much as in the ∆BuoyForc experiment (figs. S12 and S13, middle panels). In CESM, the wind stress changes contribute to strong eastward flow anomalies in the equatorial regions of the Pacific and the Atlantic Ocean and westward flow at ~10°N and 10°S in the Pacific Ocean, in accordance with the results from MITgcm (figs. S12 and S13, bottom panels). These results together offer confidence in the MITgcm decomposition results.

A 2.5-Layer RG model

To investigate ocean dynamics, we use a 2.5-layer ocean model (68), which represents the ocean into two active layers and an infinitely deep and motionless bottom layer. The momentum and continuity equations are

∂u1∂t−fv1=−(g1′+g2′)∂h1∂x−g2′∂h2∂x−k1(u1−u2)+Am∇2u1+τxρh1 (1)
∂v1∂t+fu1=−(g1′+g2′)∂h1∂y−g2′∂h2∂y−k1(v1−v2)+Am∇2v1+τyρh1 (2)
∂u2∂t−fv2=−g2′∂∂x(h1+h2)+k1(u1−u2)−k2u2+Am∇2u2 (3)
∂v2∂t+fu2=−g2′∂∂y(h1+h2)+k1(v1−v2)−k2v2+Am∇2v2 (4)
∂h1∂t=−∂(h1u1)∂x−∂(h1v1)∂y (5)
∂h2∂t=−∂(h2u2)∂x−∂(h2v2)∂y (6)

where f is the Coriolis parameter; hi, ui, and vi are the layer thickness, zonal velocity, and meridional velocity of the _i_th layer, respectively; _k_1 and _k_2 are the vertical friction coefficients (both set to 1 × 10−9 s−1); _A_m is the momentum mixing coefficient along the isopycnal surface (set to 3 × 104 m2/s); ρ0 is density; and τ is the wind stress. The model domain is 40°S to 40°N and 180°E to 180°W with a horizontal resolution of 1° latitude ×1° longitude, and the monthly climatological wind stress data diagnosed from MITgcm are used to force the model. The model is initialized with _h_1 = 300 m and _h_2 = 400 m and has parameter values g1′=0.02m/s2 and g2′=0.01m/s2. It takes 150 years to spin up the model to reach a quasi-equilibrium state, and the last 50-year outputs are used as the reference state (CTRL_RG).

In a warmer climate, the vertical stratification increases throughout most places of the world ocean. To assess the vertical stratification impacts on the currents, we have carried out a sensitivity experiment (STRAT_RG). As the RG terms (_g_1′ and _g_2′) are directly connected to vertical ocean stratification, we investigate the impacts of enhanced vertical stratification by increasing the RG terms (g1′=0.032 m/s2 and g2′=0.012 m/s2) in the STRAT_RG. The STRAT_RG is integrated forward 150 years, and the last 50-year outputs are analyzed. The ocean responses to enhanced vertical stratification can be obtained by taking the difference between STRAT_RG and CTRL_RG.

The SH in the RG model is calculated from

SH1=(g1′+g2′)h1+g2′h2g0 (9)

where _g_0 is the gravity acceleration.

Momentum budget

To understand the mechanisms governing the ocean circulation response to surface warming near the equatorial regions, it is important to quantify the sources and sinks of momentum to the circulation. The momentum budget equation is (69)

∂u∂t=(−u∂u∂x−v∂u∂y)−w∂u∂z−1ρ∂P∂x+fv+RES (11)

where, from left to right, the terms represent the time rate of change in zonal velocity, the horizontal and vertical advective terms, the zonal pressure gradient force, the Coriolis force, and the residual terms (including horizontal and vertical friction terms).

Heat storage and heat uptake

We define the zonally integrated rate of full-depth ocean heat content at each latitude as ocean heat storage (OHS) (67)

OHS=∂∂t∬−H0ρ0CPθ∆ydzdx (12)

where C_P is the specific heat of seawater (4007 J kg−1 K−1), θ denotes seawater potential temperature, ∆_y is the meridional distance (∆y = 1 m), and −H is the total depth of seawater.

The zonally integrated OHU per latitude is calculated from

where SHF is the surface net heat flux, which is the sum of shortwave fluxes, longwave fluxes, sensible fluxes, and latent heat fluxes. ∆y here is set to 1 m.

SH and geostrophic current

We also used the MITgcm outputs and CMIP6 results to compute the SH and geostrophic current as follows

Ug(p)=Ug(p0)−gf∂SH∂y (15)
Vg(p)=Vg(p0)+gf∂SH∂X (16)

where SH is the steric height at pressure p relative to reference pressure _p_0. As temperature and salinity changes are largely confined to the upper 2000 dbar and the velocity changes at 2000 dbar are much smaller than that in the upper layer, here we selected a widely used _p_0 = 2000 dbar as the level of “no motion.” The computed Ug(p) and Vg(p) are the relative geostrophic velocities at p relative to the reference level of no motion. The practical choice of _p_0 at 2000 dbar has little impact on our main conclusions. In addition, we compute the passive “SH” (SH_passive) and “geostrophic velocity” (Ug_passive, and Vg_passive) by using the temperature-like tracer from PassiveT to address the possible dynamic effects.

Mean kinetic energy

The MKE per unit mass is calculated from

Acknowledgments

We thank H. Zhang, J. Marshall, Q. Wang, J. Xiao, and B. Qiu for the helpful discussions. We also thank the two anonymous reviewers for the insightful comments.

Funding: Q.P. is supported by the National Natural Science Foundation of China (42005035), the Science and Technology Planning Project of Guangzhou (202102020935), and the Independent Research Project Program of State Key Laboratory of Tropical Oceanography (LTOZZ2102). D.W. is supported by the National Natural Science Foundation of China (92158204), and the Innovation Group Project of Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai) (311020004). S.-P.X. is supported by the National Science Foundation (AGS-1934392). Y.S. is supported by the National Key Research and Development Program of China (2016YFC1401702). G.C. is supported by National Natural Science Foundation of China (41822602). The numerical simulation is supported by the High-Performance Computing Division and HPC managers of W. Zhou and D. Sui in the South China Sea Institute of Oceanology.

Author contributions: Q.P., S.-P.X., and D.W. conceived the study. Q.P. and D.W. performed the MITgcm experiments and the analysis. Q.P. and S.-P.X. wrote the paper with input from all the authors.

Competing interests: The authors declare that they have no competing interests.

Data and materials availability: The data are available in the following links: the SODA3.4.2 (www2.atmos.umd.edu/~ocean/index_files/soda3_readme.htm), ORAS4 (www.cen.uni-hamburg.de/en/icdc/data/ocean/easy-init-ocean/ecmwf-ocean-reanalysis-system-4-oras4.html/), GECCO3 (https://icdc.cen.uni-hamburg.de/en/gecco3.html), ECCO v4r4 (www.ecco-group.org/products-ECCO-V4r4.htm), the Ishii dataset (https://climate.mri-jma.go.jp/pub/ocean/ts/v7.3/), EN4 (www.metoffice.gov.uk/hadobs/en4/), IAP (www.ocean.iap.ac.cn/pages/dataService/dataService.html?navAnchor=dataService), ERSST (psl.noaa.gov/data/gridded/data.noaa.ersst.v5.html), HadISST (www.metoffice.gov.uk/hadobs/hadisst/data/download.html), COBE SST (https://psl.noaa.gov/data/gridded/data.cobe.html), satellite altimetry data (https://resources.marine.copernicus.eu/products), ERA5 (www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5), NCEP (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html), 20CR (https://psl.noaa.gov/data/gridded/data.20thC_ReanV3.html), ERA20C (www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era-20c), the CMIP6 CESM-HR (https://esgf-node.llnl.gov/search/cmip6/), the CESM-LENS simulations (www.earthsystemgrid.org/), and the oFAFMIP (https://gws-access.ceda.ac.uk/public/ukfafmip/). All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Supplementary Materials

This PDF file includes:

Supplementary Text

Tables S1 to S3

Figs. S1 to S13

REFERENCES AND NOTES

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Supplementary Text

Tables S1 to S3

Figs. S1 to S13