Surface heating steers planetary-scale ocean circulation

The large-scale ocean circulation derives its energy from a variety of sources including wind stress (Wunsch and Ferrari, 2004; Hughes and Wilson, 2008; Jamet et al., 2021), tidal forces (Oka and Niwa, 2013), and surface and geothermal buoyancy fluxes (Hughes et al., 2009; Hogg and Gayen, 2020). Gyres are fundamental elements of the large-scale circulation and play a crucial role in the biogeochemical and hydrological cycles in the ocean by transporting momentum, heat, nutrients and chemicals within and across ocean basins (Webb, 2017). In particular, gyres contribute to global heat transport by transferring heat poleward (Palter, 2015; Zhang et al., 2021; Li et al., 2022). Despite their importance in regulating large-scale weather and climate patterns, the interplay between the processes leading to the formation and evolution of ocean gyres are not fully understood.

Traditional oceanographic literature on the drivers of ocean circulation suggest that horizontal flows are primarily caused by mechanical forcing due to wind stress (Sverdrup, 1947) and tidal forces (Wunsch and Ferrari, 2004; Oka and Niwa, 2013), and the meridional overturning circulation (MOC) is chiefly driven by surface buoyancy fluxes (Stommel and Arons, 1959). However, this simplified viewpoint has been amended with time. Recent literature points to a definitive role played by wind stress on the MOC through isopycnal upwelling in the Southern Ocean (Abernathey and Ferreira, 2015; Hogg et al., 2017) and bottom-enhanced diapycnal mixing (Stanley and Saenko, 2014; Drake et al., 2020). Similarly, buoyancy forcing is thought to exert a significant control on horizontal circulation through the conversion from potential to kinetic energy (Tailleux, 2009; Hughes et al., 2009) and control of the stratification in the ocean (Shi et al., 2020), consistent with the fundamental dynamics of rotating horizontal convection (Gayen and Griffiths, 2021). In this paper, we aim to evaluate the interplay between wind stress and surface buoyancy forcing in driving various features of the large-scale ocean circulation, with an emphasis on basin-scale ocean gyres.

Considerable progress has been made to elucidate the impact of surface wind stress on ocean gyres (Sverdrup, 1947; Stommel, 1948; Munk, 1950; Rhines and Young, 1982; Luyten et al., 1983; Pedlosky, 1986). Sverdrup (1947) proposed what became a namesake relationship that links the curl of the wind stress and the depth-integrated meridional geostrophic transport, valid in the ocean interior away from coastlines,

where $V=\rho_{0}\int\upsilon\,\mathrm{d}z$ is the time-mean depth-integrated meridional mass transport with $\rho_{0}$ the ocean’s reference density, $\hat{\bm{z}}$ is the unit vector in the vertical, $\bm{\tau}=\tau_{x}\,\hat{\bm{x}}+\tau_{y}\,\hat{\bm{y}}$ the time-mean horizontal wind stress at the ocean’s surface, and $\beta=\partial f/\partial y$ is the meridional gradient of the Coriolis frequency $f$. The return flow occurs as an intense inertial western boundary flow (see, e.g., Hughes and Cuevas (2001)). The Sverdrup relation (1) describes the dependence of the horizontal structure of barotropic gyres on the wind stress curl, and to date, remains the corner-stone theory of wind-driven gyres.

The Sverdrup relation focuses on understanding the horizontal structure of the vertically-integrated gyre transport. The ventilated thermocline theory, developed by Luyten et al. (1983), made further strides in interpreting the vertical structure of ocean gyres forced by wind stress. They obtained a layer-wise meridional transport for gyres but restricted gyre transport to only ventilated isopycnals; that is, isopycnals outcropping to the surface of the ocean, with non-ventilated isopycnals at rest. However, several studies point to a net re-circulatory transport in the non-ventilated isopycnal layers (Rhines and Young, 1982; McDowell et al., 1982) so long as they do not interact with the ocean’s surface or topography. Therefore, ocean gyres can be viewed as a combination of ventilation and recirculation regimes, and are controlled by surface wind stress curl.

