## Abstract

Weakly stratified layers over sloping topography can support a submesoscale baroclinic instability mode, a bottom boundary layer counterpart to surface mixed layer instabilities. The instability results from the release of available potential energy, which can be generated because of the observed bottom intensification of turbulent mixing in the deep ocean, or the Ekman adjustment of a current on a slope. Linear stability analysis suggests that the growth rates of bottom boundary layer baroclinic instabilities can be comparable to those of the surface mixed layer mode and are relatively insensitive to topographic slope angle, implying the instability is robust and potentially active in many areas of the global oceans. The solutions of two separate one-dimensional theories of the bottom boundary layer are both demonstrated to be linearly unstable to baroclinic instability, and results from an example nonlinear simulation are shown. Implications of these findings for understanding bottom boundary layer dynamics and processes are discussed.

## 1. Introduction

An important development in our understanding of the upper ocean, occurring largely over the last decade, has been the recognition that submesoscale (horizontal scales of approximately 0.1–10 km) fronts, eddies, and instabilities are a prominent component of the dynamics of the ocean surface mixed layer (Boccaletti et al. 2007; Callies et al. 2015). These submesoscale processes alter the classic one-dimensional picture of boundary layer dynamics, giving rise to a host of new physical processes that have been shown to affect both large-scale ocean processes (Lévy et al. 2010; Wenegrat et al. 2018) and the turbulence properties of the boundary layer itself (Taylor and Ferrari 2010; Thomas and Taylor 2010; Taylor 2016). However, despite the large body of literature that has developed on submesoscale processes in the surface boundary layer [as reviewed in Thomas et al. (2008) and McWilliams (2016)], these types of processes have received much less attention in the bottom boundary layer (BBL),^{1} despite BBLs over sloping topography exhibiting key similarities with surface mixed layers at a front, particularly the existence of available potential energy in the form of a horizontal buoyancy gradient (Fig. 1).

Recently, however, observations and numerical modeling have begun to suggest that the BBL does indeed support an active submesoscale turbulence field, with dynamical implications potentially as broad as for submesoscale processes in the surface boundary layer. Submesoscale processes in BBLs along topography have been shown to enhance cross-shelf exchanges (Gula et al. 2015), contribute to interior water mass transformation (Ruan et al. 2017), generate long-lived submesoscale coherent vortices (Molemaker et al. 2015), and provide another route for cascading energy from the large-scale flow into unbalanced motions and eventual dissipation (Gula et al. 2016). Particular focus has been given to centrifugal (CI) and symmetric instabilities (SI) in the BBL, both of which are associated with negative potential vorticity (Hoskins 1974; Haine and Marshall 1998), which can occur in BBLs over steep topography (Dewar et al. 2015) or in the presence of interior flows driving downslope Ekman transport in the boundary layer (Allen and Newberger 1998). The role of baroclinic instability in the BBL is currently less clear, despite the fact that surface mixed layer baroclinic instabilities are one of the more thoroughly studied aspects of submesoscale dynamics, with demonstrated effects on boundary layer restratification (Boccaletti et al. 2007; Fox-Kemper et al. 2008), surface potential vorticity fluxes (Wenegrat et al. 2018), boundary layer turbulence (Taylor 2016), biological productivity (Mahadevan et al. 2012), and the mesoscale eddy field (Sasaki et al. 2014). Understanding the conditions under which baroclinic instability can be expected to be active in the BBL is therefore a key step toward understanding submesoscale turbulence in the BBL.

In this article, we extend a variety of earlier works on baroclinic instability over topography (e.g., Blumsack and Gierasch 1972; Mechoso 1980; Pedlosky 2016; Solodoch et al. 2016) to the nongeostrophic regime appropriate for the low Richardson numbers typical of the BBL. Using a linear stability analysis, we explore the parameter dependence of baroclinic instabilities in a weakly stratified lower layer over sloping topography, which can be considered as a bottom counterpart to the surface mixed layer instability (Boccaletti et al. 2007). The submesoscale BBL mode is shown to be more robust to topographic slope than mesoscale baroclinic instability, and we demonstrate that the solutions to two classic one-dimensional theories of the BBL structure are unstable to BBL baroclinic instability. An example idealized nonlinear simulation further suggests that at finite amplitude the BBL baroclinic instability mode can generate vertical buoyancy fluxes, and vertical velocities, similar in magnitude to those associated with surface mixed layer instabilities. It is therefore expected that submesoscale baroclinic instability, along with the symmetric and centrifugal modes (Allen and Newberger 1998; Molemaker et al. 2015; Dewar et al. 2015; Gula et al. 2016), plays a leading-order role in the dynamics of both coastal (Brink 2016; Brink and Seo 2016; Hetland 2017), and deep ocean BBLs, with implications for closing the upwelling limb of the abyssal circulation (Ferrari et al. 2016; Callies 2018).

The article is organized as follows. In section 2 the equations governing linear perturbations in the BBL are developed, relevant quasigeostrophic results are briefly reviewed, and the parameter dependence of the baroclinic mode in the nongeostrophic limit is explored numerically. In section 3 a steady one-dimensional solution for a BBL in the presence of bottom-intensified turbulence, with parameters similar to those observed along the mid-Atlantic ridge, is shown to support growing baroclinic instability in a weakly stratified outer layer on the order of hundreds of meters thick. In section 4 numerical solutions of the time-dependent, one-dimensional, Ekman adjustment problem for flow along a slope are likewise shown to be unstable to submesoscale baroclinic instability in a well-mixed bottom boundary layer that is on the order of tens of meters thick. An example nonlinear simulation is presented in section 5, allowing us to comment briefly on the finite-amplitude behavior, and in section 6 the baroclinic mode is discussed in relation to symmetric and centrifugal instability, both of which can also be present in BBLs over sloping topography. Major findings and broader implications are summarized in section 7.

