Coupling Between Abyssal Boundary Layers and the Interior Ocean in the Absence of Along-Slope Variations

To close the overturning circulation, dense bottom water must upwell via turbulent mixing. Recent studies have identified thin bottom boundary layers (BLs) as locations of intense upwelling, yet it remains unclear how they interact with and shape the large-scale circulation of the abyssal ocean. The current understanding of this BL--interior coupling is shaped by 1D theory, suggesting that variations in locally produced BL transport generate exchange with the interior and thus a global circulation. Until now, however, this picture has been based on a 1D theory that fails to capture the local evolution in even highly idealized 2D geometries. The present work applies BL theory to revised 1D dynamics, which more naturally generalizes to two and three dimensions. The BL is assumed to be in quasi-equilibrium between the upwelling of dense water and the convergence of downward buoyancy fluxes. The BL transport, for which explicit formulae are presented, exerts an influence on the interior by modifying the bottom boundary condition. In 1D, this BL transport is independent of the interior evolution, but in 2D the BL and interior are fully coupled. Once interior variables and the bottom slope are allowed to vary in the horizontal, the resulting convergences and divergences in the BL transport exchange mass with the interior. This framework allows for the analysis of previously inaccessible problems such as the BL--interior coupling in the presence of an exponential interior stratification, laying the foundation for developing a full theory for the abyssal circulation.


Introduction
Thin boundary layers (BLs) at the ocean's bottom have recently come into focus as the primary locations in which small-scale turbulence lightens bottom waters, thus playing a crucial role in closing the overturning circulation of the abyss (Ferrari et al. 2016;de Lavergne et al. 2016). The connection between these BLs and the large-scale abyssal circulation, however, remains to be fully explained. The cornerstone of our present understanding of the mixing-generated abyssal circulation is a 1D model of a stratified, rotating fluid overlying a sloping, insulated seafloor (e.g. , Phillips 1970;Wunsch 1970;Thorpe 1987;Garrett et al. 1993). This 1D theory helped bring bottom BLs into center stage, predicting that the local response to bottom-intensified mixing is characterized by diabatic upslope flow in the thin BL compensated in part by diabatic downslope flow spread across the interior (Garrett 1990;Ferrari et al. 2016;de Lavergne et al. 2016;McDougall and Ferrari 2017;Callies 2018). Our description of large-scale abyssal dynamics is shaped by this local theory: the natural conclusion is that variations in these locally produced flows generate exchange with the interior and producing a global circulation (e.g., Phillips et al. 1986;McDougall 1989;Garrett 1991;Dell and Pratt 2015;Holmes et al. 2018). This picture fails to consider the potential feedback of the circulation produced in the interior back onto the BL, however, suggesting that this framework is incomplete.
In addition to this lack of two-way coupling, progress has also been hampered by the canonical 1D theory failing to reproduce the local evolution in simple 2D geometries. The canonical 1D model predicts slow diffusion of the interior along-slope flow (MacCready and Rhines 1991), whereas simulations of bottom-intensified mixing over an idealized 2D mid-ocean ridge display rapid spin up of the interior (Ruan and Callies 2020). In Peterson and Callies (2022, hereafter PC22), we remedied this shortcoming by including the physics of a secondary circulation and barotropic pressure gradient. The key is to constrain the vertically integrated cross-slope transport to force upwelling flow in the BL to return in the interior. This downwelling flow is then turned in the along-slope direction by the Coriolis acceleration and balanced by a barotropic pressure gradient, leading to rapid adjustment in the interior as seen in 2D. With this more faithful 1D model, we have a reliable foundation to describe the role of abyssal BLs in the large-scale circulation. Callies and Ferrari (2018) and Drake et al. (2020) connected BL dynamics to the horizontal circulation in a 3D planetary-geostrophic (PG) model with idealized bathymetry and Rayleigh friction. Callies and Ferrari (2018) found that, for vertically constant interior stratification and on moderate slopes, local 1D theory accurately emulates the 3D model's dynamics. On the sloping sidewalls of the idealized bathymetry, upslope transport in thin bottom BLs is compensated by downwelling aloft. At the base of the slopes, however, 1D theory breaks down in favor of a basin-scale circulation that feeds the BLs on slopes. An integral of the local upslope 1D BL transport along the perimeter of the basin provides an accurate estimate of the overturning. These ideas fail, however, once the interior stratification is far from constant, because 1D theory can only consider perturbations to a constant background stratification (Drake et al. 2020). This is a severe limitation, given the real ocean's near-exponential stratification (e.g., Munk 1966). For a more realistic stratification, downwelling in the interior is weakened and BL upwelling dominates, though the vertical extent and structure of the net transport remains to be explained. In this work, we provide a framework for concretely understanding this interplay between the BL and interior.
Below, we derive self-contained equations for interior 1D and 2D PG dynamics on an -plane with effective boundary conditions that capture the effects of BLs. We accomplish this using BL theory, splitting variables into their interior and BL contributions (e.g., Bender and Orszag 1999;Chang 2007, Fig. 1). This explicitly separates the interior and BL dynamics and allows for deep physical insight into their coupling. Famously, Stommel's (1948) gyre theory can be solved with BL methods (Veronis 1966), although the coupling there is one-way: the interior solution can be calculated in isolation, and the western BL is a passive element of the theory. We find that this is different for bottom BLs on slopes. Their structure is shaped by the interior solution, but the buoyancy and mass fluxes carried in the BL feed back on the interior solution in the form of boundary conditions. A central result of this paper is an explicit expression for the cross-slope BL transport (per unit along-slope distance) in terms of interior variables and flow parameters. In 1D, the BL transport takes the form cot /(1 + ), where = / is the turbulent Prandtl number with being the turbulent viscosity and the turbulent diffusivity, and = 2 tan 2 / 2 is the slope Burger number with being the background interior buoyancy frequency, the inertial frequency, and the bottom slope angle. All variables are evaluated at the bottom (or, more generally, just above the BL). In the canonical 1D framework, a steady-state balance between cross-slope upwelling of dense water and turbulent mixing requires that the total transport tends towards ∞ cot , where ∞ is the Interior BL F . 1. Illustration of the BL correction to interior solution. Shown is a typical streamfunction , defined such that = where is the cross-slope flow, after three years of mixing-generated abyssal spin up at a slope Burger number = 10 −3 (see section 2). The solution is depicted over (a) the entire 2 km domain as well as (b) a zoom-in to the bottom 100 m, shown in (a) in gray shading. The interior solution I varies slowly compared with the scale of the BL, and the BL correction B ensures that boundary conditions are satisfied.
far-field turbulent diffusivity (Thorpe 1987;Garrett et al. 1993). Our revised result instead applies to the bottom BL transport and is valid throughout transient evolution, provided that the BL has adjusted to a quasi-steady state. Unlike the canonical result, this expression smoothly approaches zero as → 0, more harmoniously connecting the model over a slope with conventional flatbottom Ekman theory (e.g., Pedlosky 1979). The expression has the same form in 2D, but there the slope Burger number is a function of interior cross-isobath buoyancy gradients as well as the local topographic slope. Thus, in 2D, variations in interior buoyancy gradients and the topographic slope cause convergence in the BL transport, generating exchange with the interior. A similar process occurs in 3D with the added physics of along-isobath variations and a modified interior balance, but we leave the details of 3D dynamics to future work.
In section 2, we begin by reviewing the transport-constrained 1D model from PC22, followed by a derivation of the 1D BL theory. We derive the 2D BL theory in section 3, applying the framework to simulations of mixing-generated spin up under a vertically varying background stratification. In section 4, we re-derive the 1D and 2D BL equations in a more rigorous fashion, quantifying the accuracy of our claims in the previous sections and uncovering some subtleties in the dynamics.
Finally, we provide discussion and conclusions in sections 5 and 6, respectively.

