2.1. Inviscid, linear flow

The inviscid, linearized forms of the governing equations (2.2)–(2.6) describe the heaving of isopycnals up and down the slope:

(2.14)\begin{gather} \partial_t\tilde{{u}}=\frac{1}{{Ro}}\tilde{v}+{C}^{2}\tilde{b}+F(t), \end{gather}

(2.15)\begin{gather}\partial_t\tilde{{v}}={-}\frac{1}{{Ro}}\tilde{u}, \end{gather}

(2.16)\begin{gather}\partial_t\tilde{{b}}={-}\tilde{{u}}. \end{gather}

Crucially, the solutions to (2.14)–(2.16) prescribe the amplitude of the non-dimensional body force:

(2.17)\begin{equation} A=\left({C}^{2}+\frac{1}{{Ro}^{2}}-1\right),\end{equation}

where $A_{d}=AU_\infty \omega$. The solutions to (2.16)–(2.16) are

(2.18)\begin{gather} \tilde{{u}}(t)={-}\cos(t), \end{gather}

(2.19)\begin{gather}\tilde{{v}}(t)={Ro}^{{-}1}\sin(t), \end{gather}

(2.20)\begin{gather}\tilde{{b}}(t)=\sin(t). \end{gather}

2.2. Resonance

The assumption that the vertical and horizontal structure of an internal wave is constant over a length scale of ${O}(10)$ m is not valid if the wave is critical. At critical slope, the forcing model $F(t)$ in (2.1) and (2.14) is degenerate and vanishes if the critical slope condition

(2.21)\begin{equation} {C}^{2}+\frac{1}{{Ro}^{2}}-1=0 \end{equation}

is satisfied. If (2.21) is satisfied, the energy of the inviscid baroclinic tide is tightly focused into narrow beams that follow the curvature of the bathymetry (Balmforth, Ierley & Young Reference Balmforth, Ierley and Young2002).

Equation (2.21) is formally consistent with internal wave theory. In internal wave theory, critical slopes are defined by

(2.22)\begin{equation} \tan\theta_c=\sqrt{\frac{\omega^{2}-f^{2}}{N^{2}-\omega^{2}}}, \end{equation}

where $\theta _c$ is the critical slope angle. If the slope angle $\theta \neq \theta _c$, then (2.17) is satisfied, $A\neq 0$. Equation (2.17) can be rearranged to obtain

(2.23)\begin{equation} \tan\theta=\sqrt{\frac{\omega^{2}(1+A)-f^{2}}{N^{2}-\omega^{2}(1+A)}}. \end{equation}

Therefore the criticality condition in (2.21) is just a rearrangement of the criticality condition defined by the slope parameter $\epsilon$ from internal wave theory:

(2.24)\begin{equation} \epsilon=\frac{\tan\theta}{\tan\theta_c}, \end{equation}

where criticality states are defined:

(2.25)\begin{equation} \left.\begin{array}{llllll@{}} \text{if} & A<0 & \text{then} & \epsilon<1 & \rightarrow & \theta\ \text{is subcritical},\\ \text{if} & A=0 & \text{then} & \epsilon=1 & \rightarrow & \theta\ \text{is critical},\\ \text{if} & A>0 & \text{then} & \epsilon>1 & \rightarrow & \theta\ \text{is supercritical}. \end{array}\right\} \end{equation}

2.3. Boundary conditions

At the solid boundary at $z=0$, the boundary conditions on the total velocity are no-slip and impermeability

(2.26)\begin{equation} \tilde{\boldsymbol{u}}=0, \end{equation}

and the boundary conditions on the total buoyancy is the adiabatic condition:

(2.27)\begin{equation} \partial_z\tilde{b}=0. \end{equation}

As $z\rightarrow \infty$, the velocity boundary conditions are the oscillatory solutions for the inviscid flow, and zero flow in the wall-normal direction: (2.18), (2.19) and $\widetilde {{w}}=0$. The non-dimensional buoyancy field as $z\rightarrow \infty$ has two components, the inviscid oscillation and the constant background stratification:

(2.28)\begin{equation} \tilde{{b}}=x+z\cot\theta+\sin(t). \end{equation}

2.4. Variable decomposition

The total velocity and buoyancy fields are decomposed into three components that when summed together satisfy (2.2)–(2.6) and (2.26)–(2.28). To distinguish the components, let ‘$H$’ denote the hydrostatic (and possibly geostrophic) component, let ‘$S$’ denote the steady component and let ‘$O$’ denote the oscillating component:

(2.29)\begin{gather} \tilde{\boldsymbol{u}}(x,y,z,t)= \boldsymbol{u}_{{H}}+\boldsymbol{u}_{{S}}(z)+\boldsymbol{u}_{{O}}(x,y,z,t), \end{gather}

(2.30)\begin{gather}\tilde{b}(x,y,z,t)=b_{{H}}(x,z)+b_{{S}}(z)+b_{{O}}(x,y,z,t), \end{gather}

(2.31)\begin{gather}\tilde{p}(x,y,z,t)=p_{{H}}(x,z)+p_{{S}}(z)+p_{{O}}(x,y,z,t). \end{gather}

The hydrostatic component of the buoyancy field is merely the background stratification in the rotated coordinate system, $b_{{H,d}}(x_{d},z_{d}) =N^{2}(x_{d}\sin \theta +z_{d}\cos \theta )$ in dimensional form, and the quiescent hydrostatic velocity field $\boldsymbol {u}_{{H}}=0$ everywhere except in the finite-Rossby-number flow regime, in which case it is the along-slope geostrophic velocity that arises from the across-isobath pressure gradient imposed by the slope (Phillips Reference Phillips1970; Wunsch Reference Wunsch1970). The buoyancy frequency $N$ is defined in the same manner as convention:

(2.32)\begin{equation} N^{2}={-}\frac{g}{\rho_0}\partial_\eta{\rho}_{{H}}=\partial_\eta{b}_{{H,d}}, \end{equation}

where $\eta$ denotes the vertical position coordinate as shown in figure 1.

The steady and oscillating flow components are anomalies from the hydrostatic background that ensure the satisfaction of frictional and diffusive boundary conditions at the wall and inviscid oscillations far from the wall. It is convenient to solve for the anomalies,