## 2. Baroclinic instability in the bottom boundary layer

### a. Theory

We consider the stability characteristics of a linearized Boussinesq system, in a coordinate system rotated to align with a linear slope in the direction (Fig. 1), where throughout the paper we will indicate quantities in the conventional, nonrotated, coordinate system using a caret notation, background fields using uppercase variables, and perturbation quantities using lowercase variables. The slope-normal coordinate is defined positive upward. An arbitrary background flow can be expressed in the rotated frame as , where *θ* is the angle of the topographic slope (Fig. 1). The equations governing the linear evolution of the perturbations are then

where *f* is the Coriolis parameter, ν is the viscosity, κ is the diffusivity, ∇*χ* = (*χ*_{x}, *χ*_{y}, *χ*_{z}) is the gradient operator in the rotated frame, and subscripts denote differentiation. Note that, consistent with Stone (1971), we retain nonhydrostatic effects but make the traditional approximation for the Coriolis force, ignoring terms involving the horizontal components of the Coriolis parameter. This is an accurate approximation for the problems considered here, and further allows the problem to remain invariant under rotation of the slope direction. Boundary conditions depend on the particular problem configuration and will be discussed in the relevant following sections.

For simplicity, we will only consider background flows that are invariant in the rotated along- and across-slope direction, which excludes slope-normal background flows (i.e., *W* = 0). These restrictions are consistent with the conceptual model of a BBL of uniform thickness along a linear slope (sections 3 and 4) and still allow for a vertically and horizontally sheared background flow in the nonrotated coordinate system (where the horizontal shear is a consequence of the sloping topography, as shown schematically in Fig. 1b). We will also retain only the slope-normal component of the turbulent diffusion terms, which for an isotropic turbulent diffusivity can be understood as an assumption of small aspect ratio, *H*^{2}/*L*^{2} ≪ 1, where *H* and *L* scale the slope-normal and rotated horizontal dimensions, respectively. This choice does not affect the conclusions of this article; however, we emphasize that the properties of diapycnal and isopycnal turbulent mixing in the BBL are currently poorly constrained from available observations, and it is possible that under some circumstances—for example, along very steep topography—the full three-dimensional diffusion operator may be important.

Assuming perturbations of the form gives the eigenvalue problem,

Given a specified background state, the eigenvalues of (6)–(10) can be found numerically. We do this by projecting the governing equations onto a Chebyshev basis in *z* (*n* = 256 in all cases) and then finding the eigenvalues of the resulting matrix using the Dedalus equation-solving framework (Burns et al. 2016).

Later in the article, general profiles of background velocity and buoyancy will be considered using (6)–(10) directly; however, to better illustrate the problem parameter dependence, it is useful to assume a background flow in the along-slope direction, with uniform vertical shear, , and stratification, , such that *V* = Λ*z* sec*θ*, and *B* = *N*^{2}(*z* cos*θ* + *x* sin*θ*). Equations (6)–(10) can then be nondimensionalized following Stone (1971) using , , , , , , , and , with primes denoting nondimensional quantities. Dropping grave and prime accents for clarity then gives

The relevant nondimensional parameters are the Richardson number of the background flow, Ri = *N*^{2}/Λ^{2}; a nonhydrostatic parameter, *δ* = *f*/Λ; the Ekman number, *E* = *A*_{υ}/(*fH*^{2}); and a slope parameter,