One-dimensional boundary layer theory
In this section, we apply BL theory to the revised 1D model from PC22 and present results from numerical integrations of both the full and BL equations. Here and throughout the paper, we employ PG scaling, thus focusing our attention on the slow and large-scale response to mixing. The PG flow should be interpreted as the residual flow after a thickness-weighted average over transients due to turbulence, waves, and baroclinic eddies, with the effect of these transients included as parameterized Eliassen-Palm and diapycnal fluxes (?).

a. Transport-constrained one-dimensional dynamics
We first consider 1D PG dynamics along a uniform slope at an angle above the horizontal.
The 1D model is typically derived by writing the Boussinesq equations in a rotated coordinate system aligned with the slope (e.g., Garrett et al. 1993). We slightly deviate from this approach by keeping the vertical coordinate aligned with gravity, which is a more natural choice if the horizontal components of the turbulent momentum and buoyancy fluxes are neglected, but it yields equivalent dynamics (PC22). Specifically, we write the 1D model in ( , , ) coordinates defined by where ( , , ) defines the usual Cartesian coordinate system with aligned with gravity. These coordinates are analogous to terrain-following coordinates (used below) in 1D with = 0 at the bottom. Neglecting all variations in and , except for the barotropic pressure gradient (equivalently, , since is independent of ), and constraining the vertically integrated cross-In the limit 1, the gravity-aligned coordinate system employed here and the previously used fully rotated coordinate system yield the same equations. slope transport to (typically to zero), the PG equations become Here, is the cross-slope velocity and is the along-slope velocity. We have split the total buoyancy into a constant background stratification and a perturbation so that = 2 + .
The fluid satisfies no-slip and insulating boundary conditions at the bottom: = 0, = 0, and = 2 + = 0 at = 0. In the far field, we impose decay conditions on the shear and anomalous buoyancy flux: → 0, → 0, and → 0 as → ∞. The extra degree of freedom supplied by allows the transport constraint (5) to be satisfied at all times. Physically, this constraint forces cross-slope upwelling in the BL to return in the interior, where it is then turned into the along-slope direction by the Coriolis force. In the PG framework, this process is instantaneous, and the far-field along-slope flow satisfies the balance: − = − . This leads to rapid spin up of the along-slope flow throughout the water column, as seen in simulations of 2D spin up (Ruan and Callies 2020, PC22).
We employ a simple down-gradient closure for the turbulent momentum and buoyancy fluxes generated by, e.g., breaking internal waves but allow for variations in the mixing coefficients and . We assume these variations to occur on a scale larger than the BL thickness. In our examples below, and are bottom-enhanced in abyssal mixing layers a few hundred meters thick, inspired by typical observations over rough mid-ocean ridges. Our main results, however, generalize to the case in which and vary rapidly within the BL, for example going to zero in a log-layer.
As in PC22, we cast equations (2) to (5) into an inversion equation for the flow, written in terms of a streamfunction ( ) defined such that = , and an evolution equation for the buoyancy Due to our non-orthogonal coordinate system, is technically the -projection of the cross-slope velocity as it would be defined in a fully rotated coordinate system (PC22, appendix A). For simplicity, we refer to it as the "cross-slope velocity" throughout. dynamics over more complicated topography. perturbation: The boundary conditions are that = 0 and = 0 at = 0 and → as → ∞. If desired, one may infer the along-slope flow from by integrating from the bottom up, using = 0 at = 0. Equations (6) and (7) fully describe the 1D PG system and can readily be solved numerically. But insight into the BL-interior coupling is more easily gained using BL theory.