(2.33)\begin{gather} \boldsymbol{u}(x,y,z,t)=\boldsymbol{u}_{{S}}(z)+\boldsymbol{u}_{{O}}(x,y,z,t), \end{gather}

(2.34)\begin{gather}b(x,y,z,t)=b_{{S}}(z)+b_{{O}}(x,y,z,t), \end{gather}

(2.35)\begin{gather}p(x,y,z,t)=p_{{S}}(z)+p_{{O}}(x,y,z,t), \end{gather}

because the removal of the hydrostatic background permits periodic analytical and numerical solutions for $\boldsymbol {u}$ and $b$.

3.1. Necessary conditions for gravitational instability

The linear solutions for the steady flow component of both rotating and non-rotating flows is always gravitationally stable, meaning that the total vertical buoyancy gradient is never negative,

(3.1)\begin{equation} \partial_{\eta}b_{H}+\partial_{\eta}b_{S} =1-\cos^{2}\theta\mathrm{e}^{-\eta/\delta_{S}} \sqrt{2}\sin\left(\frac{\eta}{\delta_{S}}+\frac{\rm \pi}{4}\right)\geqslant 0, \end{equation}

if the oscillating component vanishes. The boundary layer thickness of the steady component (Phillips Reference Phillips1970; Wunsch Reference Wunsch1970) is

(3.2)\begin{equation} \delta_{S}=\left(\frac{f^{2}\cos^{2}\theta}{4\nu^{2}} +\textit{Pr} \frac{N^{2}\sin^{2}\theta}{4\nu^{2}}\right)^{{-}1/4}. \end{equation}

The oscillating flow component solutions exhibit transient gravitationally unstable buoyancy gradients, $\delta _\eta \tilde {b}<0$, when denser fluid is advected over lighter fluid during portions of the oscillation period. If the oscillating component is non-zero, then the minimum necessary condition for gravitational instabilities is defined by

(3.3)\begin{equation} \partial_{\eta}b_{H}+\partial_{\eta}b_{S}<{-}\partial_{\eta}b_{O}, \end{equation}

because if (3.3) is satisfied, then the total vertical buoyancy gradient is negative, $\delta _\eta \tilde {b}<0$. However, instabilities can grow only if the negative buoyancy gradient is sustained for a significant amount of time and if it is negative enough to overcome resistance from friction.

A characteristic boundary layer Rayleigh number and a ratio of the time scale of the growth of an instability to the period of the oscillation are required to estimate the minimum (quasi-steady) conditions for the growth of gravitational instabilities. However, the linear solutions do not readily yield a single boundary layer buoyancy gradient length scale. If the buoyancy gradient length scale is assumed to scale with $\delta =\sqrt {2\nu /\omega }$, then a tenable time-dependent boundary layer Rayleigh number is defined:

(3.4)\begin{equation} {Ra}(t)\sim\frac{4\textit{Pr} N^{2}}{\omega^{2}}\partial_\eta\tilde{b}(t), \end{equation}

which only applies when $\partial _\eta \tilde {b}(t)<0$. To estimate the gravitational stability of the flow without explicitly accounting for the time dependence of the basic state (the quasi-steady assumption), the basic state of the flow cannot change more rapidly than the growth rate of a gravitational instability. If the instabilities are ‘slowly modulated’ by the basic state (Davis Reference Davis1976), the quasi-steady assumption is reasonable for stability analysis. The dimensional instantaneous growth rate of a gravitational instability can be estimated as

(3.5)\begin{equation} \sigma\sim{\rm Im}\left[\sqrt{\partial_\eta \tilde{b}N^{2}}\right]. \end{equation}

If $|{\omega }/{\sigma }|\ll 1$, then the modulation by the basic state is sufficiently slow for the growth of gravitational instabilities.

Figure 2 shows the minimum values of the non-dimensional total vertical buoyancy gradient from the linear solutions for non-rotating and rotating cases as a function of slope parameter and Stokes Reynolds number. Assuming that the gravitational instability in the boundary layer is physically similar to that of Rayleigh–Bénard instability in the case of one rigid and one stress-free boundary, then the critical Rayleigh number for the boundary layer is ${Ra}_c\approx 1100$ (Chandrasekhar Reference Chandrasekhar1961), which corresponds to $|\omega /\sigma |=0.06$. For the chosen fluid properties (holding $\textit {Pr}=1$, $f=10^{-4}$, $N=10^{-3}$ and ${\omega =1.4\times 10^{-4}}$ constant), ${Ra}_c\approx 1100$ corresponds to a critical non-dimensional vertical buoyancy gradient of $\partial _\eta \tilde {b}=-5.4$, which is shown as the blue lines in figure 2. The minimum boundary layer buoyancy gradient is less than zero for all non-zero $\textit {Re}$ and $\epsilon$, and the minimum boundary layer buoyancy gradient is increasingly negative with increasing Reynolds number and with increasing slope parameter. The discontinuity at $\epsilon =1$ is an artefact of the degeneracy of linear solutions at critical slope. The $\epsilon$ axis between the non-rotating and rotating cases is different because rotation alters the angle of critical slope; both plots show the same slope angle range, $0<\theta \leqslant 16^{\circ }$.

Figure 2.

Total buoyancy gradient minima. The minimum value (in both time and space) of the non-dimensional linear solution vertical buoyancy gradient for (a) the non-rotating reference frame case ($f=0$) and (b) the rotating reference frame case.

4.1. Numerical implementation

The flow anomalies, as defined by (2.33) and (2.34), are discretized to satisfy periodic boundary conditions in the wall-parallel directions via Fourier spectral bases in the across-isobath ($x$) and along-isobath ($y$) directions. Periodicity is not merely numerically convenient; it also eliminates the need to prescribe buoyancy forcing (‘restratification’) because the oscillating flow can advect the background field to gain or lose buoyancy. In the periodic domain, the boundary layer buoyancy can only reach a homogenized steady state if the turbulence sustainably converts tidal momentum to potential energy throughout the entire period.