which gives the ratio of the topographic slope to minus the isopycnal slope.^{2} Alternately, the slope parameter can be expressed as *α* = (*S*Ri)^{1/2}, where *S* = (*N* tan*θ*/*f*)^{2} is the slope Burger number. For the analyses of this paper, we retain nonhydrostatic terms (except for the nonlinear simulation of section 5); however, the nonhydrostatic parameter is generally small, even for weakly stratified layers, and hence the inviscid stability characteristics of this system of equations are primarily determined by the Richardson number and *α*.

To motivate the following sections, we first consider the baroclinic growth rates for an idealized profile of stratification (*N*^{2}) and shear (Λ), with varying slope angle. Figure 2 shows profiles of background velocity and buoyancy, where we assume there is an inviscid along-slope flow, *V*(*z*) = Λ*z* sec*θ*, and a stratified interior overlying a weakly stratified BBL. Note that these profiles are held fixed in the nonrotated frame, hence the isopycnal slope in the rotated frame changes as a function of the slope angle. The inviscid problem only requires boundary conditions on the slope-normal velocity, which are given by *w* = 0 at *z* = 0 and *w* = −*u* tan*θ* at *z* = 1000 m, representing a rigid horizontal upper boundary 1000 m from the bottom. Growth rates of the baroclinic mode (*k* = 0) are calculated numerically and are shown as a function of increasing slope angle in Fig. 3. The full-depth baroclinic mode can be seen at wavelengths of ~30 km. This mode is strongly modulated by slope angle, with maximum growth rates reduced, and shifted to higher wavenumber, as the slope angle increases, consistent with the expectation from quasigeostrophic (QG) theory (discussed further below). The fastest-growing instability, however, is found at wavelengths near 1500 m, with perturbation structure and energetics shown in Fig. 4. This is a bottom-enhanced submesoscale baroclinic instability, resulting from the interaction of counterpropagating Rossby waves along the lower-boundary and on the interface of changing stratification between the boundary layer and the interior. These ocean bottom modes are thus similar to baroclinic instabilities in the atmospheric boundary layer (Nakamura 1988) and mirror the structure of surface mixed layer instabilities (Boccaletti et al. 2007; Callies et al. 2016). As shown in Fig. 3, the growth rates of the submesoscale mode are relatively insensitive to increasing slope angle, suggesting that submesoscale baroclinic instability may be a robust feature of the bottom boundary layer. To understand this behavior, we will first briefly review relevant results from the QG limit and then explore parameter space in the non-QG limit relevant for submesoscale instabilities.

### b. Quasigeostrophic baroclinic instability over sloping topography

An analytical solution for the growth rate of inviscid baroclinic instabilities over sloping topography in the QG limit (assuming small Rossby number, large Richardson number, and small slope angle) was first provided by Blumsack and Gierasch (1972), considering instabilities in the Martian atmosphere. This work was later extended by Mechoso (1980) to include a sloping upper boundary, introducing an additional parameter,

the ratio of the upper boundary slope *θ*^{ub}, to the negative of the isopycnal slope. The case where *α*^{ub} = *α* can therefore be considered as a simple prototype model of baroclinic instability in a bottom boundary layer that is of uniform thickness in the across-slope direction, ignoring deformation of the upper boundary of the BBL (Fig. 1b). The use of a rigid lid approximation is accurate as long as the stratification contrast between the lower-layer and interior is large (Eady 1949; Boccaletti et al. 2007; Callies et al. 2016), in which case including a deformable upper boundary adds an additional free parameter and leads only to small quantitative modifications. Later in the paper, the assumption of a rigid boundary on the lower-layer will be removed. For the case of uniform stratification and vertical shear, the instability growth rate is given by (Mechoso 1980)

where *ω*_{i} is the growth rate normalized by *f*^{−1}, and we have introduced an alternate normalization of the wavenumber using the deformation radius *NH*/*f* instead of Λ*H*/*f* (such that ), as is appropriate for quasigeostrophic dynamics. This alternate normalization will prove useful for comparing the stability properties with varying Ri (below).

These growth rates are shown as a function of *α* in Fig. 5 for the case of a horizontal upper boundary (left) and for a sloping upper boundary, *α*^{ub} = *α* (right). In both cases there is a region of instability for any *α* > −1, with increasing positive values of *α* associated with reduced growth rates, and an increasingly narrow region of instability at high wavenumbers, as seen for the low-wavenumber modes in Fig. 3. While linear theory predicts finite growth rates at arbitrarily high *α*, in reality nonlinear (Pedlosky 2016), viscous, and non-QG (section 2c) effects may provide an upper bound on physically relevant wavenumbers. With an upper boundary parallel to the topography, there is no longer a barotropic potential vorticity (PV) gradient, and hence there is no longer a long-wave cutoff for *α* < 0 (Fig. 5, right). The sloping boundaries also lead to a destabilization of low wavenumbers for *α* < 0, with maximal growth rates for *α* = −0.5. This is a consequence of parcel trajectories being forced by the boundary slope to cross isopycnals in a manner that achieves maximal extraction of available potential energy (Mechoso 1980). The increase of growth rates as can be understood as a resonance between baroclinic and barotropic Rossby waves (Ripa 1999, 2001), associated with a uniform downslope flow extracting available potential energy. These results will now be compared to the non-QG limit.

### c. Nongeostrophic baroclinic instabilities over sloping topography

To explore parameter space in the non-QG limit, using (11)–(15) in the rotated frame, we continue to assume the flow is inviscid, with uniform background vertical stratification, *N*^{2}, and along-slope flow with uniform slope-normal shear. This allows the background velocity field to be written in the rotated frame as *V*(*z*) = Λ*z* sec*θ*, and in the nonrotated frame as , which has both vertical and horizontal vorticity (section 6). We also assume a rigid upper boundary, aligned parallel to the topography (i.e., *α*^{ub} = *α*, such that the boundary conditions are *w* = 0 at *z* = 0, *H*), a configuration that approximates a boundary layer of uniform thickness across the slope (Fig. 1, lower panel). The focus here is on the baroclinic mode, and hence we assume perturbation variables are uniform in the across-slope direction (*k* = 0); other instability modes (*k* ≠ 0) are discussed in section 6. This simplified model can thus be thought of as Stone’s (1966) nongeostrophic baroclinic instability problem adapted for a slope. In the inviscid and hydrostatic limit, the stability characteristics of (11)–(15) are determined by Ri and *α*, which are varied below.

In the nondimensionalization used here, *l* acts as a Rossby number, quantifying advection, vis-à-vis Doppler shifting of the perturbations by the mean flow, relative to the Coriolis term. Non-QG effects will therefore become important when , which indicates that the flow is evolving quickly relative to the inertial time scale, violating the QG assumption of slow dynamics. This can be seen in Fig. 6 for Ri = 10 and Ri = 1. In each case the low wavenumber growth rates are similar to the expectation from QG theory, whereas wavenumbers with (or in the plots ) have growth rates reduced by non-QG effects. These effects are most relevant for positive slope angles, *α* > 0, where growth rates in the high-wavenumber tail are strongly reduced. Considering the instability energetics, topographic reduction in growth rates for *α* > 0 can result from a reduction of the vertical buoyancy production by the across-slope advection of buoyancy perturbations ( appendix A). This reduction of the vertical buoyancy production scales as *α*/Ri and is therefore negligible in the QG limit, but can become significant for small Ri.

An alternate, physically intuitive, way to consider the parameter dependence of the instabilities is to hold the slope angle and horizontal buoyancy gradient fixed (Fig. 7). Varying the stratification thus varies the Richardson number and *α*, approximating interior stratification encountering deep topography, and adjusting within a turbulent boundary layer (sections 3 and 4). At low Richardson numbers the growth rates increase, and higher wavenumbers become unstable, consistent with the decreasing stratification allowing greater vertical penetration and coupling of the boundary waves. Assuming values typical of interior flows (*V* ~ 0.1 m s^{−1}, *f* = 10^{−4} s^{−1}) the most unstable modes for Ri < 10 have a horizontal scale of *O*(5) km and growth rates on the order of 0.5–3 inertial periods.

It is important to note that for isopycnals that intersect the topography at right angles, a reasonable first approximation to the BBL structure, the slope parameter is given by *α* = tan^{2}*θ*, which for small slope angle can be approximated as *α* ≈ *θ*^{2}. It is estimated that approximately 95% of the world’s ocean has (Costello et al. 2010), implying small *α* values will be most physically relevant in the BBL. Further, for a BBL with Ri = 1, as is anticipated for a bottom boundary layer that has adjusted after turbulent mixing (Haine and Marshall 1998; Fox-Kemper et al. 2008; Taylor and Ferrari 2009), *α* = *S*^{1/2}. Observations of weakly stratified BBLs in both coastal and deep oceans suggest that *S* ≪ 1 (Armi and Millard 1976; Stahr and Sanford 1999; Moum et al. 2004), and hence *α* ≪ 1 may be common in the BBL. Finally, the change in the slope parameter with increasing slope angle is given by ∂*α*/∂*θ* = *N*Ri^{1/2} sec^{2}*θ*/*f*, which is small for BBLs with low stratification and low Richardson numbers. This suggests that the properties of the BBL instability are less sensitive to changes in slope than the full-depth instability, for which stratification and Richardson number are larger. This is consistent with the robustness of the submesoscale mode shown in Fig. 3. Thus, despite the reduction of growth rates at large *α* evident in Fig. 6, the growth rates of baroclinic instability can remain significant in weakly stratified layers over topography, where both Ri and *α* are small (sections 3 and 4), and the submesoscale mode is anticipated to be a robust feature of the ocean bottom boundary layer.

Up until this point we have considered inviscid flows with greatly simplified vertical structures for the background buoyancy and shear fields, as an aid to interpreting the problem parameter dependencies. In the following sections, we will consider the stability characteristics of the bottom boundary layer structure predicted by two somewhat more realistic dynamical theories of the ocean bottom boundary layer. The first results from bottom intensified turbulent mixing, the second from Ekman adjustment of the bottom boundary layer to an imposed interior flow.

## 3. Bottom-intensified turbulent mixing

Available potential energy near the bottom boundary can be generated by turbulent mixing associated with processes such as the breaking of internal waves over sloping topography. Observations show that turbulence over rough topography is strongly bottom enhanced, suggesting there should be a dipole of vertical velocities along topography, with downwelling in an outer layer a few hundred meters thick, where diapycnal buoyancy fluxes are divergent, and upwelling in an inner layer adjacent to the bottom, where diapycnal buoyancy fluxes are convergent (e.g., Polzin 1997; Ledwell et al. 2000; St. Laurent et al. 2001; Ferrari et al. 2016; McDougall and Ferrari 2017). However, 1D boundary layer solutions, with turbulent viscosities based on observations, imply outer-layer stratification much weaker than observed, and hence weak vertical velocity dipoles (Callies 2018). Baroclinic instabilities in the mixing layer would alter the 1D buoyancy budget, with important implications for closing the upwelling branch of the abyssal circulation, and hence here we consider the stability characteristics of 1D boundary layer solutions with turbulent diffusivities based on observations.

To do this we solve the governing equations for the steady background state over sloping topography, with uniform interior stratification, and no across-slope variations of the boundary layer quantities,

where denotes a specified interior stratification far from the boundary, and the total background buoyancy field is given by , such that represents the portion of the buoyancy not associated with the imposed interior stratification. Boundary conditions are

Practically, the upper-boundary conditions are applied at a finite height, sufficiently far from the lower boundary to not influence the dynamics, which here we set to be *z* = 2500 m. These equations are solved numerically, using Dedalus, with an idealized, bottom-enhanced profile of turbulent mixing,

Parameters are chosen to approximate observations taken during the Brazil Basin Tracer Release Experiment (Ledwell et al. 2000):

These solutions are used to define the background fields in the eigenvalue problem, (6)–(10), with perturbation boundary conditions given by

The free-slip upper-boundary conditions on momentum are necessary to satisfy a no-stress boundary condition at the surface (Thomas and Rhines 2002; Wenegrat and McPhaden 2016).

The solution for the background fields is shown in Fig. 8. There is an inner layer of thickness , where isopycnals are approximately normal to the slope (Callies 2018). This layer is approximately 6 m thick for these parameters. Above this is an outer layer, with thickness , approximately 1000 m for these parameters. Stratification in the outer layer is reduced by a factor of from the interior values, where is the interior slope Burger number (Garrett 1991). The across-slope momentum balance in the outer layer is approximately , which, when combined with the force balance, , gives , that is, the along-slope flow is in geostrophic balance due to the projection of the slope-normal gradient on the true horizontal. The lower portion of the outer layer has Ri ~ *O*(10–100) and hence is stable to symmetric instability, but, as shown in Fig. 9, it can support growing baroclinic instability. The fastest-growing instability for these parameters has a wavelength of ~5 km and an *e*-folding time scale of ~5 days. For comparison, assuming an inviscid layer with a rigid lid at *H* = 500 m, Λ*H* = −0.015 m s^{−1}, and Ri = 25, the simplified analysis of section 2c would predict a fastest-growing wavelength of ~6 km and a growth rate of ~3 days.

The perturbation structure and kinetic energy tendency of the fastest-growing mode are shown in Fig. 10. Phase lines are inclined into the mean shear (note *f* < 0), and the instabilities grow by releasing available potential energy in the outer layer, with maximum vertical buoyancy fluxes near *z* = 250 m. Perturbations quantities are enhanced near the bottom boundary, except for in the thin inner layer where dissipation dominates the energetics. At finite amplitude the eddy vertical buoyancy fluxes will provide a restratifying tendency, altering the simple 1D balance given by (19)–(21), and potentially modifying the structure of the vertical velocity dipole in the boundary layer. The net secondary circulation along sloping topography may therefore involve both eddy and Eulerian components, with implications for the dynamics of the abyssal circulation (Callies 2018).

## 4. Stability of the bottom Ekman layer

The classic analysis of the laminar bottom boundary layer holds that, in the presence of an interior flow along the slope, a bottom Ekman layer must also develop to satisfy the bottom boundary condition on the momentum (Pedlosky 1979). Interior flow thus gives rise to across-slope flow in the Ekman layer that advects the across-slope buoyancy gradient, modifying the depth and stratification of the boundary layer (MacCready and Rhines 1991, 1993). For interior flow in the direction of topographic wave propagation (topography sloping upward to the right of the direction of flow in the Northern Hemisphere) the resulting Ekman flow is downslope, advecting light water under dense, causing convective mixing and deepening of the BBL (Garrett et al. 1993; Moum et al. 2004), an analog to the surface “Ekman buoyancy flux” associated with downfront winds at a surface mixed layer front (Thomas 2005). For interior flow in the opposite direction, the Ekman flow is upslope, increasing boundary layer stability by advecting dense water under light. In each case this adjustment process generates available potential energy in the BBL (Umlauf et al. 2015), and we demonstrate here that this available potential energy can allow for growing baroclinic instability in the BBL.

Following MacCready and Rhines (1993), we consider the time-dependent adjustment of a BBL over a linear slope. The initial condition is a barotropic, geostrophically balanced interior flow and uniform stable vertical stratification . The dimensional along-slope velocity can then be written as , with governing equations,

where is defined in section 3 following (21). Boundary conditions on the background flow are given by (22) and (23), and on the perturbations by (27) and (28), with the upper boundary applied at *z* = 150 m, which is sufficiently far from the lower boundary to not modify the characteristics of the BBL solution. Solutions to (29)–(31) are found using the General Ocean Turbulence Model (GOTM; Umlauf et al. 2005), modified to solve the rotated equations of motion. The barotropic interior flow is applied as an initial condition and is maintained during integration by imposing a steady barotropic pressure gradient force in the across-slope momentum equation. Turbulence closure is provided using a second-moment *k*–*ε* scheme, with a prognostic equation for the turbulence kinetic energy, and background viscosity and diffusivity of *ν* = *κ* = 10^{−5} m^{2} s^{−1}. The bottom boundary conditions on momentum are implemented using a log-layer formulation. Detailed discussion of the use of this turbulence parameterization in a BBL over topography can be found in Umlauf and Burchard (2011).

### a. Downwelling favorable interior flow

The numerical solution for the case of downwelling-favorable interior flow with magnitude typical of a deep western boundary current (Toole et al. 2011), = 0.1 m s^{−1}, is shown in Fig. 11. A well-mixed turbulent bottom boundary layer is rapidly established by convective mixing due to the across-slope buoyancy advection. This leads to a growing BBL, which reaches ~30 m thickness after only approximately five inertial periods. As the BBL thickness increases, the along-slope flow approaches thermal wind balance, and *U*(*z*) → 0, a state known as the “arrested” Ekman layer. An estimate of the Ekman shutdown time scale can be formed as (MacCready and Rhines 1993), where *M*^{x} is the across-slope transport, which for this simulation gives *τ*^{Ek} ~ 10 inertial periods.

Instability growth rates (Fig. 12) are calculated using snapshots of the numerical solutions for *U*, *V*, *B*, *κ*, and *ν* as the background fields in the eigenvalue problem, (6)–(10). The baroclinic mode emerges after approximately five inertial periods, with growth rates increasing and wavenumbers decreasing as the boundary layer deepens in time. Growth of the baroclinic mode becomes faster than the Ekman adjustment time scale after only *t* ≈ 0.5*τ*^{Ek}, emphasizing that these instabilities grow rapidly relative to the shutdown process. High-frequency variability is associated with weak inertial oscillations modifying the boundary layer structure. The structure of the fastest-growing mode at *t* = 10 inertial periods is shown in Fig. 13. Vertical buoyancy production of kinetic energy dominates the instability growth, with dissipation reducing the kinetic energy tendency by a factor of approximately 1/2. Shear production is a small contribution to the vertically integrated kinetic energy tendency and is dominated by the lateral shear production, although there are locally significant contributions from both vertical and lateral shear production terms at the top of the boundary layer.

It is important to emphasize that, for downwelling-favorable interior flow, the BBL can become unstable to symmetric or centrifugal instability, with growth rates that exceed the baroclinic mode (Allen and Newberger 1998). Future work will therefore consider the full nonlinear Ekman adjustment process; however, the findings presented here suggest that fast-growing submesoscale baroclinic instability likely plays an important role in the 3D adjustment of the BBL to interior flows, emerging after symmetric/centrifugal instability restratifies the boundary layer to a state of marginal stability (discussed further in section 6).

### b. Upwelling favorable interior flow

When the interior flow is in the opposite direction, = −0.1 m s^{−1}, the Ekman flow is upslope, which tends to restratify the BBL (Fig. 14). Initially in the simulation, shear production of turbulent kinetic energy leads to the development of a thin turbulent BBL, with a maximum depth of about 15 m after three inertial periods. A weak, low-wavenumber baroclinic mode is evident in the stability analysis (Fig. 15), with a peak growth rate of approximately 0.03*f* at *t* = 7.5 inertial periods. The growth rate, however, decreases in time as the upslope Ekman flow restratifies the thin BBL. Thus, aside from a slight initial transient destabilization of the BBL, upwelling-favorable interior flow leads to a BBL that is stable to the baroclinic mode. The transient behavior in this simple numerical integration does, however, suggest that more complex, time-dependent interior flows may lead to baroclinically unstable BBLs, even in the case of time-mean flows that are upwelling favorable (Brink 2016; Brink and Seo 2016).

## 5. Nonlinear simulation

In this section we present the results of an example of idealized nonlinear simulation. Solutions are found using Dedalus, in a computational domain that is doubly periodic in the rotated horizontal coordinate system. This is achieved by imposing a fixed across-slope buoyancy gradient, (where is an assigned interior stratification), and solving for horizontally periodic departures from the background fields, similar to the “frontal zone” configuration often used in studies of the surface boundary layer (e.g., Taylor and Ferrari 2010). Full details of the model configuration are given in appendix B.

The model is initialized with a uniformly stratified interior (*N*_{I} = 3.4 × 10^{−3} s^{−1}) and a barotropic along-slope interior flow of 0.1 m s^{−1}. The BBL is initialized in thermal wind balance, with a thickness of 100 m, such that Ri = 1.5 in the BBL (Fig. 16). The slope angle is *θ* = 10^{−2}, giving a slope parameter in the BBL of *α* ≈ 0.15. For simplicity we consider the case of weak slope-normal viscosity and diffusivity, *ν* = *κ* = 10^{−5} m^{2} s^{−1}, and, for numerical stability, a biharmonic diffusivity in the rotated horizontal (*ν*_{h} = 2 × 10^{5} m^{4} s^{−1}) ( appendix B). A more complete exploration of parameter space using nonlinear simulations, including the use of physically realistic turbulence closures, will be the subject of future work.

Snapshots of the instability evolution are shown in Fig. 17. By day 10 the baroclinic mode has reached finite amplitude and is clearly evident in both the buoyancy and slope-normal velocity fields. The along-slope wavelength of the dominant instability, *λ* ≈ 6.4 km, is similar to that predicted by linear stability analysis of the domain-averaged buoyancy and velocity fields at day 1.5, *λ* ≈ 5.5 km. Importantly, the slope-normal velocities exceed 250 m day^{−1} (Fig. 17, bottom row), comparable to values associated with baroclinic instability in the surface boundary layer. These large-magnitude vertical velocities suggest that eddy fluxes in the BBL have the potential to affect a wide range of BBL processes, including, for example, nutrient fluxes between the BBL and interior, sediment transport, and coastal hypoxia.

The energetics of the nonlinear simulation are shown in Fig. 18. In the top panel the domain-integrated eddy kinetic energy (EKE) is shown, defined as , where primes indicate departure from the rotated horizontal average. An estimate of EKE can be calculated from linear theory as , where *t*_{o} = 1.5 days, the hat notation indicates the wavenumber spectrum (calculated from the model output at *t* = *t*_{o}), and *ω*_{i} is determined by linear stability analysis of the domain-averaged profiles at *t* = *t*_{o}. This estimate thus accounts for the distribution of energy across multiple wavenumbers with varying growth rates, and agrees closely with the growth rate of EKE in the nonlinear simulation over the first 10 days. During this time the instability growth is dominated by the buoyancy production of EKE, as indicated in Fig. 18 (bottom), typical of baroclinic instability. Later, after day 10, the available potential energy of the initial condition has been depleted, EKE is saturated, and buoyancy and shear production of EKE are largely compensated by dissipation (Fig. 18, bottom).

In this simulation, domain-averaged vertical buoyancy fluxes in the BBL reach ~ *O*(10^{−8}) m^{2} s^{−3}. This value is comparable to baroclinic instability in the surface boundary layer (Boccaletti et al. 2007), and importantly, is several orders of magnitude larger than values typically associated with the global breaking of internal waves over rough topography, where it is often assumed ~ *O*(10^{−10}) m^{2} s^{−3} (Nikurashin and Ferrari 2013; Ferrari et al. 2016). This simulation thus suggests that finite-amplitude baroclinic instability in the BBL has the potential to both modify the dynamics of the BBL, and to affect the interior dynamics of both coastal and deep oceans, through large, restratifying, eddy vertical buoyancy fluxes.

## 6. Instability regimes

In addition to the baroclinic mode, (6)–(10) also admit symmetric and centrifugal instabilities (Haine and Marshall 1998), each of which may be active in the ocean BBL under certain conditions (Allen and Newberger 1998; Brink 2012; Molemaker et al. 2015; Gula et al. 2016). While the focus of this work is primarily on baroclinic instability, it is useful to also briefly summarize where in parameter space each of these modes is anticipated to be dominant.

An illustrative example can be formed by assuming a geostrophically balanced along-slope flow in the direction of topographic wave propagation, with uniform stable stratification, and uniform slope-normal shear, such that (Fig. 1b). Note that the assumed linear background shear profile excludes the Kelvin–Helmholtz mode, otherwise expected at Ri < 0.25 (Vanneste 1993). The fastest-growing instability type will therefore depend on the sign of *f* times the Ertel PV of the background state (Haine and Marshall 1998), , which for this simplified configuration can be written nondimensionally as

For *q*′ < 0, the fastest-growing instabilities will be associated with overturning circulations in the slope-normal and across-slope plane, *l* = 0, which grow by extracting kinetic energy from the background flow (Hoskins 1974). These instabilities can be of the symmetric, centrifugal, or mixed types, with the ratio of the lateral shear production (LSP) to vertical shear production (VSP) providing a useful discriminator ( appendix A; Thomas et al. 2013),

In the case that LSP/VSP < 1 the instabilities are of the symmetric type, for LSP/VSP > 1 and *α*/Ri > 1 (negative absolute vertical vorticity) the instabilities are of the centrifugal type, and for LSP/VSP > 1 and *α*/Ri < 1 (positive absolute vertical vorticity) the instabilities are a mixed SI/CI instability (Thomas et al. 2013). For *q*′ > 0 the baroclinic mode *k* = 0 will be the fastest-growing instability, with perturbation kinetic energy increasing through the extraction of the background state available potential energy (buoyancy production).

An example regime diagram is shown in Fig. 19, as a function of the slope Burger number and the Richardson number. The *q*′ = 0 line divides the domain (dashed red line), with the entire area to the right of the line linearly unstable to baroclinic instability, although instabilities at small *α* are likely the most physically relevant (section 2c). For *q*′ < 0 the domain is further divided by the energetics, with LSP/VSP = 1 shown from theory [(33); dashed blue line in Fig. 19], assuming *δ*/Ri^{1/2} = *f*/*N* = 0.1. Below this line the symmetric mode grows by extracting mean kinetic energy through the vertical shear production. Above this line the energetics are dominated by the lateral shear production, either in a pure centrifugal instability mode, to the left of the green line indicating zero absolute vertical vorticity (1 − *α*/Ri = 0), or in a mixed symmetric/centrifugal type mode, for positive absolute vertical vorticity to the right of the green line. Numerical solutions of (6)–(10) [with *l* = 0, *k* = 20*f*/(*NH*)] were also used to calculate (solid blue line), demonstrating that the theory provides an accurate approximation to the energetics.