b. Boundary layer theory
Under steady conditions, equations (6) and (7) can be combined to form a single fourth-order ordinary differential equation for . The fourth-and zeroth-order terms in that equation balance if varies on a scale −1 defined by where = √︁ 2 / is the familiar flat-bottom Ekman layer thickness, and the mixing coefficients are evaluated at = 0. This defines the BL scale of a rotating fluid adjacent to a sloping bottom (e.g., Garrett et al. 1993). For typical abyssal parameters, −1 ∼ 10 m (Callies 2018). This thinness of the BL compared to the scale of variations in the interior ocean is what allows us to apply BL theory.
We begin by splitting solutions into interior contributions I and I , which vary slowly in , and BL corrections B and B , which ensure boundary conditions are satisfied and have appreciable magnitude in the thin BL only. A similar approach was taken in Callies (2018) with the canonical 1D model, but the analysis presented here is time-dependent and extensible to higher dimensions (section 3). If the mixing coefficients and vary on a scale much larger than −1 , the fourth-order term in (6) can be neglected in the interior: assuming = 0 (see appendix A for the ≠ 0 case). Substituted back into the buoyancy equation (7), this reduces the interior dynamics to a modified diffusion equation: This is a result familiar from Gill (1981), Garrett and Loder (1981), and Garrett (1982): advection of the background stratification by the secondary circulation becomes a horizontal diffusion term, with diffusivity 2 / 2 . The form here is the result of the sloping boundary: the vertical coordinate depends on the slope-parallel distance multiplied by tan , which explains the factor tan 2 in the additional diffusion term.
This interior evolution must be complemented by a representation of the bottom BL that supplies an effective boundary condition for the interior equation. The key assumption here is that the BL scale −1 is thin compared to interior variations. This thinness of the BL also implies that it is quasi-steady on the time scales of the interior evolution. The BL correction thus satisfies the steady buoyancy equation Since all BL variables decay into the interior, i.e., as → ∞, this balance can be integrated to This relation is all that is needed to derive a boundary condition on the interior solution. At = 0, I + B = 0, such that the full = 0 boundary condition is satisfied. So, using (10), The insulating boundary condition then becomes The BL correction thus contributes an additional term I to the boundary condition for the interior buoyancy evolution (11). The added term represents physics akin to an Ekman buoyancy flux (e.g., Marshall and Nurser 1992;Thomas and Lee 2005): the BL transport I acts on the crossslope buoyancy gradient 2 tan and produces a buoyancy sink for the interior. This boundary condition on the interior problem implies a stratification at the top of the BL that is reduced from the background by a factor /(1 + ) and a BL transport, from combining (15) and (10), as claimed in the introduction (Fig. 2a). We note that the transport-constrained system, unlike the canonical one, has no steady state in a semi-infinite domain, yet previous work on the BL-interior interaction has often begun with the canonical result that the steady transport is = ∞ cot (e.g., Woods 1991; Callies and Ferrari 2018;Drake et al. 2020). The revised expression in (16) instead applies to the transport confined to the BL and more sensibly leaves the net transport (and steady-state dynamics) to be controlled by the large-scale context.
If desired, the BL correction can easily be determined from with B = − I and B = 0 at = 0 (neglecting the much smaller interior contribution to at the bottom) and B → 0 as → ∞. This has a similar form as the steady canonical 1D problem with constant mixing coefficients (e.g., Garrett et al. 1993), but the boundary conditions and righthand side are different because the transport constraint is imposed and the interior solution has been subtracted out. The general solution takes the form of the familiar Ekman spiral: where I is evaluated at = 0 as in (16).
This analytical expression for the BL correction also allows us to directly diagnose how the far-field along-slope flow is influenced by the BL. From (8) and (10), the interior along-slope shear follows thermal wind balance, which implies, upon integration in the vertical, The integration constant I (0), the flow at the upper edge of the BL, can be determined from the BL solution (18) and (8): . This BL contribution to the interior along-slope flow has the same form as the steady-state canonical result with constant mixing coefficients (Thorpe 1987;Garrett et al. 1993), but here it is rapidly spun up and accompanied by an additional interior thermal-wind component. We will see in section 4 that this BL contribution is typically of higher asymptotic order than the thermal-wind contribution.
It should be noted that the key results (15) and (16) also apply if there are variations in the mixing coefficients within the thin BL, as may be expected as the turbulence becomes suppressed very close to the bottom. The physics that lead to (15) and (16) are that the diffusive buoyancy flux into the BL is balanced by cross-slope advection within the BL and that the interior obeys (10).
While the BL corrections are more complicated if and are not approximately constant across the BL, for example including a log-layer if the mixing coefficients go to zero near the bottom, the effective boundary condition for the interior is the same.
In summary, BL theory has enabled us to elucidate the connection between the BL and interior in 1D. The BL transport quickly adjusts to (16), regardless of the interior dynamics. This transport allows the BL to communicate with the interior by moving dense water up the slope, providing a buoyancy sink and modifying the interior bottom boundary condition (15) (Fig 2a). In 1D, the BL is thus independent of the evolution of the interior, yet the cross-slope advection by the BL transport affects the interior dynamics. As we will see in the next section, the BL-interior coupling in 2D are even richer, with the interior being able to feed back onto the BL. But first, we present some illustrative 1D examples.

