Dynamics of Eddying Abyssal Mixing Layers over Sloping Rough Topography

: The abyssal overturning circulation is thought to be primarily driven by small-scale turbulent mixing. Diagnosed water-mass transformations are dominated by rough topography “ hotspots, ” where the bottom enhancement of mixing causes the diffusive buoyancy ﬂ ux to diverge, driving widespread downwelling in the interior } only to be overwhelmed by an even stronger upwelling in a thin bottom boundary layer (BBL). These water-mass transformations are signi ﬁ cantly underestimated by one-dimensional (1D) sloping boundary layer solutions, suggesting the importance of three-dimensional physics. Here, we use a hierarchy of models to generalize this 1D boundary layer approach to three-dimensional eddying ﬂ ows over realistically rough topography. When applied to the Mid-Atlantic Ridge in the Brazil Basin, the idealized simulation results are roughly consistent with available observations. Integral buoyancy budgets isolate the physical processes that contribute to realistically strong BBL upwelling. The downward diffusion of buoyancy is primarily balanced by upwelling along the sloping canyon sidewalls and the surrounding abyssal hills. These ﬂ ows are strengthened by the restratifying effects of submesoscale baroclinic eddies and by the blocking of along-ridge thermal wind within the canyon. Major topographic sills block along-thalweg ﬂ ows from restratifying the canyon trough, resulting in the continual erosion of the trough ’ s strati ﬁ cation. We propose simple modi ﬁ cations to the 1D boundary layer model that approximate each of these three-dimensional effects. These results provide local dynamical insights into mixing-driven abyssal overturning, but a complete theory will also require the nonlocal coupling to the basin-scale circulation.


Introduction
Below the oceanic pycnocline, the vast volumes of the deep ocean are ventilated by two interconnected cells of a global meridional overturning circulation (Gordon 1986).The lower cell of this circulation is sourced along the coast of Antarctica, where atmospheric cooling and brine rejection transform surface waters into the dense Antarctic Bottom Waters (AABW) that fill the global abyssal ocean at a rate of approximately 30 Sv (1 Sv ≡ 10 6 m 3 s 21 ) (Talley 2013).Since the buoyancy surface bounding AABW from above does not outcrop elsewhere in the ocean, conservation of mass implies that in steady state an equal amount of AABW must upwell across buoyancy surfaces (diabatically) from the abyss.Waters below about 2000-m depth (corresponding to the crests of major topographic features, such as midocean ridges) can upwell diabatically only in the presence of interior watermass transformations (e.g., small-scale turbulent mixing) or fluxes across the seafloor boundary (geothermal heating) (Munk 1966;Walin 1982;Tziperman 1986;Ferrari 2014).
These basic inferences of a global diabatic upwelling from the abyss (e.g., Sverdrup et al. 1942) are also consistent with more detailed inverse modeling at regional scales (e.g., Talley et al. 2003).Most notably, Hogg et al. (1982) consider the fate of 4 Sv of AABW (colder than 08C) that enters the Brazil Basin from the Southern Ocean through the Vema Channel; since there are no other exits from the basin and since geothermal fluxes are relatively weak, they infer that turbulent mixing must diffuse heat downward at a rate of O(3) cm 2 s 21 to balance the upwelling of these waters across the 08C isotherm.
Early in situ turbulence measurements in the upper 1000 m of the interior ocean suggested turbulent diffusivities an order of magnitude smaller than those predicted by the large-scale abyssal tracer budgets described above (Gregg 1987;Ledwell et al. 1993).A subsequent celebrated field campaign in the abyssal waters of the Brazil Basin reported similarly weak background diffusivities over the smooth topography of the abyssal plains, but revealed diffusivities that increased with depth by several orders of magnitude over the rough topography of the Mid-Atlantic Ridge (Polzin et al. 1997;Ledwell et al. 2000).Using regional inverse and forward approaches, respectively, St. Laurent et al. (2001) and Huang and Jin (2002) modeled the impacts of the observed bottom-enhanced mixing on the regional circulation: bottom-enhanced mixing drove interior downwelling while upwelling was restricted to a thin layer of buoyancy convergence near the bottom boundary [in contrast to Munk's (1966) uniform upwelling model] and the basin-scale horizontal circulation was dominated by narrow mixing-driven flows along ridge flanks [in contrast to the broad interior geostrophic flow predicted by Stommel (1958)].
The development of mixing parameterizations (e.g., St. Laurent and Garrett 2002;Kunze et al. 2006;Polzin 2009;Melet et al. 2014;de Lavergne et al. 2020) allowed these Brazil Basin results to be generalized to global abyssal water-mass transformations (e.g., Nikurashin and Ferrari 2013;de Lavergne et al. 2016;Kunze 2017;Cimoli et al. 2019).Based on such estimates, Ferrari et al. (2016) and McDougall and Ferrari (2017) revised the conceptual model of the global mixing-driven abyssal upwelling: mixing-driven diabatic upwelling is confined to a thin bottom boundary layer (BBL) just above the insulated (or geothermally heated) seafloor, while bottom-enhanced mixing drives diabatic downwelling in the stratified mixing layer (SML) above; the net diabatic overturning is the small remainder of these two large opposing mixing layer flows.In this emerging framework, the global overturning circulation is modulated by the dynamics of thin BBLs (Callies and Ferrari 2018;Drake et al. 2020).However, these abyssal boundary layer flows are challenging to observe (Naveira Garabato et al. 2019;Spingys et al. 2021) and are too thin to be resolved by conventional general circulation models, so they remain poorly understood [Drake (2021) and Polzin and McDougall (2022) discuss various outstanding questions].
The interpretation of the role of "boundary mixing" in the abyssal overturning circulation [dating back to Munk (1966)] has a contentious history: on the one hand, in situ observations of weakly stratified bottom mixed layers seemed to imply the existence of vigorous boundary mixing (Armi 1978); on the other hand, it was argued that mixing of already wellmixed waters was inefficient and thus did not lead to significant water-mass transformation [see Garrett's (1979) comment and Armi's (1979b) reply].Garrett (1990) later formalized his criticism using sloping boundary layer theory (Phillips 1970;Wunsch 1970) and suggested that one-dimensional (1D) flows up the sloping bottom boundary}driven by the mixing itself}could provide sufficient restratification to resolve this conundrum.Based on observations of homogeneous layers detached from the bottom boundary (but carrying distinct levels of suspended sediments), Armi (1978Armi ( , 1979a) ) instead proposed a three-dimensional boundary-interior exchange process whereby layers are rapidly mixed due to boundary layer shear turbulence when they impinge upon topographic features (e.g., seamounts or abyssal hills) and are eventually restratified by along-isopycnal exchanges with the stratified interior.
In light of recent diagnostic evidence for boundary control on the abyssal circulation (Ferrari et al. 2016), Callies (2018) revisited these ideas to test whether sloping BBL theory is quantitatively consistent with observations.In his analysis of the sloping flank of the Mid-Atlantic Ridge in the Brazil Basin (where collocated measurements of both abyssal mixing rates and stratification are available), he found that the steady state 1D boundary layer solution forced by the observed mixing exhibits a stratification an order of magnitude weaker than observed.The water-mass transformations sustained by 1D upslope advection alone (Garrett 1990) are thus too inefficient to contribute significantly to the global abyssal overturning circulation.
To reconcile boundary layer dynamics with observations, Callies (2018) argued the stratification of abyssal mixing layers may be maintained by submesoscale baroclinic eddies, which act to slump sloping buoyancy surfaces back to the horizontal.Mixing-driven 1D boundary layer solutions are linearly unstable to submesoscale baroclinic modes (Wenegrat et al. 2018;Callies 2018), in a manner similar to the well-studied analogous problem in the surface mixed layer (Boccaletti et al. 2007;Fox-Kemper et al. 2008).Callies (2018) simulated the finite amplitude evolution of these instabilities in a three-dimensional (3D) generalization of the 1D boundary layer framework and showed the solutions converge on a substantially stronger quasi-equilibrium stratification that is more consistent with observations.
As acknowledged by Callies (2018), however, it is not clear to what extent such idealized solutions are directly applicable to the midocean ridge, which is characterized by particularly rough topography.For example, many of the observations of bottom-enhanced mixing, stratification, and diabatic upwelling from the region are confined to O (500)-m-deep fracture zone canyons that cut across the ridge (Polzin et al. 1997;Ledwell et al. 2000;St. Laurent et al. 2001;Thurnherr and Speer 2003).To account for these leading-order topographic features, Ruan and Callies (2020) ran simulations of mixingdriven flow over a sinusoidal midocean ridge incised by an idealized Gaussian fracture zone canyon.They confirm Thurnherr and Speer's (2003) speculation that the canyon sidewalls suppress cross-canyon (or along-slope) flow and thus support a vigorous up-canyon (or cross-slope) mean flow.The restratifying tendency of this up-canyon mean flow is much stronger than that of either the 1D upslope flow or the submesoscale eddies on the smooth ridge flank, implying that abyssal watermass transformations are, per unit area, 4 times larger within the canyons than on the ridge flank.Ruan and Callies (2020) found, however, that the simulated stratification in the canyon is orders of magnitude larger than observed, suggesting their simulations are still missing important physics.In addition to fracture zone canyons, midocean ridges are also characterized by smaller-scale anisotropic abyssal hills; these features have characteristic scales taller than 1D BBLs and comparable to those of the fastest growing baroclinic mode (Callies 2018;Wenegrat et al. 2018), so we would expect them to affect both mean and eddying circulations.Within the fracture zone canyons, abyssal hills often manifest as sills that substantially block or constrain the deep upslope flow (Thurnherr et al. 2005;Dell 2013;Dell and Pratt 2015); hydraulic acceleration over the sill produces relatively large velocities also associated with locally enhanced turbulence (Clément et al. 2017).
Here, we use a hierarchy of analytical and numerical solutions to bridge the gap between idealized 1D BBLs and the complexity of observed flows in a region scarred by a fracture zone canyon and dotted with abyssal hills.In section 2, we review key insights from the 1D BBL buoyancy budget and derive a generalized buoyancy budget that permits topographic variations and spatiotemporal eddy correlations.In section 3, we describe the "slope-aligned" simulation configuration which leverages a coordinate frame aligned with the mean topographic slope to allow restratification by mean upslope flow across a uniform background vertical buoyancy gradient.In section 4, we describe the simulated mixing layer flows in a simulation with realistic topography and show they are qualitatively consistent with available observations.In section 5, we present simulated buoyancy budgets, and show a balance between bottom-enhanced mixing, submesoscale eddy fluxes, and the cross-slope mean flow.By using a hierarchy of models (Held 2005), we progressively simplify the configuration to isolate the roles of individual physical processes in setting the near-boundary stratification.In section 6, we discuss how our results bridge the gap between interpretations of in situ observations (e.g., Armi 1978;Thurnherr and Speer 2003;Thurnherr et al. 2020) and 1D BBL theory (e.g., Garrett 1979;Garrett et al. 1993) and how they illustrate}at a regional scale}the control of abyssal mixing layers on an "upside-down" abyssal overturning circulation (Ferrari et al. 2016).We conclude that a combination of mixing-driven upslope flows, submesoscale baroclinic eddies, and topographic control are required to maintain a steady state near-boundary stratification consistent with in situ observations and a finite global abyssal overturning circulation.