Wind-driven theories encapsulate gyre circulation to zeroth order, however, observational studies show deviations from the Sverdrup relation. Gray and Riser (2014) examined the validity of Sverdrup dynamics in a point-wise manner using observations from Argo floats (Roemmich et al., 2004) and found that it agrees well with observations in the interior subtropical gyres, with significant deviations in subpolar regions. Colin De Verdière and Ollitrault (2016) using data from Argo (Ollitrault and Rannou, 2013) and the World Ocean Atlas (Locarnini et al., 2010) computed the global Sverdrup streamfunction and found that the strength of the gyres was underrepresented by a factor of 2. These discrepancies indicate that processes other than Sverdrup dynamics could be at play in establishing the gyre strength, such as bottom pressure torques (Hughes and Cuevas, 2001), diapycnal mixing (Lavergne et al., 2021), buoyancy forcing (Hogg and Gayen, 2020; Liu et al., 2022), and coupling with the meridional overturning (Klockmann et al., 2020; Berglund et al., 2022). A unified theory encapsulating all the previously stated mechanisms is lacking and this limits our understanding of the coupling between these processes, as well as their combined effects on ocean circulation. The present paper makes an attempt to isolate and understand the roles of wind stress and surface buoyancy forcing in shaping the planetary-scale ocean gyres.

Surface buoyancy forcing alters the ocean’s density structure through heat and freshwater fluxes (Large and Yeager, 2009; Talley et al., 2011). Together with the thermal wind relation,

where $b(\bm{x},t)=g(1-\rho(\bm{x},t)/\rho_{0})$ is the buoyancy, with $g$ the gravitational acceleration, these density structure changes can be used to understand how changes in surface buoyancy forcing might impact ocean circulation. However, the relationship between surface buoyancy forcing and the horizontal circulation is complicated. For example, variations in mixed layer depth and non-linear feedbacks with the ocean circulation both influence how the surface buoyancy forcing is “felt” within the ocean. The mixed layer ingests a fraction of the surface buoyancy flux, which is reflected in the anomalous ocean’s density within the layer. The amount of buoyancy forcing reaching the layers below is thus inversely related to the mixed layer depth, with a deeper mixed layer taking a longer time to relay the excess buoyancy forcing into the subsurface layers (Xie et al., 2010) due to its higher effective heat capacity. The effect of surface heat flux on the near-surface buoyancy also depends on the thermal expansion coefficient of seawater, which varies spatially by an order of magnitude across the globe and is important in setting the ocean’s stratification (Caneill et al., 2022). Furthermore, the circulation modifies the influence of surface buoyancy forcing on the ocean’s buoyancy structure through heat advection (Bryden et al., 1991). Advection acts to alter the buoyancy structure remote from the forcing, which, in view of (2), would also cause anomalies in the ocean circulation in that remote location. In this paper, we evaluate the variability in ocean’s stratification over time to better understand the non-linear and non-local connection between the surface buoyancy forcing and the gyres.

Past studies examined the role of surface buoyancy forcing in restructuring ocean gyres. Goldsbrough (1933) observed that freshwater fluxes can drive horizontal circulation via induced sea surface height anomalies. Luyten et al. (1985) and Pedlosky (1986) extended the ventilated thermocline theory (Luyten et al., 1983) to include an interfacial mass flux between various isopycnal layers (to represent surface buoyancy forcing), and using a simplified ocean model, demonstrated a geostrophic baroclinic flow induced by the buoyancy flux and steered by wind stress. Colin de Verdière (1989) coupled the surface buoyancy forcing to wind stress via a bulk formula and concluded that the former drives a baroclinic mode to significantly recast the horizontal and vertical structure of subtropical gyres. Model studies that forced the ocean using only surface buoyancy fluxes observed single or double gyre circulation; see for example Winton (1997); LaCasce (2004); Nilsson et al. (2005). Gjermundsen et al. (2018) also conducted simulations forced only via a meridionally-varying surface temperature restoring profile and observed a broad eastward zonal flow and a western boundary current. Hogg and Gayen (2020) conducted a series of numerical simulations with a restoring temperature profile and no wind stress in two configurations: a direct numerical simulation in a 3D box domain, and a layered general circulation model in a sector configuration, and found a double (resembling a subtropical and subpolar) gyre in both scenarios. With a multitude of contrasting viewpoints on the processes leading to the formation of gyres, we are motivated to examine further the role of surface buoyancy forcing in driving ocean gyres.