c. Examples
The following experiments depict 1D PG spin up with and without BL theory. The simulations start in a state of rest: isopycnals are flat ( = 0), and the flow is zero ( = 0). The turbulent mixing then generates a buoyancy perturbation, bending isopycnals into the slope and spinning up a circulation. The transport constraint ensures that BL transport is exactly returned in the interior, and without a source of dense bottom water, the initial stratification is mixed away with time.
To numerically solve the 1D PG equations, we use second-order finite differences as in PC22.
The model can either solve for the full flow and density profiles using equations (6) and (7) or evolve the interior variables of the BL theory with equation (11). Model parameters are adapted from Callies (2018) and roughly match those of the Brazil Basin (Table 1). Mixing is represented by a bottom-intensified profile of turbulent diffusivity, Inertial frequency −5.5 × 10 −5 s −1 Far-field buoyancy frequency 10 −3 s −1 Far-field diffusivity 0 6 × 10 −5 m 2 s −1 Bottom-enhancement of diffusivity 3 years 6 years 9 years 12 years 15 years The 1D BL model yields an excellent approximation of the full 1D PG solution (Fig. 3). The interior dynamics match the interior of the full solution, and although the BL model only explicitly computes the interior evolution, the BL correction computed offline from (18) is very accurate.
The match is trivial when = 1 and = 10 −3 , because the shallow slope leads to a relatively weak BL transport, and thus the advective modification to the buoyancy flux in (11) and (15) is negligible. The interior system is then nearly identical to the full one, with diffusion dominating the dynamics. The case where = 0.5, in contrast, is a more trying test of the 1D BL theory.
The BL transport in this case is an order of magnitude larger than before, leading to enhanced stratification in the BL. This is properly captured in the BL model, with the interior stratification reaching about 0.4 × 10 −6 s −2 at the bottom and the BL correction bringing it smoothly to zero.
The assumption of 1D dynamics breaks down as soon as lateral variations in the slope are allowed, but we can anticipate the upcoming 2D results using intuition derived from the above 1D theory. Equation (16) gives an explicit expression for the BL transport in 1D depending on the local slope angle and buoyancy gradient across the slope 2 tan . In 2D, these inputs are spatially dependent, with horizontal buoyancy gradients also varying in time as part of the interior dynamics. Local 1D theory would thus predict convergences and divergences in BL transport, generating BL-interior mass exchange (Fig. 2b). This leads to a more complex picture in 2D, with interior dynamics feeding back onto the BL, as we will see in the following section.

Two-dimensional boundary layer theory
In this section, we extend the 1D BL theory to the 2D PG equations in terrain-following coordinates.
We first derive the 2D BL equations and then apply them to idealized numerical simulations.

a. Boundary layer theory
In 2D, the interaction between the BL and interior is more interesting because, in addition to the BL advection imposing a buoyancy flux on the interior, variations in the BL transport produce mass exchange with the interior (e.g., Phillips et al. 1986;McDougall 1989;Kunze et al. 2012;Dell and Pratt 2015;Ledwell 2018;Holmes et al. 2018). The BL theory generalizes from 1D to 2D and brings these physics into clearer focus.
Applying the BL theory to the 2D PG equations is most easily done in terrain-following coordinates: where ( ) is the fluid depth (Fig. 4). Under this transformation, derivatives in ( , ) space become and the contravariant velocity components are assuming no variations in (see appendix B of Callies and Ferrari (2018) for more details). The 2D PG equations in terrain-following coordinates are then where is the pressure divided by a reference density. The boundary conditions are again an insulating and no-slip bottom, = 0 and = = 0 at = −1; a constant-flux and free-slip top −1 = 2 and = = 0 at = 0; and no normal flow across both boundaries, = 0 at = −1 and = 0. We neglect horizontal turbulent fluxes, consistent with the assumption of a small aspect ratio if the turbulence is close to isotropic. This is in contrast with some other PG models, which employed horizontal diffusion terms to satisfy the no-normal-flow condition at vertical side-walls (e.g., ??).
As before, we express the momentum equations (25) to (28) as one streamfunction inversion.
We define ( , ) such that the continuity equation (28) is automatically satisfied: Integrating (26) from some level to = 0, we obtain as in equation (8). Here, is the vertically integrated transport, a constant in by continuity. Combining −1 of (25) and of (27) The boundary conditions are similar to the 1D case but for a finite domain: = 0 and = 0 at = −1 and = and 2 = 0 at = 0.
Splitting and into BL and interior contributions and neglecting the fourth-order term in (32) in the interior as before, the interior inversion reads The BL physics appear in the boundary condition on the interior buoyancy field. The BL buoyancy budget, assuming a quasi-steady state and a slowly varying interior buoyancy field, is with I evaluated at = −1. The neglected advection terms are smaller by a factor ( ) −1 1 than the terms retained in (35). This is because the boundary conditions enforce that B ∼ I and B ∼ I , such that B ∼ ( ) I and B ∼ ( ) −1 I (see section 4 for more detail). Vertically integrating (35) Substituting this bottom boundary condition for the interior into the interior inversion (33), we again arrive at an explicit formula for this BL transport: This is the generalization of the 1D result (16): − is analogous to the local slope tan and I now takes the place of the previously constant cross-slope buoyancy gradient 2 tan . Note that this expression is again well-behaved in the limit of small slopes ( → 0) and thus gives a globally valid expression for the BL transport and of the mass exchange I = − I at = −1 between the BL and the interior.
As in 1D, we can now explicitly describe contributions to the interior along-slope flow from thermal wind in the interior and a contribution from shear in the BL. Combining (31) and (33) yields the thermal-wind balance which, upon integration in the vertical, becomes The first term again represents the BL contribution I = − B at = −1, which may be computed directly from the BL solution similar to (18). Here can still be written in the same form as in (9) but with a generalized slope Burger number = − I (−1)/ 2 , which varies in the horizontal. Equation (39) has the same form as (20), except that cross-slope buoyancy gradients can now contribute to the thermal-wind term.
In 2D, we again find that the interior solution experiences a buoyancy flux due to the crossslope advection by the BL transport. In contrast to the 1D case, however, both the BL transport given by (37) and the cross-slope buoyancy gradient I may vary in time and space (Fig. 2b). Convergence in the BL transport then drives mass injection into the interior, further altering I and continuing the feedback process.
It is worth noting that BL theory can also be applied to a passive tracer, not just buoyancy. The interior tracer concentration would have a similar effective boundary condition capturing transport by BL flow. The interior tracer equation should also include a representation of along-isopycnal stirring (Redi 1982).