Theory
We review the derivation and results of sloping boundary layer theory in sections 2a and 2b in anticipation of our generalization to 3D flows over rough sloping topography in section 2c.

a. Slope-aligned equations
In sloping boundary layer theory (Wunsch 1970;Phillips 1970;Garrett et al. 1993;Thompson and Johnson 1996;Callies 2018;Holmes and McDougall 2020), analytical progress is achieved by modeling the system in a coordinate frame aligned with the mean topographic slope, rather than the typical coordinate frame (x, ŷ, ẑ) with ẑ aligned with gravity.
It is useful to decompose the buoyancy B N 2 ẑ 1 b into a background component N 2 ẑ, where N 2 is a constant vertical buoyancy gradient, and a perturbation component b(x, ŷ, ẑ, t); the background buoyancy is assumed to be in hydrostatic balance with a background pressure and we similarly decompose P (1/2)N 2 ẑ2 1 p.Then, we rotate the coordinate system to a coordinate frame aligned with the mean-slope (x, y, z) ≡ (x cos u 1 ẑ sin u, ŷ, ẑ cos u 2 x sin u), where u is the region's average slope angle in the x direction.For small slopes1 (tanu , , 1), the hydrostatic Boussinesq equations in the mean-slope coordinates are, at leading order, given by where subscripts represent partial derivatives, ∇ is the gradient operator, u is the along-canyon (or cross-slope) velocity, y is the cross-canyon (or along-slope) velocity, w is the slope-normal velocity, f is a constant Coriolis parameter, k is an isotropic eddy diffusivity, and n = sk is an isotropic eddy viscosity determined by the turbulent Prandtl number s.
The rotated along-canyon x-momentum equation is identical in form to the zonal x-momentum equation with the exception of the small but dynamically significant projection of the perturbation buoyancy force bẑ on x.
We define the anomalous seafloor depth, relative to the mean slope, as For convenient intermodel comparison later on, we set z = 0 along the sloping plane that intersects the point with the greatest anomalous seafloor depth, max(d).Boundary conditions at the seafloor, z = max(d) 2 d(x, y), are u = 0 (no-slip 2 and no-normal-flow) and n • (k∇B) = 0 (insulating 3 ), where n is a unit vector normal to the boundary.

b. Smooth planar slopes and steady 1D dynamics
Assuming a constant topographic slope (d ≡ 0) and mixing rates that vary only in the slope-normal direction, the equilibrium problem reduces to 2f y cos u b sin u 1 z (nu z ), (7) p z b cos u, (9) where the continuity equation w z = 0 combines with the nonormal-flow bottom boundary condition at z = 0 to require w ≡ 0 everywhere (no slope-normal exchange).These equations can be solved analytically in the case of constant parameter values (Wunsch 1970;Phillips 1970;Thorpe 1987;Garrett 1990), or approximately for varying parameters in some asymptotic limits (Salmon et al. 1991;Callies 2018).In either case, the slope Burger number S ≡ N 2 tan 2 u/f 2 and the BBL thickness emerge as key parameters.We recognize d as the Ekman layer thickness d E ≡ 2n/f √ , modified by buoyancy effects at the sloping boundary; for typical abyssal values, S , , 1 and s O(1) such that buoyancy effects are weak (Thurnherr and Speer 2003).
Recalling the crucial assumption of a constant background vertical stratification N 2 , the slope-aligned buoyancy Eq. ( 10) describes a direct balance between slope-normal diffusion of heat downward toward the boundary and cross-slope advection against the constant background buoyancy gradient; this balance is a near-boundary analog of Munk's (1966) classic interior ocean vertical balance.This is best illustrated by integrating (10) in the slope-normal direction, where c is the upslope transport (per along-slope unit length) and we have invoked the insulating bottom boundary condition on the full stratification, B z = 0 at z = 0. Consider the case of exponentially bottom-enhanced mixing, k(z) = k BG 1 k BOT e 2z/h with k BOT /k BG . . 1. Equation ( 12) reveals two key insights: 1) The net upslope transport, integrated over both the upwelling BBL and the downwelling SML, converges to the negligibly small value,4 since far from the boundary b z → 0 (Thorpe 1987) and k(z) → k BG is small.
2) Maximal upslope transport in the BBL is achieved when both k is large (i.e., near the boundary) and B z is large (strong restratification).If the stratification is maintained near the background value N 2 cosu where the diffusivity is large (i.e., z , , h) then the upslope transport in the BBL reaches an upper bound of max{c} k BOT cot u (k BOT /k BG )c ' . .c ' : Callies (2018) derives approximate but analytical boundary layer solutions to the steady 1D system [Eqs.( 7)-( 10)] for bottom-enhanced mixing.In the abyssal ocean regime with typical values of Ss , , 1, the equilibrium stratification B z is approximately inversely proportional to k in the SML [his Eq. (10); Fig. 1, solid line] such that the diffusive buoyancy flux kB z k BG N 2 cosu is constant}finite buoyancy flux convergence occurs only within the thin BBL.Since the BBL stratification is reduced to roughly B z ≈ (k BG /k BOT )N 2 cosu (Fig. 1, solid line) and near-boundary mixing is thus inefficient, the transport equation ( 12) implies that upslope BBL transport is roughly equal to the negligibly small integral constraint (13), max{c} k BG cotu.This weak BBL upwelling and negligible SML downwelling contrasts with the strong bidirectional flows inferred from water-mass transformation analyses (Ferrari et al. 2016;McDougall and Ferrari 2017).