It is worth emphasizing here that ocean gyres do not exist in isolation – they interact with other fundamental aspects of ocean circulation, such as the MOC. The MOC complements gyre circulation in transporting tracers across ocean basins. It has been argued that thermohaline forcing is responsible for driving the MOC (Stommel and Arons, 1959; Stommel, 1961), with gyres being driven by wind stress (Sverdrup, 1947). Luyten et al. (1985) contradicted this simplified view by establishing a link between wind stress and heat gain for the North Atlantic subpolar gyre. Yeager and Danabasoglu (2014) and Yeager (2015) used a coupled ocean–sea ice configuration of the Community Earth System Model with perturbed forcing, and found that most of the decadal variability in both the Atlantic MOC (AMOC) and the North Atlantic subpolar gyre was due to variability in surface buoyancy forcing, while interannual variability in these circulatory features were attributed to wind stress anomalies. Modeling studies (Yeager, 2015; Larson et al., 2020) identify a direct relationship between the mid-depth overturning cell and the North Atlantic subtropical gyre transport on the basis that the two circulatory features are linked through the northward flowing Gulf Stream. Thus, ocean gyres and the MOC are coupled dynamical features (Klockmann et al., 2020; Berglund et al., 2022), and a thorough analysis of the drivers behind the formation of ocean gyres requires a quantitative understanding of other processes of ocean circulation.

The objectives of our paper are: (i) to quantify how the surface buoyancy forcing affects the structure of ocean gyres, and (ii) to understand how wind stress and surface buoyancy fluxes act in concert to drive large-scale ocean circulation. Section 2 outlines the simulation setup and a gyre metric used to analyze the ocean’s circulation. Section 3 examines the sensitivity of gyre circulation to changes in surface wind stress, followed by a brief discussion of the coupling between gyres and other large-scale circulation features in the ocean. Section 4 further investigates the role of surface buoyancy forcing gradients in steering the ocean circulation. In section 5 we look at a uniform warming experiment to illustrate that changes in ocean circulation can also occur in the absence of a meridional gradient in buoyancy forcing anomaly. In section 6, we conclude by emphasizing the connected roles of wind stress and surface buoyancy fluxes in driving ocean circulation, the importance of surface buoyancy forcing in driving ocean gyres, and future directions to advance our understanding of ocean circulation.

In this section, we investigate two perturbation experiments wherein we change the magnitude of wind stress by 0.5 and 1.5 times the control experiment, which alters the time-mean vorticity input due to the wind stress curl by the same factor. We analyze short-term and long-term variations in the subtropical gyre transport for Atlantic and Pacific Oceans, along with a brief discussion of their coupling with the AMOC. Since our flux-forced experiments do not incorporate sea-ice dynamics, changes in Weddell and Ross gyres are not reported in the paper, as sea-ice can significantly alter gyre dynamics in polar regions.

The dashed lines in the time series in Fig. 2 show the expected transport as predicted by the Sverdrup linear scaling of the average gyre transport in the control experiment for the last $100$ years of the simulation. Subtropical gyres follow Sverdrup scaling to a large extent, consistent with the ventilated thermocline theory (Luyten et al., 1983); deviations are observed in the North Atlantic subtropical gyre in both wind perturbation simulations, and in the 0.5$\times$W simulation in the North Pacific subtropical gyre. These deviations from wind stress curl scaling are intriguing and may imply that the gyre is also driven by other mechanisms, including, but not limited to surface buoyancy fluxes; this is discussed in detail in the next two sections.

The gyre strengths adjust quickly to changes in wind forcing (solid curves in the time series in Fig. 2). This is likely due to a quick adjustment of fast-traveling waves (e.g., large-scale barotropic Rossby waves (Anderson and Gill, 1975) or nearly-barotropic topographic waves (LaCasce, 2017)). Subsequent smaller-magnitude adjustments in the gyre transports may be attributed to slower wave motions (e.g., baroclinic Rossby waves or surface modes (LaCasce, 2017)). The time series also show that after the initial response, gyres in the Pacific Ocean are more stable than in the Atlantic Ocean.