Although the planar extent of the computational domain is less than the excursion length of the tide, the domain size (table 2) is justifiably sufficient because the largest eddies in oscillating boundary layers are those associated with the transverse (wall-normal) length scale, which is much less than the excursion length. Indeed, at higher Reynolds number ($\textit {Re}=1790$), Gayen & Sarkar (Reference Gayen and Sarkar2010b) found the turbulent boundary layer thickness, $\delta _l$, was $\delta _l=15\delta$ for the unstratified problem and $\delta _l=17\delta$ for flat-plate stratified oscillating boundary layers at the same Reynolds number. The grid resolution parameters for the two Reynolds numbers examined are shown in table 2, where ($L_x,L_y,H$) are the domain dimensions in ($x,y,z$), $l_{K}$ and $\tau _{K}$ are the Kolmogorov length and time scales, respectively, and wall units (denoted by superscript $+$) are scaled by the viscous length scale $\delta _v=\nu /U_*$, where $U_*$ is the a priori estimate of the friction velocity, which is approximated as $U_*=\sqrt {\nu \partial _z\bar {u}}\sim \sqrt {\nu {U_\infty }/{\delta }}$.

Table 2.

Resolution parameters. The grid is identical in the $x$ and $y$ directions. Kolmogorov scales are estimated assuming that the law of the wall holds; therefore, characteristic dissipation rate is estimated a priori by $\varepsilon \sim U_*^{3}/(\delta \kappa _*)$, where the von Kármán constant is $\kappa _*=0.41$ and the friction velocity is estimated by $U_*\sim \sqrt {\nu U_\infty /\delta }$.

Our decomposition of the flow is the same as that of Phillips (Reference Phillips1970) and Wunsch (Reference Wunsch1970): we decompose the flow into a quiescent, hydrostatic background flow (denoted by the subscript ‘$H$’ variables in (2.29), (2.30) and (2.31), where $\boldsymbol {u}_{H}=0$) and solve for anomalies from the background flow (the left-hand sides of (2.33), (2.34) and (2.35)) that, when summed with the background flow component, satisfy the no-slip (2.26) and adiabatic (2.27) boundary conditions at the wall imposed on the total flow. The nonlinear governing equations for the anomalies are

(4.1)\begin{gather} \partial_x{u}+\partial_y{v}+\partial_z{w}=0, \end{gather}

(4.2)\begin{gather}{d}_t{{u}}=\frac{1}{{Ro}}{v}-\partial_x{p}+ \frac{1}{\textit{Re}_L}\left(\partial_{xx}+\partial_{yy}+ \partial_{zz}\right){u}+{C}^{2}{b}+F(t), \end{gather}

(4.3)\begin{gather}{d}_t{{v}}={-}\frac{1}{{Ro}}{u}-\partial_y{p}+ \frac{1}{\textit{Re}_L}\left(\partial_{xx}+\partial_{yy}+\partial_{zz}\right){v}, \end{gather}

(4.4)\begin{gather}{d}_t{{w}}={-}\partial_z{p}+ \frac{1}{\textit{Re}_L}\left(\partial_{xx}+\partial_{yy}+\partial_{zz}\right){w} +{C}^{2}{b}\cot\theta, \end{gather}

(4.5)\begin{gather}{d}_t{{b}}=\frac{1}{\textit{Pr} \,\textit{Re}_L} \left(\partial_{xx}+\partial_{yy}+\partial_{zz}\right){b}, \end{gather}

where the anomalies ${u},v,w,p$ and $b$ are defined by (2.33)–(2.35). The material derivative on the left-hand side of (4.5) contains nonlinear buoyancy anomaly advection terms as well as the terms describing the time rate of change of anomalous buoyancy, $-(u+w\cot \theta )$, by the advection of background buoyancy. Since the linear solutions for the oscillating laminar flow (Thorpe Reference Thorpe1987; Baidulov Reference Baidulov2010) contain only the advective terms corresponding to the background buoyancy advection, the laminar stratification oscillates from negative to positive and back again in a linear advection–diffusion balance. However, the boundary layer stratification can weaken, relative to the oscillating laminar stratification, if the vertical gradient of the nonlinear anomalous buoyancy advection in (4.5) is negative and dominates the vertical gradient of background buoyancy advection. Therefore, gradients of nonlinear anomalous buoyancy advection must continually overcome the restratifying tendency of background buoyancy advection (e.g. $u>0$ corresponds to a local decrease in buoyancy as relatively heavy fluid is advected upslope) for weakened boundary layer stratification to persist.

The anomalies were computed using the MPI-parallel pseudo-spectral partial differential equation solver Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2019) using $128^{3}$ modes. A third-order, four-stage, implicit–explicit Runge–Kutta method derived by Ascher, Ruuth & Spiteri (Reference Ascher, Ruuth and Spiteri1997) was used for temporal integration. Chebyshev polynomial bases of the first kind were employed for spatial discretization on a cosine grid in the wall-normal direction. Chebyshev polynomials permit the exact enforcement of the adiabatic wall-boundary condition ((2.27) minus the background component) on the buoyancy field and no-slip/impermeability wall-boundary conditions on the velocities ((2.26) minus the background component). The $3/2$ rule dealiasing scheme is used not only for dealiasing the spatial modes online but also for dealiasing post-processed flow statistics.

At the maximum wall-normal extent of the domain, the boundary conditions at infinity (2.18), (2.19) and (2.28) were approximated for the anomalies as free-slip, impermeable conditions:

(4.6)\begin{equation} \partial_zu=\partial_zv=w=0, \end{equation}

and an adiabatic condition on just the anomaly:

(4.7)\begin{equation} \partial_zb=0, \end{equation}

such that the total flow buoyancy gradient at $z=H$ is the background buoyancy gradient in that direction.

Although the impermeability condition causes the reflection of internal waves that reach the upper boundary, the effects are assumed to be negligible because of the negligible amount of energy propagated by such high-wavenumber waves in low-Reynolds-number flow. Gayen & Sarkar (Reference Gayen and Sarkar2010a) found that for flat-bottomed stratified oscillatory flow at larger-Reynolds-number flow ($\textit {Re}=1790$), the vertical wave energy flux is less than 1 % of the boundary layer dissipation and production rates. Indeed, small but non-zero dissipation rates of turbulent kinetic energy (TKE) were found near the upper boundary in some of the simulations, presumably from subharmonic parametric instability or other wave–wave instabilities because of the free-slip reflective upper boundary condition. However, 99.9 % of the shear production rate and dissipation rate occurred within one Ozmidov length of the wall ($L_O=\sqrt {\varepsilon /N^{3}}$, using the magnitude of the background stratification for $N$) at the lower boundary for all simulations.