b. Examples
We now illustrate these theoretical results using numerical simulations over idealized topographies.
We solve the full 2D PG system (29) and (32) and the 2D BL PG system (33) and (34) using numerical methods and model parameters similar to the 1D case described above. The mixing profile is now written as Similar to the analysis in Ledwell (2018), we consider an azimuthally symmetric Gaussian seamount in axisymmetric coordinates (Fig. 5). On an -plane, the flow is invariant under rotation about the center of the seamount, allowing us to fully describe the flow using 2D theory (see appendix B). The depth of the seafloor as a function of distance from the symmetry axis is given by where the maximum depth is 0 = 5.5 km, the height of the seamount is = 3 km, the width of the seamount is ℓ = 50 km, and the width of the domain is = 200 km. We assume no flow at = 0 and allow the flow to evolve freely at = , consistent with our assumption that horizontal diffusion may be neglected. In the horizontal, the grid has an even spacing of about 0.8 km. As in the 1D models, we use Chebyshev nodes in the vertical when solving the full 2D PG equations (with a near-bottom resolution of about 10 −5 in -space) and uniform grid spacing for the 2D BL equations (with a resolution of about 10 −3 in -space). We initialize the model at rest with a constant stratification = 2 and use a mixed implicit-explicit time integration scheme with a timestep of one day.
At the steepest point on the seamount ( = 50 km, red lines in Fig. 5), the slope Burger number is order unity. The 1D BL solution applied at this position over-predicts the stratification in the bottom 500 m and under-predicts it above (Fig. 6). This leads to errors in the predicted interior along-slope flow, which can be understood from (20) and (39): even subtle changes in the buoyancy field can lead to substantial impacts on I after being integrated throughout the column. The 1D BL solution's buoyancy field differs from that of the 2D solution because its secondary circulation, enforced simply by a transport constraint, is stronger. This is due to the lack of a two-way feedback in 1D; the BL cannot exchange mass with the interior and the induced changes in the interior do not reduce the BL transport. The 2D BL theory, in contrast, captures these physics and agrees well with the full 2D model. This confirms that the 2D BL equations are capable of fully capturing 2D PG spin up, even in regimes with relatively large variations in local slope.

2) E
The simulations presented so far were initialized with a constant background stratification. In the real ocean, the stratification varies significantly in the vertical, often decreasing close to exponentially with depth (e.g., Munk 1966). A number of studies have attempted to discern how this may shape the abyssal circulation, often qualitatively arguing that variations in stratification across slopes must lead to gradients in BL transports, inducing BL-interior exchange (e.g., Phillips et al. 1986;Salmun et al. 1991). Quantitative explanations of this process, however, have remained complicated and opaque at best. A major benefit of the BL theory framework built up here is that it provides concise expressions for the BL transport in terms of interior variables, allowing us to reason about how varying background stratification might impact the abyss with minimal mathematical gymnastics.
Let us consider an idealized mid-Atlantic ridge, following previous studies of mixing-generated spin up in the abyss (e.g., Ruan and Callies 2020;Drake et al. 2020, PC22). The depth of the 2D ridge is given by where the mean depth is 0 = 2 km, the amplitude is = 800 m, and the width is = 2000 km ( Fig. 7). At the steepest point on the ridge, the slope Burger number is approximately 2 × 10 −3 , typical of the mid-Atlantic ridge. We apply periodic boundary conditions at = 0 and = and use a constant horizontal grid spacing of about 8 km. The vertical grid spacing is as before. We run one simulation with constant initial stratification as before and one initialized with an exponential stratification: ∝ / . We set the decay scale to = 1000 m and choose the proportionality constant such that the bottom stratification at the center of the ridge flank matches that of the simulation with constant 2 = 10 −6 s −2 . We again use a mixed implicit-explicit timestepping scheme, this time with a timestep of 10 days, enabled by the much weaker advective terms.
The circulation in the case with exponential initial stratification is stronger and more confined to the peak of the ridge compared to the case with constant initial stratification (Fig. 7a,b). This is better understood by the explicit formula for 2D BL transport derived in the previous subsection.
Evaluating equation (37) for these simulations, we see that the BL transport is enhanced at the peak of the ridge with exponential background stratification (Fig. 7c). For the small slopes in this simulation, equation (37) reduces to In the case with constant stratification, the initial cross-slope buoyancy gradient is proportional to − and does not change appreciably with time, explaining the sinusoidal BL transport. For exponential stratification, in contrast, we have I ∝ − − / , which is enhanced at shallower depths. As a result, the exchange velocity is also enhanced for the case with exponential stratification (Fig. 7d). In both cases, I does not evolve much in the first three years, so the exchange does not either. The BL theory enables us to easily and quantitatively understand this behavior.