c. Rough topography and eddy fluxes
We now derive the 3D integral buoyancy budget, which allows for topographic and flow variations along the plane of the slope.Consider the buoyancy budget for a volume V within a height z above the mean slope (Fig. 2): FIG. 1. Height-above-bottom stratification profiles at steady state for 1D BBL models with the same external parameters as the BBTRE simulations (solid), without rotation ( f = 0; dotted), and with enhanced vertical diffusion of along-slope momentum, s y (z) . . 1 (dash-dotted; see appendix).The gray line marks the background stratification N 2 .
where we invoked the divergence theorem to rewrite the righthand-side terms in terms of fluxes normal to the bounding surface V (Fig. 2).Fluxes through the seafloor at z = max(d) 2 d(x, y) vanish due to the no-flow and insulating bottom boundary conditions.Motivated by the simulations in section 3, we assume fluxes through cross-slope and along-slope boundaries cancel due to periodicity in the perturbations [i.e., b(x) = b(x 1 L (x) )], except for the upslope flow across the background buoyancy gradient (recall u dy dz : (15) This, combined with the slope-normal component of the flux through the z = z surface, gives where we define f ≡ A(z) f dx dy as the slope-integral operation and C(z) ≡ A(x;z # z) u dy dz as the upslope transport across the periodic boundary (Fig. 2a).At equilibrium, the form of the generalized volume-integral buoyancy equation ( 16) is similar to the 1D transport equation ( 12), although there is now an additional eddy flux of buoyancy toward or away from the boundary, and the turbulent buoyancy flux may be modified by along-and cross-slope correlations between k and B z .Assuming a steady state and integrating up far into the interior, where k → k BG and the perturbations vanish, we recover the integral constraint (13) on the net upslope transport (per unit length) from the 1D solution, C ' /L (y) = k BG cotu.Callies (2018) proposes a simple representation of restratification by 3D submesoscale baroclinic eddies as a way to account for these missing physics in the 1D boundary layer solution.The main effect of baroclinic eddies is to extract available potential energy from the mean flow by slumping sloping buoyancy surfaces back toward the horizontal, thereby maintaining a realistically large near-bottom stratification; this adiabatic process is most conventionally parameterized as an eddy overturning circulation (Gent and McWilliams 1990;Fox-Kemper et al. 2008).Taking advantage of thermal wind balance (f y ẑ b x ), the slumping of isopycnals by baroclinic instability}which decreases horizontal buoyancy gradients b x }can equivalently be parameterized as a reduction in the vertical shear y ẑ , e.g., by enhanced vertical momentum diffusion (Rhines and Young 1982;Greatbatch and Lamb 1990;Young 2012).We provide a derivation of this closure in the appendix, in which we apply Andrews and McIntyre's (1976) transformed Eulerian mean and Gent and McWilliams's (1990) baroclinic eddy parameterization scheme to the slopealigned framework.
Following Callies (2018), we thus represent submesoscale eddy restratification by artificially increasing the vertical eddy viscosity n = sk.Unlike Callies ( 2018), who chose s = 230 to match the mean behavior of his 3D model, however, we: 1) only enhance the viscosity n y = s y k acting on the alongslope thermal wind (as in Holmes et al. 2019) since the available potential energy that fuels the instabilities is stored in cross-slope buoyancy gradients; 2) we allow the eddy viscosity to have slope-normal structure, s y = s y (z); and 3) we estimate the magnitude and structure of s y (z) from the eddy fluxes resolved by a 3D simulation (Fig. A1b; see appendix).We refer to C 1 ( wb /N 2 sinuL (x) ) as the cross-slope residual transport (analogous to that of the Southern Ocean, e.g., Marshall and Radko 2003), since the additional eddy flux term is approximately equal to the eddy-induced overturning streamfunction defined in the appendix for the smooth planar-slope case.
Applying this simple closure to the 1D model results in a reduction of the slope-normal shear of the along-slope flow and, because of the approximate thermal wind balance fy z 2b z sinu that holds in the SML [(Eq.( 7)], in a corresponding reduction of the negative perturbation stratification b z (equivalent to a strengthening of the total stratification B z ; compare dash-dotted and solid lines in Fig. 1).In this context, the 1D model's upslope transport c is reinterpreted as the residual transport, since it implicitly includes the eddy-induced overturning.At equilibrium, this representation of eddy restratification triples B z and thus also kB z and the residual flow c at the top of the 1D solution's BBL (Fig. 1), bringing the water-mass transformations of the 1D BBL more in line with the basin-scale overturning (Morris et al. 2001;Callies 2018).

Numerical models
We simulate 3D mixing-driven flows using the hydrostatic Boussinesq equations in the MIT General Circulation Model (MITgcm; Marshall et al. 1997).For simplicity, we assume a linear equation of state; because temperature units are more intuitive, we use temperature T and buoyancy b ≡ gaT interchangeably throughout, where g = 9.81 m s 22 is the gravitational acceleration and a = 2 3 10 24 8C 21 is a constant thermal expansion coefficient.

a. Realistic bathymetry
Most of the results describe a core realistic-bathymetry simulation of the Brazil Basin subregion sampled by both the Brazil Basin Tracer Release Experiment (BBTRE; Ledwell et al. 2000) and the Dynamics of the Mid-Ocean Ridge Experiment (DoMORE; Clément et al. 2017), located on the western flank of the Mid-Atlantic Ridge.We extract the Brazil Basin's seafloor topography from the Global Bathymetry and Topography at 15 dataset (SRTM151; Tozer et al. 2019), which includes many more multibeam measurements than previous products (e.g., Smith and Sandwell 1997) and thus better resolves both the BBTRE fracture zone canyon at 21830 S and the smallerscale abyssal hills characteristic of midocean ridges (Fig. 3a).High-resolution multibeam measurements are not available for shallower parts of the hilly ridge flank surrounding the canyon (Fig. 3a) so we underestimate the impacts of small-scale roughness there.We interpolate the bathymetry onto a locally tangent Cartesian grid (x, ŷ, ẑ) aligned with the BBTRE canyon, where x denotes the along-canyon dimension and ŷ denotes the cross-canyon dimension (Fig. 3a), and produce a gridded bathymetry field d(x, ŷ).The simulation domain stretches from a few km west of the Tracer Release Experiment site around 18.58W (Ledwell et al. 2000) to a few kilometers east of the DoMORE sill that dramatically constrains the up-canyon flow at 14.58W (Clément et al. 2017).

b. Implementing the perturbation Boussinesq equations in the mean-slope coordinate frame
Following section 2a, we solve Eqs. ( 1)-( 5) in a coordinate frame aligned with the domain's mean slope.Equations ( 1)-( 5) are solved in terms of the perturbation variables, with the background buoyancy field N 2 ẑ entering only indirectly via linear and inhomogeneous terms in the perturbation buoyancy equation, implemented as additional explicit tendency terms in the MITgcm.To stabilize the numerical solution without damping submesoscale eddies, we additionally implement horizontal (in the rotated frame) biharmonic hyperdiffusion of momentum and buoyancy which acts only at scales close to the grid resolution.Horizontal hyperdiffusive tendencies vanish in the budgets presented here, so we omit them in all of our analyses.We enforce an insulating boundary condition on the full buoyancy at the seafloor: n • (k∇B) = 0.
Relative to the mean slope, the anomalous seafloor topography d(x, y) ≡ d(x, ŷ) 2 x tan u is nearly continuous across the periodic boundaries in the along-canyon direction x and in the cross-canyon direction y; however, to eliminate any remaining discontinuities across these boundaries, we join the two boundaries smoothly by linear interpolation in both x and y over a few grid cells.
By 1) removing the uniformly stratified background state from the prognostic variables, 2) formulating the model in the slope coordinate frame, and 3) making the boundary conditions and forcing terms periodic in the (x, y) plane, we are free to apply periodic boundary conditions to the perturbation state variables u, y, b, and p in both x and y.