For small Ri, the dominant instability types in the BBL are therefore likely to be of the centrifugal, symmetric, or mixed types, with the centrifugal mode becoming dominant for large *S*, where the topography begins to act as an effectively vertical boundary (Dewar et al. 2015). For and small *S*, the symmetric mode is likely most relevant (Allen and Newberger 1998). In the case that , the baroclinic mode will dominate. Further, when considering the nonlinear dynamics of the BBL, it is expected that the time-dependent adjustment of the background fields by finite-amplitude instabilities will lead to a traversal of parameter space as the shear and stratification are adjusted by eddy fluxes. Thus, even for initial background conditions that have fastest-growing modes of the overturning type, the baroclinic mode may emerge later as the instabilities restratify the boundary layer. For example, during the Ekman adjustment to a downwelling-favorable interior flow, potential vorticity is removed at the boundary (Benthuysen and Thomas 2012), and in the BBL the potential vorticity can become negative. The fastest-growing modes of the 1D solutions of section 4a (Fig. 11) will therefore likely be of the symmetric or centrifugal type. However, when these instabilities reach finite amplitude they will adjust the boundary layer stratification such that Ri ≈ 1, at which point the baroclinic mode will dominate (Haine and Marshall 1998; Brink and Cherian 2013). Likewise, even in the case where the boundary layer is forced to remain near *q* = 0, as, for instance, might occur in the case of continuing downslope Ekman transport, growing baroclinic instability can still be present (Callies and Ferrari 2018).