The AMOC (evaluated by integrating horizontal transport between $\sigma_{2}\in[1035.5,1038]\;\mathrm{kg}\,\mathrm{m}^{-3}$ density classes) is only marginally impacted by changes in wind stress. For the first 10 years, the AMOC strength is inversely related to the wind stress magnitude (Fig. 3a), which is likely linked to a change in northward Gulf Stream barotropic transport (as is also observed by Hazeleger and Drijfhout (2006) and Yang et al. (2016)). After the initial transient response, we observe a slight decline in AMOC strength for the 0.5$\times$W experiment compared with the control, but no robust amplification in the 1.5$\times$W experiment, which could be associated with partial eddy compensation in the Southern Ocean (Morrison and Hogg, 2013), or other non-linear feedback in the ocean circulation.

We observe a near-perfect linear scaling of the globally integrated kinetic energy with the wind stress magnitude (Fig. 3b). Winds supply energy to the ocean through generation of eddies and enhanced circulation, which leads to an increase in the global kinetic energy (Wunsch and Ferrari, 2004). We discern that the change in kinetic energy due to wind stress is majorly due to a change in the gyre circulation (Fig. 2) as well as mesoscale eddies. A complete energy budget calculation is beyond the scope of the study.

In the previous section, we showed that variations in surface heat fluxes, with opposite signed anomalies applied over subpolar and subtropical regions, can produce anomalies in the ocean circulation. In this section we ask; is the same true if the surface heat flux anomalies are globally uniform? We expect changes in the ocean circulation due to a spatially uniform surface heat flux due to several processes. Firstly, lateral variations in mixed layer depth imply that buoyancy anomalies induced by the uniform heating will be non-uniform. In addition, changes in circulation continuously alter the buoyancy structure of the ocean through advection. Finally, a spatially varying thermal expansion coefficient implies that a constant surface heat flux manifests as a non-uniform surface buoyancy flux, which can cause non-uniform buoyancy anomalies in the ocean. To understand the combined effects of (i) spatial variations in the mixed layer depth and the thermal expansion coefficient, and (ii) advective feedbacks on the ocean circulation, we analyze a uniform warming experiment, where a globally constant heat flux of $+5\;\mathrm{W}\,\mathrm{m}^{-2}$ is applied at the ocean’s surface.

Changes in the strength of each gyre do occur under the uniform warming perturbation (Fig. 12a-c). Focusing first on the North Pacific basin (Fig. 12b), we observe $\sim$22% intensification of the subtropical gyre, consistent with the results of Sakamoto et al. (2005) and Chen et al. (2019). This increase can be attributed to spatial variations in mixed layer depth near the western boundary region.The mixed layer traps the excess heat received from the ocean’s surface and distributes only a fraction of this heat to the layer below. Furthermore, a deeper mixed layer has a higher heat capacity due to its ability to store more heat. Deep mixed layers at the western boundary of the North Pacific subtropical gyre moderately shield the region from developing stratification. Conversely, shallower mixed layers to the east of the western boundary lead to more stratification in that region. The spatially uneven growth of stratification in the subtropical gyre strengthens zonal buoyancy gradients near the western boundary region, which in view of the thermal wind relation, intensifies the meridional gyre flow.

We record a strong ($\sim$50%) intensification of the South Atlantic subtropical gyre (Fig. 12c) due to a similar laterally-varying mixed layer depth as observed in the North Pacific basin, which augments the zonal buoyancy gradients near the western boundary in response to surface heating. The South Pacific gyre anomalies show minor oscillations (with an amplitude of $\sim 1.5\,\mathrm{Sv}$) due to uniform surface heating (not shown).