c. Forcing by observation-inspired bottom-enhanced turbulent mixing
Following the classic 1D boundary layer configuration (Wunsch 1970), we represent small-scale turbulent mixing as a slope-normal5 diffusive buoyancy flux 2k z Bz.We use Callies' (2018) self-similar height-above-bottom (hab) profile with k BOT = 1.8 3 10 23 m 2 s 21 , k BG = 5.3 3 10 25 m 2 s 21 , and h = 230 m; these parameter values are chosen by performing a least squares fit to the hab average of 126 microstructure profiles in the BBTRE region.The sparsity and noisiness of individual mixing profiles, and disagreements in the literature about where mixing in strongest (Polzin et al. 1997;St. Laurent et al. 2001;Polzin 2009;Clément et al. 2017;Thurnherr et al. 2020), prohibit the formulation of a robust parameterization with a richer spatial structure.We imagine this imposed bottomenhanced mixing to represent a variety of turbulent ocean processes (see Thorpe 2005), especially the breaking of internal waves (Whalen et al. 2020) but also including unspecified boundary mixing processes (Armi 1978;Armi and D'Asaro 1980;Polzin et al. 2021).

d. Numerics
The horizontal grid spacing of Dx = D y = 600 m is fine enough to permit the anticipated submesoscale baroclinic turbulence, which for the 1D sloping BBL problem has a maximum linear growth rate near the local deformation radius L ∼ NH ML /f 6 km (Stone 1966;Wenegrat et al. 2018), where H ML ≈ 250 m is the thickness of the weakly stratified bottom layer (Callies 2018).Yet, the grid spacing is also coarse enough for a 3D simulation of the entire 480 km 3 60 km region to be computationally feasible.We set the hyperdiffusivities k 4 ≡ n 4 = 2 3 10 4 m 4 s 21 , the smallest value that maintains a stable solution, so that hyperdiffusion interferes minimally with diapycnal buoyancy fluxes and the growth of submesoscale instabilities (Callies 2018;Ruan and Callies 2020).In the vertical, a cell thickness of Dz = 6 m (with partial cells down to 1.2 m) marginally resolves the predicted O(10)-m-thick BBL.A 1D spinup experiment confirmed this vertical resolution is sufficient to accurately reproduce all features of the analytical solution.Above z = 1000 m, the cell thickness Dz is increasingly stretched (up to Dz = 50 m at the top of the domain) to efficiently fit both the h log(k BOT /k BG ) ≈ 1300 m vertical scale of abyssal mixing layers (Callies 2018) and the O(800) m topographic extent into a domain that spans a height H = 2700 m above the deepest point (relative to the mean slope).Averaging the bathymetry d(x, ŷ) in ŷ and applying a linear fit in x yields the domain's average topographic slope angle u = 1.26 3 10 23 .We assume that small-scale turbulent mixing acts similarly to mix buoyancy and momentum, i.e., we assume a turbulent Prandtl number of s ≡ n/k 1. (Because we resolve submesoscale instabilities, we do not need to approximate their restratification by increasing s.) Mixing layers are thus characterized by weak stratification and gentle large-scale slopes, equating to a small slope Burger number, S ≡ N 2 tan 2 u/f 2 = 10 23 , , 1 and BBL thickness d ≈ 8 m [Eq.( 11)].
We spin up the simulations from a uniformly stratified rest state (b = 0, p = 0, u = 0).The BBL adjusts rapidly on a time scale t BBL = d 2 /k BOT = 10 h.While the full equilibration of the solution occurs over a prohibitively long diffusive time scale characteristic of the abyssal ocean interior, t INT = H 2 /k BG ≈ 4000 years, buoyancy tendencies are small enough by t = 13 years in the bottom 1000 m (see section 5) that we consider the solution sufficiently equilibrated for the analyses presented here.

f. Hierarchy of progressively idealized simulations
The simulations in our model hierarchy differ only in their seafloor topography, domain length, and dimensionality.We progressively idealize the BBTRE canyon configuration (Fig. 3f): first, we remove the abyssal hills along the ridge flank and idealize the geometry of the remaining canyon and sill features (Canyon1Sill; Fig. 3e); second, we remove the sill (Canyon; Fig. 3d); third, we remove the canyon entirely (Smooth3D; Fig. 3c); and finally, we eliminate flow variations along the plane of the slope, collapsing the solution onto a single slope-normal dimension as in classical BBL theory (1D).For reference, we also include some additional variants on the 1D model where we vary one parameter at a time: nonrotating (1D f=0 ), nonsloping (1D u=0 ), and with a crude closure for restratification by submesoscale eddies (1D sy(z) ; see appendix).Unless we specify otherwise, results refer to the realistic-topography BBTRE simulation.

Mixing-driven up-canyon flow, submesoscale turbulence, and stratification
At quasi-equilibrium, the time-mean flow (averaged over days 5000-5500) is dominated by a vigorous up-canyon jet along the canyon thalweg, banked along the steeper southern sidewall of the canyon (as in Dell 2013; Ruan and Callies 2020).The up-canyon jet exhibits a maximum along-canyonaveraged velocity of u x 0:75 cm s 21 about 400 m above the seafloor (Fig. 4a).This upslope jet is nonuniform and partially compensated by a downslope jet on the gentler northern sidewall, such that the maximum cross-canyon-averaged upcanyon velocity is reduced to u y O(0:1) cm s 21 (Figs.4a,b).The upslope jet accelerates as it spills over two major crosscanyon sills: the BBTRE sill at x = 110 km and the DoMORE sill at x = 420 km (Figs.4a,b); this acceleration and the spilling over of isopycnals at both sills is suggestive of hydraulic control6 (Pratt and Whitehead 2008).The vertically integrated cross-slope transport H z 0 u dz is dominated by standing eddy features above the canyon (Fig. 4c, recall z = 0 at the deepest point relative to the mean slope), but prominently features meandering up-and down-canyon jets when integration is restricted to just the canyon itself, 800 m z 0 u dz (Fig. 4d).These simulated mixing-driven mean flows can be compared against two in situ mooring observations: the BBTRE mooring at x = 110 km, several kilometers upstream of the BBTRE sill (Toole 2007; also analyzed by Thurnherr et al. 2005), and a DoMORE mooring at x = 420 km, atop the DoMORE sill (Clément et al. 2017).At the DoMORE sill, horizontal and vertical constrictions accelerate the simulated up-canyon flow to 5 cm s 21 over a layer dz = 150 m thick and dx = 2.5 km wide (Fig. 5a).The resulting velocities are roughly constant in time, also suggestive of hydraulic control (Pratt and Whitehead 2008), and are about 25% those measured by the mooring (half as fast and half as thick; Fig. 5b).By contrast, the simulated up-canyon flow is much weaker at the BBTRE mooring (u ≈ 0.75 cm s 21 ) but spread over a thicker (dz ≈ 600 m) and wider (dx ≈ 5 km) layer, such that the total up-canyon transports at the two sections are similar (Fig. 5c).It is impossible to compare against observed transports because single mooring velocity profiles (e.g., Thurnherr et al. 2005) cannot be reliably extrapolated across the canyon, although such errors may be smaller at constrictions considerably narrower than the local deformation radius (Thurnherr 2000), as at the DoMORE sill.The simulated flow at the BBTRE mooring has roughly the same vertical structure as in the moored current meter velocities, but about half their magnitude (Fig. 5d).The relative weakness of the simulated flows suggests that either the imposed microstructure-based mixing rates are bi- dotted red lines).Most of this discrepancy is resolved by sampling the simulation at the exact locations of the observational profiles (Fig. 6b), and comparing their sample mean to that of the observations (Fig. 6a, red lines).Since the BBTRE sampling strategy was to find as much tracer as possible, the field campaign specifically focused on sampling the deep depressions in the BBTRE canyon, which appear to exhibit unusually weak stratification compared to the canyon sidewalls and sills, and the surrounding ridge flank.However, several microstructure profiles from the 1996 cruise are available along the canyon crests}just north of the domain}and on average exhibit strong near-bottom stratification (Fig. 6a, dashed blue line) similar to the simulation's domain average.This conditional averaging exercise clarifies the significant disagreements in reported estimates of the BBTRE region's average stratification (Polzin et al. 1997;St. Laurent et al. 2001;Polzin 2009).But even accounting for sampling bias, the simulated canyon is more stratified by about a factor of 2 relative to the observations (see section 6).
The time-mean view of the up-canyon circulation described above filters out a rich field of submesoscale eddies which have radii comparable to the deformation radius and are trapped within a few hundred meters of the seafloor, including within the O(10)-km-wide canyon (Fig. 7; Supplemental Movie M1).These eddies manifest themselves as spatial and temporal meanders of the mean up-canyon jet (Supplemental Movie M2), which in the following section are shown to contribute significantly to the simulation's buoyancy budget and to maintaining its strong near-bottom stratification.