## 7. Summary

In this article, we considered baroclinic instability over sloping topography, extending earlier QG results to the nongeostrophic regime appropriate for the ocean BBL. Importantly, weakly stratified BBLs can support a submesoscale baroclinic instability, a BBL counterpart to the surface baroclinic mixed layer instability (Boccaletti et al. 2007). In the BBL these submesoscale instabilities are relatively insensitive to topographic slope, suggesting they are likely a robust part of the dynamics of the ocean BBL, along with the symmetric (Allen and Newberger 1998) and centrifugal (Molemaker et al. 2015) modes.

Two 1D theories of the BBL structure over topography were also shown to result in solutions that are susceptible to growing baroclinic instability, with important implications for our understanding of the dynamics of both the BBL and the interior. In the case of a BBL generated by turbulent mixing, a thick outer layer supports a rapidly growing baroclinic mode that releases available potential energy and provides a restratifying tendency that is not accounted for in the 1D formulation. This suggests that deep submesoscale baroclinic instabilities could alter the stratification and across-slope flow along topographic features such as the mid-Atlantic ridge, with important implications for closure of the abyssal circulation (Callies 2018). Likewise, interior flow along a slope generates available potential energy in the BBL through an across-slope Ekman flow, leading to baroclinic growth rates that quickly exceed the rate of Ekman adjustment. These instabilities thus have the potential to affect both the BBL dynamics and the exchange between the boundary layer and interior, in both the coastal and deep ocean.