We observe no significant anomalies in the North Atlantic subtropical gyre strength in the first 20 years of the uniform warming simulation (Fig. 12a), followed by a systematic slowdown over the next 75 years. Several modeling studies have reported a correlation between the North Atlantic subtropical gyre strength and the AMOC (Yeager, 2015; Larson et al., 2020), and we observe a reduction in the AMOC over the same time period (Fig. 12d). This reduction is explained by two processes: (i) in the first 20 years, the Gulf Stream transports slightly larger volume of warm water poleward and (ii) the uniform warming applied at the ocean’s surface promotes the generation of lighter waters in the subpolar regions at the ocean’s surface. These two processes limit the North Atlantic deep water formation, causing an AMOC slowdown (Lohmann et al., 2008; Cheng et al., 2013).

Finally, the globally integrated kinetic energy increases by almost $50\%$ over the full experiment (Fig. 12e). We can ascribe the resulting kinetic energy increase to mean circulation (for example, the gyres) as well as mesoscale eddies, and is consistent with the energy conversion argument put forth by Hughes et al. (2009) that surface buoyancy forcing could induce kinetic energy in the system through a conversion from available potential energy.

We conducted a series of perturbed forcing simulations using a partially eddy-resolving ocean model (at 0.25^∘ lateral resolution) to understand the importance of wind stress and surface buoyancy forcing in steering planetary-scale ocean circulation. Our perturbation experiments (listed in Table 1) are forced using surface boundary fluxes (and are thus called “flux-forced simulations”) to separate the contribution of winds and surface buoyancy in driving the circulation, and are classified into three categories: (i) wind perturbation experiments, (ii) surface buoyancy flux contrast perturbations, and (iii) a spatially uniform warming perturbation.

The flux-forced simulations illustrate that both wind stress (Fig. 2) and surface buoyancy forcing (Figs. 4, 8, and 12a-c) are crucial in shaping the planetary-scale subtropical gyres. We find that perturbations in surface buoyancy flux gradients modify the ocean’s buoyancy structure (Fig. 5, 6, 9, and 10), and thus the circulation through the thermal wind relation (2). In addition, spatially uniform surface heat fluxes still induce anomalies in horizontal buoyancy gradients (and hence, the circulation; Fig. 12) due to lateral differences in mixed layer depth and thermal expansion coefficient, and heat advection by the circulation.

On short timescales ($<1$ decade), the anomalies in horizontal density gradients are proportional to the surface buoyancy flux gradient perturbations (compare panels a and d in Figs. 5, 6, 9, and 10). Through the thermal wind relation (2), we diagnose a linear relationship ($\mathcal{R}^{2}>0.9$ in Table 2) between the anomalous gyre circulation and the magnitude of the surface buoyancy flux gradient perturbation on short timescales. Over this period, the Atlantic gyres are observed to be 2-4 times more susceptible to changes in surface buoyancy flux gradients than the Pacific gyres, with as much as a $0.15\,\mathrm{Sv}$ anomaly per $\mathrm{W}\,\mathrm{m}^{-2}$ change in the subtropical/subpolar surface heat flux in the North Atlantic subtropical gyre.

Over time, the induced lateral buoyancy gradients become less proportional to the surface buoyancy flux perturbations (compare panels c and f in Figs. 5, 6, 9, and  10). This divergence from a linear relationship over longer times scales could be attributed to several factors including lateral variations in mixed layer depth and thermal expansion coefficient, and advective feedbacks between the surface buoyancy forcing anomalies and the ocean circulation. For example, in the $-15\,\mathrm{W}\,\mathrm{m}^{-2}$ simulation, this non-linear connection can be observed in the spin-up of the North Atlantic subtropical gyre in the last 20 years (Fig. 4a) and a surge in the AMOC in the last 40 years (Fig. 7a). The ocean circulation in the increased buoyancy flux contrast simulations is more non-linearly related to surface buoyancy flux gradient perturbation than the reduced buoyancy flux contrast simulations. For example, we observe a slowdown in North Atlantic, South Atlantic, and South Pacific gyre strength in the last 40 years of the $+15\,\mathrm{W}\,\mathrm{m}^{-2}$ simulation (Fig. 8) and an oscillatory AMOC in both $+15\,\mathrm{W}\,\mathrm{m}^{-2}$ and $+30\,\mathrm{W}\,\mathrm{m}^{-2}$ simulations (Fig. 11a).