Asymptotic theory
In the previous sections, we derived the BL equations somewhat heuristically, glossing over some detail of the underlying asymptotics. In this section, we present a more rigorous derivation of the BL theory that justifies the claims in the previous sections and sheds light on the asymptotic orders of the various components of the flow. The casual reader should note that the contents of this section are not required to understand the main results of the paper.
We show below that, in both 1D and 2D, the cross-slope flow is of lower order than the alongslope flow in the interior, aligning with our intuition from the examples above. The interior flow evolves on a slow timescale driven by diffusion and second-order advection of the leading-order buoyancy in the interior. The BL flow is of first order, in between the orders of the interior alongand cross-slope flows. If the transport is constrained to zero, this implies that the leading-order interior flow vanishes at the bottom. These results do not generally hold in 3D, but we leave this generalization to future work.

a. One-dimensional asymptotics
To begin the formal derivation of the 1D BL equations, we first nondimensionalize the 1D equations (2)-(5) in order to isolate the key parameters in the problem. We define characteristic scales for the vertical coordinate, velocities, and mixing coefficients such that where 0 and 0 are characteristic values of and . We assume that the pressure and buoyancy terms in (2) scale with the Coriolis term and that the buoyancy perturbation scales with the background buoyancy scale: Assuming an advective timescale, so that then yields the nondimensional 1D equations where all variables are redefined to their scaled versions. The nondimensional parameters for the 1D problem are thus the Ekman number 2 = 0 / 2 0 , the Prandtl number = 0 / 0 , and the slope Burger number = 2 tan 2 / 2 , although and only appear as a product, so can be considered a single parameter. The reason for defining the Ekman number as 2 will become clear in the BL analysis below.
To develop the asymptotic theory, we assume the scaling 1 and ∼ 1. While the Burger number is typically small in the abyss, the turbulent Prandtl number may be large if momentum fluxes by baroclinic eddies are taken into account. If instead 1, buoyancy advection is negligible in the BL, and the theory developed with ∼ 1 remains accurate (Fig. 3a).
We begin with the interior and expand all variables in 2 : I = I0 + 2 I2 + . . . , etc. This expansion into even powers of is sufficient because only appears as 2 in the interior equations.
The (1) interior flow then satisfies At this order, the interior along-slope flow is in balance with the barotropic pressure gradient and the projection of the buoyancy perturbation, and the interior cross-slope flow is zero. The (1) buoyancy equation is trivial, implying that the interior buoyancy evolution is slow compared to the advective timescale assumed in the scaling.
To obtain the evolution of the (1) interior buoyancy, we need to go to ( 2 ) and also expand the time coordinate, = 0 + 2 2 + . . . Higher-order buoyancy terms inherit the slow evolution from the low orders, so 0 I2 = 0. The buoyancy equation (51) at ( 2 ) is then This implies that advection and turbulent diffusion operate on a slow time 2 . Since the (1) and ( ) interior cross-slope flows are zero, the dominant buoyancy advection is by the second-order flow in the interior, given by (50) at ( 2 ): Equations (53), (56), and (57) comprise the leading-order interior dynamics. They can be expressed in terms of the streamfunction I , whose leading non-zero component is I2 , recovering (10) and (11) above (assuming = 0). The interior along-slope flow can be obtained by integrating the thermal-wind balance I0 = − I0 , which follows from a -derivative of (53): The integration constant I0 (0) must be determined from the BL correction. If the transport constraint is = 0, one finds that I0 (0) = 0. In the thin bottom BL, -derivatives are enhanced, elevating the diffusion terms in (49)-(51) to (1). Given that the BL thickness scales with , we assume the BL variables to depend on the re-scaled vertical coordinate¯= / , with which = −1¯. The nondimensional BL equations are then Crucially, the insulating bottom boundary condition picks up a factor of −1 after this re-scaling: This factor of −1 means that we need an ( ) BL buoyancy to absorb the (1) interior buoyancy flux into the BL. We thus expand the BL variables in terms of rather than 2 . We immediately find that the (1) BL buoyancy flux must vanish at the bottom:¯ B0 = 0. In the case with zero net transport ( = 0), this condition, along with the boundary conditions on the flow B0 = 0 and B0 = − I0 at¯= 0, forces the (1) BL flow to vanish and the (1) interior along-slope flow to go to zero at the bottom, consistent with the examples shown in Fig. 3 (see appendix A for the ≠ 0 case). The BL flow instead comes in at ( ), in between the orders of the interior alongand cross-slope flows. This ( ) BL flow satisfies with the bottom boundary conditions¯ B1 = −(1 + I0 ), B1 = 0, and B1 = 0. The tendency term 0 B1 is dropped because the interior does not evolve on this timescale, so the BL will not either. These BL equations are equivalent to (12) and (17).
This more rigorous derivation of the 1D BL equations clarifies the asymptotic orders of the various components of the flow. The leading-order contributions are ( 2 ) for the interior crossslope flow, (1) for the interior along-slope flow, and ( ) for both components of the BL flow.
Buoyancy does not have an (1) BL correction-only its derivative does.

b. Two-dimensional asymptotics
The 2D asymptotics follow in much the same way as in 1D. We again nondimensionalize the equations of motion (25)-(29), setting characteristic scales equivalent to (46)-(48): We then arrive at the nondimensional 2D PG equations where = 2 2 0 / 2 2 is now the conventional Burger number. Again assuming the scaling 1 and ∼ 1, expanding interior variables in 2 , and matching orders as before, we arrive at the complete set of interior equations We again find that the interior along-slope flow is of lower order than the interior cross-slope flow, and the interior buoyancy evolution is again slow. In 2D, the interior slope-normal flow I2 comes in, contributing a second-order advective flux in the vertical, along with the cross-slope advection. The BL contribution can again be assessed after a re-scaling of the vertical coordinate such that¯= / . We again find that the (1) BL flow, along with the interior along-slope flow I0 at the bottom, vanishes when = 0. The BL flow is instead of ( ), satisfying with hydrostatic balance and continuity implying that B1 = 0 and B1 = 0, respectively. The BL is again characterized by a balance between cross-slope advection and down-gradient diffusion of buoyancy, with the BL buoyancy flux due to B1 balancing the interior buoyancy flux due to I0 at the bottom as before: 1 + I0 = −¯B 1 at = −1. The tendency term in (79) is again dropped because the interior evolution is slow, so the BL evolution must be slow as well. Expressing B1 =¯B 2 , vertically integrating (79), and enforcing I2 + B2 = 0 at = −1 yields an effective boundary condition on the interior. The BL-interior exchange velocity I2 = − B2 at = −1 may be obtained by vertically integrating The leading-order equations obtained using this more rigorous approach again match the expressions derived heuristically above. The asymptotic orders revealed by this approach are the same as in the 1D case.