Extensive work on the surface boundary layer has demonstrated the leading-order role that submesoscale baroclinic instability plays in boundary layer restratification (Fox-Kemper et al. 2008; Johnson et al. 2016; Su et al. 2018), vertical exchange between the boundary layer and interior (Mahadevan and Tandon 2006), and the surface flux of potential vorticity (Wenegrat et al. 2018). The similarities between the BBL and surface mode (section 2), and the large-amplitude eddy fluxes found in the nonlinear simulation of section 5, suggest that many of the same results may hold in the BBL. However, there are also important differences between the surface and bottom boundary layer that require further consideration. For example, as noted in section 2, topography shapes parcel trajectories, modifying the instability energetics and potentially the associated eddy fluxes of buoyancy and tracers (Spall 2004; Isachsen 2011; Brink 2013). Further, variations in topography, not considered here, can impose an external length scale on instabilities (de Szoeke 1983), generate slope-normal flow and secondary circulations (Benthuysen et al. 2015), and modify instability growth rates (de Szoeke 1983; Solodoch et al. 2016), all of which may alter the role of baroclinic instability in BBL dynamics. Relaxing the assumption of across-slope uniformity in boundary layer height or structure may also lead to an arrest of the inverse cascade due to the topographic beta effect (Rhines 1975; Brink 2012). Full nonlinear simulations to explore the finite-amplitude behavior, and parameter dependence, of the submesoscale BBL baroclinic mode will therefore be the subject of a future paper.