A small amount of grid-scale noise was observed primarily in the $x,y$ directions for $\textit {Re}=840$ cases. However, the wall-normal integrated TKE balances (figure 4) the residual of the right-hand-side TKE equation terms and the calculation of $\partial _tK$ matched to graphical accuracy, suggesting the grid-scale noise did not significantly alter energetics. The resolution of $\Delta x/l_K=\Delta y/l_K=13.8,18.0$ (table 2) was chosen by estimating the characteristic scale of $\varepsilon ={O}(10^{-10}) \ \textrm {m}^{2}\ \textrm {s}^{-3}$ a priori assuming fully developed steady turbulence that satisfies the law of the wall, even though law-of-the-wall turbulence was not anticipated because of the transitional Reynolds numbers prescribed. The a priori estimate $\varepsilon ={O}(10^{-10})\ \textrm {m}^{2}\ \textrm {s}^{-3}$ roughly agrees with the majority of the corresponding a posteriori average dissipation rates shown in table 4, despite the false assumptions of the a priori estimate. Since the TKE equation closes to graphical accuracy and the instabilities investigated are laminar in origin, we anticipate that the results of this study will remain the same at higher resolution. However, the small amount of signal attenuation at the Nyquist wavenumbers suggests that the shear production and dissipation rates may be slightly underestimated.

The initial conditions were specified as the sum of the steady component (Phillips Reference Phillips1970; Wunsch Reference Wunsch1970), the oscillating component (Thorpe Reference Thorpe1987) at time $t=0$ and uniformly distributed white noise corresponding to buoyancy anomaly perturbations of magnitude $10^{-10}\ \textrm {m}\ \textrm {s}^{-2}$. All of the simulations that exhibited turbulence (defined as wall-normal integrated production rates of TKE greater than $10^{-10}\ \textrm {m}^{3}\ \textrm {s}^{-3}$) did so within two oscillations.

The parameter regimes shown in table 1 qualitatively describe flows forced by the $M_2$ tide frequency, which is specified as $\omega =1.4\times 10^{-4}$ ($\textrm {rad}\ \textrm {s}^{-1}$, for a 12.4 h tide period) and the Coriolis parameter is specified as $f=10^{-4}$ ($\textrm {rad}\ \textrm {s}^{-1}$). Much of the abyssal ocean is filled with Antarctic Bottom Water (AABW), characterized by temperatures near $0\,^{\circ }\textrm {C}$ and practical salinities of approximately 35 psu. At $0\,^{\circ }\textrm {C}$ and 35 psu, the kinematic viscosity is $1.8\times 10^{-6}\ \textrm {m}^{2} \ \textrm {s}^{-1}$ (Chen et al. Reference Chen, Chan, Read and Bromley1973; Talley Reference Talley2011) and the thermal diffusivity is $1.4\times 10^{-7}\ \textrm {m}^{2}\ \textrm {s}^{-1}$ (Thorpe Reference Thorpe2007); thus $Pr\approx 13$ for AABW. We specify the kinematic viscosity as $\nu =2\times 10^{-6}\ \textrm {m}^{2}\ \textrm {s}^{-1}$ in the relevant parameters in table 1 and $Pr=1$ for simplicity. We approximate the background buoyancy frequency at mid-latitude abyssal depths as $N=10^{-3}\ \textrm {rad}\ \textrm {s}^{-1}$ (Thurnherr & Speer Reference Thurnherr and Speer2003). Baroclinic tide amplitudes of $U_{\infty }\leqslant 0.01\ \textrm {m}\ \textrm {s}^{-1}$ in tidal models (Turnewitsch et al. Reference Turnewitsch, Falahat, Nycander, Dale, Scott and Furnival2013) agree with observations at the northwestern edge of the Sargasso Sea at $34^{\circ }\textrm {N}\ 70^{\circ }\textrm {W}$ (Tarbell et al. Reference Tarbell, Montgomery and Briscoe1985; Nash et al. Reference Nash, Kunze, Toole and Schmitt2004). Global assessments of slope criticality suggest that the same region is characterized by subcritical, low-angle slopes. We also examine supercritical slopes for completeness. The prescribed non-dimensional parameters in table 1 were varied by only varying the velocity, $U_\infty =0.01,0.005\ \textrm {m}\ \textrm {s}^{-1}$ corresponding to the two Reynolds numbers, and the slope angle (shown in the rightmost column of table 1).

4.2. Intermittent turbulent bursts

The integrated TKE budget of each simulation was computed in order to distinguish the laminar and turbulent regimes and to quantify turbulence production mechanisms. The planar mean TKE is defined as