The flux-forced simulations allowed us to conduct perturbation simulations in which each surface forcing could be altered independently. However, the present study has numerous caveats. Firstly, in reality the wind and buoyancy forcing are strongly coupled. In earlier experiments that used bulk formula for the surface heat fluxes (not shown), we found that decreasing the wind forcing strongly reduced the surface buoyancy fluxes as well due to the reduction in poleward heat transport. Next, the flux-forced simulations are conducted at 0.25^∘ resolution and can only partially capture the mesoscale eddies. Moreover, the flux-forced control simulation is not fully equilibrated (see Fig. 1c) due to the dynamic frazil formation at high latitudes that continuously adds heat in the polar regions. The frazil formation limits the magnitude of surface buoyancy flux perturbation we can apply in the polar regions. We have partially muted the frazil heat gain in the increased buoyancy flux contrast experiments by adding a globally uniform heat loss. In regions of extreme buoyancy anomalies, the isopycnal outcropping method (section 2b) is prone to capturing other elements of ocean circulation especially in regions of heat gain, such as the deep cell of the AMOC, which may produce erroneous results. Finally, the surface buoyancy flux perturbation experiments have not equilibrated even after 100 years, and hence, should not be misunderstood as the final response.

Our study reinforces recent evidence for a buoyancy-driven component of the ocean gyres (Gjermundsen et al., 2018; Hogg and Gayen, 2020; Liu et al., 2022). We envisage that a complete theory describing the formation of ocean gyres should incorporate the effects of surface buoyancy forcing, in addition to surface wind stress (Sverdrup, 1947). However, the influence of surface buoyancy forcing on gyres depends on the ocean state, and this induces non-linear and non-local feedbacks with the ocean circulation that obscure the formulation of a simple theory. These feedbacks include the role of the mixed layer in capturing excess heat, the non-uniform ingestion of surface heat fluxes by the ocean due to spatially varying thermal expansion coefficient, and the horizontal transport of heat by the circulation, all of which influence the background stratification.

where, e.g., $b(-h)\equiv b(x,y,z=-h,t)$. In (4), the numerator is the mean stratification averaged over depth $h$ and in the denominator, $|\bm{u}(0)-\bm{u}(-h)|/h$ is the magnitude of the resolved velocity shear averaged over depth $h$ while term $u_{\rm turb}/h$ quantifies the unresolved/turbulent velocity shear. The parameterization determines the mixing layer depth $h$ such that $\mathrm{Ri}_{b}(h)$ is equal to a critical Richardson number (taken to be $0.3$ in our simulations).

Prior to creating the flux-forced simulations, we conducted wind sensitivity experiments using ACCESS-OM2-025 (not presented here) along with the traditional K-profile parameterization reported in Large et al. (1994). However, the surface buoyancy forcing was inadvertently modified in these sensitivity experiments through anomalies in the mixing layer depth and ocean circulation. In an attempt to minimize surface buoyancy flux variations in the ACCESS-OM2-025 wind stress sensitivity experiments, we reconstructed the resolved velocity shear term $|\bm{u}(0)-\bm{u}(-h)|/h$ in (4) in the K-profile parameterization, which was found to primarily cause mixing layer depth anomalies in the sensitivity experiments. We parameterized $|\bm{u}(0)-\bm{u}(-h)|$ as a function of the friction velocity, $u_{*}=(\bm{|\tau|}/\rho_{0})^{1/2}$ and depth $h$:

where $c_{a}$, $c_{b}$, and $c_{e}$ are coefficients determined via multi-variate linear regression. Table 3 lists typical ranges of the three parameters for an optimal solution. The parameterization (5) performs well in the tropical and subtropical regions. Errors in the polar regions are expected because our parameterization (5) does not account for sea ice and marginal ice zone so as to stay consistent with the flux-forced experiments, which could alter resolved velocity shear in these regions.

The parameterization prevented discrepancies in the mixing layer depth due to alterations in the wind stress, and was implemented in the ACCESS-OM2-025 control experiment (Fig. 1c). For consistency with the ACCESS-OM2-025 control simulation, we retained the resolved shear parameterization in the flux-forced simulations.