The large-scale circulation of the abyssal ocean is enabled by small-scale diapycnal mixing, which observations suggest is strongly enhanced toward the ocean bottom, where the breaking of internal tides and lee waves is most vigorous. As discussed recently, bottom-intensified mixing induces a pattern of near-bottom up- and downwelling that is quite different from the traditionally assumed widespread upwelling. Here the consequences of bottom-intensified mixing for the horizontal circulation of the abyssal ocean are explored by considering planetary geostrophic dynamics in an idealized “bathtub geometry.” Up- and downwelling layers develop on bottom slopes as expected, and these layers are well described by boundary layer theory. The basin-scale circulation is driven by flows in and out of these boundary layers at the base of the sloping topography, which creates primarily zonal currents in the interior and a net meridional exchange along western boundaries. The rate of the net overturning is controlled by the up- and downslope transports in boundary layers on slopes and can be predicted with boundary layer theory.
The large-scale circulation of the abyssal ocean—below about 2000-m depth—is enabled by diabatic water-mass transformation. Antarctic Bottom Water sinks to the ocean bottom at the Antarctic margin, and it must cross density surfaces to come back up toward the surface (e.g., Lumpkin and Speer 2007; Talley 2013; Ferrari 2014). In steady state, this net upwelling across density surfaces must be balanced by diabatic transformation. This transformation is achieved by turbulence on scales smaller than about 100 m, produced primarily by breaking internal waves, which mix light water down into dense water (e.g., MacKinnon et al. 2013).1
Observations of the distribution of small-scale turbulence in the water column, however, suggest that mixing induces sinking instead of upwelling in the interior ocean (Polzin et al. 1997; Ledwell et al. 2000; St. Laurent et al. 2001; Ferrari et al. 2016; de Lavergne et al. 2016). This pattern emerged through microstructure measurements and tracer release experiments, which repeatedly showed that small-scale turbulence is strongly enhanced over rough topographic features and decays upward on a scale of a few hundred meters (e.g., Polzin et al. 1997; Ledwell et al. 2000; Waterhouse et al. 2014). Within this layer of bottom-intensified mixing, water parcels mix more with the dense water below than with the light water above, so they lose buoyancy and sink. This transformation has the wrong sign for allowing Antarctic Bottom Water to come back toward the surface.
This apparent conundrum is resolved if the ocean’s bathymetry is taken into account (Ferrari et al. 2016; de Lavergne et al. 2016; McDougall and Ferrari 2017). Water parcels adjacent to the bottom only mix with the light water above and thus gain buoyancy. If the bottom is sloping, this allows upwelling along the bottom and across density surfaces. The bathymetry also determines where in the water column the mixing occurs—and the geometries of density surfaces and the ocean bottom conspire to let bottom-trapped upwelling outweigh interior downwelling.
These arguments have so far only considered the water-mass budget. The implied pattern of up- and downwelling along topographic slopes, however, suggests that the dynamics of the abyssal circulation, also known as the lower cell of the meridional overturning circulation, are drastically different than previously assumed. We here explore the dynamical implications of bottom-intensified mixing on sloping bathymetry using a planetary geostrophic model in a simple “bathtub geometry.”
The theory of the abyssal circulation goes back to a set of classic papers by Stommel and Arons (1960a,b). Lacking firm observational constraints on small-scale turbulence at the time, they postulated uniform upwelling out of the abyssal ocean, as suggested by diffusive thermocline theory (Robinson and Stommel 1959). The horizontal circulation that develops in response to the upwelling can be inferred using planetary geostrophic vorticity dynamics (e.g., Pedlosky 1998):
where uy and uz are the meridional and vertical components of the velocity vector, f is the Coriolis parameter, and β is the meridional gradient of the Coriolis parameter.2 For simplicity, we assume a β-plane geometry here and throughout the paper. A spherical geometry would make a quantitative difference, but it is not expected to change the qualitative description of the abyssal circulation, which is what we are concerned with here.
Stommel and Arons (1960a,b) inferred from (1) that if there is upwelling at the top of the abyssal layer, uz > 0, and if the bottom is flat and thus uz = 0 there, fluid columns are stretched and move poleward. This poleward flow in the interior of the basin is fed by western boundary currents that connect the interior flow with the sources of abyssal water at high latitudes.
The predicted poleward interior flow is very weak and has never been observed. The great success of the Stommel–Arons theory was instead its prediction of deep western boundary currents, observations of which were emerging around the same time (Wüst 1955; Swallow 1957; Warren 1981). Deep western boundary currents, however, are not unique to the uniform upwelling envisioned by Stommel and Arons (1960a,b) and should not be taken as its confirmation. We will show below that bottom-intensified mixing on slopes drives a basin-scale circulation that similarly includes deep western boundary currents but no uniform upwelling and no meridional flows away from the ocean boundaries.
It is unlikely that the real ocean has widespread upwelling because there is not enough mixing to allow for it. Upwelling across density surfaces must be sustained by the diabatic transformation achieved by small-scale turbulence. Assuming the dominant balance in the buoyancy equation to be between vertical advection and turbulent diffusion,
Munk (1966) argued that the observed stratification below 1000 m in the central Pacific was “not inconsistent” with a constant upwelling and a constant turbulent diffusivity κ of the order 10−4 m2 s−1. But subsequent estimates of the turbulent diffusivity from observations of small-scale turbulence in the ocean interior turned out to be an order of magnitude smaller than Munk’s estimate (e.g., Gregg 1987; Ledwell et al. 1993; Munk and Wunsch 1998).
It was instead discovered that mixing is strongly enhanced over rough topographic features (e.g., Polzin et al. 1997; Ledwell et al. 2000; Waterhouse et al. 2014). This enhancement of mixing appears to be due to the nature of the physical processes that produce it. Tidal and geostrophic flows passing over rough topography displace isopycnals and thus generate internal waves (e.g., Garrett and Kunze 2007). These waves, when they are of large enough amplitude, induce convective and shear instabilities and thus produce turbulence, which is strongest near the bottom and rapidly decays away from it (e.g., Muller and Bühler 2009; Nikurashin and Ferrari 2010; Nikurashin and Legg 2011).
This bottom enhancement of small-scale turbulence has important consequences for the circulation (Ledwell et al. 2000; St. Laurent et al. 2001; Ferrari et al. 2016). Assuming mixing to be balanced by vertical advection as in (2), but allowing for a bottom-intensified turbulent diffusivity κ(z),
the observed stratification and estimates of mixing can be used to infer vertical velocities. If the turbulent diffusivity κ decreases with z more rapidly than the stratification ∂b/∂z increases, the turbulent buoyancy flux (which is equal to −κ∂b/∂z) is divergent, the right-hand side of (3) is negative, and water parcels become denser through the effect of mixing. In stable stratification, this means that there must be downwelling to balance the mixing.3 Only very close to the sea floor, where the turbulent buoyancy flux must go to zero, do water parcels gain buoyancy through mixing and does upwelling occur.
These dynamics of bottom-intensified mixing on slopes can be described by boundary layer theory, which considers the local response of a stratified ocean to mixing and to the no-flux condition on buoyancy along a sloping boundary (e.g., Phillips 1970; Wunsch 1970; Garrett et al. 1993). The theory describes how turbulent diffusion can be balanced by across-slope flow and, for bottom-intensified mixing, indeed predicts upwelling along the bottom, downwelling right above, and no up- or downwelling in the far field.
The bottom-intensified nature of mixing thus suggests that all up- and downwelling occurs in boundary layers on slopes and that there is little vortex stretching and hence meridional flow in the interior. The goal here is to understand how the bottom boundary layer flows drive a basin-scale circulation and produce net upwelling and overturning, as required to maintain an abyssal stratification.
It should be noted that the case of bottom-intensified mixing considered here produces very different dynamics than the previously considered case of enhanced mixing near the vertical sidewalls of a rectangular ocean basin (Marotzke 1997; Samelson 1998; Callies and Marotzke 2012). While the interior is largely adiabatic and interior meridional flow is suppressed in both cases, the way a basin-scale circulation is forced is very different. With vertical sidewalls and laterally intensified mixing, the circulation consists of upwelling along the zonal sidewall, balanced by downwelling associated with upright convection along high-latitude sidewalls. With a sloping bottom and bottom-intensified mixing, the circulation consists of large up- and downwelling along slopes, whose small residual balances the convective sinking at high latitudes (Ferrari et al. 2016; McDougall and Ferrari 2017). It is the vertical structure of mixing that produces this pattern of up- and downwelling, not its horizontal structure.
In the following, we explore the dynamics of an abyssal circulation driven by bottom-intensified mixing in two steps. First, we discuss the transient response of a uniformly stratified ocean to the insulating bottom boundary condition, which bends isopycnals and induces flow if the bottom is sloping (section 4). The bottom boundary layers that emerge exhibit the up- and downwelling layers anticipated from (3) for bottom-intensified mixing and are well described by boundary layer theory (section 5). The transient solutions eventually tend to a homogeneous ocean with no flow. To get a steady circulation, we subsequently force the production of dense water in the southern high latitudes of a closed basin. In addition to the bottom boundary layers, the circulation then develops basin-scale flows that feed the boundary layers, as well as western boundary currents that exchange fluid meridionally (section 6). The net upwelling and meridional exchange is tightly controlled by the bottom boundary layers on slopes, and boundary layer theory is used to predict the net overturning (section 7). The modeling approach and setup are discussed in the next section (section 3) and in an appendix with details of the implementation ( appendix B).
a. Planetary geostrophic equations
To understand the essence of how bottom-intensified mixing on slopes drives the abyssal circulation, we consider the planetary geostrophic equations,
in an ocean of idealized shape, where f = βy is the Coriolis parameter, u is the velocity vector, p is pressure divided by a reference density, b is buoyancy, F is a vector field representing friction, z is the vertical unit vector, and t is time. This system, with appropriate boundary conditions, yields a self-consistent and fully predictive description of the circulation and stratification. It satisfies a simple energy budget, as discussed in appendix A.
The planetary geostrophic equations (4)–(6) are an approximation to the Boussinesq equations under the assumption that the Rossby number is small, an assumption that is well justified for the large-scale circulation of the abyssal ocean. This approximation prevents the development of mesoscale eddies, and it filters out both inertia–gravity waves and small-scale turbulence. The buoyancy transfer by small-scale turbulence is represented by diffusion. The planetary geostrophic equations have been used widely as a starting point in the analytical and numerical study of the large-scale ocean circulation (e.g., Robinson and Stommel 1959; Welander 1959; de Verdière 1988; Samelson and Vallis 1997a; Salmon 1998a; Pedlosky 1998). Notably, the neglect of inertia in (4) has no effect on the steady flow in boundary layers on slopes (cf. Garrett et al. 1993).
We consider the planetary geostrophic equations with a continuous vertical coordinate because that allows the representation of bottom-intensified mixing with a bottom-intensified diffusivity. Representing small-scale mixing in layered models is more difficult; the effect of mixing is often parameterized by restoring or the prescription of water-mass transformation (e.g., Tziperman 1986; Kawase 1987; Spall 2001). The continuous system instead allows a self-consistent balance between advective and diffusive terms in the buoyancy budget, and a buoyancy flux boundary condition can easily be prescribed at a sloping bottom.
For the planetary geostrophic flow described by (4)–(6) to satisfy a no-normal flow condition on all boundaries, a form of friction F must be specified. We follow Salmon (1990, 1992, 1994) in resorting to Rayleigh friction,
instead of the more common Fickian diffusion of momentum. While Rayleigh friction implies an unphysical loss of momentum in the interior of the fluid, it has the practical advantage of leading to mathematically simpler boundary layers, both on the western boundary (cf. Stommel 1948; Munk 1950) and in boundary layers on slopes (section 5). These drag-controlled boundary layers, however, should be regarded as nothing more than simple stand-ins for the more complicated and incompletely understood turbulent dynamics of the real boundary regions. The friction parameter r is chosen to be small enough for the momentum balance to be dominantly geostrophic away from boundary layers and the equator. The sensitivity of our solutions to the friction parameter is discussed in section 8.
Modern coarse-resolution ocean models employ a Gent and McWilliams (1990) parameterization of mesoscale eddies, which introduces a tendency to flatten out isopycnals in the buoyancy budget. We forgo such a parameterization here for the sake of simplicity. Implementing such a mesoscale parameterization in a thickness-weighted average framework would move the additional tendency term to the momentum equations, where it would appear as the divergences of Eliassen–Palm vectors (Young 2012) and replace our Rayleigh drag parameterization.
The friction parameterization in (7) imposes drag only in the horizontal because the depth of the ocean basin we will consider goes to zeros continuously at its margins (Salmon 1992). If vertical sidewalls were present, we would need additional terms to accommodate thermal-wind shear there. This is typically done either through friction in the vertical, which generates nonhydrostatic upwelling layers (Salmon 1986, 1990) or through hyperdiffusion of buoyancy, which generates complicated thermal boundary layers (Samelson and Vallis 1997b). Neither of these is necessary in our case, an additional indication that bottom slopes are a crucial element of the low-Rossby-number dynamics.
c. Boundary conditions
With Rayleigh friction (7), the flow can satisfy a no-normal flow boundary condition at the bottom, at , where is the bottom-normal and is the basin depth, as well as a rigid-lid condition at the surface, at . We neglect wind forcing because we are concerned with the abyssal ocean, but in principle this could be included through the prescription of an Ekman pumping velocity at the surface.
For buoyancy, we impose an insulating boundary condition at the bottom, at , neglecting any geothermal heat flux. Given that our focus is on the abyssal circulation, the surface of the domain at can be thought of as being located at 2000-m depth in the real ocean. Since isopycnals are approximately flat in the deep ocean basins, we prescribe a constant buoyancy at . This is not appropriate for the Southern Ocean, however, where isopycnal slopes are large. In reality, these sloping isopycnals provide an adiabatic pathway to the surface for the deep waters flowing out of the oceans basins. Talley (2013) showed that the deep waters that enter the Southern Ocean at a depth below 2000 m come to the surface in the Weddell and Ross Seas, where strong surface buoyancy forcing through brine rejection by freezing seawater and atmospheric cooling transforms them back into denser waters, which sink to the bottom. This transformation of deep to bottom waters is in balance with the diapycnal transformation of bottom to deep waters through mixing that characterizes the lower cell of the meridional overturning circulation. Without the surface transformation in our model (i.e., with at ), all solutions eventually tend to a homogeneous ocean. The transients nevertheless provide crucial insight into the boundary layer dynamics and are discussed in section 4.
To allow for a steady overturning and stable stratification, we subsequently include a crude representation of dense water formation. In the real ocean, winds and mesoscale eddies in the Antarctic Circumpolar Channel are thought to set the isopycnal slope, mapping the meridional distribution of surface buoyancy to a vertical stratification at the northern edge of the Southern Ocean (e.g., Marshall and Radko 2003; Nikurashin and Vallis 2011). This process operates on a time scale of decades at most, so it is much faster than the diffusive dynamics of the abyssal overturning, which operates on centennial to millennial time scales. In section 6, the effect of this fast process is represented through restoring of buoyancy to a prescribed stratification in the southern part of the basin. This restoring acts to transform deep to bottom waters, and solutions can reach a steady state with realistic stratification and overturning in the basin.
We consider a “bathtub geometry” with the bottom sloping like half-Gaussians on all sides:
The basin extends zonally from to and meridionally from to . The maximum depth is scaled by H, and the slope width is controlled by . We use throughout (Fig. 1), which gives slopes comparable to those on the flanks of deep-ocean ridges, where most of the abyssal mixing is thought to occur.
d. Nondimensionalization and parameters
We nondimensionalize the system with the coordinate transformation
and the substitutions
Note that the factor in the velocity substitution has units of inverse time because the coordinate transformation takes care of the length units (, , ):
The distinction between the coordinate transformation (9) and the substitutions (10) is important when employing tensor calculus notation (e.g., Grinfeld 2013), and it becomes useful when we use terrain-following coordinates and slope-aligned coordinates ( appendix B).
In the substitutions (10), we assumed to have available a stratification scale N2. In the solutions below, this will be the initial stratification or the stratification we restore to in the southern part of the domain. In both cases, the actual bulk stratification remains fairly close to this scale.
The equations in nondimensional Cartesian coordinates then read
and they contain the parameters
These parameters have straightforward interpretations: is the aspect ratio of the basin, is the ratio of the time it takes a long Rossby wave to propagate across the basin (using ) to the diffusive time (if κ is constant), and is the ratio of the Stommel boundary layer width to the basin width L. If present, any spatial dependence in κ is inherited by .
To get a sense for the size of these parameters, let us consider typical scales:
The aspect ratio of the ocean is very small, . This renders horizontal diffusion very small compared to vertical diffusion, which is challenging to capture in a numerical model. We will thus artificially increase the aspect ratio to . Such an increase is equivalent to increasing the horizontal diffusion coefficient, which is not expected to qualitatively affect the solution as long as (cf. Veronis 1975). Note that an increase in has no effect on the dynamics or geometry of the problem other than an increase in horizontal diffusion— appears nowhere in the equations but in front of the horizontal diffusion term in (16), and the nondimensionalized bottom depth in (8) becomes independent of as well. Using Munk’s (1966) estimate of an average diffusivity, , gives , suggesting that the ocean is generally in a regime of weak diffusion (i.e., ). We will discuss different values and spatial distributions of in the following sections. The friction parameter is somewhat arbitrary; we will use , such that the Stommel boundary layer width is a tenth of the domain width. The sensitivity to this parameter is investigated in section 8.
e. Numerical model
(e.g., Phillips 1957; Salmon 1998b; Haidvogel et al. 2000; appendix B). These coordinates simplify the bottom boundary conditions: and at become and at , where upper indices denote contravariant components.4 The price to pay is that buoyancy forces enter the horizontal momentum equations because -coordinate lines are not horizontal and that cross terms appear in the diffusion operator because the coordinates are not orthogonal. These extra terms are easier to keep track of if tensor calculus notation is employed. The pressure gradient error arising from the cancellation of two large terms in the pressure gradient (Haney 1991) is alleviated by a relatively high horizontal resolution.
The system (12)–(16) consists of a conservation equation for buoyancy and a diagnostic relation between buoyancy and the flow. We thus implement the dynamics by computing the flow from a given buoyancy field and then using that flow field to time step buoyancy. The flow field is obtained by first solving a two-dimensional elliptic problem for the streamfunction of the vertically integrated flow. Subsequently, the vertical shear of the horizontal flow is obtained from the frictionally modified thermal-wind balance and is added to the depth-average flow. The velocity component is obtained by integrating the continuity equation in . See section b of appendix B for details.
For the discretization, we use standard centered finite differences on a grid that is equally spaced in all three coordinates. One-sided second-order differences are used in the cross terms of the diffusive fluxes near boundaries, where centered differences are not possible. All solutions are obtained with a grid spacing of and , a choice motivated by the aspect ratio and a grid fine enough for the discretization error to be small. The grid is staggered, with located on faces, on faces, on faces, at cell centers, and the streamfunction for the vertically integrated flow at (horizontal) cell nodes. Among a number of advantages, this staggering allows us to sidestep problems arising from the singularity of the coordinate system at the coasts, where . Time stepping with is explicit for the advection step, the horizontal diffusion step, and the step involving the cross terms in the vertical diffusion term. Vertical diffusion (excluding cross terms) is treated implicitly because the vertical grid spacing becomes very fine where the depth goes to zero, such that an explicit treatment would be overly restrictive on the time step. The restoring step in the equilibrating solutions is performed through an exact integration of the continuous equation over a time . The Julia implementation (Bezanson et al. 2017) is available online (at https://github.com/joernc/pgcm).
4. Transient solutions
We begin by considering the transient adjustment of a stratified ocean to the insulating bottom boundary condition on buoyancy. The solutions are initialized with a constant and uniform stratification, (dimensionally ), a buoyancy field that in the absence of wind forcing is associated with no flow. This initial condition satisfies at , but not at .
While these transient solutions eventually tend to a uniform buoyancy field because there is no source of dense water, the initial phase reveals key characteristics of the diffusively driven abyssal flow. Diffusion sets up bottom boundary layers, in which the bulk of the water-mass transformation occurs, and across- and along-slope flows develop in response. We will recognize these bottom boundary layers in the forced steady solutions discussed in section 6, and we will see that they exert important control on the overall basin circulation. We discuss the transient solutions at , which is too short for diffusion to significantly alter the bulk interior stratification, but it is long enough for the circulation in the boundary layers on slopes to develop.
a. Flat bottom
To understand the role of bottom slopes in shaping the circulation, it is instructive to briefly consider what happens in their absence. If the bottom is flat, , no horizontal buoyancy gradients develop, and the ocean remains motionless.5 The problem reduces to a one-dimensional diffusion equation in , which for a uniform diffusivity has the solution
Buoyancy converges near the bottom, resulting in a growing well-mixed boundary layer there. Eventually, after a diffusive time , buoyancy becomes fully homogenized.
b. Uniform mixing
This trivial behavior on a flat bottom contrasts with the case with the “bathtub bathymetry” given by (8). On slopes, buoyancy diffusion tilts isopycnals to satisfy the insulating boundary condition—and thus induces flow.
To illustrate this, we begin with the case of uniform diffusivity . We solve the equations numerically, as described in the previous section. The relatively large value of (corresponding roughly to ) is chosen mostly because of numerical constraints; a smaller would require higher resolution and thus increase the computational demand beyond what is easily achievable without parallelization. With this value of , all circulation features appear to be reasonably well resolved. A smaller value of would not change the results qualitatively but simply reduce the thickness of boundary layers that develop along the sloping boundaries (cf. section 5).
As in the flat-bottom case, the diffusive buoyancy fluxes converge near the bottom, where the no-flux boundary condition must be satisfied, and boundary layers develop rapidly (Figs. 2a–c). These boundary layers are of different character on slopes, however, where isopycnals tilt. This tilting causes a flow that quickly arrests the growth of the boundary layers: the convergence of buoyancy on the slope is balanced by an upslope advection of dense water, which tends to flatten isopycnals and thus maintains stratification. The tilted isopycnals on slopes are associated with along-slope geostrophic flow in the direction opposite that of Kelvin wave propagation (cf. Benthuysen and Thomas 2012): clockwise in the Northern Hemisphere and counterclockwise in the Southern Hemisphere. There are weak interior flows, as mandated by mass balance.
The adjustment to the insulating boundary condition on slopes thus induces a boundary layer flow that was not anticipated by Stommel and Arons (1960b). Even though the mixing is uniform, vigorous upwelling occurs on the slopes, with broad compensating downwelling in the interior. There is no net upwelling in this transient case because there is no site of dense water formation. But the flow’s tendency to localize upwelling in boundary layers on slopes carries over to the equilibrium case and can be understood with the boundary layer theory of Phillips (1970), Wunsch (1970), and Garrett et al. (1993), as discussed in section 5.
c. Bottom-intensified mixing
Before discussing boundary layer theory, we turn to the more realistic case with bottom-intensified mixing, in which the boundary layers on slopes acquire a layer of downslope flow above the bottom-trapped upslope flow. We set
with and . The volume-average diffusivity here is roughly the same as in the case with uniform mixing. If the stratification remains similar between the two cases, the sources of potential energy will have similar magnitudes ( appendix A), and a circulation of similar strength is expected. The relatively large value of is chosen for numerical convenience because we are after the qualitative character of the induced circulation. An analysis of the dependence of the solution on the bottom value is performed in section 7.
Flow quickly sets up in boundary layers on slopes (Figs. 2d–f). The flow right above the bottom is similar to the boundary layer flow of the uniform-mixing case, though somewhat stronger. In this bottom-trapped layer, the insulating bottom boundary condition again renders the diffusive buoyancy flux convergent. Isopycnals dip down, inducing upslope flow advecting dense water from below and an along-slope flow in the direction opposite that of Kelvin wave propagation.
The bottom intensification of , however, renders the diffusive buoyancy flux divergent above the bottom layer. In the layer of diffusive buoyancy flux divergence, which Ferrari et al. (2016) called the “stratified mixing layer,” a negative buoyancy anomaly prevails, such that isopycnals slope up slightly before they dip down in the bottom layer. In response to the negative buoyancy anomaly, a downslope flow develops. The resulting downslope advection of buoyant water from above balances the diffusive buoyancy flux divergence, such that the boundary region is close to a steady state at the time of the diagnostics. The downslope flow is weaker but broader than the upslope flow below, carrying about the same transport. The along-slope flow also reverses direction above the bottom layer, becoming aligned with the direction of Kelvin wave propagation. The interior flow is weak.
5. Boundary layer theory
The dynamics of the near-bottom flow in the transient solutions, both for uniform and bottom-intensified mixing, can be understood with boundary layer theory. As discussed in the next section, these boundary layers also emerge in equilibrated solutions of the full system—and in fact exert a controlling influence on the entire abyssal circulation.
The approach taken here is analogous to that of Phillips (1970), Wunsch (1970), and Garrett et al. (1993). We consider the local adjustment to the insulating bottom boundary condition of a fluid with prescribed stratification in the far field. Assuming there are no variations in the across- and along-slope directions, the dynamics reduce to a one-dimensional problem in the slope-normal direction. For ease of exposition, we work in dimensional quantities but perform all calculations in nondimensional quantities. The nondimensional equations are given in appendix C.
In contrast to previous studies of the dynamics of boundary layers on slopes, we use Rayleigh friction instead of Fickian diffusion of momentum. Rayleigh friction is crude but yields simple solutions. It should be regarded as a stand-in for the insufficiently understood physics of turbulent boundary layers on a rough and sloping bottom.
We begin by transforming the planetary geostrophic equations (4)–(6) to a coordinate system that is aligned with an infinitely extending bottom at , with fixed slope angle θ: , , and . The across-slope, along-slope, and slope-normal velocity components are then , , and . We split the buoyancy and pressure fields into an imposed background, which has constant stratification and is in hydrostatic balance, plus a perturbation: , and . As the slope is infinitely extending in and , we assume the perturbations to depend on the slope-normal distance only. For consistency, we also neglect variations of the Coriolis parameter, , which is justified as long as the horizontal width of the diffusive boundary layer on the slope is smaller than the Stommel boundary layer width . In the absence of and variations, the continuity equation can only be satisfied if . Then,
The only terms left in the buoyancy equation are the tendency term, the advection of the mean buoyancy field by the across-slope flow , and the divergence of the diffusive flux (including that due to the background stratification). The insulating boundary condition at is
Together with the boundary condition (25), this determines the time evolution of ; the up- and along-slope velocities follow from the momentum equations. The system can reach a steady state by balancing diffusion with across-slope advection.
For constant κ, there is a simple analytical solution. In steady state,
The solution has a bottom boundary layer of thickness , in which buoyancy is increased () to cancel the contribution of the background stratification and satisfy the insulating bottom-boundary condition (Fig. 3). The across-slope flow is in the upslope direction ( for ), as required to balance the convergence of the diffusive flux with advection of dense water from below. The along-slope flow is in the direction opposite that of Kelvin wave propagation ( for ) and is roughly in geostrophic balance away from the equator, where or . This flow has the same structure as that found in the full transient solution with uniform mixing (Figs. 2a–c).
A more quantitative comparison between boundary layer theory and numerical solutions is obtained by considering the approach of the steady state, that is, the transient solution to (26) with initial condition . This corresponds to an initially constant stratification, the same case as considered in the full numerical solution described in section 4. The time-dependent solution to (26) then is
where erfc is the complementary error function. This shows that the steady state (27) is approached on a time scale of , the time it takes to diffuse across the steady boundary layer of thickness .
To see to what degree this local analytical solution is in agreement with the full numerical solution of the previous section, we compare (28) pointwise to the full solution. At every point on the ocean bottom, we pick the local slope and Coriolis parameter, and we assign (28) in the column above. The buoyancy and velocity fields so constructed only approximately satisfy the full equations of motion because we disregard variations in slope and Coriolis parameter as well as surface boundary conditions. Nevertheless, this heuristic solution captures the shape of the isopycnals in the bottom boundary layer, the direction and magnitudes of the flow, and the horizontal variation of these boundary layer properties (Figs. 2a–c and 4c–e). This suggests that we were justified in neglecting slow variations in slope and Coriolis parameter in the spirit of a WKB approach.
On the slopes, the solution in (27) predicts a steady boundary layer that has a thickness much less than the domain depth, (Fig. 4a; cf. appendix C). The time it takes to set up these steady boundary layers is short compared to a Rossby wave transit time, or (Fig. 4b). This time scale is thus much shorter than the time that we picked for Figs. 2a–c and 4, so the boundary layers there have likely reached steady state.
The thickness of the steady boundary layer and the time scale over which the steady solution is approached go to infinity where the bottom becomes flat (). On a perfectly flat bottom (), there is no steady-state solution. The problem in (26) reduces to a pure diffusion equation, with solution
which describes a boundary layer that keeps growing with time. This is the same scenario as the transient solution with a flat bottom discussed in section 4a, except that here we have assumed a semi-infinite ocean, so (29) is required to decay as instead of having to satisfy at the surface, as in (20).
Unlike on the slopes, where downward buoyancy diffusion can be balanced by across-slope upwelling of dense water, the growth of bottom boundary layers on the flat part of the basin cannot be arrested locally. The time scale of the boundary layer evolution becomes so long, however, that basin-scale processes can enter the budget. In the steady solutions that include dense water formation discussed in the next section, lateral advection of dense water arrests the growth of boundary layers on flat bathymetry. On the slopes, on the other hand, the growth of boundary layers is arrested locally by upwelling along the seafloor.
The simple diffusive boundary layer solution on slopes also breaks down when the β effect becomes important in the vorticity budget of the horizontal flow. This occurs when the horizontal width of the diffusive boundary layer becomes comparable to the Stommel boundary layer width . This breakdown again occurs on gentle slopes. With the parameters chosen for our solution, the β effect is of leading order in the section between and (Fig. 4).
If mixing is bottom intensified, the structure of the boundary layer solution changes, consistent with the results from the full transient solution in the previous section. With the diffusivity decaying exponentially with height above the bottom,
the boundary layer equation in (26) is most easily solved numerically. As in the case with Fickian diffusion of momentum (Garrett 1990), an analytical solution is possible, but the result is unwieldy. The steady solutions on slopes now have a two-layer structure (Fig. 3). Isopycnals bend up () before dipping down and intersecting the bottom (). The across-slope velocity is upslope at the bottom and downslope above, balancing the diffusive flux convergence near the bottom and its divergence above. The along-slope flow also reverses sign, as required by the reversing isopycnal slopes. The up- and downslope transports balance exactly—the total transport vanishes, as required by the integral of (26) in the slope-normal coordinate and the fact that as (Phillips et al. 1986; Garrett et al. 1993).
On the slopes, this steady boundary layer solution matches the boundary layers of the full transient solution (Figs. 2d–f). On the flat parts of the bathymetry, on the other hand, the boundary layers again are predicted to grow without limit, and lateral advection by a basin-scale circulation must enter the budget to reach a steady state.
6. Steady solutions
The transient solutions discussed above eventually tend to a homogeneous ocean () with no flow. In the real ocean, a deep stratification is believed to be maintained by diabatic transformation in the Southern Ocean. Through the joint effect of winds and mesoscale eddies, deep waters are brought to the surface, where surface processes turn them into denser waters that sink to the ocean bottom: salinification through brine rejection of freezing seawater and atmospheric cooling. This surface transformation and the circulation along strongly sloping isopycnals in the Southern Ocean is what prevents the abyss from becoming a stagnant homogeneous pool.
Here we are interested in the diabatic transformations away from the surface, which to a large degree occur in the ocean basins. To keep the problem as tractable as possible, we dramatically simplify the processes in the Southern Ocean. Instead of explicitly representing the wind-driven circulation in a circumpolar channel, we pretend the effect of this circulation is simply to restore buoyancy to a particular stratification at the northern edge of the Southern Ocean and to provide any mass flux required by the basin circulation. As a substitute for the effects of wind forcing, eddies, and surface transformations, we thus simply restore buoyancy to a given stratification N2, which is taken to be constant:
The restoring constant λ is chosen such that it transitions from λ0 to zero over a distance lt around the latitude ys:
We make the restoring fast compared to diffusion by picking the nondimensional such that it is about 100 times the average (Table 1). We choose and .
All solutions converge to a steady state. There is no time dependence in the equilibrated states, consistent with what Salmon (1990) found across a wide range of solutions to the planetary geostrophic equations. In steady state, the water-mass transformation in the basin is balanced by the restoring in the southern part of the domain. The restoring adjusts to balance the diabatic transformation in the basins and thus accepts the net overturning that is induced by mixing. One could also consider the opposite experiment: impose in the south a buoyancy sink and thus net overturning and let the stratification adjust until mixing balances that transformation. Such an experiment is left to future work.
a. Uniform mixing
Before considering the more realistic case with bottom-intensified mixing, we discuss uniform mixing with (and ). We time step the transient system, starting from , until the solution has converged to its steady state. This has occurred after one diffusive time scale , the longest time scale in the system, as confirmed by inspecting the time evolution of the volume-averaged buoyancy and the cross-equatorial overturning (discussed below).
Throughout the basin, the steady-state solution has a stratification that is not too far from the stratification prescribed in the south (Fig. 5a). Exceptions are a basinwide benthic layer overlying the flat bottom, in which isopycnals bend down, and thin boundary layers on slopes.
The boundary layers on slopes have a structure familiar from the transient solutions and from the boundary layer theory of section 5 (Figs. 6a–c). As before, these boundary layers are associated with up- and along-slope flow. At the base of the sloping topography, where the bathymetry flattens, the upwelling in the boundary layers has to be fed by dense water. This is achieved by a basin-scale circulation that is now present. This circulation carries dense water north in a boundary current flowing along the western slope (Figs. 6b and 7a), directly feeds the upwelling on the western boundary, and connects to the upwelling on the eastern boundary through a zonal current in the benthic layer (Fig. 6a). The water upwelled on the eastern boundary is then returned westward by a zonal current above the benthic layer (Fig. 6a). A southward western boundary current above the northward one returns the upwelled water to the south (Figs. 6b and 7a), where the loop is closed by the transformation to dense water achieved by the restoring. The western boundary currents are classic Stommel boundary currents and have a width scaling with .
The upwelling is concentrated in the boundary layers on slopes, but there is weak widespread interior upwelling as well (Fig. 6c). This occurs where changes enough in the transition from the weakly stratified benthic layer to the stratified interior, such that diffusive fluxes converge, and this does drive a weak Stommel–Arons-like poleward flow there (Fig. 6b). But the full circulation has an important contribution controlled by the boundary layers on slopes, even in this case with uniform mixing.
A striking feature of this solution—and all steady solutions discussed in the following—is that isopycnals are to leading order flat in the interior of the domain. This is consistent with observations of the real ocean’s deep hydrography—as seen, for example, in WOCE sections (Talley 2007)—but it is far from obvious. The explanation typically given for this observation is that at a vertical eastern boundary, the no-normal flow boundary condition through thermal wind requires that meridional buoyancy gradients vanish and that Rossby wave radiation carries this signal into the interior. This argument, however, fails if the eastern boundary is sloping instead of being vertical. There is then no a priori restriction on meridional buoyancy gradients because even with finite zonal shear all the way to the coast, the flow itself smoothly goes to zero as the depth goes to zero.
That the stratification is close to that prescribed in the south—or in the real ocean to that at the northern edge of the Southern Ocean—is then more usefully understood as a consequence of the weakness of mixing. If there were no mixing (), a steady solution to the equations of motion would be , that is, a uniform stratification matching that prescribed in the south and no flow. The nontrivial solution arising in the presence of mixing can then be considered a perturbation to this state. This produces deviations from the uniform stratification, particularly near the boundaries, and thus drives a circulation. As mixing is small, , the perturbation to the state is expected to be small.6
The fact that means that Rossby wave radiation across the ocean basin is much faster than diffusion across the full depth of the domain. Rossby waves can flatten out isopycnals much more effectively than diffusion can tilt them (except in boundary layers). That Rossby waves tend to flatten isopycnals is expected from the energy budget in (A1) and (A2) because all planetary geostrophic flow releases potential energy: . This process appears to be sufficient to yield approximately flat isopycnals—the boundary condition at a fictitious vertical wall in the east is not required.
The approximate flatness of isopycnals also means that there are strong meridional gradients of potential vorticity . In the absence of potential vorticity sources or sinks in the interior, water parcels thus travel zonally. To travel meridionally in the overturning circulation, water parcels must adjust their potential vorticity. Since in our solution meridional flow occurs primarily in boundary layers, this adjustment is easily achieved by frictional potential vorticity fluxes. In the real ocean, lateral stirring likely plays a role in the exchange of potential vorticity with the boundary (Edwards and Pedlosky 1998), but how this occurs on a sloping bottom remains insufficiently understood. Since potential vorticity dynamics do not appear to yield any deeper understanding of the circulation dominated by boundary layer flows, we do not discuss them any further.
b. Bottom-intensified mixing
When mixing is bottom intensified, the circulation becomes even more strongly controlled by the boundary layers on slopes than already was the case with uniform mixing. We begin with the case and , whose transient evolution we considered above. We integrate the system (with restoring ) to , at which point the solution has converged to its steady state.
The solution has a bulk interior stratification close to that prescribed in the south, except in a weakly stratified benthic layer on the flat part of the bathymetry and in thin boundary layers on the slopes (Fig. 5c). The structure of the boundary layers on slopes in this equilibrated solution is familiar from the transient solution and boundary layer theory: isopycnals slightly slope up before dipping down, there is strong narrow upwelling on the slope and weaker but broader downwelling above, and the along-slope current similarly shows a two-layer structure (Figs. 6d–f). As with uniform mixing, there is a basin-scale circulation connecting the eastern boundary layers on slopes to the western boundary as well as western boundary currents connecting the circulation to the dense-water formation in the south (Figs. 6e and 7c). Western boundary currents are required in order to close the mass budget because there is little net meridional transport in boundary layers on slopes—the transports approximately cancel between the upwelling and downwelling layers, as expected from boundary layer theory.
The circulation thus consists of a deep northward Stommel boundary current on the western slope, upwelling on the slopes (Figs. 6e and 7c), downwelling above (Fig. 6f), and a return southward Stommel boundary current at middepth (Figs. 6e and 7c). The upwelling on the eastern boundary is supplied with dense water from the western boundary by an eastward zonal current in the benthic layer (Fig. 6d). Water downwelled in the east is brought back west by a westward zonal current just above the benthic layer (Fig. 6d). A simple schematic of the these flows is shown in Fig. 8. In addition, there are along-slope currents in the diffusive boundary layers on slopes (Fig. 6e).
From boundary layer theory, as discussed in the previous section, the transports up and down the slope are expected to be equal because decays to zero away from the boundary (Fig. 3b). This prediction is roughly consistent with the full solution—we will see in the water-mass transformation below that there is a large degree of compensation between up- and downwelling. Because of this vanishing net boundary layer transport, changes in slope do not induce any convergences or divergences of the net across-slope flow, which would drive exchanges with the interior (Phillips et al. 1986). This argument for no exchange with the interior, however, fails at the base of the sloping topography. There, boundary layer theory itself fails because no steady boundary layer solution is possible. Instead, the up- and downwelling on the slopes must be accommodated by a basin-scale circulation. On the western boundary, Stommel boundary currents achieve this: they import dense water from the south at the bottom, and they export more buoyant water above. For the eastern boundary, zonal currents in the benthic layer and just above connect the up- and downwelling to the western boundary currents. It thus appears that the in- and outflows at the base of the sloping topography, which are driven by the boundary layers on the slopes, require a basin-scale flow. Because these in- and outflows at the base of the sloping topography occur at different density classes, they drive a meridional overturning.
To diagnose more quantitatively the compensation between up- and downwelling in the boundary layers on slopes, we perform a water-mass transformation analysis of the full steady solution with bottom-intensified mixing (Walin 1982; Garrett et al. 1995; Marshall et al. 1999; Ferrari et al. 2016). The water-mass transformation in the Northern Hemisphere is defined as
where is the (nondimensional) sum of the diabatic terms on the right of (31), and
is the volume between the isopycnals and , bounded to the south by the equator. The restoring term very nearly vanishes in any , so the only term that significantly contributes to in (33) is the diffusive term. Note that corresponds to the negative of the streamfunction of the zonally integrated meridional circulation at buoyancy , evaluated at the equator, as required by continuity. To illustrate the compensation between positive and negative transformations, we split (33) into positive and negative contributions: , where are defined by replacing by in (33), with and .
The water-mass transformation in the Northern Hemisphere shows that there is significant compensation between positive and negative contributions (Fig. 9b). The two contributions nearly cancel in the upper part of the water column. This is expected from boundary layer theory—the transformation is largely confined to boundary layers on slopes. The compensation is instead weaker deeper down, where boundary layer solutions start breaking down at the base of the sloping topography, and a basin-scale circulation feeds the boundary layers (Figs. 6d–f). As the inflow into the boundary layer occurs at a lower buoyancy than the outflow, there is now a net positive transformation. This is the net transformation that leads to the cross-equatorial flow and thus the net overturning. We will see in the next section that the boundary layers on slopes control the magnitude of the in- and outflow at the base of the sloping topography. The boundary layer solutions are shown to yield a prediction for the overturning in our simple bathymetry.
It should be noted that there is some water-mass transformation also in the benthic layer: bottom water experiences buoyancy flux convergence and becomes lighter as it travels east across the flat part of the basin, and the overlying water experiences buoyancy flux divergence and becomes denser as it travels back west. These transformations in the flat part of the basin very nearly cancel and therefore do not contribute to the net overturning strength (Fig. 9b). Instead, they act to shift the net overturning upward somewhat in buoyancy space.
The net water-mass transformation and thus overturning is about 0.01 (Fig. 9b). Given the parameters in (18) and , this corresponds to about 5 Sv (1 Sv ≡ 106 m3 s−1), which is somewhat small compared to the Lumpkin and Speer (2007) estimate of about 15 Sv in the abyssal Pacific, considering that (corresponding to ) is relatively large. The overturning has the right order of magnitude, however, which is all we can expect given the idealized nature of this study.
It is customary to display the net meridional overturning using an overturning streamfunction defined by zonal averages in buoyancy space. The structure of this streamfunction, mapped back into physical space, resembles observational estimates (Fig. 10; cf. Lumpkin and Speer 2007). When displaying such zonal averages, however, it should be kept in mind that the net upwelling is a residual of large up- and downwelling flows on the slopes. The circulation sketched in Fig. 8 is only partially visible in the overturning streamfunction.
7. Predicting the overturning
The phenomenology of the circulation that arises in response to bottom-intensified mixing suggests that the boundary layers on slopes exert a strong control on the circulation. The net overturning circulation appears to result from the net transformation in the boundary layers. This net transformation occurs at the base of the sloping topography, where the boundary layers are fed with inflow that is denser than the outflow above. If the magnitude of these in- and outflows is determined by the boundary layer solutions on the slopes, we can integrate these solutions along the perimeter of the basin and get an estimate of the overturning.
We thus obtain a prediction for the cross-equatorial overturning:
where is a particular depth contour, confined to the Northern Hemisphere, is the slope-normal coordinate, is where changes sign, and is an along-contour coordinate. The boundary layer solution is obtained by numerically solving for the steady state of the boundary layer equations (22)–(24) at every point along the depth contour, using the local slope and local latitude. The boundary layer solutions and thus can be calculated without any input from the full solutions because the interior stratification remains close to that prescribed in the south, so the far-field stratification is an external parameter.
To test this prediction, we obtain a range of steady numerical solutions by varying from 0.025 to 0.8, keeping fixed (Table 1). The time it takes for the numerical solutions to converge to a steady state roughly varies with . All solutions have converged by the time chosen for diagnosis.
Across the explored range of values, the solutions retain the structure described in section 6b. The smaller the , the thinner the weakly stratified benthic layer (Figs. 5b–d). The water-mass transformation is then increasingly compensated and the net transformation moves to a narrow density range near the bottom (Fig. 9). Consistently, the cross-equatorial flow becomes progressively confined to a thin layer above the bottom (Figs. 5b–d). The boundary layers on slopes are also thinner, the smaller the , as expected from boundary layer theory (Figs. 5b–d and 11a). The interior retains a stratification close to that prescribed in the south (Figs. 5b–d).
The prediction in (35) compares well with the diagnosed cross-equatorial overturning from (33), both in magnitude and in the scaling with (Fig. 11b). There is only a weak dependence of the prediction in (35) on what depth contour is chosen, as long as that depth contour is above the weakly stratified benthic layer. The overturning roughly scales with , but there is some curvature.
The success in predicting the net transformation with boundary layer theory confirms that the boundary layers on slopes control much of the circulation. While the net transformation occurs at the base of the sloping topography, where the boundary layer solutions break down, the flow there appears to be slaved to the boundary layer flows on the slopes above.
This also suggests that the transformation in the benthic layer, even if mixing is (unrealistically) strong there, has little effect on the net overturning. To test this assertion, we performed an additional simulation with mixing coefficients reduced by an order of magnitude over the flat part of the bathymetry, with the rest of the setup identical to the reference case with bottom-intensified mixing. The scaling argument predicts the same overturning strength, because it depends on the mixing on the slopes only, which is unchanged. This prediction is borne out in the additional simulation: the benthic layer becomes more stratified and the overturning circulation shifts downward somewhat in buoyancy space (Fig. 9d), but neither the shape nor the strength of the overturning circulation changes substantially compared to the case with enhanced mixing everywhere. This again confirms the tight control of the overturning circulation by the boundary layers on slopes.
8. Friction dependence
The control of the overturning by boundary layers on slopes and the crucial role that friction plays in these boundary layers raises the question of how the circulation depends on friction. The control of the circulation by boundary layers is different from the situation in Stommel–Arons theory or in linear gyre theory, where boundary currents are passive and simply close the mass budget as required by interior transport.
Unlike the diffusivity , which directly relates to observable buoyancy flux divergences achieved by small-scale turbulence, the friction parameter is quite arbitrary and relies on the heuristic representation of turbulent processes by Rayleigh friction. We discuss the dependence of the circulation on this parameter to get a sense for the sensitivity of our solutions.
We obtain two additional full steady solutions with changed from 0.1 to 0.05 and 0.2, keeping all other parameters the same as in the bottom-intensified mixing case with discussed in section 6b. The qualitative structure of the circulation remains unchanged across these solutions. The boundary layers on slopes respond to the change in friction in the way expected from boundary layer theory (Figs. 12a,b). Away from the equator, where f ≫ r or , boundary layers on slopes become thinner and upslope velocities are enhanced as the friction parameter is increased. Near the equator, where f ≪ r or , boundary layers on slopes become thicker and upslope velocities are reduced. This compensation between boundary layer thickness and upslope velocities appears to result in the upslope transports being relatively insensitive to changes in . Across these solutions, western boundary currents carry about the same transport, and they have widths that scale with the Stommel boundary layer width , as expected.
As a consequence, the cross-equatorial overturning is also relatively insensitive to the value of friction (Fig. 12c). This weak dependence of the overturning is expected from the behavior of upslope transports in boundary layers on slopes, as confirmed by the contour integral in (35) (Fig. 12c). The friction dependence in the full solutions is even smaller than that predicted by boundary layer theory.
The results of this section suggest that our solutions are relatively insensitive to the value of friction, both in the qualitative structure of the circulation and in the net overturning. It should be kept in mind, however, that this may be particular to the Rayleigh friction used in these solutions. The control of the circulation by boundary layers on slopes suggests that a better understanding of the actual turbulent momentum transport in these boundary layers should be sought in future work. Preliminary results suggest that the qualitative nature of boundary layers on slopes is likely robust, but the magnitude of the up- and downwelling dipole can be quite sensitive to the choice of momentum closure (J. Callies 2018, unpublished manuscript). One-dimensional boundary layer solutions with a Fickian momentum flux closure produce much too weak a stratification compared to observations, which would give much too weak a dipole in buoyancy flux divergence and thus in up- and downwelling. This inconsistency argues for more complicated three-dimensional dynamics in the balances of boundary layers on slopes, which should be explored in future work.
The solutions to the planetary geostrophic equations in an idealized “bathtub geometry” with an idealized distribution of bottom-intensified mixing exhibit an abyssal circulation that is tightly controlled by diffusive boundary layers on slopes. The up- and downwelling in these boundary layers on slopes drives a basin-scale circulation and overturning because the inflow at the base of the sloping topography is denser than the outflow above. Zonal currents connect these in- and outflows at the base of the sloping topography to the western boundary, where Stommel boundary currents transport the water meridionally. A simple schematic of this circulation is shown in Fig. 8. Boundary layer theory captures the structure of the diffusive boundary layers on slopes, and integrating the upslope transport of these boundary layer solutions along the perimeter of the basin yields a prediction for the overturning. Our idealized solutions are instructive for understanding the abyssal circulation in the real ocean, but a number of complications must be considered to bridge that gap.
First, the representation of the Southern Ocean is obviously simplistic. This choice was motivated by the intention to keep the system as simple as possible and to direct the focus on the dynamics in the basin. This simplification defers the question of what sets the stratification at the northern edge of the Southern Ocean, as well as the question of whether the basin dynamics might affect that stratification. The point of view taken here is that this stratification is independent of the basin dynamics because Southern Ocean dynamics are much faster than the diffusive dynamics in the basins. But in reality, the water-mass transformation in the basin must be matched by surface transformation in the Southern Ocean, and the circulation and stratification must arrange themselves to satisfy that balance. These dynamics are excluded from the setup studied here.
Second, the real ocean’s bathymetry is much more complicated than the “bathtub geometry” considered here. There are midocean ridges, deep basins, fracture zones, and other geologic features. This complex bathymetry complicates the circulation considerably by steering the flow topographically and by changing where mixing is enhanced. For example, midocean ridges allow boundary currents along their western flanks; midocean ridges and seamounts can enhance mixing to middepth, and they may also affect the circulation differently than the slopes on the sides that come all the way to the surface. It seems likely that the circulation on many of these bathymetric features has characteristics captured by boundary layer theory, such that elements of the circulation described here carry over to the more complicated case. But whether a simple integration along depth contours can predict the overturning in more complex geometry remains unclear. The shape of bathymetry probably also affects the degree to which there is compensation between transformation to light and dense water (cf. McDougall 1989; McDougall and Ferrari 2017; Holmes et al. 2018), which may account for the larger degree of compensation estimated for the real ocean by Ferrari et al. (2016) and de Lavergne et al. (2017).
Third, the stratification we restore to in the south is unrealistically constant. The real stratification in the abyss is closer to exponential, which changes how much buoyancy flux convergence or divergence there is in the real ocean’s interior. That said, observations suggest that the buoyancy flux is divergent in the interior (Polzin et al. 1997; St. Laurent et al. 2001; Waterhouse et al. 2014), such that the dynamics should qualitatively be similar to the bottom-intensified case considered here.
Fourth, diabatic transformation in the real ocean appears to be confined to regions with strong flow over a rough bottom. It is thus expected that boundary layer flows are significant in such locations only, and that the basin-scale circulation is forced more heterogeneously than in the case with uniformly bottom-intensified mixing discussed here. The uniformly strong bottom value employed here is certainly unrealistic.
Fifth, the neglect of inertia, while drastically simplifying the dynamics, comes at the cost of preventing a whole host of physics, from recirculation gyres to time-dependent eddy dynamics. Particularly near the equator and in the western boundary currents, inertial effects may be important. The western boundary currents obtained here are forced to close the mass balance, similar to those postulated by Stommel and Arons (1960b) or those in linear gyre theory. But these boundary currents may be nothing more than caricatures of the boundary currents of the real ocean. Eddy stirring of the boundary layers on slopes may also be important.
Sixth, the physics of the boundary layers on slopes depend on the chosen form of friction and, more generally, on the physics of the small-scale flows. The boundary layers in our solutions assume a particularly simple form because Rayleigh friction is used, but they would be different in character if a Fickian closure for the momentum flux was used instead. The boundary layers are constrained to balancing the diffusive water-mass transformation with across-slope advection, but the degree to which the stratification and thus the diffusive fluxes are altered, depends on friction. The net overturning in our solutions is relatively insensitive to the value of friction, but that may be specific to Rayleigh friction. How the true time-dependent flows in bottom layers on rough slopes transport buoyancy and momentum is largely unclear. The control these bottom layers appear to exert on the abyssal circulation—they are not passive like western boundary currents in Stommel–Arons or linear gyre theory—suggests that more detailed understanding of the small-scale near-bottom dynamics is key to making progress.
We hope that despite these gaps to reality, the dynamics described here are instructive for understanding the abyssal circulation of the real ocean. It seems likely that elements of the solutions—in particular the control by boundary layers on slopes—survive in the full system. At the very least, the solutions illustrate that relaxing the assumptions made by Stommel and Arons (1960b) dramatically changes the abyssal circulation.
We thank Chris Garrett and Carl Wunsch for stimulating discussions. Support from the U.S. National Science Foundation through Grants OCE-1233832 and OCE-1736109 and from the National Aeronautics and Space Administration through Grant NNX16AH77G is gratefully acknowledged.
If restoring to a prescribed stratification in the south is present, an additional term representing the loss of potential energy in the restoring region appears on the right of (A2). In planetary geostrophic dynamics, any production of kinetic energy through is instantaneously balanced by friction; this balance is nonlocal, however, because pressure work redistributes kinetic energy. Potential energy is created by the diabatic term, lost to kinetic energy through buoyancy production, and redistributed advectively. Upon integration over the domain, indicated by angle brackets, we find if the bottom is insulating and the top is at . In the absence of wind forcing and restoring, this term represents the only source of energy for the system. It is always a source because in stable stratification. Even though diffusion does not generate kinetic energy directly, we thus speak of a diffusively driven circulation.
It should be noted that the fact that does not necessarily imply that a circulation is generated: if isopycnals are flat and remain flat, the potential energy generation is balanced by the tendency term—none of the generated potential energy is available for release (see section 4a). This situation, however, cannot occur in a steady state, and it cannot occur if the ocean bottom is sloping (see section 4b).
Planetary Geostrophic Equations in Terrain-Following Coordinates
a. Transformation to terrain-following coordinates
This appendix expresses (4)–(6) in terrain-following coordinates. We employ the notation of tensor calculus (e.g., Grinfeld 2013), which minimizes the amount of algebraic manipulation, helps identify the correct transformed components of the velocity vector, and easily yields expressions for buoyancy diffusion in flux form for arbitrary coordinates.
We begin by transforming the momentum, continuity, and buoyancy equations to the dimensional terrain-following coordinates:
In a second step, we will transform to the nondimensional coordinates in (19). One could do these transformations in one step, but the exposition is a bit clearer if this is done one step at a time.
The covariant metric tensor for terrain-following coordinates is
where hx and hy are the bottom slopes; everywhere else partial derivatives are written out. The contravariant metric tensor is the matrix inverse:
The volume element is the square root of the determinant of .
We now write the momentum equation (4) in covariant form. Noting that the contravariant components of the vertical unit vector are , , and yields the covariant components of the cross product :
where ε is the Levi–Civita symbol. The contravariant velocity components transform as
and the covariant components of friction become
Furthermore using , , and for the buoyancy term, the covariant components of the momentum equation (4) become
where we also used that the components of the covariant derivative of a scalar, in this case pressure, are simply its partial derivatives. These equations could also have been derived from the Cartesian form of the equations using the chain rule. But the true power of tensor calculus becomes apparent when writing the continuity equation (5) and the buoyancy equation (6) in terrain-following coordinates. We make use of the Voss–Weyl formula for the divergence of a vector field A:
The continuity equation (5) thus simply reads
In the buoyancy equation (6), the advective term, written as a flux divergence, becomes
The diffusive term has the same form of a flux divergence, with the diffusive fluxes being
which follows from .
The main drawbacks of the terrain-following coordinates are that buoyancy terms appear in the momentum equations in (B6) and (B7) because and are nonzero and that cross terms are introduced in the diffusive fluxes in (B12)–(B14) because the metric tensor is not diagonal. These drawbacks are balanced by the simplicity of the boundary conditions: the no-normal flow condition simply becomes at and , and the no-flux condition on buoyancy becomes at . Both of these conditions are easily implemented in our finite-difference scheme because these fluxes are located at σ faces.
b. Depth-integrated flow and thermal wind
To calculate the depth-integrated flow, we solve a two-dimensional elliptic problem for the streamfunction of the integrated flow ψ:
The right-hand side of (B16) is the joint effect of baroclinicity and relief (JEBAR) term. Given a buoyancy field b and on all coasts, (B16) can be solved for ψ. If there was wind forcing or other imposed upwelling, additional forcing terms would appear in (B16). If there were no stratification () or no bathymetry (), the JEBAR term would vanish and the depth-integrated flow could be determined from (B16) alone (Welander 1968). In the more general case considered here, the buoyancy field must be known.
We must add the baroclinic component to the depth-independent flow obtained from (B16). The baroclinic component is obtained from the frictionally modified thermal wind equations:
c. Nondimensional equations in terrain-following coordinates
For completeness, we now list the equations in the nondimensional terrain-following coordinates defined by the transformation
With and the substitutions in (10), we have
and define consistently
The depth-integrated flow
is then obtained from
Thermal wind reads
and the continuity equation becomes
The nondimensional buoyancy equation is
where the advective term is simply
and the diffusive fluxes are
Nondimensional Boundary Layer Equations
We include some detail of the transformation to a slope-aligned coordinate system in nondimensional form because there are a few quirks introduced by the squeezing of the system in the horizontal. This means that angles are not preserved, and a geometric factor arises in the nondimensional equations.
The transformation is
and its inverse is
where is the slope angle in the nondimensional coordinates and is related to the true slope angle θ through . This transformation yields the covariant metric tensor
and the contravariant metric tensor
The volume element is . The vertical unit vector now has contravariant components
and covariant components
The covariant components of friction are
and the momentum equations, making use of the substitutions (10), become
Splitting buoyancy and pressure into a background state and a perturbation,
and assuming no variations along the slope, , , and then simplifies the across- and along-slope momentum equations to
and the buoyancy equation becomes
Substituting from the momentum equations yields the nondimensional version of (26):
The insulating boundary condition at becomes
Note the appearance of the geometric factor
which results from the horizontal squeezing of the domain when going to nondimensional coordinates.
The analytical solution for uniform mixing now reads
For bottom-intensified mixing,
and numerical solutions are obtained by solving the steady-state version of (C18) using finite differences.
The unconventional notation for the velocity components is common in tensor calculus, which we will employ in earnest for coordinate transformations below.
In the case of microstructure measurements, what is inferred from the shear measurements is the kinetic energy dissipation rate ε, which is related to the turbulent buoyancy flux through a flux ratio Γ. The balance in (3) then becomes The available observations of ε and estimates of the flux ratio (and its variations; e.g., Gregg et al. 2018) strongly suggest the right-hand side is negative over the layer of enhanced mixing near the rough ocean floor.
As discussed above, the set (4)–(6) with horizontal Rayleigh friction (7) generally cannot satisfy no-normal flow boundary conditions at vertical sidewalls because the system cannot prevent thermal-wind shear at the walls. This violation of boundary conditions does not occur in the special case considered here because there is no flow whatsoever.
This also suggests that a linear solution may be a reasonable approximation to the full nonlinear solution. The extent to which this is true is not further explored here.