(4.8)\begin{equation} K(z,t)\equiv\tfrac{1}{2}\left(\overline{u'^{2}}+\overline{v'^{2}}+\overline{w'^{2}}\right), \end{equation}

where the planar mean operator and variable decomposition are defined:

(4.9)\begin{gather} \bar{\phi}(z,t)=\frac{1}{L_xL_y}\int_{{-}L_x/2}^{L_x/2}\int_{{-}L_y/2}^{L_y/2} \phi(y,z,t)\,{\rm d} y\,{\rm d}x, \end{gather}

(4.10)\begin{gather}\phi(x,y,z,t)=\bar{\phi}(z,t)+\phi'(x,y,z,t), \end{gather}

and $\phi$ is any of the anomalous variables defined by (2.33)–(2.35). The planar mean TKE evolution equation is

(4.11)\begin{equation} {\partial_t{K}} +\partial_{z}\mathcal{T} =\mathcal{P} +\mathcal{B} -\varepsilon. \end{equation}

The TKE transport term $\partial _{z}\mathcal {T}$ includes all TKE flux divergences (mean, turbulent, pressure, diffusion), which vanish upon wall-normal integration of (4.11). The rate production of TKE by mean shear is $\mathcal {P}$ (production in the sense that, generally, $\mathcal {P}>0$), and it is defined as

(4.12)\begin{gather} \mathcal{P}(z,t)={-}\overline{u'w'}\partial_z\bar{u}-\overline{v'w'}\partial_z\bar{v}, \end{gather}

(4.13)\begin{gather}\mathcal{P}_{13}={-}\overline{u'w'}\partial_z\bar{u}, \end{gather}

(4.14)\begin{gather}\mathcal{P}_{23}={-}\overline{v'w'}\partial_z\bar{v}. \end{gather}

The buoyancy flux $\mathcal {B}$ is typically downgradient ($\mathcal {B}<0$ amidst $\partial _z\bar {b}>0$ or $\mathcal {B}>0$ amidst $\partial _z\bar {b}<0$), in which case it represents the conversion of TKE into potential energy, and it is defined as

(4.15)\begin{gather} \mathcal{B}(z,t)=\overline{w_\eta'b'}=\overline{u'b'}\sin\theta+\overline{w'b'}\cos\theta, \end{gather}

(4.16)\begin{gather}\mathcal{B}_1=\overline{u'b'}\sin\theta, \end{gather}

(4.17)\begin{gather}\mathcal{B}_3=\overline{w'b'}\cos\theta \end{gather}

in the rotated reference frame (where $w_\eta ={d}_t\eta$ is the velocity in the vertical, not the wall-normal velocity $w$). Defined in this manner a downgradient buoyancy flux may be reversible. A reversible buoyancy flux may be thought of as a buoyancy flux that converts TKE into potential energy through stirring alone. Finally, the dissipation rate of TKE

(4.18)\begin{align} \varepsilon(z,t)&=\nu\left( \overline{(\partial_x{u'})^{2}} +\overline{(\partial_x{v'})^{2}} +\overline{(\partial_x{w'})^{2}} +\overline{(\partial_y{u'})^{2}} +\overline{(\partial_y{v'})^{2}}\right. \nonumber\\ &\quad \left.+\overline{(\partial_y{w'})^{2}} +\overline{(\partial_z{u'})^{2}} +\overline{(\partial_z{v'})^{2}} +\overline{(\partial_z{w'})^{2}}\right) \end{align}

is positive definite and therefore the last term of (4.11) is always a sink of TKE.

The laminar oscillating boundary layer buoyancy gradient is transient and evanescent; therefore, we seek an integral quantity to measure the stabilizing/destabilizing effects of the time-periodic laminar buoyancy gradient. We borrow the concept of boundary layer displacement thickness (Monin & Yaglom Reference Monin and Yaglom1971) and apply it to transient bulk buoyancy gradients rather than momentum. We refer to this measure as the boundary layer stratification thickness, $\delta _s$. If the total buoyancy field is not constant over small distance in the wall-normal direction $z_1$, then it can be approximated as constant over some distance $z_0$ from the wall, where $z_0$ is defined by

(4.19)\begin{equation} N^{2}z_0=\int_0^{z_1}\partial_z\tilde{b}\,\mathrm{d}z. \end{equation}

It follows that

(4.20)\begin{equation} \delta_s=z_1-z_0=\int_0^{z_1}\left(1-\frac{\partial_z\tilde{b}}{N^{2}}\right)\mathrm{d}z. \end{equation}

Since $\partial _z\tilde {b}\rightarrow N^{2}$ as $z_1\rightarrow \infty$ the integrand of (4.20) vanishes as $z_1\rightarrow \infty$.

Figure 3 shows the geometric interpretation of (4.19). Here $\delta _s>0$ and indicates that, in bulk, the oscillating laminar boundary layer stratification is less than the background stratification, and vice versa as the sign of the laminar boundary layer buoyancy gradient oscillates from negative to positive in the laminar flow solutions provided in Appendix A. The stratification thickness is positive when heavier fluid is advected over lighter fluid that has been impeded from flowing upslope by molecular friction at the boundary. Note that $\delta _s<0$ for the laminar steady flow component solutions because the analytical solutions (Phillips Reference Phillips1970; Wunsch Reference Wunsch1970) require steady positive bulk stratification near the boundary, where isopycnals curve downwards. The stratification thickness was calculated for each flow investigated here, using the linear analytical solutions for the laminar oscillating flows and setting $z_1\gg 1$.

Figure 3.

(a–d) Stratification thickness concept.

The wave- and planar-averaged, wall-normal integrated TKE budget statistics for $\textit {Re}=840$ are shown in figure 4. The statistics were wave-averaged over 5–10 oscillations. All of the integrated TKE budgets at $\textit {Re}=840$, with the exceptions of case 5 and arguably case 7, possess a single burst of chaotic three-dimensional motion that is characterized by a rapid increase in the production rate of the TKE from the across-slope shear, $\mathcal {P}_{13}$, the component of shear parallel to the direction of the oscillating body force. The turbulent bursts, which occur shortly after $t/T\approx 0.5$, preferentially select the phase regime during which the velocity is upslope but decelerating, the sign of the oscillating buoyancy changes from positive to negative and the stratification thickness is negative (as indicated by the dark grey shading in figure 4). The negative stratification thickness preference of the bursts contrasts the low-Reynolds-number, intermittent turbulence regime of Stokes’ second problem, in which a single burst occurs per oscillation, corresponding to the random selection of one of two shear maxima that occur within one period (Spalart & Baldwin Reference Spalart and Baldwin1987).

Figure 4.

(a–h) Wall-normal integrated, planar mean TKE budgets. The grey shading corresponds to the sign of the stratification thickness (negative (positive) represents enhanced (weak or negative) bulk boundary layer stratification). The dashed lines correspond to the time of the minimum total vertical buoyancy gradient in the linear solutions. The far-field velocity oscillates with $-\cos (t)$.

To investigate the role of the linear buoyancy dynamics in the formation of the turbulent bursts in figure 4, the time of the minimum total vertical buoyancy gradient in the linear solutions, which the reader may recall is negative for all of the considered parameter space as shown in figure 2, is depicted as the vertical dashed black line. The maximum TKE production rate by the mean shear approximately coincides in the time of the minimum total vertical buoyancy gradient for cases 1 and 6, the smallest intensity turbulent bursts shown in figure 4. The shear production is predominantly by $\mathcal {P}_{13}$, consistent with the notion that $w'>0$ disturbances initiated by gravitational instabilities impede the mean shear and trigger $u'<0$. The negative Reynolds stresses combined with $\partial _z\bar {u}>0$ during the upslope flow at $t\approx 0.5$ subsequently produce $\mathcal {P}_{13}>0$ as plotted in figure 4. In the rotating cases, $\mathcal {P}_{23}$ is non-zero as mean shear in the along-isobath direction provides a source of TKE which is present in the non-rotating regime. This result suggests that the bursts of TKE production rate by the mean shear are triggered by buoyant ejections of low-momentum fluid upward. However, the triggering of buoyancy ejections is brief, weak and not sustained, because the buoyancy fluxes in figure 4 are negligible prior to bursts in the along-isobath shear production. Thus the gravitational instabilities appear to initiate, but not drive, bursts of chaotic three-dimensional motion.

It is well known that boundary layer turbulence is inherently anisotropic. However, the majority of ocean turbulence measurements measure the fluctuations of the vertical shear of the horizontal velocities and then assume homogeneous isotropic turbulent motion to subsequently estimate the dissipation rate (Polzin & Montgomery Reference Polzin and Montgomery1996; St. Laurent, Toole & Schmitt Reference St. Laurent, Toole and Schmitt2001). The isotropic, homogeneous turbulence dissipation rate of TKE (Taylor Reference Taylor1935) is defined in the sloped coordinate frame as

(4.21)\begin{equation} \varepsilon_{hi}= \tfrac{15}{4}\nu\cos^{2}\theta\left(\overline{(\partial_zu')^{2}} +\overline{(\partial_zv')^{2}}\right). \end{equation}

The wall-normal integrated forms of $\varepsilon _{hi}$ are plotted for the rotating reference frame cases in figure 4. Comparison of the dissipation rate assuming homogeneous isotropic turbulence (4.21) and the full dissipation rate (4.18), shown for the rotating cases, indicates that the assumption of homogeneous isotropic turbulence leads to overpredictions of the dissipation rate by a factor of approximately two. At high Reynolds number in unstratified boundary layers it is well known that the boundary layer dissipation rate is highly anisotropic in the viscous sublayer (Pope Reference Pope2000); therefore we speculate that, if the Prandtl number remains unity, then the anisotropy of the dissipation rate will be similarly significant in the viscous sublayer of the high-Reynolds-number regimes of the flow investigated here.

4.3. Gravitationally unstable rolls

Except for case 6, all of the simulations at $\textit {Re}=840$ exhibited rolls characterized by growing streamwise vorticity. Figure 5 shows the instantaneous vertical velocity of case 2 to illustrate the life cycle of the rolls. In figure 5, red corresponds to upward motions and blue approximately corresponds to downward motions. The generation of two-dimensional convective rolls in the along-isobath/wall-normal ($y$–$z$) plane is visible just prior to the beginning of a burst. At time $t=0.51$ the rolls appear and by $t=0.55$ the rolls have begun to shear apart, erupting into the three-dimensional turbulence at the time of increase in TKE production by mean shear at $t=0.55$.

Figure 5.

(a–f) Contours of the vertical velocity, $w$, for case 2. The contour plots show the vertical velocity at a fixed distance (roughly $\delta$) in the wall-normal direction at six consecutive times. Here $w>0$ is coloured red, while $w<0$ is coloured blue. At $t=0.5$ the across-isobath velocity is positive but begins to decelerate. Simultaneously, (b) two-dimensional rolls form in the $y$–$z$ plane, as heavier fluid is advected over lighter fluid trapped near the wall by the friction.

The rolls formed by gravitational instabilities in figure 5 are qualitatively consistent with rolls observed in oscillating sloping stratified boundary layer experiments by Hart (Reference Hart1971). Hart (Reference Hart1971) identified spanwise plumes and rolls associated with the periodic reversals of the density gradient that qualitatively resembled the rolls that appeared in high-Rayleigh-number Couette flow experiments by Bénard & Avsec (Reference Bénard and Avsec1938), Chandra (Reference Chandra1938) and Brunt (Reference Brunt1951). Perhaps due to the similarity to the convection experiments, the rolls observed by Hart (Reference Hart1971) were referred to as ‘convective rolls’. Linear stability analyses by Deardorff (Reference Deardorff1965), Gallagher & Mercer (Reference Gallagher and Mercer1965) and Ingersoll (Reference Ingersoll1966) revealed that the growth of gravitationally unstable disturbances in high-Rayleigh-number Couette flows is suppressed in the plane of the shear (the streamwise–vertical plane) by the shear (i.e. the suppression of the spanwise vorticity disturbances). They also found that the growth of disturbances in the spanwise–vertical plane (steamwise vorticity disturbances) is unimpeded by the shear and grows in the same manner as pure convection. It has since been established that streamwise (the across-isobath direction) vortices with axes in the direction of a mean shear flow (or ‘rolls’) can arise due to heating or centrifugal effects (Hu & Kelly Reference Hu and Kelly1997). Linear stability analyses by Kaiser (Reference Kaiser2020) of the same regimes indicate that the streamwise (across-isobath) vorticity, which describes the rolls in figure 5, is linearly gravitationally unstable at the low Reynolds numbers investigated here. The initial growth of the rolls in figure 5 appears to have similar attributes.

To verify the hypothesis that gravitational instabilities spawn the rolls, which in turn spawn the turbulent burst, an additional simulation with the same parameters as those of case 2, but with the nonlinear advective terms in the buoyancy equation removed (i.e. replacing the material derivative ${d}_t\tilde {b}$ with $\partial _t\tilde {b}$ in (2.6)), was executed. The simulation of the non-advective buoyancy equation version of case 2 had no turbulent bursts over 10 cycles (all other simulations with bursts developed a burst within 2 cycles). The rolls are a bypass transition mechanism, lifting low-momentum fluid up and bringing high-momentum fluid down into the near-wall flow, and so they transiently destabilize the shear. The transient gravitationally unstable buoyancy gradients, discussed previously, can trigger oscillating boundary layer turbulent bursts even if the buoyancy fluxes are a negligible source of TKE.

Large-eddy simulation of a non-rotating internal wave beam at critical slope and at $Re=10\,500$ ($U_\infty =0.125\ \textrm {m}\ \textrm {s}^{-1}$) by Gayen & Sarkar (Reference Gayen and Sarkar2011) revealed turbulent bursts at the same time in the phase as those shown in figure 4: the bursts occur as the far-field velocity changes from downslope to upslope flow. However, the wall-normal structure of the prescribed internal wave beam investigated in Gayen & Sarkar (Reference Gayen and Sarkar2011) produced turbulence dominated by buoyancy fluxes as the beam generated a much larger supply of potential energy through isopycnal tucking (dense water advected over light water) than is possible by molecular diffusion. Given the differences in the magnitudes of the buoyancy fluxes in Gayen & Sarkar (Reference Gayen and Sarkar2011) and those shown in figure 4, we speculate that the low- to transitional-Reynolds-number regime is characterized by shear-driven turbulent bursts triggered by gravitational instability, while the higher Reynolds number/smaller internal tide vertical wavelength regime is characterized by large buoyancy fluxes due to significant isopycnal tucking.

Although the streamwise rolls are initially two-dimensional, they produce a three-dimensional vorticity field. The inherent three-dimensionality of the gravitational instability is evident in the Boussinesq baroclinic production of vorticity term ($\boldsymbol {\nabla }\times \tilde {b}$) in the absolute vorticity budget for rotating and non-rotating oscillating boundary layers:

(4.22)\begin{equation} \text{baroclinic production} ={C}^{2}\big(\hspace{-4mm}\overbrace{\partial_y\tilde{b}\cot\theta\boldsymbol{i}}^{\substack{\text{nonlinear 2-D rolls} \\ \text{and 3-D bursts}}}\hspace{-1.5mm} +\hspace{2.5mm}(\overbrace{\partial_x\tilde{b}\cot\theta}^{\substack{\text{nonlinear} \\ \text{3-D bursts}}}\hspace{2mm}-\hspace{-2.5mm}\overbrace{\partial_z\tilde{b}}^{\substack{\text{linear flow} \\ \text{and 3-D bursts}}}\hspace{-4mm})\boldsymbol{j}\hspace{4mm} -\hspace{-4mm}\overbrace{\partial_y\tilde{b}\boldsymbol{k}}^{\substack{\text{nonlinear 2-D rolls} \\ \text{and 3-D bursts}}}\hspace{-7.5mm}\big).\end{equation}

The linear oscillating boundary layer vorticity field has only one vorticity component, the spanwise vorticity in the $y$ direction, and the linear rotating oscillating boundary layer vorticity field is comprised of the spanwise vorticity and the streamwise vorticity in the $x$ direction. In either case, only the $\partial _z\tilde {b}\boldsymbol {j}$ term in (4.22) is non-zero. However, the rolls produce gradients in the buoyancy field in the $y$ direction. The first and last terms on the right-hand side of (4.22) indicate that the rolls in the $y$–$z$ plane will inevitably generate vorticity in the streamwise and wall-normal directions; therefore, the rolls may induce three-dimensional motion in the oscillating boundary layers if viscosity and/or other mechanisms that suppress secondary instabilities are overcome. Therefore the coherent structures shown in figure 5, which facilitate the transition to turbulence, may initiate secondary instabilities through three-dimensional baroclinic production of vorticity, a phenomenon that is widely observed in other stratified shear flow instabilities (e.g. Peltier & Caulfield Reference Peltier and Caulfield2003).

4.4. Relaminarization and instability suppression

During the phase of enhanced boundary layer stratification relative to the background, $\delta _s<0$ (the white regions in figure 4), the boundary shear instabilities must overcome not only the stabilizing effect of increased stratification but also the stabilizing effect of the wall that is present regardless of phase. Linear stability analysis by Schlichting (Reference Schlichting1935) yielded a critical gradient Richardson number of 1/24 for a stratified Blasius boundary layer, notably lower than the Miles–Howard theorem threshold for inviscid stratified shear. This suggests that the flow is stabilized with respect to shear perturbations during $\delta _s<0$, although it is difficult to derive a suitable gradient Richardson number from the time-periodic boundary layer velocity and buoyancy solutions (Appendix A) because of inner and outer boundary layer gradients with opposite sign as well as shear and buoyancy gradient sign changes (figure 6). However, if the flow is turbulent during the phase of $\delta _s<0$, turbulence model simulations at high Reynolds number by Umlauf & Burchard (Reference Umlauf and Burchard2011) suggest mixing is much more efficient, presumably because the turbulence intensity required to overcome the strengthened boundary layer stratification must be considerable.

Figure 6.

(a–h) Hövmuller plots of mean stratification for $\textit {Re}=840$. The total wall-normal buoyancy gradients are non-dimensionalized by $N^{2}$. The colour bar axes show that the boundary layer stratification maxima/minima increase/decrease with increasing slope.

Three other mechanisms contribute to the relaminarization of the turbulent bursts in figure 4. First, the turbulent burst diffuses the mean shear and thus its primary energy source. Second, the tidal acceleration opposes the mean shear during the second half of the phase, so the decay and reversal of the shear amplitude imply that less mean flow kinetic energy is available. Third, once the flow reverses the outer boundary layer becomes increasingly stratified, as mentioned previously, when $\delta _s<0$. For a burst to persist across the entire period it must have a constant source of mean shear of large enough magnitude to sustain production of TKE throughout flow reversals and increased stratification.

4.5. The effect of rotation

Only case 5 features sustained turbulence throughout the period (figure 4) and sustained boundary layer buoyancy homogenization (figure 6), because projection of the Coriolis force onto the wall-normal direction increases with decreasing slope, and thus oscillating Stokes–Ekman-like boundary layers develop at low angle. As the slope angle decreases, the boundary layer velocity spiral broadens in the along-isobath direction and mean shear is produced in the along-isobath direction, providing an additional source of TKE through a component of shear production, $\mathcal {P}_{23}$, that is not active in the non-rotating cases (figure 4). The Burger number is the ratio of the squares of the inertial period to the time scale associated with the buoyant restoring force. The slope Burger number accounts for the rotated reference frame, defined as

(4.23)\begin{equation} {Bu}={Ro}^{2}{C}^{2}=\frac{N^{2}\tan^{2}\theta}{f^{2}}. \end{equation}

If $Bu <1$, the buoyant restoring force acts more slowly than the Coriolis force and rotation is significant. Table 3 shows the Burger number for the rotating boundary layer flows, cases 5–8 and 13–16. Cases of the lowest slope Burger number, cases 5 and 13, are influenced the most by rotation. In case 5, the highest-Reynolds-number and lowest-Burger-number flow, the turbulence is sustained throughout the period and the mean velocity field oscillates in a tidal ellipse. These results suggest that, within the investigated parameter regime, low-Burger-number flows are more likely to sustain turbulence throughout the entire period because the tidal ellipse of the boundary layer mean shear can supply TKE even at times in the oscillation when the component of mean shear in the across-isobath direction is weak.

Table 3.

Rotating boundary layer simulation slope Burger numbers.

4.6. Barotropic tide dissipation

In oceanography, the rate of energy loss to drag per square metre of the barotropic flow by tidal bottom boundary layers is often estimated using the quasi-empirical model $D\approx \rho _0c_D|U|U^{2}$ (Hoerner Reference Hoerner1965), where $c_D$ is the dimensionless drag coefficient and $U$ is an estimate of the bulk velocity (Jayne & St. Laurent Reference Jayne and St. Laurent2001; St. Laurent & Garrett Reference St. Laurent and Garrett2002). The rate of energy dissipated by the tide can be represented in terms of watts per metre squared by

(4.24)\begin{equation} D\approx\rho_0\int_0^{30\delta}\bar{\varepsilon}\,\mathrm{d}z, \end{equation}

where $\bar {\varepsilon }$ is the dimensional time mean dissipation rate of TKE (units $\textrm {m}^{2}\ \textrm {s}^{-3}$). Lacking a time-mean length scale for the boundary layer thickness that applies to bursting flows, we spatially average the dissipation rate of TKE over the region of 5 m near the wall, equivalent to averaging over approximately $0< z<30\delta$. For flat-plate boundary layers, the transitional flow regime is characterized by drag coefficients in the range $0.001\leqslant c_D\leqslant 0.005$, for $1<\textit {Re}<10^{3}$ or equivalently $1<\textit {Re}_L<10^{6}$ (Hoerner Reference Hoerner1965). The drag coefficients for the low-Reynolds-number boundary layers in this study were calculated as

(4.25)\begin{equation} c_D=\frac{D}{\rho_0U_\infty^{3}}. \end{equation}

The drag coefficients are shown in table 4. The drag coefficient values shown in table 4 mostly fall within the expected range for flat-plate boundary layers, with the exception of the steepest slope case at $\textit {Re}=840$, case 8. The drag coefficients are small at $\textit {Re}=420$ for all but the steepest slope angles in the rotating reference frame (cases 15 and 16). In table 4, $\sim$0 represents negligible drag. The drag coefficients increase with slope at constant Reynolds number, and they effectively vanish in a portion of the range $420<\textit {Re}<840$ on lower slopes where the flow is in the laminar or disturbed laminar abyssal slope regime.

Table 4.

Drag coefficients, averaged dissipation rates and Ozmidov–Kolmogorov length-scale ratios estimated from time-mean simulation results.

The average boundary layer dissipation rates $\varepsilon _{{BL}}$ and the ratio of Ozmidov to Kolmogorov length scales $(L_O/L_K)^{4/3}$ are also shown in table 4, where

(4.26)\begin{gather} \varepsilon_{{BL}}\approx\frac{1}{30\delta}\int_0^{30\delta}\bar{\varepsilon}\,\mathrm{d}z, \end{gather}

(4.27)\begin{gather}\left(\frac{L_O}{L_K}\right)^{4/3}=\frac{\varepsilon_{{BL}}}{\nu N^{2}}. \end{gather}

The use of the ratio of Ozmidov to Kolmogorov length scales $L_O/L_K$ to characterize the bulk effect of stratification on turbulence can be traced to the isotropy parameter of Gargett, Osborn & Nasmyth (Reference Gargett, Osborn and Nasmyth1984). More recently, Gargett's ratio raised to the 4/3 power has been referred to as the isotropy parameter $I$ (Thorpe Reference Thorpe2007; Ivey, Winters & Koseff Reference Ivey, Winters and Koseff2008) and as the buoyancy Reynolds number $Re_b$ (Ivey et al. Reference Ivey, Bluteau, Gayen, Jones and Sohail2021; Mashayek et al. Reference Mashayek, Baker, Cael and Caulfield2022). To avoid confusion, we refer to Gargett's ratio raised to the 4/3 power as the Ozmidov–Kolmogorov length-scale ratio.

Sign changes in the oscillating boundary layer buoyancy gradient (see figure 6 and the buoyancy solutions for the laminar flow provided in Appendix A) and therefore it does not lend itself to a simple time-mean estimate of the buoyancy gradient; therefore, the background buoyancy gradient $N^{2}$ is used in (4.27). Since the Ozmidov–Kolmogorov length-scale ratio does not distinguish between weak stratification and strong turbulence, the ${O}(100)$ values presented in table 4 should be interpreted as weak turbulence, or chaotic motions with little scale separation that overcome weaker stratification. The dissipation rates presented in table 4 are consistent with observed minimum abyssal dissipation rates of $\varepsilon \sim 10^{-10}\ \textrm {m}^{2}\ \textrm {s}^{-3}$ (Thorpe Reference Thorpe2007). The low levels of dissipation rate and the relaminarization of the turbulent bursts imply that in all but the sustained turbulence case (case 5) the irreversible buoyancy flux is negligible. However, $I>100$ in stratified homogeneous turbulence is associated with significant irreversible buoyancy fluxes in laboratory experiments (Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005; Ivey et al. Reference Ivey, Winters and Koseff2008), suggesting (1) that a better estimation of the characteristic boundary layer buoyancy gradient is required and (2) that the interpretation of the Ozmidov–Kolmogorov length-scale ratio near adiabatic boundaries is ambiguous because $(L_O/L_K)^{4/3}\rightarrow \infty$ on a steady adiabat.