Buoyancy budgets: Mixing, eddies, and mean flow
In this section, we use a hierarchy of models to elucidate the complicated dynamics that support the vigorous up-canyon mean flows described in the previous section.Volume-integrated buoyancy budgets [Eq.( 16)] provide the major insights and are presented in Fig. 8 for each model in the hierarchy.We further separate the eddy term wb into contributions from standing eddies (spatial correlations of the time-mean w and b) and transient eddies (the remainder).All of the solutions exhibit substantial residual tendencies several hundred meters above the topography; however, within a few hundred meters of the seafloor, tendencies are an order of magnitude smaller than other terms in the budgets because the dynamics (vigorous mixing and submesoscale processes) within the bottom few hundred meters are much faster than the weak diffusion in the interior (Fig. 8, black).The 1D simulations are computationally inexpensive, so we also provide their fully equilibrated solutions for context (Figs.8a,b, dotted).
In the classical 1D solution, a weak upslope transport in the BBL (Fig. 8a, blue line) maintains a weak near-boundary stratification, although it is already much stronger than in the nonsloping case after 5000 days of spinup (Fig. 9a).The evolution of the Smooth3D solution follows the 1D solution closely until about 800 days (not shown), at which point the laminar solution becomes unstable to submesoscale baroclinic modes which rapidly grow and equilibrate at finite amplitude (Callies 2018;Wenegrat et al. 2018).At quasi-equilibrium, these transient eddies advect denser waters from the SML back into the BBL (Fig. 8b, orange), effectively restratifying the BBL (Fig. 9b) and thus strengthening the maximum diffusive buoyancy flux (Fig. 8b, red).It is helpful to interpret the  16)], averaged over days 5000 to 5200, for a layer bounded by a given HAMS.We interpret the sum of the mean flow and eddy terms as a residual flow.The left-hand-side tendencies (LHS) are equal to the remainder of the approximate balance (RHS) between slope-normal turbulent diffusion and the cross-slope residual circulation, which includes both mean and eddy components.We divide [Eq.( 16)] by the factor 2N 2 L (x) sinu to conveniently express the budget in terms of the quantity of interest, the upslope volume transport C with units of mSv ≡ 10 3 m 3 s 21 .Dotted lines in (a) and (c) show 1D steady state solutions; dashdotted lines in (a) show the 1D f=0 steady state solution; and the dashed red line shows the integral constraint [Eq.( 13 combination of the mean flow and the eddy fluxes as the residual circulation that advects tracers (Ferrari and Plumb 2003; see appendix).In this framing, the slope-normal eddy flux nearly doubles the residual upwelling in the BBL (Figs. 8b,a, olive lines).The crude closure for eddies in 1D sy(z) qualitatively captures this restratifying effect (Fig. 9b, compare dash-dotted and blue against solid gray) and enhances the residual BBL upwelling by a factor of 2-3 relative to the 1D model, both transiently and at equilibrium (Figs.8b,c, solid and dotted green lines, respectively).The 1D sy(z) model likely overestimates eddy restratification relative to the Smooth3D model because our simple closure (A14) for s y (z) is tuned based on a snapshot of the solution at 5000 days, when the eddies are the strongest.Even at 5000 days, the Smooth3D solution exhibits substantial buoyancy tendencies in the interior so it is possible that s y (z) will change substantially as the solution approaches steady state, invalidating our tuning of the 1D sy(z) model.Proper parameterization of bottom mixing layer eddies would avoid these issues through explicit parametric dependence on the time-varying large-scale model state, but is well beyond the scope of this paper.
The volume-integrated buoyancy budget is more complicated to interpret in the presence of variable topography.In the Canyon solution, a substantial diffusive buoyancy flux convergence drives a vigorous upslope mean flow within the bottom 200 m along the narrow trough of the canyon, producing a transport of 5 mSv (Fig. 8d, blue) which is already larger than the total BBL transport in the 1D model (Fig. 8a, blue).This strong mean flow maintains a stratification near the large background value within the canyon trough (Figs.9c, orange line; 10b).Thurnherr and Speer (2003) hypothesizes this efficient restratification is due to the canyon sidewalls blocking the along-slope thermal wind, such that the momentum is redirected into the cross-slope flow.The Canyon simulation's excellent agreement with the 1D f=0 model, in which rotation is turned off and thus the along-slope thermal wind is suppressed by construction, supports this hypothesis (Fig. 9c, orange and dotted lines).Ruan and Callies (2020) hypothesize that flow across the steep canyon sidewalls with S O(1) also contributes significantly to the strong stratification in the canyon.However, this hypothesis does not explain the strong stratification along the canyon thalweg, where the cross-canyon slope goes to zero and local dynamics cannot sustain a finite stratification at equilibrium in the absence of an along-canyon topographic slope.
The turbulent buoyancy flux also converges around a height above the mean slope (HAMS) of z = 800 m, driving an additional residual upwelling of about 13 mSv from z = 600 to 800 m dominated by the BBLs on the upper canyon sidewalls and on the smooth ridge flank surrounding the canyon (Fig. 8d, green line).The upwelling along the smooth ridge flank of the Canyon simulation is about twice as large as that of the Smooth3D simulation, despite covering a smaller area, because along-slope buoyancy gradients above the canyon sidewalls provide an additional energy source for submesoscale instabilities (Fig. 10d), driving an isopycnal thickness flux between the canyon and surrounding ridge flank and thus maintaining a much larger stratification on the flank (Fig. 9b).In the Canyon simulation's quasi-equilibrium state, much of the turbulent buoyancy flux divergence in the upper SML (far above the seafloor) is not yet equilibrated: the bottom-enhanced diffusion of buoyancy toward the boundary slowly cools the interior (Fig. 8d, red and black lines).
In the Canyon1Sill simulation, the sill blocks upslope flow within the trough of the canyon (Figs.8e,d).This is expected, since the up-canyon flows of O(1) cm s 21 only carry sufficient kinetic energy to lift a parcel across a stratification of N ∼ O(10 24 -10 23 ) s 21 by a height d Fr = U/N ∼ 20-200 m (based on a topographic Froude number of Fr ≡ Nd Fr /U ∼ 2), much smaller than the sill height of h sill = 800 m and resulting in a blocked flow layer of thickness h sill 2 d Fr (Baines 1979;Winters and Armi 2012), both up-and downstream of the sill  3c-f).Arrows show how the stratification profiles evolve when processes are added: 1) adding a mean slope, 2) allowing 3D eddies, 3) introducing a cross-slope canyon, 4) blocking the canyon with a sill, and 5) adding realistic hills (i.e., the BBTRE topography).
(recall the cross-slope periodicity).No upslope mean flow is available to restratify the trough of the canyon, so it slowly homogenizes due to mixing (Fig. 10e; as in Dell 2013) as in the nonsloping limit 1D u=0 (Figs.9c, green and dotted gray lines).Within a thin layer extending d Fr below the sill,7 however, vigorous mean flows along the upper canyon sidewalls benefit from the along-slope blocking by the canyon sidewalls and yet are also able to make it over the sill, resulting in strong near-boundary stratification (Figs.8e and 10e,f).
The structure of the stratification in the BBTRE simulation is qualitatively similar to that of the Canyon1Sill simulation, although the rougher abyssal hill topography acts to thicken the layer of enhanced stratification near the DoMORE sill height and supports a large near-bottom stratification on the hilly ridge flank surrounding the canyon (Figs.10g,h and 9b).The slope-normal structure of the BBTRE canyon's buoyancy budgets (Fig. 8f) is remarkably similar to that of the Canyon1Sill simulation and can thus be explained as the combination of the processes described above}only slightly blurred in the slopenormal direction by the additional topographic roughness.