## Acknowledgments

The authors thank the contributors to the Dedalus Project (www.dedalus-project.org) and GOTM (www.gotm.net), both of which were extremely useful for the analyses in this paper. Authors JOW and LNT were supported during this work by NSF Grant OCE-1459677. Comments from three anonymous reviewers are also gratefully acknowledged.

### APPENDIX A

#### Energetics

In this appendix we derive several results related to the instability energetics. For simplicity we assume hydrostatic dynamics, with background fields that are uniform in the along-slope direction, and diffusion only in the slope-normal direction. The eddy potential energy equation can be formed by taking *b*/*N*^{2} × (5), giving

where

is the conversion of mean potential energy to eddy potential energy through the horizontal buoyancy flux,

is the conversion of eddy potential energy to EKE through the vertical buoyancy production, and

is the dissipation of eddy potential energy. In this appendix, angle brackets indicate spatial averaging, both vertically and horizontally over a wavelength.

where

and

give the conversion of mean kinetic energy to EKE through vertical and lateral shear production, respectively. Dissipation of EKE is given by

For the baroclinic mode the EKE source is the vertical buoyancy production (VBP), which consists of two terms related to the slope-normal and across-slope advection of buoyancy. The ratio of these terms scale as

Noting that PE_{t} > 0 requires 〈*ub*〉 < 0, the across-slope advection of buoyancy can be seen to reduce (increase) VBP for *α* > 0 (*α* < 0).

Finally, for the purposes of discriminating between centrifugal- and symmetric-type instabilities (section 6), it is useful to find the ratio of the lateral and vertical shear production terms. Following Thomas et al. (2013), and assuming a geostrophically balanced along-slope background flow, this ratio can be written as

where

Using , this can be written as

and

### APPENDIX B

#### Configuration of the Nonlinear Model

The nonlinear simulation of section 5 is performed using Dedalus (Burns et al. 2016). The domain is 32 km in each rotated horizontal direction and 1000 m in the slope-normal direction and is decomposed using 64 Fourier components in each horizontal direction and 256 Chebyshev modes in the vertical, giving horizontal and vertical resolutions of Δ*x* = Δ*y* = 500 m, and Δ*z* = 0.075–6 m. The rotated equations are solved in a horizontally periodic domain by assuming a fixed mean across-slope gradient in buoyancy, such that , where is an imposed value of the interior stratification. This configuration is thus similar to the frontal-zone setup commonly used for studying surface mixed layer fronts in periodic domains (e.g., Taylor and Ferrari 2010). The initial condition is assumed to be an along-slope flow in hydrostatic and geostrophic balance that varies only in the slope-normal direction *V*(*z*). The rotated equations of motion can then be written in the hydrostatic limit, that is, and , as

In the above, uppercase quantities indicate the imposed initial conditions, and time-dependent, horizontally periodic quantities are denoted with the tilde notation. Specifically, the total velocity is decomposed as , the total pressure divided by the reference density as , and the buoyancy as . Horizontal mixing is parameterized using a biharmonic diffusion operator for numerical stability, . Boundary conditions are given by

## REFERENCES

*Geophysical Fluid Dynamics*. Springer, 710 pp.

*Ocean Modeling in an Eddying Regime, Geophys. Monogr.*, Vol. 177, Amer. Geophys. Union, 17–38, https://doi.org/10.1029/177GM04.

## Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

Throughout this paper, we will use the term bottom boundary layer to denote any weakly stratified lower layer, a definition that includes both very well mixed layers adjacent to topography and the weakly stratified overlying layer that can appear in some solutions to the 1D dynamics (e.g., section 3).

^{2}

Note that we choose to follow the notation of Stone (1966), and hence we denote the slope parameter using *α* and the nonhydrostatic parameter using *δ* (cf. Hetland 2017). As in Hetland (2017) our definition of the slope parameter is of the opposite sign to that used in Blumsack and Gierasch (1972) and other preceding work, such that *α* > 0 implies isopycnal and topographic slopes of opposite sign (see, e.g., Fig. 5).