Discussion
Callies and Ferrari (2018) studied the mixing-generated abyssal circulation in an idealized global basin using PG dynamics, but their model employed Rayleigh drag rather than a Fickian friction.
The models and theory presented here make use of a down-gradient turbulence closure of the momentum fluxes, allowing them to produce more realistic BLs and avoid unphysical interior momentum sinks. Still, the results presented here provide some insight into the conclusions from this previous study. With Rayleigh drag, Callies and Ferrari (2018) found that the canonical 1D model was a reasonably accurate emulator for the full dynamics over slopes with a constant initial stratification. This may have been somewhat of a coincidence, as in their case the steady state canonical transport ∞ cot was zero everywhere, adding a transport constraint to the canonical 1D model. With Fickian friction, setting ∞ = 0 does not immediately make the canonical 1D model equivalent to the transport-constrained 1D model because it still evolves diffusively and with nonzero transport, taking thousands of years to equilibrate (PC22). Rayleigh drag, in contrast, damps flow in the interior, allowing for fast adjustment (in a matter of years, not shown) to the = 0 steady state. The combination of ∞ = 0 and Rayleigh drag thus conspired to let Callies and Ferrari (2018) get the right answer from the canonical model, but modifying either of these choices would have made the argument fall apart.
Furthermore, Callies and Ferrari's (2018) application of BL theory was somewhat ad hoc. For slopes steep enough for the canonical BL theory to apply, the steady-state transport was exactly zero, meaning that all upslope transport was exactly balanced by downslope transport above. The BL theory broke down at the base of the slopes, allowing the BLs to be fed by dense water from the south and the less dense downwelled water to return south, forming a basin-wide circulation that constituted an overturning. The overturning transport could thus be estimated with an isobath integral of the upslope transport in BLs on the slopes. As Drake et al. (2020) pointed out, however, this approach is not successful if the interior stratification is far from constant and canonical BL theory does not apply. The theory presented here supplies a globally valid expression for the BL transport that allows for variations in the interior stratification. At this point, this expression is only a diagnostic tool, itself depending on the interior dynamics, but it unambiguously describes how the interior can exert control on the BL, and vice versa, ultimately generating a basin-wide circulation that involves both BL and interior pathways-and mass exchange between them. This sharpens our view of the abyssal overturning, with no confusion about the roles of the BL and interior.
The framework presented here can also help understand the results from Drake et al. (2020) regarding how water mass transformations are affected by changes in the interior stratification.
Using the same 3D PG model with Rayleigh friction as in Callies and Ferrari (2018), they found that the degree of compensation between BL upwelling and interior downwelling is strongly dependent on vertical variations in the initial stratification. With only the canonical 1D theory as a starting point, they were unable to explain the vertical extent and structure of water mass transformations. The BL theory presented here would enable us to understand these physics more clearly, because it explicitly separates the BL and interior components of the flow. This allows us to describe the abyssal circulation in terms of flows into and out of the BL, rather than simply bulk diapycnal motion throughout the water column. In section 3, we demonstrated the power of this framework in describing abyssal spin up in 2D with exponential initial stratification. Applied to 3D simulations such as those in Drake et al. (2020), this approach would undoubtedly shed light on what shapes the vertical structure of water mass transformations in the abyss.
Here, we have only presented results in 1D and 2D. We leave the 3D case to a future paper, but preliminary work indicates that much of the theory developed here carries over, although there are some key differences. In 3D, the interior dynamics satisfies geostrophic balance in both the and directions. Because of this, the asymptotics in 3D are qualitatively different from those presented in section 4 of this paper: instead of evolving on a slower timescale, the leading-order 3D interior buoyancy field is advected by the geostrophic velocities, with diffusion only playing a role at higher order. We anticipate that this qualitative difference between 2D and 3D may be crucial in explaining the full 3D abyssal circulation. In 3D, it is also no longer possible to write the PG inversion in terms of a scalar streamfunction. This makes the mathematics more complicated, but it is still possible to write down an expression for the 3D BL transport in terms of interior variables evaluated at the bottom. As in 2D, the 3D BL mass and buoyancy transports feed back on the interior, now with gradients in the direction shaping the flow field. A future extension to 3D will allow us to explain the dynamics of abyssal circulations in more complicated and realistic geometries, including cases with variations in .
Our BL theory results are not only theoretically useful but could also lighten the computational demand of simulating the abyssal circulation. The interior solution can be computed without the need to resolve the thin BL, allowing numerical models to have coarser grids and larger timesteps. This is crucial when studying the 3D system over long abyssal timescales of thousands or tens of thousands of years (e.g., Wunsch and Heimbach 2008;Liu et al. 2009;Jansen et al. 2018).
This framework could even be used to analyze tracer transport without explicitly resolving the BL, allowing us to better understand carbon and heat storage (e.g., Sarmiento and Toggweiler 1984) Variations in will allow for vortex stretching in the absence of friction: = , where = . In the -plane solutions considered here, a non-zero interior vertical velocity only appears at second order (see section 4). and Lagrangian pathways (e.g., Rousselet et al. 2021) in the abyss. If needed, the BL correction can be computed after the fact on a finer grid as was done for Figs. 3 and 6.
Although the results presented here are derived in the context of PG dynamics, they might also point the way towards a parameterization of the effects of BLs over a sloping seafloor in primitiveequation models. Applying effective boundary conditions on the interior evolution, following the BL framework, should most easily be accomplished in models with terrain-following coordinates But a translation to -coordinates also appears feasible, which would alleviate not only the need to resolve thin boundary layers in the vertical but also the need to capture BL flow across the artificial steps in the topography in such models. An extension of the BL theory to 3D is needed, however, to produce expressions directly useful for such a parameterization effort.
The circulation in the examples presented in this paper depend on the particular, simple closure of turbulent momentum and buoyancy fluxes employed in all of them. Although Fickian friction is much more physical than Rayleigh drag, our use of it with a simple profile for still glosses over the true complexity of turbulence in the abyss. Without a more faithful representation of the internal-wave field and baroclinic eddies in abyssal mixing layers, we cannot claim to be accurately simulating the dynamics of the real ocean. The BL framework, however, is robust to the choice of turbulence parameterization-as long as the vertical scale of the turbulent mixing in the interior is larger than the thickness of the BL, our approach should require minimal modification. The results presented here are in terms of a particular choice of parameterization, but the general themes describing how the BL and interior communicate will carry over to more complex closures. This flexibility makes BL theory an attractive tool for understanding the mixing-generated abyss over a hierarchy of complexities.