Conclusions and discussion
By generalizing the methods of classical 1D sloping bottom boundary layer (BBL) theory (Garrett et al. 1993), we construct a hierarchy of mixing-driven flow simulations that bridge the gap between 3D (Armi 1978) and 1D (Garrett 1979) conceptual models of abyssal mixing layer restratification.Our choice to represent small-scale turbulence as a bottomenhanced turbulent diffusivity}inspired by local microstructure measurements}considerably simplifies the analysis but may not adequately represent the underlying small-scale physics (see Polzin and McDougall 2022).Nevertheless, in this conventional prescribed-diffusivity framework we demonstrate that the homogenizing tendency due to bottom-enhanced small-scale mixing is balanced by the restratifying effects of the residual overturning circulation, which is a combination of mean and submesoscale eddy flows [Eq.( 16)].At equilibrium, the slow interior diffusion of heat into the abyss is balanced by a weak net upwelling [Eq.( 13)], the result of substantial cancellation of up-and downslope flows.
The simulations' steady states are never achieved here due to the prohibitively slow diffusive adjustment in the interior (MacCready and Rhines 1991); in more realistic contexts, cross-slope pressure gradients due to coupling with the nonlocal circulation would support a much more rapid adjustment process (Peterson and Callies 2022).Despite the nonequilibrated nature of our solutions, the slope-aligned framework permits simplified buoyancy budgets which facilitate our dynamical interpretation and the derivation of a simple closure for eddy fluxes (see appendix).Another advantage of the slope-aligned framework is that the solutions are less ambiguous than previous approaches, which either require ad hoc sponge layers at distant horizontal boundaries (Dell 2013) or can only be analyzed transiently before mixing completely homogenizes buoyancy (Ruan and Callies 2020;Peterson and Callies 2022), challenging interpretations of their water-mass transformations.The slope-aligned framework also permits a consistent exploration of ever more realistic configurations: from a constant topographic slope}well described by 1D BBL models (Garrett et al. 1993)}to the complex geometry of the region surrounding the BBTRE canyon.While the local nature of the sloping BBL framework is conceptually convenient for all of the above reasons, several important nonlocal factors have been ignored.For example, the inclusion of cross-slope pressure gradients (Peterson and Callies 2022) or large-scale boundary currents (MacCready and Rhines 1991;Naveira Garabato et al. 2019) would fundamentally alter the transient spinup problem.The periodic nature of the simulation may also overemphasize topographic blocking effects since upstream topographic sills unrealistically reappear downstream.
The results of our quasi-realistic simulation of the Brazil Basin Tracer Release Experiment (BBTRE) reconciles two dominant boundary mixing paradigms: yes, bottom-enhanced mixing drives a restratifying upslope flow in the BBL (Garrett 1979(Garrett , 1990)), but this flow is much stronger than predicted by 1D theory due to net restratification by 3D transient baroclinic eddies and topographic steering/blocking effects (Armi 1978(Armi , 1979a;;Thurnherr and Speer 2003;Callies 2018;Ruan and Callies 2020).The net restratifying effect can to a large extent be attributed to three distinct physical restratification/ destratification processes: 1) slumping of isopycnals by finite-amplitude submesoscale baroclinic instabilities (Wenegrat et al. 2018;Callies 2018), 2) the blocking of cross-canyon thermal winds within narrow fracture zone canyons (Thurnherr and Speer 2003;Dell 2013;Ruan and Callies 2020), and 3) the effect of sills in blocking up-canyon mean flows and homogenizing depressions well below the sill height (Baines 1979;Winters and Armi 2012;Dell 2013).
We propose a simple closure for the restratifying effects of submesoscale baroclinic eddies in terms of a vertically varying enhancement of vertical momentum diffusion in the 1D model (see appendix).The blocking of along-slope flow by canyon walls can be captured in the 1D model by inhibiting the development of along-slope thermal wind, such as by setting f = 0.
The blocking of up-canyon flow by sills results in the endless homogenization of depressions, comparable to the diffusive growth of a bottom mixed layer over a flat bottom, i.e., setting u = 0. Above the ridge flank that surrounds the canyon, the solution is best represented by the 1D model with the eddy closure.
Applied to the BBTRE model, the slope-averaged buoyancy budget (16) confirms Thurnherr et al.'s (2020) hypothesis that spatial averaging reconciles the thin local BBL transformations implied by vertical microstructure profiles with the thicker bulk BBL transformations implied by a decreasing topographic perimeter}or mixing area}with depth (Polzin 2009;Kunze et al. 2012;Holmes et al. 2018): water below the canyon crest upwells in the net, while water above downwells (Fig. 8f).The spatial heterogeneity of the simulated up-canyon flow (Figs. 5 and 6) may explain why the buoyancy fluxes estimated from microstructure profiles are much too weak to balance the upwelling transports inferred by uniformly extrapolated moored velocity estimates (Thurnherr et al. 2005).
Our quasi-realistic simulations provide the first BBL-and submesoscale-resolving simulations of the mixing-driven abyssal overturning in the Brazil Basin, complementing Huang and Jin (2002) and Ogden and Ferrari's (K.Ogden and R. Ferrari 2022, unpublished manuscript).Despite the idealization of our numerical setup, we qualitatively reproduce key features of the observations: broad upslope flow and near-boundary stratification of B z ≈ O(10 27 ) s 22 along the canyon trough (Toole 2007;Ledwell et al. 2000), with relatively stronger nearbottom stratification along the hills surrounding the canyon (Polzin 2009); hydraulically accelerated flow over blocking sills (Clément et al. 2017); and the mean diapycnal downwelling and spreading of a tracer released in the SML (Ledwell et al. 2000;Drake et al. 2022).Despite this qualitative agreement, the simulated diapycnal transports within the canyon are too weak}and the stratification too strong}by roughly a factor of 2. These remaining discrepancies could be explained by the previously mentioned limitations of the inherently local slope-aligned modeling framework and the self-similar representation of small-scale mixing.The lack of full equilibration of the simulations could explain the too-strong stratification}the 1D models become about half as stratified at equilibrium}but not the tooweak up-canyon flow.Too-weak mixing, on the other hand, could potentially explain both biases; we speculate that microstructure-based estimates of the turbulent diffusivity may be biased low due to sampling biases (Watson and Ledwell 1988;Voet et al. 2015;Cael and Mashayek 2021;Whalen 2021) or biases in mixing parameterizations (Ijichi et al. 2020;Spingys et al. 2021).Based on and basin-scale simulations of tracer spreading, respectively, Ledwell and Ogden and Ferrari similarly conclude that tracer observations are more consistent with diffusivities about 2 times larger than those inferred from microstructure (J.Ledwell 2022, unpublished manuscript;K. Ogden and R. Ferrari 2022, unpublished manuscript).Given that estimates span several orders of magnitude spatiotemporally and are inherently uncertain, agreement within a factor of 2 is generally considered to be acceptable (e.g., Gregg et al. 2018).
The characteristic topographic features in the BBTRE (largescale slope, canyon, and hills) are typical of midocean ridges, such that the dynamics described here can be thought to apply to the global midocean ridge system (with the steepness of slopes and hills modulated by the age of the rift valley and the Coriolis parameter by its latitude).The BBTRE simulation exhibits an instantaneous diapycnal upwelling transport in the BBL of E BBL 60 mSv, where is the average water-mass transformation rate within a volume V for a layer of thickness Db; for E BBL , we confine this integral to the near-bottom regions of buoyancy flux convergence only (see Drake et al. 2022).The upwelling transport suggested by the bulk buoyancy budget presented here (Fig. 8f) is smaller than E BBL by a factor of 3 due to substantial cancellation from temporal averaging and opposing cross-slope flows at the same height above the mean slope (e.g., Fig. 4a).Extrapolating these BBL water-mass transformations to the length of the Mid-Atlantic Ridge in the Brazil Basin (about 55 times the domain width L (y) = 60 km), this 3.3 Sv of BBL upwelling8 would alone balance much of the 3.7-4.0-Svnet inflow of Antarctic Bottom Water in the Brazil Basin (Hogg et al. 1982;Morris et al. 2001).Extrapolating even further to a global midocean ridge system of length 80 3 10 3 km (including both flanks of the ridge; Thurnherr et al. 2005) leads to a global BBL upwelling of 80 Sv due to upwelling along midocean ridges, roughly consistent with global diagnostic estimates of BBL upwelling (Ferrari et al. 2016;McDougall and Ferrari 2017;Drake et al. 2020).
Global extrapolations of localized estimates of BBL upwelling, such as the above, have been used to attribute significant fractions of the net abyssal overturning to localized mixing hotspots (e.g., Ferron et al. 1998;Voet et al. 2015;Thurnherr et al. 2020;Spingys et al. 2021).These observations, however, generally also imply significant downwelling in adjacent buoyancy classes that may offset the local upwelling when integrated over a larger area.For example, Thurnherr et al. ( 2020) argue that the observed turbulent buoyancy flux convergence in the BBTRE canyon, extrapolated to all of the fracture zone canyons in the Brazil Basin, is sufficient to transform "the total inflow of AABW."Above the canyon, however, their own observations imply an opposing buoyancy flux divergence of comparable magnitude}upwelling within the canyon is thus only half of the story.Consider the following heuristic argument which applies the slope-aligned buoyancy budgets derived in section 2c in buoyancy coordinates.Following the g n ∈ {28.1, 28.15} kg m 23 neutral density class in Thurnherr et al.'s (2020) Fig. 3, we can apply Eq. ( 16) to the integrated buoyancy fluxes in their Fig.7  within the canyon at the DoMORE site.9A few hundred kilometers down-canyon, however, this same density class rests above the canyon and experiences a net buoyancy flux divergence, driving a downwelling of C(z crest 1 500 m) 2 C(z crest ) ≈ 28 mSv that partially compensates for the upwelling in the canyon and suggests a significantly weaker net upwelling of 12 mSv for the BBTRE canyon.This heuristic exercise serves as a cautionary tale for attributing abyssal upwelling to individual regions or processes: both strictly positive and strictly negative components of watermass transformations along a buoyancy surface must be accounted for to robustly characterize the net overturning circulation.At a global scale, diagnostic estimates of water-mass transformations suggest significant compensation is the norm, exhibiting typical amplification factors of A ≡ E BBL /E of 2-5, where E E BBL 1 E SML is the net diapyncal transport and E SML is the downwelling in the stratified mixing layer (Ferrari et al. 2016;McDougall and Ferrari 2017;Cimoli et al. 2019).However, these diagnostic exercises do not provide any insight into the physics underlying the observed density structure that supports these transformations.More problematically, these results seem to contradict the weak upwelling with A 1 implied by 1D boundary layer dynamics (section 2b).Building upon Callies (2018) and Ruan and Callies (2020), our prognostic modeling approach demonstrates how three-dimensional eddy and topographic effects conspire to provide sufficient restratification to support a significant upwelling/downwelling dipole, i.e., A . . 1 (Figs.8a,f).Our results inspire two open questions: 1) which topographic regimes (e.g., ridges, slopes, plains) or topographic roughness features (e.g., hills, canyons, channels, sills, or seamounts) contribute the most to abyssal water-mass transformations (e.g., Armi and D'Asaro 1980;Bryden and Nurser 2003;Thurnherr et al. 2005;Legg et al. 2009;Nazarian et al. 2021;Mashayek et al. 2021) and 2) what are the equilibrium dynamics that support finite water-mass transformations in these regions (Garrett 1979(Garrett , 1990;;Callies 2018;Drake et al. 2020)?
Our combined assumptions of constant background stratification and zero barotropic cross-slope pressure gradient assert that the net upwelling scales with the background diffusivity [Eq.( 13)] and thus that the net upwelling C ' E is very small.While our local model helps explain the magnitude of bottom boundary layer upwelling E BBL , it does not meaningfully constrain E SML or A. Salmon et al. (1991) use asymptotic analysis to show that small perturbations away from a constant interior stratification drive an exchange flow between the boundary and the interior, which then feeds back on the interior stratification.In the context of the abyssal ocean, vertical variations in the basin-scale interior stratification are relatively large, such that they enter as leading-order terms in water-mass transformations (Spingys et al. 2021) and drive substantial exchange between the mixing layers and the interior (Holmes et al. 2018).In Drake et al.'s (2020) idealized basin-scale simulations, this boundary-interior coupling results in a substantial reduction of E SML , permitting an amplification factor of A 1:5 much smaller than the A . . 1 governed by local dynamics.These idealized prognostic model results are qualitatively consistent with the diagnostic approaches described above, but quantitative understanding of E BBL , E SML , and A remains incomplete.
Understanding of bottom-enhanced mixing has advanced considerably in recent years due to a combination of breakthroughs in observation (e.g., Polzin et al. 1997;Ledwell et al. 2000), theory (e.g., Polzin 2009), and modeling (e.g., Nikurashin and Legg 2011).The interpretation of these results in terms of broad diapycnal downwelling in the SML atop vigorous diapycnal upwelling in a BBL (Ferrari et al. 2016), however, is challenged by higher-resolution observations (van Haren 2018;Polzin et al. 2021) and simulations (Gayen and Sarkar 2011;Kaiser 2020) of mixing processes within the bottom few dozen meters of the ocean.This picture is further complicated by the realization that submesoscale instabilities may be common within sloping bottom boundary layers (Ruan et al. 2017;Naveira Garabato et al. 2019;Wenegrat et al. 2018); our results here demonstrate their importance in a realistic context and suggests a path toward parameterization.In addition to the debate on the nature of boundary mixing itself (see Polzin and McDougall 2022), the role of the resulting boundary layer flows in the global overturning circulation remains shrouded by poor understanding of the coupling between bottom mixing layers and the far-field interior (Drake et al. 2020;Peterson and Callies 2022).
FIG. 2. A generalized slope-normal buoyancy budget (16), derived by integrating the buoyancy equation below a given height above the mean slope z (volume shown in light blue); at equilibrium, the mean upslope transport (across the background stratification N 2 ) into the box (blue lines) is given by the net flux of buoyancy into the box from above (red line), C∝ 2 2kB z 2 wb .We assume no buoyancy flux across the seafloor (black line) at z = max(d) 2 d(x, y).

FIG. 3 .
FIG. 3. Numerical model domains.(a) Seafloor elevation 2 d(x, ŷ), including the doubly periodic simulation domain centered on the Brazil Basin Tracer Release Experiment (BBTRE) canyon.Red markers show the locations of moorings from Clément et al. (2017) (CTS17) and Thurnherr et al. (2005) (T05).The inset highlights the DoMORE sill that dramatically constrains up-canyon flow.White markers mark the injection location from the BBTRE (Ledwell et al. 2000).(b) Imposed slope-normal diffusivity field, the result of applying a self-similar exponential profile as a function of the height above bottom [Eq.(17)] to variable topography.Arrows show the original along-canyon ŷ and cross-canyon x directions in (a) as well as the transformed slope-normal z and along-canyon x coordinate vectors in (b), which appear distorted because the vertical dimension is exaggerated.(c)-(f) A hierarchy of simulations with progressively complex seafloor bathymetry geometries [relative to a constant mean slope of angle u; see dashed lines in (b)].Thin black lines distinguish three subregions: the canyon trough, the canyon's sloping sidewalls, and the ridge flank surrounding the canyon.

FIG. 4 .
FIG. 4. Structure of up-canyon mean flow in the BBTRE canyon.(a) Along-canyon-averaged up-canyon flow u x , with the mean canyon seafloor outlined in transparent gray shading and cross-canyon thalweg shown in the dark gray shading.(b) Cross-canyon-averaged up-canyon flow u y in the original coordinate frame (x, ẑ) Gray lines represent equally spaced buoyancy surfaces.The much gentler isopycnal slopes suggested by some visualizations of hydrographic section in canyons, as in Thurnherr et al. (2020), are largely an artifact of their much lower horizontal resolution, as evidenced by the favorable comparison of vertical density profiles in Fig. 6.The black line marks the mean seafloor depth of the half of the domain furthest from the canyon thalweg and acts as a proxy for the crests of the canyon sidewalls.(c) Up-canyon flux, integrated in the slope-normal direction z.Black and gray contours show a deep and shallow isobath, respectively, to highlight the canyon topography that shallows to the right.(d) As in (c), but integrated only from z = 0 to 800 m to highlight the core up-canyon jet within the canyon.
FIG. 5. Structure of up-canyon flow at two mooring sites.(a),(c) Cross-canyon sections of the up-canyon flow at the locations of the DoMORE sill mooring (Clément et al. 2017) (CTS17-P1) and the BBTRE mooring (Thurnherr et al. 2005) (T05).Light gray shading shows the local seafloor depth while the dark gray shading in (a) shows the mean height of the canyon floor above the mean slope, highlighting the significant vertical and cross-canyon constriction introduced by the sill.Gray lines denote equally spaced buoyancy surfaces.(b),(d) Height-above-bottom profiles of the up-canyon flow at the locations of the two moorings (light gray lines) and shifted a few grid columns over to improve capture the core of the jet (black lines), which is somewhat displaced due to the coarse model bathymetry.Red crosses denote Aquadopp current meter measurements, and the red (dotted) line denotes (extrapolated) McLane Moored Profiler (MMP) measurements.

FIG. 6 .
FIG. 6.Comparison between observed and simulated stratification in the BBTRE canyon region.(a) Height-above-bottom averaged profile of stratification for the full simulation domain (solid blue), for the sample-mean of 9 collocated CTD casts (dotted red; Ledwell et al. 2000), free-falling HRP-microstructure profiles (dashed red; Polzin et al. 1997), and simulated CTD casts (solid red).The dashed blue line shows the sample-mean of 10 HRP profiles that follow the hilly ridge flank just north of the domain.(b) Observed (solid) and simulated (dashed) density profiles at the 9 locations sampled by the BBTRE observational campaign, overlaid on a map of the simulation's seafloor elevation.An additional simulated profile typical of the hilly ridge flank surrounding the canyon is also shown, revealing an apparent sampling bias due to the strategy of measuring weakly stratified deep depressions along the trough of the canyon in search of the released tracer(Ledwell et al. 2000).

FIG. 8 .
FIG.8.Generalized integral buoyancy budget in a hierarchy of increasingly complex simulations of mixing-driven flows up a mean slope of angle u: (a) 1D, (b) Smooth3D, (c) 1D sy(z) , (d) Canyon, (e) Canyon1Sill, (f) BBTRE.Solid lines show terms of the volume-integrated buoyancy budget [Eq.(16)], averaged over days 5000 to 5200, for a layer bounded by a given HAMS.We interpret the sum of the mean flow and eddy terms as a residual flow.The left-hand-side tendencies (LHS) are equal to the remainder of the approximate balance (RHS) between slope-normal turbulent diffusion and the cross-slope residual circulation, which includes both mean and eddy components.We divide [Eq.(16)] by the factor 2N 2 L (x) sinu to conveniently express the budget in terms of the quantity of interest, the upslope volume transport C with units of mSv ≡ 10 3 m 3 s 21 .Dotted lines in (a) and (c) show 1D steady state solutions; dashdotted lines in (a) show the 1D f=0 steady state solution; and the dashed red line shows the integral constraint [Eq.(13)].[In (a) and (b), some of the dotted lines appear missing because they overlap with others.]Gray shading shows the HAMS range spanned by the canyon, if present.
FIG.8.Generalized integral buoyancy budget in a hierarchy of increasingly complex simulations of mixing-driven flows up a mean slope of angle u: (a) 1D, (b) Smooth3D, (c) 1D sy(z) , (d) Canyon, (e) Canyon1Sill, (f) BBTRE.Solid lines show terms of the volume-integrated buoyancy budget [Eq.(16)], averaged over days 5000 to 5200, for a layer bounded by a given HAMS.We interpret the sum of the mean flow and eddy terms as a residual flow.The left-hand-side tendencies (LHS) are equal to the remainder of the approximate balance (RHS) between slope-normal turbulent diffusion and the cross-slope residual circulation, which includes both mean and eddy components.We divide [Eq.(16)] by the factor 2N 2 L (x) sinu to conveniently express the budget in terms of the quantity of interest, the upslope volume transport C with units of mSv ≡ 10 3 m 3 s 21 .Dotted lines in (a) and (c) show 1D steady state solutions; dashdotted lines in (a) show the 1D f=0 steady state solution; and the dashed red line shows the integral constraint [Eq.(13)].[In (a) and (b), some of the dotted lines appear missing because they overlap with others.]Gray shading shows the HAMS range spanned by the canyon, if present.

FIG. 9
FIG. 9. (b),(c) Height above bottom-averaged stratification profiles at t = 5000 days, as a function of model complexity (lines) and domain subregion.(a) 1D solutions with the same parameters as the BBTRE simulations (solid); without a mean-slope (u = 0; dashed), without rotation ( f = 0; dotted); and with an enhanced along-slope turbulent Prandtl number s y (z), a crude proxy for restratification by submesoscale baroclinic eddies (dash-dotted) [also shown in gray lines in (b) and (c)].Colored lines show a hierarchy of 3D simulations with increasingly complex topographies (see Figs.3c-f).Arrows show how the stratification profiles evolve when processes are added: 1) adding a mean slope, 2) allowing 3D eddies, 3) introducing a cross-slope canyon, 4) blocking the canyon with a sill, and 5) adding realistic hills (i.e., the BBTRE topography).

FIG. 10
FIG. 10. (a),(c),(e),(g) Cross-slope and (b),(d),(f),(h) along-slope sections of the stratification along the trough of a canyon in a hierarchy of numerical simulations (Smooth3D has no canyon, so the section is arbitrary).Solid gray lines in the left column show the approximate elevation of the ridge flank surrounding the canyon while in the right column they show HAMS of the topographic sill (if present).Dashed gray lines show the locations of the respective sections.Black lines in (d) represent equally spaced buoyancy surfaces.
to infer a bulk upwelling ofC(z crest ) kB z /L (x)