Conclusions
Motivated by observations of bottom-enhanced mixing, recent work on the abyssal circulation has focused on the role of thin bottom BLs (Ferrari et al. 2016;de Lavergne et al. 2016;McDougall and Ferrari 2017;Holmes et al. 2018;Callies and Ferrari 2018;Drake et al. 2020). Until now, the coupling between these BLs and the interior circulation remained opaque, with most of our understanding coming from somewhat heuristic arguments using 1D theory. The framework presented in this work uses BL theory to paint a clear picture of the interior-BL interaction of the mixing-generated abyssal circulation. By explicitly defining BL and interior contributions to the flow, we obtain expressions for the BL transport in 1D and 2D that are bounded for all bottom slopes, solving the old 1D conundrum of the steady total transport ∞ cot being set by the farfield mixing and diverging for small slopes. In the revised theory, the BL transport is set by local flow parameters and interior variables evaluated at the bottom, with the total transport allowed to evolve according to the global context. The interior dynamics are themselves modified by this BL transport, which advects dense water up-slope and thus modifies the interior bottom boundary condition. This two-way coupling provides a complex yet transparent story of how BLs influence the abyssal circulation, and this framework makes previously unwieldy problems, such as determining the response to vertically varying initial stratification, comparatively simple. With these promising results, we anticipate that BL theory will play a crucial role in the development of a more complete understanding of the abyssal circulation in the real ocean.
Data availability statement. The numerical models for all the simulations presented here are hosted at https://github.com/hgpeterson/nuPGCM.

BL Theory when
≠ 0 For completeness, we here show how the BL theory derivations in sections 2 and 3 are slightly altered when the transport is non-zero. In both 1D and 2D, the interior inversion is modified to include the added transport term. The BL accounts for 1/(1 + ) of the total transport, leading to a modified interior bottom boundary condition compared to before. The 2D case is special in that the total transport is itself a function of the flow and geometry of the domain (see appendix B of PC22), allowing us to derive an explicit equation for in that case.

a. One-dimensional theory
The 1D interior inversion for general is 2 ( I − ) = − I tan .
This does not affect I , leaving the interior evolution equation (11) unchanged. This new interior balance results in a modified bottom boundary condition compared with equation (15): 2 + (1 + ) I = 2 tan at = 0.
The added flux on the right-hand side represents the integrated buoyancy supplied to the column by the net transport . The BL transport (cf. 16) now takes the form supplying a fraction of the total transport. For 1, the BL absorbs the majority of the added transport. Note that the bottom boundary condition may be written as ( 2 + I ) = I 2 tan regardless of whether is nonzero. The BL correction B remains the same as in (18), with I at = 0 now coming from (A3).
The asymptotic order of must match that of I , so it must be restricted to be ( 2 ). It is then simple to incorporate ≠ 0 into the theory presented in section 4.

b. Two-dimensional theory
In 2D, the general interior inversion is 2 and again the BL absorbs a fraction of the added transport so that (37) becomes where = − I / 2 at = −1. For symmetric topography, = 0, but this is not the case in general. We can infer for asymmetric geometries with knowledge of the interior buoyancy distribution. Evaluating (39) at = 0 and taking the mean in , denoted by · , we have where, crucially, the BL transport from equation (A5) now depends on . We have assumed that the domain is tall enough such that gradients in buoyancy at = 0 are small and therefore I (0) = 0. Solving for yields where all variables are evaluated at = −1 unless otherwise noted. Simulations of an asymmetric ridge, similar to that in appendix B of PC22, confirm the accuracy of this formula (not shown).
Again, we restrict ourselves to cases where the non-dimensional is ( 2 ), the same order as I . This is true when the second term on the right in equation (A6) of lower order than the first. This is always the case after a fast initial adjustment.

Axisymmetric Coordinates
For simulations of an idealized seamount, we transform to axisymmetric coordinates, assuming rotational symmetry. The depth is then a function of the radial distance and invariant under rotation about the origin by some angle , leading to effectively 2D flow. Defining = and = / , we have