1. Introduction
In both the atmosphere and the ocean, a significant portion of energy is stored in and transported by internal gravity wave motions. The means by which these waves break and dissipate their energy via small-scale turbulence is a classical and important problem with an extremely broad scope, encompassing a variety of linear and nonlinear interactions and instabilities over a range of scales. Wave breaking is often inferred observationally by considering one-dimensional vertical spectra of shear and isopycnal strain. At vertical scales larger than O(10 m), there is generally a close match between observed spectra and the empirical Garrett–Munk spectrum modelling dynamics dominated by weakly nonlinear internal waves (Garrett & Munk Reference Garrett and Munk1975). On the other hand, energetic turbulence typically exhibits an inertial subrange over scales smaller than the Ozmidov length scale (below which the vertical extents of turbulent motions are not affected at leading order by the stratification). It is routinely observed that the smallest scales of internal waves are larger than the largest scales of isotropic turbulence (Kunze Reference Kunze2019; D’Asaro Reference D’Asaro2022). The dynamics across the resulting intermediate range of scales are in general not well understood, though evidence suggests they may be of leading-order importance in determining the properties of turbulence and mixing at dissipative scales. For example, it is becoming increasingly apparent from observations in a variety of environments that the mixing efficiency of stratified turbulence – often assumed to be constant throughout the ocean interior – exhibits up to an order of magnitude variability over spatial and temporal scales significantly larger than those associated with turbulence; for recent examples, see Ijichi & Hibiya (Reference Ijichi and Hibiya2018), Alford et al. (Reference Alford, Wynne-Cattanach, Le Boyer, Couto, Voet, Spingys, Castro, Forryan, Garabato and van Haren2025), Lewin et al. (Reference Lewin, Kaminski, McSweeney and Waterhouse2025). The implication is therefore that the dynamics at the scales of wave breaking are of leading-order importance for the subsequent properties of mixing, a feature framed by Caulfield (Reference Caulfield2021) in the terms ‘history matters’.
To capture the influence of dynamics at scales over which internal waves dominate energetics and to model downscale energy transfer in a mathematically tractable way, our study will consider the evolution and eventual instability of linear monochromatic internal waves under the Boussinesq approximation and in the limit of weak rotation
$f\ll N$
, where
$f$
and
$N$
are the Coriolis and buoyancy frequency, respectively. Floquet theory may be used to show that even plane periodic waves of this kind are unstable at arbitrarily small amplitudes to a parametric instability driven by wave self-interaction (Mied Reference Mied1976; Klostermeyer Reference Klostermeyer1982; Lombard & Riley Reference Lombard and Riley1996a
). However, the growth rate of this instability is typically expected to be small, in particular compared with the rate of modulation of the local wave amplitude due to wave–mean flow interactions (see, e.g. Sutherland Reference Sutherland2006a
). In either case, when the wave amplitude is sufficiently large (such that instabilities grow quickly) or the frequency is sufficiently small, a local linear stability analysis may be performed on vertical velocity and buoyancy profiles derived from or approximating the wave polarisation relations (Dunkerton Reference Dunkerton1997; Fritts & Alexander Reference Fritts and Alexander2003). This provides potentially useful information about, for example, the roles of convection and shear in the resulting turbulent state (Lelong & Dunkerton Reference Lelong and Dunkerton1998a
,
Reference Lelong and Dunkertonb
), which are found to vary significantly with the properties of the background wave. In this framework, there is therefore motivation to explore the range of possible unstable configurations that are dynamically accessible during the lifecycle of an internal wave.
A classical paradigm for the modulation and eventual breaking of small-scale, vertically propagating wave packets is the rapid increase in vertical wavenumber due to strong refraction near critical levels in a vertical shear flow. Such critical levels, originally studied by Bretherton (Reference Bretherton1966) and Booker & Bretherton (Reference Booker and Bretherton1967), are locations where the horizontal phase velocity of a wave group approaches the speed of the horizontal background flow. A Wentzel–Kramers–Brillouin–Jeffreys (WKBJ) analysis may be used to demonstrate that the intrinsic frequency (the frequency measured by an observer travelling at the speed of the background flow) decreases and approaches zero at the critical level (in the absence of rotation). The wave becomes increasingly steep leading to the possibility of breaking, as has been captured experimentally by Koop (Reference Koop1981), Thorpe (Reference Thorpe1981) and in turbulence-resolving numerical simulations by Winters & D’Asaro (Reference Winters and D’Asaro1994), Dörnbrack (Reference Dörnbrack1998), Howland, Taylor & Caulfield (Reference Howland, Taylor and Caulfield2021). Approximate local forms for the vertical velocity and density profiles indicate that strong vertical shear and statically unstable density profiles are both possible (Winters & Riley Reference Winters and Riley1992); correspondingly, the emerging modes of instability are dynamically rich and inherently three dimensional, with important consequences for the resulting turbulent energetics (Winters & D’Asaro Reference Winters and D’Asaro1994). A notable success of this paradigm has been the eikonal theory of ray tracing, which, when applied to test waves in a measured background shear flow, may be used to estimate the fraction of energy that is cascaded to small scales and eventually dissipated (Henyey, Wright & Flatté Reference Henyey, Wright and Flatté1986). This method forms the basis of the fine-scale parameterisation approach for estimating energy dissipation from measurements of the larger-scale internal wave field (Polzin et al. Reference Polzin, Garabato, Alberto, Huussen, Sloyan and Waterman2014).
The problem of internal wave modulation by a horizontal (as opposed to vertical) shear has received relatively less attention in the literature. In this case, WKBJ analysis indicates the existence of a new type of critical level where the intrinsic wave frequency approaches the buoyancy frequency (Olbers Reference Olbers1981; Basovich & Tsimring Reference Basovich and Tsimring1984; Badulin, Shrira & Tsimring Reference Badulin, Shrira and Tsimring1985), resulting in dynamics that are characteristically distinct from the vertical shear scenario described above. These critical levels are sometimes referred to as baroclinic critical layers and have been shown to play a significant role in the dynamics of rotating stratified shear flows (Marcus et al. Reference Marcus, Pei, Jiang and Hassanzadeh2013) and tilted stratified vortices (Boulanger et al. Reference Boulanger, Meunier and Le Dizès2007). Nonlinear dynamics within baroclinic critical levels under a steady wave forcing in a stratified uniform horizontal shear flow have been described by Wang & Balmforth (Reference Wang and Balmforth2020, Reference Wang and Balmforth2021). More relevant to this study is the work of Staquet & Huerre (Reference Staquet and Huerre2002), who investigated the transient problem of an inertia-gravity wave packet approaching a critical level in a horizontal shear layer using numerical simulations. Wave breaking and significant subsequent turbulence production was observed due to a transient growth in wave energy extracted from the shear in the vicinity of the critical level. Building on this result, Bakas & Farrell (Reference Bakas and Farrell2009a , Reference Bakas and Farrellb ) analysed the process by which linear internal waves may extract energy from a background unbounded uniform horizontal shear flow, thus potentially leading to instability and wave breaking. They showed that, in fact, two distinct mechanisms for energy perturbation growth exist. The first corresponds to wave amplitude growth due to extraction of energy from the background shear through downgradient Reynolds stresses, consistent with the dynamics near critical levels observed by Staquet & Huerre (Reference Staquet and Huerre2002). The second results from wave advection of the mean flow momentum that may produce significant vertical shear, thus providing an alternative route to wave breaking. This latter mechanism operates in the same way as the well-known lift-up mechanism for unstratified flows (Ellingsen & Palm Reference Ellingsen and Palm1975).
Motivated by the discussion above, in this study we revisit the initial value problem considered by Staquet & Huerre (Reference Staquet and Huerre2002) consisting of a monochromatic internal wave packet in a background horizontal shear layer
\begin{align} \boldsymbol{U}^*(y^*) = \Delta u^* \tanh (y^*/h^*)\boldsymbol{e}_x, \end{align}
where
$\Delta u^*$
and
$h^*$
are constants characterising the velocity and length scales of the shear, and stars denote dimensional variables. The linear stability of this base flow in a uniform stratification was studied by Deloncle, Chomaz & Billant (Reference Deloncle, Chomaz and Billant2007), who showed that the most unstable mode is two-dimensional with a zero vertical wavenumber. However, as the (horizontal) Froude number
$ \textit{Fr}_h = \Delta u^* /N^*h^*$
decreases, the growth rates of modes with vertical wavenumber
$k_z\gt 0$
increase, thereby increasing the tendency for a vertical structure to emerge in practice. Indeed, the fully nonlinear simulations of Basak & Sarkar (Reference Basak and Sarkar2006) highlight such three-dimensional behaviour, though turbulence was not sustained for Froude numbers
$ \textit{Fr}_h\leqslant 1$
. Recently, Lewin & Caulfield (Reference Lewin and Caulfield2024) demonstrated the existence of energetic turbulence in this strongly stratified regime, provided an additional small-amplitude wave-like perturbation was supplied as an initial condition. Waves were able to rapidly extract energy from the horizontal shear and amplify vertical shear in the flow by multiple orders of magnitude through the lift-up mechanism described by Bakas & Farrell (Reference Bakas and Farrell2009a
), resulting in a significant departure from the dynamics observed by Basak & Sarkar (Reference Basak and Sarkar2006). In this study, the linear instability of the stratified horizontal shear flow is suppressed by choosing a sufficiently large shear layer width
$h^*$
relative to the domain size in the
$x$
direction (precisely, unstable perturbations are restricted to wavelengths
$\lambda _x^*\gt 2\pi h^*$
; see Deloncle et al. Reference Deloncle, Chomaz and Billant2007).
We investigate the local structure and amplitude of wave-like perturbations in the shear flow (1.1) using ray-tracing theory derived under the WKBJ approximation. A principle aim is to explore whether this method, similar to that of Winters & Riley (Reference Winters and Riley1992) for the vertical shear case, can make useful predictions regarding the local flow fields of the waves from which instabilities may develop and lead to breaking. In particular, we wish to determine whether the effects of the two mechanisms for wave energy growth in a horizontal shear described by Bakas & Farrell (Reference Bakas and Farrell2009a ) may be captured by the local theory. Using this theory as a guide, we then proceed to construct a parameter space (based on the scales associated with the wave and shear flow) for investigation using fully nonlinear direct numerical simulations (DNS). In cases where wave breaking leads to turbulence, a key question to be addressed is: How do the eventual momentum and energy exchanges between the mean flow and perturbations depend on the route to wave breaking? An obvious point of reference is the vertical critical level problem, where it has been typically observed that a significant fraction of the incoming wave energy is absorbed by the mean flow near the critical level, even in cases where wave breaking results in considerable turbulent dissipation (Howland et al. Reference Howland, Taylor and Caulfield2021; Winters & D’Asaro Reference Winters and D’Asaro1994; Thorpe Reference Thorpe1981).
The remainder of the paper proceeds as follows. In § 2 we describe the basic problem set-up and detail the ray-tracing theory used in the analysis. In § 3 we outline a suite of DNS with parameters chosen based on predictions from the theoretical analysis. A qualitative overview of the DNS results and a comparison with the theory are presented in § 4, followed by a detailed description of momentum and energy budgets in §§ 5 and 6. Finally, we discuss some broader implications of the results and conclude in § 7.
2. Problem set-up and ray-tracing theory
2.1. Initial value problem
We consider a Boussinesq fluid with uniform density stratification characterised by buoyancy frequency
$N^{*2}$
. Precisely, the buoyancy field is decomposed as
$b^* = N^{*2}z^* + \theta ^*$
, where
$\theta ^* = -g^*\rho ^*/\rho _0^*$
is the buoyancy perturbation corresponding to density perturbation
$\rho ^*$
away from hydrostatic balance. Here,
$g^*$
is the acceleration due to gravity in the vertical (
$z$
) direction and
$\rho _0^*$
is a constant reference density. The equations of motion are non-dimensionalised with time and length scales obtained from the buoyancy frequency
$N^*$
and a characteristic wavenumber (or inverse length scale)
$K^*$
to be associated with the initial wave. In this way, dimensionless variables are defined as
\begin{align} t = N^*t^*, \quad x_i = K^* x_i^*, \quad u_i = \frac {K^*}{N^*} u_i^*, \quad \theta = \frac {K^*}{N^{*2}}\theta ^*, \quad p = \frac {K^{*2}}{N^{*2}} \frac {p^*}{\rho _0^*}, \end{align}
so that the equations of motion are
\begin{align} \displaystyle \frac {\mathrm{\partial }\boldsymbol{u}}{\mathrm{\partial }t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} = -\boldsymbol{\nabla }p + \theta \boldsymbol{e}_z+ \frac {1}{Re_w} {\nabla} ^2\boldsymbol{u}, \\[-28pt] \nonumber \end{align}
\begin{align} \displaystyle \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0, \\[-28pt] \nonumber \end{align}
\begin{align} \displaystyle \frac {\mathrm{\partial }\theta }{\mathrm{\partial }t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\theta + w = \frac {1}{\textit{Re}_w\textit{Pr}}{\nabla} ^2 \theta . \\[3pt] \nonumber \end{align}
The wave Reynolds number
$Re_w \equiv N^*/(K^{*2}\nu ^*)$
(cf. Fringer & Street Reference Fringer and Street2003) and Prandtl number
$ \textit{Pr} \equiv \nu ^*/\kappa ^*$
are defined in terms of the molecular kinematic viscosity
$\nu ^*$
and buoyancy diffusivity
$\kappa ^*$
. Note that, under this choice of non-dimensionalisation, the coefficient in front of the buoyancy term in the momentum equation (2.2a
) is unity.
In the absence of any background flow, monochromatic periodic linear internal waves with wavenumber
$\boldsymbol{k}_0 = (k_{x0},k_{y0},k_{z0})$
and frequency
$\varOmega _0$
are given by
\begin{eqnarray} [u,v,w,\theta ] = [U_0,V_0,W_0,\varTheta _0]s_0\exp (i\boldsymbol{k}_0\boldsymbol{\cdot }\boldsymbol{x} - i\varOmega _0 t), \end{eqnarray}
where the amplitude coefficients are defined by the polarisation relations
\begin{align} \begin{split} \varTheta _0 = \frac {1}{k_{z0}} ,\quad U_0 = -\frac {\mathrm{i} k_{x0} }{\varOmega _0 |\boldsymbol{k}_0|^2}, \\[6pt]V_0 = -\frac {\mathrm{i} k_{y0}}{\varOmega _0|\boldsymbol{k}_0|^2}, \quad W_0 = \mathrm{i}\frac {\varOmega _0}{k_{z0}}, \end{split} \end{align}
and the frequency
$\varOmega _0$
is given by the dispersion relation
\begin{eqnarray} \varOmega _0^2 = \frac {k_{x0}^2+k_{y0}^2}{|\boldsymbol{k}|^2}. \end{eqnarray}
The wave steepness
$s_0$
appearing in (2.3) is a dimensionless amplitude defined such that
$s\gt 1$
corresponds to the existence of regions of local static instability in the total buoyancy field with
$\partial b/\partial z \lt 0$
.
Schematic showing a representative initial condition for the problem under consideration. The shape and location of the incident wave packet is shown by plotting isocontours of the
$y$
component of the velocity field (
$v$
). Notably, the wave packet is attenuated in
$y$
based on the exponential envelope defined in (2.6). The tilt of the isocontours in the
$x-z$
plane is governed by the wavenumber ratio
$k_x/k_z$
.
We consider a wave packet bounded by a Gaussian envelope in
$y$
of width
$a$
, centred at an initial location
$y_0$
. At time
$t=0$
, the velocity and buoyancy perturbation fields are given by
$\boldsymbol{u}_{{wave}}(\boldsymbol{x})$
and
$\theta _{{wave}}(\boldsymbol{x})$
, where
\begin{eqnarray} [u_{{wave}},v_{{wave}},w_{{wave}},\theta _{{wave}}] = [U_0,V_0,W_0,\varTheta _0]s_0\exp \left (\mathrm{i}\boldsymbol{k}_0\boldsymbol{\cdot }\boldsymbol{x}-\frac {(y-y_0)^2}{a^2}\right )\!, \end{eqnarray}
with amplitudes given by (2.4). Taking
$y_0\lt 0$
, the sign of
$k_{y0}$
is chosen to be positive such that the wave packet propagates in the positive
$y$
direction. We note that wave packets of this form have also been studied by both Staquet & Huerre (Reference Staquet and Huerre2002) and Bakas & Farrell (Reference Bakas and Farrell2009a
). In dimensionless form, the background shear flow into which the wave will propagate is
\begin{align} \boldsymbol{U}(y) = \Delta u \tanh (y/h)\boldsymbol{e}_x, \end{align}
where
$\Delta u = \Delta u^* K^*/N^*$
and
$h = h^*K^*$
are dimensionless constants. Putting everything together, the initial velocity and buoyancy fields are
$\boldsymbol{u}(\boldsymbol{x},t=0) = \boldsymbol{U}(y) + \boldsymbol{u}_{{wave}}(\boldsymbol{x})$
and
$ \theta (\boldsymbol{x},t=0) = \theta _{{wave}}(\boldsymbol{x})$
. A schematic of the flow described by this initial condition is shown in figure 1.
2.2. Linear dynamics
We wish to describe the propagation of the wave packet given by (2.6) in the background flow (2.7), with an eye towards predicting properties of the local wave state that determine the dynamics of wave breaking. A natural approach is to start from the inviscid equations of motion linearised about the mean shear
$U(y)$
, i.e.
\begin{align} \displaystyle \frac {\partial u}{\partial t} + U(y)\frac {\partial u}{\partial x} + v\frac {\partial U}{\partial y} = -\frac {\partial p}{\partial x}, \\[-28pt] \nonumber \end{align}
\begin{align} \displaystyle \frac {\partial v}{\partial t} + U(y)\frac {\partial v}{\partial x} = -\frac {\partial p}{\partial y}, \\[-28pt] \nonumber \end{align}
\begin{align} \displaystyle \frac {\partial w}{\partial t} + U(y)\frac {\partial w}{\partial x} - \theta = -\frac {\partial p}{\partial z}, \\[-28pt] \nonumber \end{align}
\begin{align} \displaystyle \frac {\partial \theta }{\partial t} + U(y)\frac {\partial \theta }{\partial x} +w= 0, \\[0pt] \nonumber \end{align}
seeking WKBJ solutions of the form
\begin{eqnarray} [u,v,w,\theta , p] = [\tilde {U}, \tilde {V}, \tilde {W}, \tilde {\varTheta }, \tilde {P}](\boldsymbol{x},t)s(\boldsymbol{x},t) \exp (\mathrm{i}\varphi (\boldsymbol{x},t)/\epsilon ) + O(\epsilon ), \end{eqnarray}
where
$\epsilon$
is a small parameter characterising the relative scale of phase variations. In the problem at hand, the WKBJ framework describes waves with wavelengths that are asymptotically small relative to the width of the wave envelope
$a$
and the width of shear layer
$h$
. Precisely, a local frequency
$\omega (\boldsymbol{x},t)$
and wavenumber
$\boldsymbol{k}(\boldsymbol{x},t)$
may be defined in terms of the phase function
$\varphi (\boldsymbol{x},t)$
as
\begin{eqnarray} \omega = -\frac {1}{\epsilon }\partial _t \varphi , \quad \boldsymbol{k} = \frac {1}{\epsilon }\boldsymbol{\nabla }\varphi . \end{eqnarray}
Substituting into (2.8), one obtains locally a set of algebraic equations for
$U_0$
,
$V_0$
,
$W_0$
,
$\varTheta _0$
and
$P_0$
. Eliminating variables, at leading order the incompressibility condition gives
\begin{eqnarray} |\boldsymbol{k}|^2\varOmega ^3 - (k_x^2+k_y^2)\varOmega -\mathrm{i}\frac {\partial U}{\partial y}k_xk_y(\varOmega ^2-1) = 0, \end{eqnarray}
where
$\varOmega \equiv \omega - Uk_x$
is the intrinsic frequency measured by an observer travelling at the speed of the background flow. The frequency
$\omega$
measured by a stationary observer is often referred to as the absolute (or extrinsic) frequency. The two are related by the Doppler shift term
$U(y)k_x$
. For a non-trivial balance in (2.11), we must have
$\varOmega = O(1)$
. It follows that the first two terms are
$O(\epsilon ^{-2})$
, so that the standard dispersion relation (2.5) is recovered when
$(\partial U/\partial y )k_xk_y = O(\epsilon ^{-1})$
. Because
$k_x$
is a priori assumed to be
$O(1/\epsilon )$
, the usual approach in the literature is to assume that
$U$
(and hence
$\partial U/\partial y$
) is
$O(\epsilon )$
(see, e.g. Olbers Reference Olbers1981). However, it will be shown below that, due to the fact that the background flow does not vary in the
$x$
direction,
$k_x$
is constant along rays. Therefore, the assumption that
$k_x=1/\epsilon$
and, hence,
$U = O(\epsilon )$
may be relaxed (provided the wave amplitude varies sufficiently slowly in
$x$
), as Bretherton (Reference Bretherton1966) himself pointed out. Thus, we may instead consider the case
$k_x = O(1)$
with
$U=O(1)$
, in which case local polarisation relations are given by
\begin{align} \begin{array}{cc} \tilde {\varTheta } = \frac {1}{k_{z}}, & \tilde {U} = -\frac {\mathrm{i} k_{x}}{\varOmega |\boldsymbol{k}|^2} - \frac {k_y}{\varOmega ^2 |\boldsymbol{k}|^2} \frac {\partial U}{\partial y}, \\[12pt] \tilde {V} = -\frac {\mathrm{i} k_{y}}{\varOmega |\boldsymbol{k}|^2} , & \tilde {W} = \mathrm{i} \frac { \varOmega }{k_{z}}. \end{array} \end{align}
The existence of the second term in the expression for
$\tilde {U}$
in (2.12) at leading order whilst the local dispersion relation remains equivalent to (2.5) is novel to the horizontally sheared system in the case when
$U=O(1)$
. It is important to keep in mind here that
$s$
,
$\boldsymbol{k}$
and
$\varOmega$
are functions of
$\boldsymbol{x}$
and t.
Changes in
$\boldsymbol{k}$
(and hence
$\varOmega$
, via the dispersion relation) may be determined using the standard kinematic theory of ray tracing, which involves computing wave properties along ray paths
$\boldsymbol{x}(t)$
that are everywhere tangent to the group velocity field measured by a stationary observer
$\boldsymbol{c}_g = \partial \omega / \partial \boldsymbol{k} = U\boldsymbol{e}_x + \partial \varOmega /\partial \boldsymbol{k}$
. Precisely,
\begin{eqnarray} \frac {{\rm d} x}{{\rm d} t} = U + \frac {\partial \varOmega }{\partial k_x}, \ \ \frac {{\rm d} y}{{\rm d} t} = \frac {\partial \varOmega }{\partial k_y}, \ \ \frac {{\rm d} z}{{\rm d} t} = \frac {\partial \varOmega }{\partial k_z}, \end{eqnarray}
where
${\rm d}/{\rm d}t = \partial /\partial t + \boldsymbol{c}_g\boldsymbol{\cdot }\boldsymbol{\nabla}$
. The change in
$\boldsymbol{k}$
along rays is given by
\begin{eqnarray} \frac {{\rm d}k_x}{{\rm d}t} = 0 , \ \ \frac {{\rm d}k_y}{{\rm d}t} = -k_x\frac {\partial U}{\partial y}, \ \ \frac {{\rm d}k_z}{{\rm d}t} = 0. \end{eqnarray}
This justifies the relaxation of the assumption that
$k_x=O(1/\varepsilon )$
above. Note in particular that when
$k_x=0$
initially such that the group and phase velocities are orthogonal to the direction of the background flow,
$k_y$
and hence
$\varOmega$
, remain constant, that is, the intrinsic wave properties are entirely unaffected by the shear.
2.2.1. Trapping and turning levels
Given
$U(y)=\Delta u \tanh (y/h)$
, it follows from (2.14) that
$k_y$
increases or decreases monotonically for
$k_x\lt 0$
or
$k_x\gt 0$
, respectively. We assume henceforth that
$k_y\gt 0$
initially, so that
$k_x\lt 0$
means that
$k_y\gt 0$
for all times. In particular, the situation may arise where
$k_y\to \infty$
, in which case the intrinsic frequency approaches the buoyancy frequency
$\varOmega \to 1$
(corresponding to
$\varOmega ^*\to N^*$
in dimensional variables). The location
$y=y_t$
where
$\varOmega = 1$
(when it exists) is referred to as a trapping level, where the group velocity in the
$y$
direction approaches zero asymptotically. (Formally however, the WKBJ theory becomes invalid before the trapping level is reached.) It follows from the Hamiltonian structure of (2.13) and (2.14) that
$\omega$
is constant along rays (e.g. Sutherland Reference Sutherland2010; Bühler Reference Bühler2014). Then, equating the initial value of
$\omega =\omega _0$
(for which the intrinsic frequency
$\varOmega = \varOmega _0)$
with the value at the trapping level determines
$y_t$
implicitly as
\begin{eqnarray} \tanh (y_t/h) = \frac {\varOmega _0-1}{\Delta u k_x}+\tanh (y_0/h), \end{eqnarray}
where
$y_0$
is the initial location of the wave packet. A detailed description of the asymptotic approach to
$y_t$
is provided by Badulin et al. (Reference Badulin, Shrira and Tsimring1985). The trapping level exists only when the right-hand side of (2.15) falls between
$-1$
and
$1$
. Assuming without loss of generality that the wave packet propagates in the positive
$y$
direction and is also located initially outside of the shear layer such that
$\tanh (y_0/h)\approx -1$
, it follows that we must have
$0\leqslant (\varOmega _0-1)/(\Delta u k_x)\leqslant 2$
for the trapping level to exist.
When
$k_x\gt 0$
,
$k_y$
decreases and internal reflection may occur where
$k_y$
passes through
$0$
and becomes negative. However, the WKBJ analysis fails for small
$k_y$
and a separate solution must be obtained in the vicinity of the turning level or an alternative approach used (e.g. Bakas & Farrell Reference Bakas and Farrell2009a
). The location of the turning level where
$k_y=0$
(when it exists) may be obtained from an equation analogous to (2.15).
2.2.2. Kinetic and potential energy partition
The local wave energy density
$E=E_K + E_P = |\boldsymbol{u}|^2/2+|\theta |^2/2$
, where
$E_K$
and
$E_P$
are the kinetic and potential energy densities, is
\begin{eqnarray} E = \frac {s^2}{2}\Bigl [ \overbrace {\frac {1}{k_z^2}}^{O(\epsilon ^2)} +\overbrace {\left (\frac {k_y}{|\boldsymbol{k}|^2\varOmega ^2}\frac {\partial U}{\partial y}\right )^2}^{O(\epsilon ^2 (\partial U/\partial y)^2)} \Bigr ]. \end{eqnarray}
It follows from (2.12) and (2.16) that wave advection of momentum can result in ‘excess’ kinetic energy contained purely in the streamwise velocity component, here referring to energy not accounted for by the usual wave energy density
$s^2/2k_z^2$
. The presence of this additional energy has no effect on the dynamics of wave propagation (for example, the group velocity and dispersion relation) at leading order because
$Uk_x\sim O(1)$
. In particular, this applies to the conservation equation for wave action
$\mathcal{A}$
(e.g. Bretherton & Garrett Reference Bretherton and Garrett1968), which for the case of a wave packet periodic in
$x$
and
$z$
reads
\begin{eqnarray} \frac {\partial \mathcal{A}}{\partial t} + \frac {\partial }{\partial y}(c_{g,y}\mathcal{A}) = -\frac {1}{Re}|\boldsymbol{k}|^2\mathcal{A} \iff \frac {{\rm d}\mathcal{A}}{{\rm d}t} = -\mathcal{A}\frac {\partial c_{g,y}}{\partial y} - \frac {1}{Re}|\boldsymbol{k}|^2\mathcal{A}, \end{eqnarray}
where
$\mathcal{A} \equiv s^2/(2 k_z^2)/\varOmega$
. Here we have also included a viscous term on the right-hand side due to Grimshaw (Reference Grimshaw1974) (see also Koop Reference Koop1981; Howland et al. Reference Howland, Taylor and Caulfield2021).
The discussion above suggests that perturbation energy may change due to two distinct mechanisms. Wave advection of momentum leads to a change in the kinetic energy contained entirely in the streamwise velocity component according to the local polarisation relations (2.12). On the other hand, wave refraction by the horizontal shear modifies the intrinsic frequency and group velocity and, hence, changes the local wave steepness
$s$
via wave action conservation. In the context of the ray-tracing equations, these behaviours may be interpreted as local manifestations of the global transient energy growth mechanisms referred to by Bakas & Farrell (Reference Bakas and Farrell2009a
,Reference Bakas and Farrell
b
) as growth due to the roll (lift-up) mechanism and to downgradient Reynolds stresses, respectively. Note that these mechanisms operate independently of one another when
$Uk_x\sim O(1)$
, but may be simultaneously active.
2.2.3. Limiting dynamics and instability
As noted above, when
$k_x=0$
, it is clear from (2.14) that the local wave state is unaffected by the shear flow. Thus, in this case, energy growth occurs purely through wave advection of mean flow momentum. As long as perturbations remain independent of
$x$
, the total streamwise momentum
$U(y)+u(y,z,t)$
behaves simply as a passive scalar in the system. This behaviour is often referred to as the ‘lift-up’ mechanism for shear flow perturbations dating back to the works of Moffat (Reference Moffat1967), Ellingsen & Palm (Reference Ellingsen and Palm1975), Landahl (Reference Landahl1980) and leads to the generation of
$y$
vorticity and, hence, vertical shear via vortex tilting (Butler & Farrell Reference Butler and Farrell1992). Note that
$\varOmega ^2 = k_y^2/|\boldsymbol{k}|^2$
in this case so that
\begin{eqnarray} E = \frac {s^2}{2}\left [ \frac {1}{k_z^2} +\left (\frac {1}{k_y}\frac {\partial U}{\partial y}\right )^2 \right ]\!. \end{eqnarray}
Hence, the ratio of the bracketed terms in (2.18) is
$1+ Fr_h^2 k_z^2/k_y^2$
, where
$ \textit{Fr}_h=\partial U/\partial y = (\partial U^*/\partial y^*)/N^*$
is the local horizontal Froude number associated with the background shear flow. Though the analysis here is strictly local, it is worth pointing out that this is an equivalent scaling to that derived for the optimal maximum total energy growth for
$x$
-independent perturbations in an unbounded uniform horizontal shear flow by Bakas & Farrell (Reference Bakas and Farrell2009a
). The greatest perturbation energies result from waves with a large horizontal to vertical aspect ratio
$k_z/k_y$
in a flow with large
$ \textit{Fr}_h$
.
More generally, the ratio of the two components of the wave energy density in (2.16) forms a dimensionless parameter
$F$
defined by
\begin{eqnarray} F^2 = \frac {k_y^2k_z^2}{\varOmega ^4|\boldsymbol{k}|^4}\left (\frac {\partial U}{\partial y}\right )^2 =\frac {k_z^2}{k_y^2(1+k_x^2/k_y^2)^2}\left (\frac {\partial U}{\partial y}\right )^2 , \end{eqnarray}
describing the instantaneous proportion of excess kinetic energy due to wave advection of momentum. Note that
$F$
varies along rays due to changes in
$k_y$
and
$\partial U/\partial y$
. The evolution may be readily calculated by solving the ray-tracing equations (2.14) numerically. Badulin et al. (Reference Badulin, Shrira and Tsimring1985) demonstrated that the wave steepness
$s$
increases due to energy extraction from the shear near trapping levels. As suggested by Bakas & Farrell (Reference Bakas and Farrell2009a
), when
$k_x\lt 0$
(so that
$k_y$
is increasing) and
$F$
is small so that wave advection of momentum is relatively unimportant, corresponding amplitude growth may lead to the formation of regions of unstable buoyancy gradients and, hence, the potential for convective instability. As
$F$
increases, the importance of wave advection of momentum increases, which leads to a potentially significant increase in the vertical shear. This can be seen in our analysis by using the local wave polarisation relations (2.12) to compute the fact that
$|\partial u/\partial z|^2 +|\partial v/\partial z|^2 = k_z^2/|\boldsymbol{k}|^2 +F^2$
.
Dynamically, the tendency for vertical shear to destabilise the flow is characterised by a local gradient Richardson number
\begin{eqnarray} \textit{Ri}_g(\boldsymbol{x},t) = \frac {1 + \partial \theta /\partial z}{S^2} = \frac {1 + \partial \theta /\partial z}{(\partial u_{\beta }/\partial z)^2}, \end{eqnarray}
defined here based on a locally parallel shear flow
$u_{\beta }= u\sin \beta +v\cos \beta$
obtained from the projection of the horizontal velocity field onto a chosen azimuthal direction
$\beta$
in which unstable perturbations are assumed to grow. A criterion of
$ \textit{Ri}_g\lt 1/4$
somewhere in the flow is often adopted as an indicator of the potential for shear instability based on the classical Miles–Howard theorem (Miles Reference Miles1961; Howard Reference Howard1961), though the formal assumptions required for the theorem to be applicable are invalid even for the simple linear wave state considered here. Nonetheless,
$ \textit{Ri}_g\lt 1/4$
has proven useful as a qualitative indicator of flow instability in a wide range of flows. For a linear internal gravity wave given by the polarisation relations (2.4), it may be shown (e.g. Thorpe Reference Thorpe1999) that, for statically stable configurations
$s\lt 1$
, the minimum value of
$ \textit{Ri}_g$
is greater than
$1/2$
, and so shear instability is not expected. However, when the wave interacts with a horizontal shear flow such that the streamwise velocity field is modified according to (2.12), significant changes to
$ \textit{Ri}_g$
may be expected.
In the limit
$k_x=0$
, the horizontal velocity components
$u$
and
$v$
differ by a phase of
$\pi /2$
. In this case, an analogy with inertia-gravity waves in rotating stratified flows provides an interesting point of comparison. In particular, the polarisation between these two components is the same as a linear inertia-gravity wave propagating in the
$y$
–
$z$
plane under the association
$\partial U/\partial y \leftrightarrow -f$
(importantly however, the same is not true of the dispersion relation). Unlike those in systems with
$f=0$
, such waves may exhibit
$ \textit{Ri}_g\lt 1/4$
for
$s\lt 1$
, with the resulting dynamics closely resembling classical parallel shear instability (Lelong & Dunkerton Reference Lelong and Dunkerton1998a
). The differences between the resulting energetics and those associated with convective instabilities may be significant (Thorpe Reference Thorpe2018). For
$F\neq 0$
, the local gradient Richardson number for the horizontally sheared system considered here is calculated to be
\begin{equation} Ri_g(\boldsymbol{x},t) = \frac {1+s\cos \phi }{F^2s^2 [\cos \phi \sin \beta +(k_x/k_y)(\sin \phi \sin \beta )/R +(\sin \phi \cos \beta )/R]^2}, \end{equation}
where
$\phi = \boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x} - \varOmega t$
is the local wave phase and
$R= (\partial U/\partial y)/\varOmega$
. It is demonstrated in Appendix A that the minimum value of
$ \textit{Ri}_g$
is very close to (and in some cases exactly equal to)
$(1-s)/F^2s^2$
when this quantity is smaller than
$0.5$
. This value is achieved at
$\beta =\pi /2$
and
$\phi = \pi$
, corresponding to perturbations transverse to the direction of wave propagation (cf. Lelong & Dunkerton Reference Lelong and Dunkerton1998a
). It is therefore apparent that
$F$
may be useful as a proxy for shear instability. Indeed, because
$k_z/k_y$
represents a horizontal to vertical wave aspect ratio, the expression (2.19) takes the form of an inverse Richardson number.
The physical meaning of
$F$
breaks down when
$k_y$
becomes small such that the WKBJ approximation is no longer valid. However, for rapid growth in the streamwise velocity field (i.e. large
$\partial U/\partial y$
), the small-amplitude linear approximation may break down before this point. In the special case
$k_x=0$
, Lagrangian conservation of total streamwise momentum limits the maximum amplitude of the streamwise velocity to be at most
$\Delta u$
and the expression for
$F$
is certainly meaningless beyond this limit. With these restrictions in mind, we consider the ray-tracing equations under some representative cases to illustrate the variability in dynamics that can occur due to the influence of horizontal shear. We fix
$k_z=1$
and choose initial values of
$k_x$
,
$k_y$
,
$\Delta u$
and
$h$
to produce scenarios where the wave is trapped (
$k_x\lt 0$
), reflected (
$k_x\gt 0$
) and unaffected by the shear (
$k_x=0$
), consistent with the DNS described in § 3. The ray-tracing equations (2.13) and (2.14) are solved numerically alongside (2.17) governing the evolution of wave action
$\mathcal{A}$
. Local values of the wave steepness
$s$
are computed from
$s^2 = 2\mathcal{A}\varOmega k_z^2$
.
Plots of the wave steepness
$s$
derived from the local wave action
$\mathcal{A}$
(dark blue, left axes) and the parameter
$F$
defined in (2.19) (dark red, right axes) along rays in the
$y$
direction for parameter values corresponding to selected simulations from table 1. Results are shown for
$(a)$
F0.1s2.28(T) and
$(b)$
F1.0s5.0(T) (wave trapping),
$(c)$
F6.0s0.75_kx0 (no refraction) and
$(d)$
F6.0s0.75 (wave reflection). Lines become dashed in panel (c) where the perturbation velocity exceeds
$\Delta u$
and dotted in panel
$(d)$
where
$k_y\lt 0.1$
. Vertical dashed lines indicate a trapping level in panels (a,b) and a turning level in panel
$(d)$
. Horizontal dot-dashed lines represent the value of
$F$
corresponding to a minimum
$ \textit{Ri}_g$
smaller than
$1/4$
(assuming that
$s=0.75)$
.
Values of
$s$
and
$F$
along rays are shown in figure 2 for four representative cases. The first two cases have
$k_x\lt 0$
and both exhibit wave trapping, but differ in the value of
$\Delta u/h$
by approximately an order of magnitude. The third case is the special case
$k_x=0$
(no refraction) and the fourth case has
$k_x\gt 0$
such that wave reflection is predicted. It is important to keep in mind that the WKBJ theory is formally invalid near trapping and turning levels.
We first consider the behaviour of
$s$
. When
$k_x\lt 0$
,
$k_y$
increases monotonically along rays, leading to a decrease in the group velocity as the trapping level at
$y\approx 0$
is approached. This leads to the significant increase in the wave steepness observed in figure 2(a,b) before viscosity acts to damp further energy growth, consistent with the results of Badulin et al. (Reference Badulin, Shrira and Tsimring1985). When
$k_x=0$
, no refraction takes place and wave energy is conserved along rays in the inviscid limit. This corresponds to the approximately constant value of
$s$
seen in panel
$(c)$
, which decreases only due to viscosity. When
$k_x=0.05$
in panel
$(d)$
such that a turning level exists at
$y\approx 0$
,
$s$
similarly remains approximately constant as the turning level is approached, until the linear WKBJ approximation becomes invalid at
$y\approx -2$
. The behaviour of
$F$
is distinctly different to that of
$s$
. When
$k_x\lt 0$
, the growth in
$F$
depends on
$\partial U/\partial y$
. When
$\Delta u/h = 0.125$
is small, as in panel
$(a)$
, wave refraction dominates and
$F$
remains close to
$0$
. For larger
$\Delta u/h=1$
(panel
$b$
),
$F$
is able to grow to
$O(1)$
before growth in
$k_y$
eventually dominates. When
$k_x=0$
(panel
$c$
), energy growth is due entirely to the lift-up mechanism, demonstrated by the increase in
$F$
as the wave approaches the centre of the shear layer. A similar behaviour is observed for
$k_x\gt 0$
(panel
$d$
). Note that, for a given
$k_{x0}$
,
$k_{z0}$
and
$\Delta u/h$
, the maximum value of
$F$
is obtained when
$k_y^2=k_{x0}^2\lt k_{y0}^2$
, explaining why a larger maximum is predicted for the reflected case where
$k_y$
may vary.
2.2.4. Summary
Before continuing, it is worth briefly summarising the use of the ray-tracing theory described above for predicting properties of the local wave state prior to breaking. For a monochromatic wave packet with wavevector
$\boldsymbol{k}_0$
and steepness
$s_0$
initially centred at
$y=y_0$
in the horizontal shear flow
$U(y)$
given by (2.7), the ray-tracing equations (2.13) and (2.14) are solved numerically alongside (2.17) for the wave action
$\mathcal{A}$
. The resulting local values of
$k_y$
and
$\mathcal{A}$
along rays are then used to compute corresponding values of
$s$
and
$F$
describing the relative importance of buoyancy overturns versus enhanced vertical shear in driving instability, as inferred from the local wave state defined by the polarisation relations (2.12).
When
$k_x\lt 0$
and the wave is trapped or refracted strongly by the horizontal shear flow such that
$s$
becomes larger than 1, if
$F\ll 1$
remains small then mixed shear/convective instability associated with a statically unstable local buoyancy profile is expected. If
$F\gtrsim 1$
, wave advection of momentum leads to increased vertical shear that may modify the properties of instability. On the other hand, when
$k_x\geqslant 0$
, trapping does not occur (so that the wave passes through the shear layer or is reflected) and
$s$
remains smaller than
$1$
. In such cases,
$F\gtrsim 1$
may lead to growth of local horizontal velocity perturbations that are sufficient to destabilise the flow to vertical shear instability alone.
3. Numerical simulations
The linear theory in the previous section suggests that a range of wave-breaking behaviours are possible in a horizontal shear layer, depending on the initial flow parameters. However, this analysis is at best heuristic, being strictly valid only for small-amplitude wave perturbations in the WKBJ approximation. To further understand the flow evolution, we perform a series of complementary fully nonlinear DNS intended to span the range of predicted breaking behaviours. Simulations are performed using the multi-parallel pseudo-spectral channel flow solver Diablo (Taylor Reference Taylor2008), which employs a second-order finite difference discretisation in the wall-bounded direction, here taken to be in the direction of horizontal shear. The code uses a third-order mixed implicit/explicit Runge–Kutta/Crank–Nicolson scheme for time marching, with a 2/3 dealiasing rule applied to the nonlinear terms.
3.1. Overview
Under the expectation that
$k_z$
remains constant throughout the wave–shear interaction, it is convenient to take
$K^* = k^*_{z0}$
in the non-dimensionalisation outlined in § 2.1. In this way, the dimensionless initial wavenumber vector
$\boldsymbol{k}_0 = (k_{x0}, k_{y0}, 1)$
has corresponding frequency
$\varOmega _0$
satisfying
$\varOmega _0^2 = (k_{x0}^2+k_{y0}^2)/(1+k_{x0}^2+ k_{y0}^2)$
, and the initial (maximum) wave energy density
$E_0 = s_0^2/(2k_{z0}^2)=s_0^2/2$
. Most of the simulations presented here have
$k_{y0}=0.25$
, with values of
$k_{x0}$
then selected to span a range of predicted trapping and reflecting scenarios. The dimensionless initial frequency
$0.125\lt \varOmega _0\lt 0.38$
(corresponding to propagation angles of
$13^{\circ }$
–
$21^{\circ }$
between
$\boldsymbol{k}_0$
and the horizontal) is large enough such that the neglect of planetary rotation is reasonable (under the assumption that
$f\ll N$
) without being so large that waves quickly become unstable to a parametric instability (Lombard & Riley Reference Lombard and Riley1996b
). We note that these values are within the estimated range of high frequency waves observed in the ocean thermocline by Alford & Pinkel (Reference Alford and Pinkel2000). The packet width scale is chosen as
$a=6\pi$
, with the envelope centred at
$y_0 = -2a$
. Simulations are initialised with
$\boldsymbol{u}(\boldsymbol{x},t=0)=\boldsymbol{U}(y)+\boldsymbol{u}_{{wave}}(\boldsymbol{x})$
and
$b=b_{{wave}}(\boldsymbol{x})$
from (2.6) and (2.7). Additionally, white noise of amplitude
$10^{-5}$
is added to each component of the velocity field in order to trigger instabilities leading to turbulence.
3.2. Inclined computational domain
In order to allow for waves with
$k_{x0}\ll 1$
, it is convenient to perform simulations in an inclined reference frame
$(x,y,z)\to (\hat {x}, y, \hat {z}) \equiv (x\cos \alpha +z\sin \alpha , y, z\cos \alpha - x\sin \alpha )$
obtained by rotation in the
$x$
–
$z$
plane by an angle
$\alpha = -\arctan (k_{x0}/k_{z0})$
. In this frame,
$k_{x0}x+k_{z0}z = \tilde {k}_{z0}\tilde {z}$
, where
$\tilde {k}_{z0} = k_z/\cos \alpha$
. As a result, the initial wave fields
$[\hat {u}, \hat {v}, \hat {w}, \hat {b}]$
are proportional to
$\exp (\mathrm{i}k_{y0} y + \mathrm{i}\tilde {k}_{z0} \hat {z})$
and the initial condition is independent of
$\hat {x}$
. This removes the restriction that the domain length in the
$x$
direction must be an integer multiple of
$2\pi /k_x$
and allows the possibility for arbitrarily small
$k_{x0}$
to be investigated. Further details may be found in Appendix B. A similar approach for studying the instability of plane periodic internal waves is adopted by Fritts & Wang (Reference Fritts and Wang2013), Fritts et al. (Reference Fritts, Bizon, Werne and Meyer2003, Reference Fritts, Wang and Werne2013). As explained by these authors, it is important to note that due to the numerical discretisation and finite domain size, the permitted modes of instability are not the same between the original and inclined coordinate frames. In practice however, provided the largest-scale modes are smaller than the domain size and the resolution is sufficiently high, the subsequent dynamics is expected to be similar between the two. Indeed, in general we find that the most energetic modes of instability in the
$\hat {x}$
direction have wavenumber
$k_x \gtrsim 1$
. Additionally, some preliminary simulations with
$k_{x0}=-1/3$
(not shown) were performed in physical coordinates and compared with calculations performed in the inclined frame, with instability and turbulent breakdown being qualitatively and quantitatively very similar between the two.
The domain length in the periodic
$\hat {x}$
and
$\hat {z}$
directions is taken to be
$L_x=L_z=8\pi$
, whilst a larger extent
$y\in (-100,100)$
is used in the direction of shear to accommodate the initially separated wave packet and shear layer. The grid is evenly spaced in the
$\hat {x}$
and
$\hat {z}$
directions with
$N_x=N_z=648$
grid points in physical space. Polynomial grid stretching is used in the finite difference
$y$
direction so that the finest resolution is achieved at the location
$y=0$
of maximum shear in the centre of the domain and is equivalent to the grid spacing in the
$\hat {x}$
and
$\hat {z}$
directions. A total of
$N_y=1177$
grid points are used in the
$y$
direction, with
$577$
of these points concentrated in the central region
$-12.5 \leqslant y \leqslant 12.5$
. The grid resolution is selected such that grid scale remains at most 2.5 times larger than the Kolmogorov length scale
$L_K = \max (\overline {\varepsilon } Re_w^3)^{-1/4}$
throughout each simulation, where
$\overline {\varepsilon }$
is the kinetic energy dissipation rate averaged over
$\hat {x}$
and
$\hat {z}$
(cf. Kaminski & Smyth Reference Kaminski and Smyth2019). Free-slip/no-flux boundary conditions are applied at the walls
$y=\pm 100$
(
$v=0$
,
$\partial \phi /\partial y=0$
for
$\phi =u,w,b,p$
). A sponge layer is employed in regions of width
$5$
dimensionless units of length adjacent to the boundaries, within which perturbations are damped to zero quadratically in order to prevent wave reflection back into the centre of the domain.
Importantly, though the simulations presented here are performed in the inclined reference frame, output fields shown in what follows have been transformed back into the physical coordinate system.
3.3. Parameter space
Values of
$k_{x0}$
,
$\Delta u$
,
$h$
and
$s_0$
are selected to span a range of wave–shear interaction dynamics. When negative, the value of
$k_{x0}$
is chosen such that a trapping level exists at
$y=y_t=0$
according to (2.15). Analogously, when positive,
$k_{x0}$
is chosen such that a turning level exists at the same location. The horizontal Froude number
$ \textit{Fr}_h = \Delta u/h$
of the initial shear flow (2.7) varies from
$0.125$
to
$1$
. Because the wave extracts kinetic energy from the background shear prior to breaking, a Reynolds number
$Re_I$
more appropriate for describing the resulting turbulence may be defined using the shear flow velocity scale
$\Delta u$
rather than the wave velocity scale. As a result,
\begin{eqnarray} \textit{Re}_I \equiv \frac {\Delta u^*}{k_{z0}^* \nu ^*} = \textit{Re}_w\Delta u \end{eqnarray}
is chosen as
$Re_I=1200$
for all simulations. Most simulations use an initial wave steepness of
$s_0=0.75$
, though, similar to Winters & D’Asaro (Reference Winters and D’Asaro1994) and Howland et al. (Reference Howland, Taylor and Caulfield2021), we also perform some simulations at a smaller amplitude of
$s_0=0.25$
to investigate the effect of wave amplitude on the resulting instability. Additional simulations with
$\Delta u =0$
(not shown) were performed to verify that the propagating initial wave packet remains stable in the absence of the shear flow. Thus, the dynamics of wave instability studied below are instigated by wave–shear interaction and not due solely to parametric self-interaction instability (Lombard & Riley Reference Lombard and Riley1996a
) or the interaction of the wave with its amplitude envelope (Sutherland Reference Sutherland2006b
).
Table of simulations with parameters (as defined in the text) listed, where the naming convention specifies
$F_{\textit{max}}$
and
$s_{\textit{max}}$
, with (T) referring to simulations where the incident wave is trapped. All simulations have
$a=6\pi$
,
$y_0=-2a$
and
$ \textit{Pr}=1$
. Note all simulations have
$k_{y0}=0.25$
except for simulation F6.0s0.75_kx0, which has
$k_{y0}=0.125$
. The shear layer Froude number
$ \textit{Fr}_h$
, wave Reynolds number
$Re_w$
, maximum perturbation energy ratio
$F_{\textit{max}}$
and maximum wave steepness
$s_{\textit{max}}$
are defined in the text.
A list of DNS is displayed in table 1. Included are the maximum values of
$F_{\textit{max}}$
and
$s_{\textit{max}}$
of the parameter
$F$
defined in (2.19) and local steepness
$s$
computed along rays using the ray-tracing equations outlined above. Larger values of
$F_{\textit{max}}$
correspond to an increase in the importance of wave advection of streamwise momentum resulting in excess perturbation kinetic energy, whilst values of
$s_{\textit{max}}\gt 1$
indicate the potential for convective instability. We stress that the precise magnitudes of these parameters are valid only for linear wave perturbations that satisfy the assumptions required for validity of the WKBJ approximation. However, it will be shown that
$F_{\textit{max}}$
remains useful as a predictor for the partition between kinetic and potential energy in the system even for the strongly nonlinear waves considered here. For the purposes of comparison, we note that in our non-dimensionalisation, simulations 1–4 of Staquet & Huerre (Reference Staquet and Huerre2002) have
$\Delta u = 4$
,
$h=4$
,
$k_x=-0.25$
,
$k_y=0.25$
,
$a=4$
,
$y_0=-8$
,
$s_0=0.3$
and
$Re_w\approx 2500$
.
4. Results
We now present results obtained from the DNS outlined in the previous section, starting with a qualitative picture of the wave–shear interaction. For ease of communication, the results are split into simulations with
$k_x\lt 0$
in which the wave is trapped, and those with
$k_x\geqslant 0$
, where the wave either reflects or passes through the shear layer. As detailed above, the dynamics is expected to be characteristically different between these cases.
Snapshots in a vertical plane of the central region
$-15\lt y\lt 15$
taken at successive time points during the flow evolution for simulations (a–d) F0.1s0.76(T), (e–h) F0.1s2.28(T), (i–l) F1.0s5.0(T). Colours illustrate the streamwise vorticity field
$\zeta _x$
. Triangle markers along the horizontal
$y$
axis are located at
$y = \pm h$
, indicating the approximate width of the shear layer.
4.1. Evolution for
$k_x\lt 0$
(wave trapping)
The evolution of the wave packet in the
$(y,z)$
plane is shown in figure 3 by plotting contours of the streamwise vorticity field
$\zeta _x = \partial w/\partial y - \partial v/\partial z$
. Results from three representative simulations with
$k_x\lt 0$
are presented to illustrate similarities and differences in the dynamics between cases with varying
$s_{\textit{max}}$
and
$F_{\textit{max}}$
. Snapshots are taken at times showing the steepening of the incident wave due to refraction in the shear layer, the emergence of instabilities, the subsequent development of small-scale turbulent structures with enhanced vorticity and the eventual decay of these motions. Note that simulations shown in the second and third columns have the same initial parameter values as the results from the linear calculations shown in figures 2
$(a)$
and 2
$(b)$
.
When
$s_0=0.25$
(first column), distortion of the refracted wave packet is visible in figure 2
$(a)$
at time
$t=108$
(approximately 17 dimensionless buoyancy periods). Shortly afterwards, in panel
$(b)$
, small-scale vorticity structures emerge and propagate rapidly across the trapping level. Disordered small-scale motions develop in the centre of the shear layer but are quickly damped by viscosity by time
$t=150$
in panel
$(d)$
at the chosen Reynolds number owing to their relatively small amplitude. Nonetheless, small-scale internal waves, first seen in panel
$(c)$
, are generated at the edge of the central turbulent region and radiate outward in both directions in a manner qualitatively similar to the results of Maffioli et al. (Reference Maffioli, Davidson, Dalziel and Swaminathan2014). For a larger initial wave amplitude
$s_0=0.75$
but keeping all other parameters the same (second column), instability develops much earlier and at a considerable distance from the predicted trapping level. Without appealing to the linear theory, there is no obvious indication as to the location of the trapping level from the vorticity snapshots, with a small amount of wave energy clearly able to penetrate through
$y=0$
. Despite these differences, the subsequent evolution of the instability itself is somewhat similar between the two cases, with the formation of localised vorticity structures that propagate towards and across the centre of the shear layer and the generation of internal waves by turbulence. Likewise, the picture of instability remains similar when the shear layer width is decreased by a factor of four and the magnitude of the velocity is increased by a factor of two (third column), though turbulence is clearly most energetic in this case. The qualitative behaviour in this case is also consistent with results obtained by Staquet & Huerre (Reference Staquet and Huerre2002) for comparable initial values of
$\Delta u$
,
$h$
and
$k_{x0}$
.
Panels (a,d,g) show vertical slices of the perturbation streamwise velocity field for simulations F0.1s0.76(T), F0.1s2.28(T) and F1.0s5.0(T) taken at times corresponding to the panels in the top row of figure 3. Solid black lines are evenly spaced isopycnals (surfaces of constant buoyancy
$b=z+\theta$
). Vertical profiles of the local buoyancy frequency
$N^2 = 1 +\partial \theta /\partial z$
and vertical shear
$S$
are shown in panels (b,c), (e, f) and (h,i). The locations at which these profiles are taken are indicated by the dashed magenta lines. Hatched regions show vertical locations where
$N^2\lt 0$
.
Three-dimensional renderings of the
$b = 6.5$
isopycnal from simulation F0.1s2.28(T) at successive times during the flow evolution. Colours indicate the magnitude of the vertical vorticity field
$\zeta _z$
.
Snapshots showing the local wave state prior to breaking are shown in figure 4. Colours show the perturbation streamwise velocity field
$ u' = u - \overline {u{}}$
, where an overline represents an average in
$x$
and
$z$
. Isopycnals are included to highlight overturning regions. Local vertical profiles of the buoyancy gradient
$N^2(z)$
and vertical shear (of the streamwise velocity)
$S(z) = |\partial u/\partial z|$
taken at the
$y$
location of the maximum velocity amplitude are plotted in the right-hand panels. The results are qualitatively consistent with the predictions from figure 2. In particular, the refraction of the wave by the shear leads to an increase in steepness
$s$
, as seen from the isopycnals. When the initial amplitude is sufficiently large, local regions of overturning develop (hatched regions in the profiles). It is worth noting that the linear theory predicts values of
$s\gt 1$
only for parameter values corresponding to the bottom two rows, which is seen to be consistent with the DNS. Larger values of
$s$
are naturally associated with increased vertical shear. More important for the shear however is the advection of momentum when
$\partial U/\partial y$
is large, indicated by the large values of
$u'$
and corresponding enhanced vertical shear seen in the bottom row of figure 4. This enhanced vertical shear may be associated with the larger value of
$F_{\textit{max}}$
in simulation F1.0s5.0(T). Notably however,
$ \textit{Ri}_g\lt 1/4$
only when
$N^2\lt 0$
: there are no regions of the flow susceptible to pure shear instability. This is also consistent with figure 2(b), which shows that although
$F$
increases, it never becomes large enough to support a corresponding minimum
$ \textit{Ri}_g\lt 1/4$
whilst
$s\lt 1$
.
The three-dimensional behaviour of the wave-breaking instability is shown in figure 5 for simulation F0.1s2.28(T). The development of convectively unstable overturns is clear in panel (a). At low frequencies, convectively unstable waves break via the production of vorticity that emerges in the form of roll-like structures in the direction parallel to wave crests (Howland et al. Reference Howland, Taylor and Caulfield2021). This behaviour has been explained by analysing perturbations to an assumed quasi-steady background state representing the wave at the point of overturning (Winters & Riley Reference Winters and Riley1992; Parker et al. Reference Parker, Howland, Caulfield and Kerswell2021). In general, perturbation energy is extracted from the wave by a mixture of vertical shear production and convection via the buoyancy flux. In this study however, it is less clear that the locally large frequency
$\varOmega \sim N$
expected in the case of wave trapping satisfies the quasi-steady assumption. Indeed, the vortex structures seen in figure 3(e,i) and figure 5(b) appear unique to the horizontal shear scenario. Nonetheless, though the primary instability is fundamentally different to the classical quasi-steady convective wave instability, it is expected that the secondary instability of the streamwise vortices exhibits mixed shear/convection properties (Caulfield & Peltier Reference Caulfield and Peltier2000). The resulting three-dimensionalisation results in the rapid collapse of the convectively unstable region and development of small-scale turbulent structures seen in figure 5(c).
The same three-dimensional rendering as in figure 5 but for simulation F0.1s0.76(T), at two successive time points during the flow evolution.
It is tempting to conclude that turbulent breakdown in the case of wave trapping is ultimately due to a primarily convective instability. However, the case of a smaller-amplitude wave is not consistent with this picture. Figure 6 shows a three-dimensional visualisation of the breaking process for simulation F0.1s0.76(T) with a smaller initial wave amplitude, where it is seen in panel
$(a)$
that a three-dimensional structure develops prior to the overturning of isopycnals. It is shown in § 5 that wave trapping results in the formation of a narrow velocity jet in the streamwise
$x$
direction that is colocated with the region of large isopycnal displacements. Perturbations may then extract energy from the horizontal shear associated with the jet, potentially leading to the generation of vertical vorticity seen in panel
$(a)$
and the subsequent breakdown in panel
$(b)$
. It is worth noting that this behaviour is similar in character to the instability mechanism proposed by Wang & Balmforth (Reference Wang and Balmforth2021) for forced critical layers in a stratified uniform horizontal shear flow.
4.2. Evolution for
$k_x\geqslant 0$
The same as figure 3, for simulations (a–d) F6.0s0.75_kx0, (e–f) F6.0s0.25, (i–l) F6.0s0.75. Note the difference in the colour bar scale.
The time evolution of
$\zeta _x$
for three representative simulations in the case
$k_x\geqslant 0$
is shown in figure 7 (simulations F3.0s0.25 and F3.0s0.75 are qualitatively similar to F6.0s0.25 and F6.0s0.75). Note that the left and right columns have parameter values corresponding to the linear calculations shown in figures 2
$(c)$
and 2
$(d)$
. As expected, the dynamics is fundamentally distinct from the wave trapping scenario explored in the previous section. In the special case
$k_x=0$
shown in the left-hand column of figure 7, the shear flow does not affect the propagation of the wave so that the slope of phase lines is unchanged in
$y$
. However, as will shortly be explained, wave advection of mean momentum across the centre of the shear layer leads to the local generation of vertical shear. The resulting shear instability creates velocity perturbations that correspond to the small-scale vorticity structures that emerge in the centre of the shear layer in panel
$(c)$
. These local shear instabilities lead to thin layers of turbulence that persist even after the initial wave has passed through the shear layer completely. Curiously, there are noticeable variations in the development of shear instabilities at different vertical heights (most clearly seen around
$z\approx 20$
). Vorticity perturbations have a sufficiently large vertical extent such that it seems reasonable to suppose that the instabilities developing in one layer influence the dynamics of adjacent layers. We speculate that interactions between layers could quickly magnify initially small differences in the timing of instability onset in neighbouring layers that occur due to random fluctuations in the background noise field.
When
$k_x$
is strictly positive, refraction of the incoming wave is observed, though phase lines are now tilted towards the horizontal instead of the vertical direction, as seen in figures 7
$(e)$
and 7
$(i)$
. Wave reflection caused by the shear leads to the generation of an interference pattern (seen most clearly in panel e) between the incident and reflected waves, complicating the subsequent dynamics. Note that the linear theory predicts a turning level at
$y=0$
for both cases shown. The flow remains laminar for the smaller-amplitude wave scenario; however, for a larger initial wave amplitude, instabilities start to develop in the region where incident and reflected waves interact (panel
$j$
). Panel (k) shows that subsequent instabilities in the form of classical Kelvin–Helmholtz billows develop towards the centre of the shear layer, leading to greatly enhanced vorticity and eventually fully developed turbulence.
The influence of the advection of background shear by the wave is made abundantly clear by plotting the streamwise perturbation velocity
$u'$
in figure 8. Horizontal layers of positive and negative
$u'$
are seen in all three cases, with shear magnitudes
$O(1)$
and above. In particular, for the large-amplitude initial conditions,
$u'$
may be sufficiently large such that regions with
$0\lt Ri_g\lt 0.25$
develop (top and bottom rows), whose locations correspond to the layers of shear instabilities seen developing in figure 7
$(c{,}k)$
. This behaviour is in accordance with the ray-tracing predictions of figure 2
$(c{,}d)$
. Specifically, the linear theory correctly predicts values of
$F$
large enough such that
$ \textit{Ri}_g\lt 0.25$
in simulations F6.0s0.75 and F6.0s0.75_kx0, whereas for simulation F6.0s0.25 (which does not become turbulent),
$\min (Ri_g)\approx (1-s_{\textit{max}})/F_{\textit{max}}^2s_{\textit{max}}^2 \approx 0.42 \gt 0.25$
. On the other hand, a notable difference from the theory arises from the interference of the incident and reflected waves in simulation F6.0s0.75, leading to narrow regions of local static instability with
$N^2\lt 0$
in locations where the instabilities seen in figure 7
$(j)$
develop. Clearly in this case, convective instabilities develop on a shorter time scale than shear instabilities in the adjacent regions with
$0\lt Ri_g\lt 0.25$
. However, though regions with
$N^2\lt 0$
play a dominant role in the triggering of secondary instabilities, it is shown that the subsequent development of turbulence and mixing remains very much driven by the enhanced vertical shear.
The same as figure 4 for simulations (a–c) F6.0s0.75_kx0, (d–f) F6.0s0.25, (g–i) F6.0s0.75, where snapshots are taken at times corresponding to the top row panels in figure 7. In the right-hand panels, dotted hatching indicates shear-unstable regions where
$0\leqslant Ri_g \lt 1/4$
, and line hatching shows vertical locations where
$N^2\lt 0$
, as before.
The same as figure 5 for simulation F6.0s0.75.
The same as figure 5 for simulation F6.0s0.75_kx0.
A three-dimensional visualisation of simulation F6.0s0.75 is shown in figure 9. The interaction of the reflected and incident waves at the centre of the shear layer leads to the isopycnal overturning seen in panel
$(a)$
. Panel
$(b)$
shows this overturning region spreading out across the shear layer and, in particular, highlights the wave advection of background flow vertical vorticity that leads to the vertical shear generation. At the same time, streak-like
$\zeta _z$
perturbations have emerged in the centre of the overturning region. These instabilities evolve to have billow-like structures seen in panel
$(c)$
closely resembling those of classical shear instability, where vortices are oriented in the horizontal direction at an oblique angle to the
$x$
direction of the background flow. Note that this orientation is not consistent with the prediction
$\beta =\pi /2$
from (2.21), that is, that the most unstable direction is the streamwise direction (and, thus, the billow axis should be in the
$y$
direction). This is likely due to the continuous shearing of the growing perturbations by the background flow. In fact, a similar phenomenon was also reported by Lelong & Dunkerton (Reference Lelong and Dunkerton1998a
) for shear instabilities in convectively stable inertia-gravity waves, where the discrepancy between theory and simulations was attributed to the lack of scale separation between the time scales of wave propagation and the growth rate of instability assumed when appealing to a local gradient Richardson number criterion for instability. The development of small-scale turbulent structures seen in panel
$(c)$
is due to a combination of convective instability of the overturning region seen in panel
$(b)$
and secondary instabilities of the developing Kelvin–Helmholtz billows. As such, the precise realisation of the turbulent transition is likely to be sensitive to
$Re_w$
(and indeed
$ \textit{Pr}$
).
Further evidence for behaviour dominated by vertical shear instability is provided by simulation F6.0s0.75_kx0, which does not exhibit the wave reflection and subsequent isopycnal overturns seen in simulation F6.0s0.75. A three-dimensional picture of developing secondary instabilities for F6.0s0.75_kx0 is shown in figure 10. The dynamics is similar to those in figure 9, with perturbations oriented at an oblique angle to the
$x$
direction. In this case, the angle
$\beta$
is even further from the predicted value of
$\pi /2$
, possibly due to the slower growth rate of instability. As in figure 9, turbulence in figure 10 follows from a subsequent instability of these perturbations. In this case the secondary instabilities seen in figure 10
$(b)$
resemble the instability due to misaligned vortex cores studied by Fritts et al. (Reference Fritts, Wang, Lund and Thorpe2022).
5. Momentum transport
The evolution of the horizontally averaged background flow
$\overline {u}$
is given by
\begin{eqnarray} \frac {\partial \overline {u}}{\partial t} = -\frac {\partial }{\partial y}(\overline {u'v' }) + \frac {1}{\textit{Re}_w}\frac {\partial ^2 \overline {u}}{\partial y^2}. \end{eqnarray}
Note that, neglecting the viscous term and using the wave polarisation relations (2.4),
$ \overline {u'v'}$
is initially negative (positive) when
$k_x\lt 0$
(
$k_x\gt 0$
) and the right-hand side is thus negative (positive) at the leading edge of the wave packet. It follows that the wave packet carries with it a negative mean flow anomaly when
$k_x\lt 0$
and a positive mean flow anomaly when
$k_x\gt 0$
. This corresponds to a unidirectional Bretherton flow derived in Appendix C. It is also convenient to quantify bulk changes in the momentum of the background flow in terms of the shear layer thickness
$\delta _u$
(cf. Staquet & Huerre Reference Staquet and Huerre2002; Pham, Sarkar & Brucker Reference Pham, Sarkar and Brucker2009), defined here as
\begin{eqnarray} \delta _u = \frac {1}{\Delta u^2}\int _{-15}^{15}\Delta u^2 - \overline {u}^2(y)\, {\rm d}y. \end{eqnarray}
Henceforth, for convenience, we denote an integral over
$-15\lt y\lt 15$
with angle brackets
$\langle \boldsymbol{\cdot }\rangle$
. The evolution of
$\delta _u$
is given by
\begin{align} \frac {1}{2}\frac {{\rm d} \delta _u}{{\rm d}t} = -\frac 12 \frac {1}{\Delta u^2}\frac {{\rm d}\langle \overline {u}^2\rangle }{{\rm d}t} = \left \langle \frac {\partial \overline {u}}{\partial y} \overline {u'v'}\right \rangle -\left [ \overline {u}(\overline {u'v'}) - \frac {1}{Re_w}\overline {u}\frac {\partial \overline {u}}{\partial y}\right ]_{-15}^{15} -\frac {1}{Re_w}\left \langle \left (\frac {\partial \overline {u}}{\partial y}\right )^2\right \rangle , \end{align}
where the first and second terms on the right-hand side represent the exchange of energy between the mean flow and perturbations (shear production) and the mean momentum flux into and out of the region. The final two terms are due to diffusion of the mean flow and are negligible at the values of
$Re_w$
considered here.
5.1. Cases with
$k_x\lt 0$
(wave trapping)
The evolution of
$\overline {u}(y)$
relative to the initial shear profile
$U(y)$
is shown in figures 11
$(a)$
and 11
$(b)$
for simulations F0.1s2.28(T) and F1.0s5.0(T). The negative anomaly propagating with the wave packet is clearly seen. As the wave packet is refracted and becomes locally steep in panel
$(a)$
, the anomaly is focused into a narrow jet and increases in amplitude in a manner closely resembling the self-acceleration mechanism described by Sutherland (Reference Sutherland2006b
) for internal waves propagating in a vertical plane. For the narrower shear layer in panel
$(b)$
, a thin jet also develops, though in this case the acceleration does not appear to coincide with a focusing of the wave-induced mean flow but rather emerges within the critical layer located at
$y=0$
prior to the arrival of the wave packet centre. These dynamics are more similar to the situation described by Wang & Balmforth (Reference Wang and Balmforth2020) with a steady small-amplitude wave forcing. The cumulative contribution of the wave-induced mean flow to the shear layer thickness is clear from panels
$(c)$
and
$(d)$
, which reveal the competition between momentum flux from the wave that acts to decrease the shear layer width and the extraction of momentum by shear production that has the opposite effect. In both cases, the development of turbulence after the onset of wave breaking results in a drag on the mean flow leading to a reduction in speed. This drag may be significant and in fact overcome the prior acceleration, resulting in the net increase in the shear layer thickness seen in panel
$(d)$
. Despite the differences in the manner of formation of the jet, the similarities in its structure and the subsequent wave breaking suggests underlying similarities in the dynamics driving instability. A key dynamical feature of the nonlinear critical layer studied by Wang & Balmforth (Reference Wang and Balmforth2020) is the local acceleration of the mean flow by the accumulation of wave pseudo-momentum (see Bühler Reference Bühler2014) within the critical level. Mean profiles of the root-mean-square buoyancy perturbation are shown in figure 12 alongside the adjustment to the mean flow, which show that formation of the mean flow jet defect is associated with a similarly sharp increase in the local buoyancy perturbation amplitude in all cases. For case F0.1s0.76(T) with a small amplitude, local perturbation and mean flow growth are associated with a focusing of the initial wave energy seen in panel
$(a)$
. On the other hand, in case F1.0s5.0(T) it is clear from panel
$(c)$
that amplification occurs rapidly within the predicted critical layer location significantly prior to the arrival of the centre of the wave packet. Despite these differences, the increase in both local mean flow and perturbation energy growth described by Wang & Balmforth (Reference Wang and Balmforth2020) for a steady wave forcing in unbounded stratified horizontal shear flow appears to be quite generic.
Panels (a,b) show time-evolving,
$x$
- and
$z$
-averaged plots of the change in the mean streamwise velocity field from the initial profile
$U(y)$
, for simulations
$(a)$
F0.1s2.28(T),
$(b)$
F1.0s5.0(T). Dotted vertical lines indicate the onset of turbulence. Arrows point to the development of a mean flow jet prior to breaking. Panels (c,d) show corresponding cumulative contributions to the change in the shear layer thickness
$\Delta \delta _u$
given by the (time-integrated) terms in (5.3), averaged over the region indicated by the dashed horizontal lines in the panels above.
Horizontally averaged profiles of the root-mean-square buoyancy perturbation field (top row) and root-mean-square wave-induced mean streamwise flow (bottom row) at successive time points indicated by the dotted to solid lines. Panels (a,d) correspond to simulation F0.1s0.76(T), (b,e) to F0.1s2.28(T), (c, f) to F1.0s5.0(T).
5.2. Cases with
$k_x\geqslant 0$
Figures 13
$(a)$
and 13
$(b)$
show the evolution of
$\overline {u}-U(y)$
for cases with
$k_x=0$
and
$k_x=0.05\gt 0$
. The signature of the Bretherton flow in the latter case is not visible due to the comparatively small magnitude of
$k_x$
. In both cases, significant mean flow deceleration occurs prior to the development of turbulence as momentum is transferred into the streamwise perturbation component, as the local ray-tracing theory suggests. When
$k_x=0$
(panel
$a$
) and the wave passes through the shear layer, some of this energy is transferred back to the mean flow as the wave starts to exit the shear layer at around
$t=50$
, before the development of turbulence at around
$t=60$
. When
$k_x=0.05$
and turbulence is energetic, significant wave perturbations radiate outward from the shear layer and decelerate the flow at considerable distances beyond the central region
$y=\pm 15$
, as seen in panel
$(b)$
. The asymmetry between the regions to the left and right of the shear layer is due to the reflection of the incident wave. Panel
$(c)$
shows that when
$k_x=0$
, changes to the shear layer thickness are driven entirely by shear production. When
$k_x\gt 0$
(panel
$d$
), the wave-induced mean flow, being proportional to
$k_x$
, makes a small but positive contribution to the shear layer thickness. The increase in
$\delta _u$
is eventually attenuated by mean momentum fluxes out of the domain at both
$y=-15$
and
$y=15$
due to outward propagating waves, comprised of the reflected incident wave and waves generated by the interaction with the shear.
The same as figure 11 for simulations F6.0s0.75_kx0 (a,c) and F6.0s0.75 (b,d).
6. Energy transport
We define the total energy density in the region
$-15\lt y\lt 15$
as
$E_{\textit{tot}}= \overline {E}+E'$
, where
\begin{align} \overline {E} = \overline {E}_K+\overline {E}_P=\frac 12 \left \langle \overline {u}_i\overline {u}_i\right \rangle +\frac 12 \big \langle \overline {\theta }^2\big \rangle , \\[-28pt] \nonumber \end{align}
\begin{align} E' = E'_K +E'_P = \frac 12 \big \langle \overline {u_i'u_i'}\big \rangle +\frac 12 \big \langle \overline {\theta '^2}\big \rangle . \\[0pt] \nonumber \end{align}
The time evolution of the mean and perturbation components
$\overline {E}$
and
$E'$
is given by
\begin{align} \frac {\partial \overline {E}}{\partial t} &= -\langle \mathcal{P}\rangle - \langle \mathcal{D}_{mean}\rangle - [T_{mean}]^{15}_{-15}- [ T_{D,mean}]^{15}_{-15}, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial E'}{\partial t} &= \langle \mathcal{P}\rangle - \langle \mathcal{D}\rangle - [T]^{15}_{-15} - [T_D]^{15}_{-15}. \\[0pt] \nonumber \end{align}
Here, the turbulent production
$\mathcal{P}$
represents energy exchanges between the mean flow and perturbations, with dissipation rates for each given by
$\mathcal{D}_{mean}$
and
$\mathcal{D}$
, respectively. Energy fluxes are denoted by
$T$
and are separated into diffusive (subscript
$D$
) and non-diffusive contributions. As in (5.3), at sufficiently high
$Re_w$
diffusive contributions are typically small. The remaining five terms are defined explicitly as
\begin{align} \mathcal{P} = - \overline {u'v'}\frac {\partial \overline {u}}{\partial y} - \overline {w'v'}\frac {\partial \overline {w}}{\partial y} - \overline {\theta 'v'}\frac {\partial \overline {\theta }}{\partial y}, \\[-28pt] \nonumber \end{align}
\begin{align} \mathcal{D} = \varepsilon ' + \chi ', \\[-28pt] \nonumber \end{align}
\begin{align} T = \overline {p'v'} + \frac 12 (\overline {u_i'u_i'v'} + \overline {\theta '\theta 'v'}), \\[-28pt] \nonumber \end{align}
\begin{align} T_{mean} = \frac {1}{2}\overline {u}(\overline {u'v'}) +\frac {1}{2}\overline {w}(\overline {w'v'}), \\[-28pt] \nonumber \end{align}
\begin{align} \mathcal{D}_{mean} = \frac {1}{Re_w}\left [\left (\frac {\partial \overline {u}}{\partial y}\right )^2 + \left (\frac {\partial \overline {w}}{\partial y}\right )^2 \right ] + \frac {1}{\textit{Re}_w \textit{Pr}}\left (\frac {\partial \overline {\theta }}{\partial y}\right )^2, \\[0pt] \nonumber \end{align}
where
$\varepsilon '$
and
$\chi '$
are the perturbation dissipation rates of kinetic and (available) potential energy:
\begin{eqnarray} \varepsilon ' = \frac {1}{Re_w}\left (\frac {\partial u'_i}{\partial x_j}\frac {\partial u'_i}{\partial x_j}\right ), \quad \chi ' = \frac {1}{Re_w Pr}\left (\frac {\partial \theta ' }{\partial x_j}\frac {\partial \theta ' }{\partial x_j}\right ). \end{eqnarray}
Panels (a–d) show the evolution of the relative total, mean and perturbation energies alongside the cumulative energy dissipation, whilst panels (e–h) show the cumulative perturbation energy fluxes into and out of the region
$-15\lt y\lt 15$
for the following simulations: (a,e) F0.1s2.28(T), (b,f) F1.0s5.0(T), (c,g) F6.0s0.75_kx0 and (d,h) F6.0s0.75. Dashed lines in the bottom panels show the cumulative energy fluxes computed from simulations with equivalent initial wave parameters but without the background shear flow. All quantities have been normalised by the initial wave energy.
The evolutions of
$ E_{\textit{tot}}$
,
$\overline {E}$
and
$E'$
within the region
$-15\lt y\lt 15$
are plotted in figure 14 alongside the cumulative energy loss due to viscous dissipation
$ \int _0^t \langle \mathcal{D}+ \mathcal{D}_{mean}\rangle {\rm d}t$
. The parameter values for the simulations shown correspond to the cases in figure 2. For simulations F0.1s2.28(T) and F1.0s5.0(T) with
$k_x\lt 0$
shown in figure 14(a,b), the total energy increases as the wave propagates into the shear layer, with positive contributions from both mean and perturbation energy. This is consistent with the expected local acceleration due to the wave-induced mean flow described above. The perturbation energy growth reaches a peak at the point where the wave starts to break, which is followed by a rapid loss of mean flow and perturbation energy due to turbulent dissipation. For weaker shear (panel
$a$
), the cumulative dissipation over the wave-breaking event approaches a value approximately equal to the maximum energy gain due to the wave–shear interaction, which is approximately twice the initial wave energy. However, for stronger shear (panel
$b$
), the cumulative dissipation is 7.5 times larger than the initial wave energy. This greatly exceeds the earlier energy growth, suggesting that the turbulent state in this case might more appropriately be thought of as being associated with a subcritical transition of the otherwise linearly stable horizontal shear flow, rather than simply balancing the total increase in energy due to the wave at the point prior to breaking. In particular, the inertial instability suggested by Staquet & Huerre (Reference Staquet and Huerre2002) (in which rotation provides a means for additional perturbation energy to be supplied to turbulence) does not appear to be necessary to drive significantly enhanced dissipation.
The behaviour in simulations F6.0s0.75 and F6.0s0.75_kx0, shown in figures 14
$(c)$
and 14
$(d)$
, is quite different, with wave advection of momentum leading to an immediate and dramatic increase in perturbation energy sourced from the mean flow. Once again, the peak in energy growth corresponds to wave breaking and the onset of enhanced turbulent dissipation. In simulation F6.0s0.75_kx0 (panel
$c$
), the perturbation energy starts to decrease as the wave packet exits the shear layer. While some of this energy is returned to the mean flow, a steady but significant decrease in total energy occurs due to turbulent dissipation associated with vertical shear instability. In contrast, for simulation F6.0s0.75 (panel
$d$
), there is a much more rapid decrease in the total energy due to enhanced dissipation, indicating a more transient burst of turbulence. The total energy dissipation is significant in both scenarios at between 5 and 10 times the initial wave energy.
A significant flux of energy out of the region
$-15\lt y\lt 15$
may be expected when the primary wave either passes through or is reflected, or when the wave–shear interaction and/or resulting turbulence excite further waves that radiate outwards. These effects are shown in figure 14(e–h). To provide reference values for wave transmission, four additional simulations were performed with the same wave parameters but no background shear flow (cf. Winters & D’Asaro Reference Winters and D’Asaro1994), indicated by the dashed lines in panels (e–h). Cumulative energy fluxes through the
$y=\pm 15$
boundaries from each simulation where the background shear is present are compared with those from the reference no-shear cases. Note that in all cases the total flux of energy is smaller than the initial wave energy due to losses associated with (non-turbulent) viscous dissipation and dispersive effects, the former of which dominates. In simulations F0.1s2.28(T) and F1.0s5.0(T), wave trapping initially prevents the expected transmission of energy, evident from the solid blue lines in panels
$(e)$
and
$(f)$
. It is notable however that significant generation of internal waves in F1.0s5.0(T) results in outgoing wave radiation, thereby enabling energy transmission. The decrease in the cumulative energy flux is approximately the same through the left and right boundaries suggesting wave generation is roughly symmetric. The cumulative outgoing energy flux (approximately twice the net flux through
$y=15$
) is almost equal to the energy of the initial wave packet and, comparing to panel
$(b)$
, approximately 15 % of the total energy dissipated at
$t=165$
(see also figure 16 and the accompanying discussion below). In simulation F6.0s0.75_kx0 (panel
$g$
), the wave energy transmission through the shear layer very closely matches the case with no shear, indicating that despite the significant total dissipation, there is very little energy loss from the incident wave and little to no generation of internal waves by turbulence. In simulation F6.0s0.75 (panel
$h$
), the reflection of the initial wave is clear, with a negative cumulative energy flux through
$y=-15$
. A substantial generation of waves also occurs in this case, adding to the reflected outgoing wave energy and providing significant energy flux through
$y=15$
. Excluding the reflected incident wave, the total outgoing energy is approximately 20 % of the total dissipation at
$t=165$
.
6.1. Energy partition
The ray-tracing predictions in § 2 suggest that the interaction between the incident wave and horizontal shear can, at least locally, lead to characteristic differences in the partition between perturbation kinetic and potential energy, and more specifically, between the individual components of kinetic energy. Plots of the total perturbation energy
$ E'$
and its partition into separate components (6.2) are shown in figure 15 where, once again, panels correspond to the cases from figure 2. Despite the limitations of the analysis describing local wave energy densities, the bulk energy partition is remarkably consistent with the predictions from the theory. Recall that
$F$
represents the excess energy in the streamwise velocity perturbation component
$u'$
due to wave advection of momentum, with the maximum instantaneous value
$F_{\textit{max}}$
from the ray-tracing theory used to characterise simulations in table 1. It is immediately clear that perturbation energy growth is dominated by
$ (1/2) \langle u'^2\rangle$
in panels
$(c)$
and
$(d)$
, for which
$F$
is predicted to increase significantly without corresponding growth in the wave steepness (and, hence, the other components of the total energy). On the other hand, wave steepness growth is seen to dominate proceedings in panels
$(a)$
and
$(b)$
, as is clear from looking at, for example, the larger contribution of the potential energy component
$ (1/2) \langle \theta '^2\rangle$
to the total perturbation energy. However, as was suggested by the growth of
$F$
in figure 2
$(b)$
and the snapshots in figure 4, advection of momentum may still play a significant role in the instability dynamics in such cases. This prediction is further confirmed here by observing that the (maximum) relative magnitude of
$ (1/2) \langle u'^2\rangle$
is significantly larger in simulation F1.0s5.0(T) with
$F_{\textit{max}}=1$
than in F0.1s2.28(T) with
$F_{\textit{max}}=0.1$
. It is important to keep in mind that the linear theory only describes the dynamics up to the point of wave breaking (signified by the peaks in perturbation energy in figure 15).
Stacked area plots illustrating the split of the total perturbation energy in the region
$-15\lt y\lt 15$
into contributions from individual velocity and buoyancy components for simulations
$(a)$
F0.1s2.28(T),
$(b)$
F1.0s5.0(T),
$(c)$
F6.0s0.75_kx0 and
$(d)$
F6.0s0.75.
A more detailed picture of the exchange of energy between kinetic and potential reservoirs arises from the evolution equations for
$E'_K$
and
$E'_P$
, given by
\begin{align} \frac {\partial E'_K}{\partial t} &= \langle \mathcal{P}_{\mathcal{S}} \rangle + \langle \mathcal{B}\rangle - \langle \varepsilon '\rangle - \left [\overline {p'v'} + \frac 12 (\overline {u_i'u_i'v'})\right ]_{-15}^{15} , \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial E'_P}{\partial t} &= -\langle \mathcal{B}\rangle - \langle \chi '\rangle -\left \langle \overline {\theta 'v'}\frac {\partial \overline {\theta }}{\partial y}\right \rangle - \left [\frac 12(\overline {\theta '\theta 'v'})\right ]^{15}_{-15}, \\[0pt] \nonumber \end{align}
where diffusive fluxes have been omitted for brevity. Here the buoyancy flux
$ \mathcal{B}\equiv \overline {\theta 'w'}$
is the rate of conversion of potential to kinetic energy and the shear production
$\mathcal{P}_S\equiv -\overline {u'v'}({\partial \overline {u}}/{\partial y}) - \overline {w'v'}({\partial \overline {w}}/{\partial y})$
is the rate of extraction of kinetic energy from the background flow. The cumulative evolution of these terms is shown in figure 16 alongside the cumulative kinetic and available potential energy dissipation, where terms are normalised by the total cumulative turbulent dissipation late in the simulation at
$t=165$
. In the wave trapping cases shown in panels
$(a)$
and
$(b)$
, a significant fraction of the initial energy perturbation growth through the production term is transferred to potential energy via the buoyancy flux. This fraction is slightly reduced in panel
$(b)$
when
$F$
is larger. On the other hand, panels
$(c)$
and
$(d)$
confirm the significantly smaller fraction of kinetic to potential energy transfer obvious from figure 15.
A significant result from figure 16 is that the cumulative dissipation rates of turbulent kinetic and potential energy depend sensitively on the energetics of the wave prior to breaking. Put simply, small-scale turbulent statistics apparently inherit properties of the larger-scale breaking wave. For example, in panels
$(c)$
and
$(d)$
where wave advection of momentum and, hence, kinetic energy dominate at the time of breaking, the cumulative dissipation at the end of the mixing event is dominated by
$\varepsilon$
. On the other hand, when kinetic and potential energy is more equally partitioned at the time of breaking, as in panels
$(a)$
and
$(b)$
,
$\chi$
and
$\varepsilon$
contribute approximately equally to the total dissipation. Note that
$\chi$
is often interpreted as a diapycnal mixing rate (see, e.g. Howland et al. Reference Howland, Taylor and Caulfield2021), with a corresponding (cumulative) mixing efficiency given by
$\eta (t) \equiv \int _0^t \langle \chi '\rangle {\rm d}t /\int _0^t \langle \mathcal{D} \rangle {\rm d}t$
. An important observed feature of the horizontally sheared system is that
$\eta$
varies considerably between wave-breaking scenarios, asymptoting to values from
$0.08$
in panel (c) to
$0.42$
in panel (a).
Time-integrated terms from (6.11) and (6.12) for simulations (a) F0.1s2.28(T), (b) F1.0s5.0(T), (c) F6.0s0.75_kx0 and (d) F6.0s0.75. The black dashed lines in panels (b,d), for which there is significant generation of internal waves, show an estimate for the cumulative (total) energy loss due to this emission obtained from doubling the outgoing energy flux through
$y=15$
(i.e. assuming symmetric emission). All quantities have been normalised by the cumulative total energy dissipation
$\int _0^{165} \langle \mathcal{D} \rangle {\rm d}t$
.
7. Discussion and conclusions
In this study we investigated the problem of an internal gravity wave packet propagating into a horizontal shear layer. When the wave vector forms an oblique angle with the direction of the background flow (
$k_x\neq 0$
), Doppler shifting leads to local changes in the intrinsic frequency measured by a stationary observer and gradients in the flow velocity refract the wave, leading to changes in the local steepness. When
$k_x=0$
, wave propagation is unaffected by the shear layer, however significant local perturbation energy changes are possible due to advection of mean flow momentum (equivalent to the lift-up mechanism), as pointed out by Bakas & Farrell (Reference Bakas and Farrell2009a
). In general, the parameter space describing the initial condition is extremely large; thus, a reduced-order model capturing the essential dynamical variability is desirable. Here, we constructed such a model by studying locally wave-like solutions to the ray-tracing equations under the WKBJ approximation. To provide a simple estimate from the initial conditions for the subsequent relative importance of wave advection of momentum, we derived the perturbation energy ratio
$F$
, where large
$F$
corresponds to an excess of kinetic over potential energy due to wave advection of mean flow momentum.
The value of
$F$
at a given time during the wave–shear interaction is inversely proportional to the square root of the minimum gradient Richardson number and, thus, at least when the latter is positive, serves as an indicator of the susceptibility of the flow to vertical shear instability. Indeed, fully nonlinear DNS with large
$F$
and no increase in wave steepness exhibited wave-breaking dynamics driven by the emergence of billow-like vortices closely resembling those associated with a classical Kelvin–Helmholtz instability (figures 7 and 9). Notably, this picture is also consistent with simulations of plane periodic inertia-gravity waves that, unlike in the non-rotating case, may satisfy
$0\lt Ri_g\lt 0.25$
locally and exhibit similar breaking dynamics to those observed here (Lelong & Dunkerton Reference Lelong and Dunkerton1998a
). On the other hand, when
$F$
remains small but wave refraction leads to a significant increase in the local wave steepness, large overturns in the density field are expected to play a dominant role in turbulent transition (convective breaking). This prediction is broadly consistent with the dynamics observed in simulations in which waves were trapped by the shear near critical levels, though questions remain regarding the role of the significant local acceleration of the mean streamwise flow in wave breaking and, in particular, the development of three-dimensional instabilities. The development of this flow appears to be driven by the same fundamental dynamics governing local energy growth in baroclinic critical layers under a steady wave forcing, as described by Wang & Balmforth (Reference Wang and Balmforth2020).
Predicted ratio of local maximum kinetic to potential energy
$1+F_{\textit{max}}^2$
from ray-tracing theory versus the value obtained from DNS. Dark blue symbols correspond to simulations with an initial wave amplitude of
$s_0=0.75$
(where a triangle marks the special case
$k_x=0$
) and red symbols an initial wave amplitude of
$s_0=0.25$
.
The effectiveness of the linear theory for describing the local wave state prior to breaking is summarised in figure 17. In general, there is a strong positive correlation between the DNS data and the predicted local maximum ratio of kinetic to potential energy density. The agreement, though far from exact, is quite satisfactory given the assumptions of linearity and scale separation
$k_y h \gg 1$
are generally invalid. At the very least, the use of the theory for describing qualitatively the essential dynamical variability in the wave–mean flow interaction is justified. Unsurprisingly, better agreement is found for simulations with a smaller initial wave amplitude. The under-prediction for larger values of
$F_{\textit{max}}$
when
$k_x\gt 0$
in cases where
$s_0=0.75$
(compared with smaller values of
$s_0$
and the special case
$k_x=0$
) is a product of the nonlinear interaction between incident and reflected waves leading to a local increase in the wave amplitude and corresponding overturning. It is worth noting here that the surprising utility of the ray-tracing theory in capturing essential dynamics of wave breaking was also noted for the problem of internal waves in a vertical shear by Howland et al. (Reference Howland, Taylor and Caulfield2021).
It is important to note that, whilst we have shown that linear ray-tracing theory provides reasonable predictions for the dynamics of instability, energy growth may additionally result from mechanisms not captured under the WKBJ approximation. Firstly, the WKBJ theory is invalid when
$k_y h \ll 1$
so that the assumption of scale separation between the waves and the background flow no longer holds. In particular, this occurs near turning levels where
$k_y\to 0$
. In any case, the theory is invalid even in the linear regime, however, (linear) transient growth may still occur. These dynamics are studied formally in Arratia (Reference Arratia2011), where the formation of layers with strong vertical shear for
$k_x\ll 1$
is also observed. The mechanism for producing these layers, referred to as stratified lift-up, comes from wave advection of mean flow momentum, which also supplies the excess vertical shear predicted in the ray-tracing theory. Secondly, near trapping levels
$k_y\to \infty$
, so the scale separation assumption becomes increasingly valid, however, the ray-tracing theory instead becomes invalid due to growing nonlinearities as
$s$
increases. Indeed, as discussed in § 5.1, evidence for nonlinear adjustments to the mean flow (through wave self-acceleration or within linear critical levels) has been observed here and warrants attention in future studies.
In almost all of the DNS, energy perturbation growth eventually led to wave breaking, turbulent dissipation and mixing. As summarised in figure 14, the dynamics of turbulent breakdown and the resulting energy exchanges between the mean flow, incident wave and turbulence depend sensitively on the route to wave breaking, and hence, on
$F$
. In both wave trapping and wave reflection scenarios, a large fraction of the total dissipation takes place during a transient burst of turbulence, much like the behaviour observed for wave breaking in a vertical shear flow (Winters & D’Asaro Reference Winters and D’Asaro1994; Howland et al. Reference Howland, Taylor and Caulfield2021). In notable contrast however, for the horizontal shear flow studied here, larger values of
$F$
are associated with a total dissipation that may be several times larger than the energy of the primary wave, thus requiring a significant input of energy to turbulence from the background flow, whose linear stability is enforced here. The corresponding horizontal mixing of momentum results in significant deceleration of the background shear flow, as discussed in § 5. This may be compared with Winters & D’Asaro (Reference Winters and D’Asaro1994) for example, who found that the total dissipation is only around 20 %–25 % of the initial wave energy in a vertical shear with a minimum gradient Richardson number equal to
$0.5$
.
In some cases, the wave–shear interaction led to the generation of internal waves whose total energy was seen to be up to 25 % of the total energy dissipated by turbulence and more than twice that contained in the initial wave packet (figure 16). We note that even in the absence of turbulent breaking, strong generation of internal waves for
$F=O(1)$
was predicted by Bakas & Farrell (Reference Bakas and Farrell2009b
), who found that up to 30 % of the maximum wave energy within the shear layer may be radiated away. On the other hand, localised decaying initially homogeneous and isotropic turbulence can also produce significant internal wave emission in a background stratification. Using a combination of laboratory experiments and DNS, Maffioli et al. (Reference Maffioli, Davidson, Dalziel and Swaminathan2014) found the total outgoing wave energy to be 13 % of the cumulative dissipation for a horizontally localised patch of turbulence. To our knowledge, the substantial emission of internal waves by turbulence has not been reported for wave breaking at a vertical critical level (this being differentiated from the often-observed phenomenon of partial reflection most obviously by the symmetry of wave emission between each side the critical level). However, for turbulence produced in an unstable vertical shear layer in the presence of a strongly stratified region below, we can infer from Pham et al. (Reference Pham, Sarkar and Brucker2009) that internal wave emission constituted up to 22 % of the total energy dissipated.
An analysis of the turbulent kinetic and potential energy budgets demonstrated that the exchanges between (perturbation) kinetic and potential energy in the lead up to wave breaking had a substantial influence on the subsequent turbulent state, most notably being closely related to the cumulative dissipation rates of kinetic and available potential energy (figure 16). As a result, waves that broke more convectively, here characterised by approximate equipartition between turbulent kinetic and potential energy at the point of breaking, exhibited corresponding values of the (cumulative) mixing efficiency up to
$0.42$
, considerably larger than the typical value of
$0.17$
assumed in the literature. When instability was driven by vertical shear, mixing efficiencies as low as
$0.08$
were observed. In this study, a relatively small number of representative simulations were studied in detail to highlight the influence of the breaking mechanisms on the subsequent turbulent dynamics. The results naturally motivate further simulations spanning a wider range of (at least)
$F$
,
$\Delta u/h$
and
$Re_w$
to investigate in more detail the mixing variability and its dependence on turbulent parameters in addition to the mechanisms driven wave breaking.
Nonetheless, the close relationship between the partition of kinetic and potential energy (which is dominated by larger-scale motions) and mixing efficiency is generally consistent with results obtained from strongly stratified forced systems for which turbulence production is maintained against dissipation in a statistically stationary state by continual power input at large scales (see, for example, Maffioli, Brethouwer & Lindborg Reference Maffioli, Brethouwer and Lindborg2016). However, the precise relationship to the mixing efficiency is sensitive to other dimensionless turbulent parameters and, importantly, the type of forcing used. Despite broad similarities in the anisotropic structure of turbulence and the shapes of turbulent kinetic and potential energy spectra, considerable differences in mixing properties are observed between forcing concentrated in vertically rotational modes with
$k_z=0$
, forcing by internal waves and forcing with a mean vertical shear (Lindborg & Brethouwer Reference Lindborg and Brethouwer2007; Portwood, de Bruyn Kops & Caulfield Reference Portwood, de Bruyn Kops and Caulfield2019; Howland, Taylor & Caulfield Reference Howland, Taylor and Caulfield2020). The results of this study indicate that a wide range of possible local transient wave-breaking dynamics may contribute towards the considerable mixing variability seen in forced systems, which could in theory collectively comprise many such local wave–mean flow interactions depending on the type of forcing. Further investigation is required in order to link robustly the local properties of wave breaking studied here to their collective influence in complex multiscale systems, where the distinction between the largest scales of turbulence and the smallest scales of internal waves is unclear.
Throughout this study, we have neglected the influence of planetary rotation that will undoubtedly be important when the time scales of either the shear layer or the internal wave are sufficiently long, that is, when
$\varOmega \sim O(f)$
or
$\partial U/\partial y\sim O(f)$
. The latter condition was studied for
$N\sim f$
by Staquet & Huerre (Reference Staquet and Huerre2002), who suggested that in the case of wave trapping, inertial instability may act as an additional source of turbulence energy. The former condition may be particularly relevant for scenarios where wave advection of momentum dominates energy perturbation growth, which is favoured for small horizontal wavenumbers
$k_y,k_x \ll 1$
and, hence, small
$\varOmega /N$
. In this case, the dynamics in the presence of rotation may be considerably more complicated, as is suggested by the recent results of Hilditch, Taylor & Thomas (Reference Hilditch, Taylor and Thomas2025). The effect of changing the Reynolds number
$Re_w$
will undoubtedly be important during wave breaking and the turbulence that follows, especially at the modest values considered here. In particular, the extents to which the critical layer dynamics (e.g. Wang & Balmforth Reference Wang and Balmforth2021), secondary instabilities (e.g. Mashayek & Peltier Reference Mashayek and Peltier2013) and the dissipation and mixing rates shown in figure 16 are affected by changes in
$Re_w$
remain open problems. A further parameter that continues to present a considerable computational challenge in DNS is the Prandtl number
$ \textit{Pr}$
that we have assumed to be unity, but in reality is expected to be
$O(10)$
in thermally stratified ocean environments.
Finally, to what extent the variability in wave-breaking dynamics explored here extends to a more general class of background flows – which might include cyclostrophically balanced flows or, for non-zero
$f$
, rotationally balanced flows (see, e.g. Edwards & Staquet Reference Edwards and Staquet2005; Whitt & Thomas Reference Whitt and Thomas2013) – remains an open question. Scenarios such as these, in addition to the problem studied here, may help to provide a crucial link between wave breaking and the asymptotic regime of anistropic stratified turbulence that has been proposed as a universal description for the dynamics of turbulent patches in the ocean interior (Riley & Lindborg Reference Riley and Lindborg2008; Kunze Reference Kunze2019; D’Asaro Reference D’Asaro2022). As D’Asaro (Reference D’Asaro2022) points out, the initial condition for the classical problem of a (horizontally infinite) wave packet incident on a vertical shear flow lacks potential vorticity, thought to be a key ingredient in the energy cascade of strongly stratified turbulent flows (Lindborg Reference Lindborg2006). The inclusion of a horizontally varying mean flow provides a possible source of potential vorticity at large scales. It would be worthwhile to investigate whether the turbulence produced in the wave-breaking scenarios studied here can access this strongly stratified regime in the case when the mean flow is linearly stable (the unstable case is discussed in detail in Lewin & Caulfield Reference Lewin and Caulfield2024). It is also worth recognising that, in general (and as was the case in this study), horizontally bounded wave packets carry with them their own source of potential vorticity in the form of a Lagrangian vorticity dipole (Bühler & McIntyre Reference Bühler and McIntyre2005; Bretherton Reference Bretherton1969). This may, as suggested by D’Asaro (Reference D’Asaro2022) and by the focusing of the mean flow shown in figure 12, play an important role in the resulting energy cascade to small scales during breaking.
Acknowledgements
This work was initiated at the Stanford Center for Turbulence Research Summer Program 2024. We thank the organisers and participants for their efforts in facilitating an engaging and stimulating research environment. Discussions with Leif Thomas, Jeff Koseff and Bruce Sutherland are gratefully acknowledged.
Funding
MMPC acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC) through RGPIN-2024–06184. SFL and AKK were supported by the Ho-Shang and Mei-Li Lee Faculty Fellowship in Mechanical Engineering at UC Berkeley. AB acknowledges support from a National Science Foundation Graduate Research Fellowship (Grant Number DGE-1656518). This research used the Savio computational cluster resource provided by the Berkeley Research Computing programme at the University of California, Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor for Research, and Chief Information Officer).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Minimum gradient Richardson number
After some algebra, equation (2.21) for the gradient Richardson number can be written as
\begin{eqnarray} Ri_g = \frac {R^2}{F^2s^2}\frac {1+s\cos \phi }{\big [\big (R\cos \phi + (k_x/k_y)\sin \phi \big )^2 + \sin ^2\phi \big ]\sin ^2(\beta +\gamma )}, \end{eqnarray}
where
$\tan \gamma = \sin \phi /(R\cos \phi + k_x\sin \phi /k_y)$
. Thus,
\begin{eqnarray} \min _{\phi , \beta } (Ri_g) = \frac {R^2}{F^2s^2}\min _{\phi } \left (\frac {1+s\cos \phi }{\left (R\cos \phi +(k_x/k_y)\sin \phi \right )^2+\sin ^2\phi }\right ). \end{eqnarray}
In the following we assume that
$s\lt 1$
so that
$ \textit{Ri}_g\gt 0$
. Consider first the case
$k_x=0$
. For
$R\geqslant 1$
, the denominator and numerator of the function in brackets are respectively maximised and minimised at the same value of
$\phi =\pi$
, making this a global minimiser. For
$R\lt 1$
, let
$c=\cos \phi$
, then stationary points of the function in the brackets above are given by the roots of a quadratic equation:
\begin{eqnarray} c_{min}=\frac {1}{s}\left (-1 \pm \sqrt {1-s^2/(1-R^2)}\right ). \end{eqnarray}
There are two cases to consider. When
$1-R^2\lt s^2$
, the discriminant in the expression above is negative and there are no stationary points in the range
$-1\lt c\lt 1$
. The global minimum thus occurs at
$c=-1 \iff \phi =\pi$
. On the other hand, when
$1-R^2\gt s^2$
, there is a single global minimum in the range
$-1\lt c\lt 1$
at
$c=c_{min}$
. Substituting this into (A2) we obtain
\begin{align} \min (Ri_g)= \begin{cases} \dfrac {R^2}{2F^2s^2} \left ( \sqrt {1-\dfrac {s^2}{1-R^2}}+1 \right ), & R^2\leqslant 1-s^2,\\[2ex] \dfrac {1-s}{F^2s^2}, & R^2\gt 1-s^2.\\ \end{cases} \end{align}
Note, in particular, that
$R^2/F^2\gt 1 =$
so that, when
$R^2\leqslant 1-s^2$
,
$\min (Ri_g)\gt 1/2$
for
$s\lt 1$
.
When
$k_x\neq 0$
, the denominator of (A2) is not maximised at
$\phi =\pi /2$
for any value of
$R$
. However, the value at
$\phi =\pi /2$
will be similar to the maximum value when
$Rk_y/k_x$
is large. Because
$Rk_y/k_x = F k_h |\boldsymbol{k}|/(k_xk_z)\gt F$
, we can argue heuristically that this may be anticipated near locations of wave instability. Where dynamics are dominated by wave steepening near trapping levels,
$k_h|\boldsymbol{k}|/(k_xk_z) \gg 1$
. On the other hand, where vertical shear produced by wave advection of momentum dominates,
$F$
is necessarily large. Figure 18 demonstrates that, indeed,
$(1-s)/F^2s^2$
becomes an increasingly good approximation for
$\min (Ri_g)$
as the latter decreases towards zero. Note that, in particular, the point at which
$\min (Ri_g)$
becomes smaller than
$1/4$
is very well estimated using
$(1-s)/F^2s^2$
in all cases.
Evolution of the minimum gradient Richardson number along rays for the parameter values considered in figure 2, computed by evaluating (A2) numerically. Dotted curves show the estimate
$\min (Ri_g)=(1-s)/F^2s^2$
. The heuristic local stability threshold
$\min (Ri_g)=1/4$
is also indicated.
Appendix B. Inclined coordinates
A schematic of the inclined coordinate system introduced in § 3.2 is shown in figure 19. Under rotation by an angle
$\alpha = -\arctan (k_{x0}/k_{z0})$
in the
$x$
–
$z$
plane, the velocities in the
$\tilde {x}$
and
$\tilde {z}$
directions are
$(\hat {u}, \hat {w}) = (u\cos \alpha + w\sin \alpha , -u\sin \alpha +w\cos \alpha$
). The equations of motion are
\begin{align} \displaystyle \frac {\mathrm{D}\hat {\boldsymbol{u}}}{\mathrm{D}t} = \boldsymbol{\nabla }\hat {p} + \hat {\theta }(\hat {\boldsymbol{e}}_z\cos \alpha + \hat {\boldsymbol{e}}_x\sin \alpha ) + \frac {1}{Re_w} {\nabla} ^2\hat {\boldsymbol{u}}, \\[-28pt] \nonumber \end{align}
\begin{align} \displaystyle \boldsymbol{\nabla }\boldsymbol{\cdot }\hat {\boldsymbol{u}} = 0, \\[-28pt] \nonumber \end{align}
\begin{align} \displaystyle \frac {\mathrm{D}\hat {\theta }}{\mathrm{D}t} = -\hat {w}\cos \alpha - \hat {u}\sin \alpha + \frac {1}{Re_wPr}{\nabla} ^2 \hat {\theta }, \\[0pt] \nonumber \end{align}
with initial conditions
\begin{align} \hat {\boldsymbol{u}}(\hat {\boldsymbol{x}},t=0) = \Delta u \cos \alpha \tanh (y/h) \hat {\boldsymbol{e}}_x - \Delta u \sin \alpha \tanh (y/h) \hat {\boldsymbol{e}}_z +\hat {\boldsymbol{u}}_{{wave}}, \\[-28pt] \nonumber \end{align}
\begin{align} \hat {\theta }(\hat {\boldsymbol{x}},t=0)=\hat {\theta }_{{wave}}, \\[0pt] \nonumber \end{align}
where
$\hat {\boldsymbol{u}}_{{wave}}$
and
$\hat {\theta }_{{wave}}$
are the plane wave solutions from (5) that have been transformed into the inclined reference frame and expressed in terms of
$\hat {\boldsymbol{x}}$
.
Schematic illustrating the coordinate transformation used in simulations. A translucent slice of the initial vertical velocity field is shown to demonstrate that wave phase lines in the
$\hat {x}$
–
$\hat {z}$
plane are parallel to the
$\hat {x}$
direction.
Appendix C. Wave-induced mean flow
As originally demonstrated by Bretherton (Reference Bretherton1969), the stress exerted by an internal wave packet on the fluid through which it propagates induces a mean flow at second order in the wave amplitude. The nature of this mean flow depends on the geometry of the wave packet, in particular, the aspect ratio of its envelope. Van den Bremer & Sutherland (Reference Van den Bremer and Sutherland2018) provide a general expression for the wave-induced mean flow in a wave packet with arbitrary amplitude envelope
$A(x,y,z,t)$
, writing the result as the sum of the horizontal Bretherton flow forced by vertical vorticity and a long wave flow forced by horizontal vorticity. Because the wave packets we consider here are unbounded in the vertical, the long wave flow is zero at leading order. Extending the analysis of Van den Bremer & Sutherland (Reference Van den Bremer and Sutherland2018) to a three-dimensional wave packet, horizontal divergence of the momentum flux drives a mean flow
$\boldsymbol{u}_{\textit{DF}}$
given by
\begin{eqnarray} \boldsymbol{u}_{\textit{DF}} = u_{\textit{DF}}\left (1,\frac {c_{gy}}{c_{gx}},\frac {c_{gz}}{c_{gx}}\right ), \end{eqnarray}
where
$\boldsymbol{c}_g = (c_{gx}, c_{gy}, c_{gz})$
is the linear group velocity of the waves and
\begin{eqnarray} u_{\textit{DF}} =\frac {k_x}{\omega } \frac {N^2}{2k_z^2}s^2A^2, \end{eqnarray}
taking the form of the pseudo-momentum per unit mass in the
$x$
direction. This so-called divergent-flux induced flow must induce a secondary response flow to ensure incompressibility; the sum of these two flows gives the Bretherton flow
$\boldsymbol{u}_{BF}$
. As shown in Van den Bremer & Sutherland (Reference Van den Bremer and Sutherland2018),
$w_{BF}=0$
and we must have additionally that
$(\boldsymbol{\nabla }\times \boldsymbol{u}_{\textit{DF}})\boldsymbol{\cdot }\hat {\boldsymbol{k}} = (\boldsymbol{\nabla }\times \boldsymbol{u}_{BF})\boldsymbol{\cdot }\hat {\boldsymbol{k}}$
. In combination with the incompressibility condition
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{BF}=0$
, it follows that, in Fourier space with coordinates
$(\kappa ,\lambda ,\mu )$
,
\begin{eqnarray} \hat {\boldsymbol{u}}_{BF}= \hat {u}_{\textit{DF}}\begin{bmatrix} [\lambda ^2 - (k_y/k_x)\mu \kappa ]/\kappa _H^2\\[4pt] [(k_y/k_x)\kappa ^2 - \kappa \lambda ]/\kappa _H^2 \\[4pt] 0 \end{bmatrix}, \end{eqnarray}
where a hat denotes the three-dimensional Fourier transform and
$\kappa _H^2 = \kappa ^2+\lambda ^2$
. Note that, with
$A \sim \exp (-y^2/a^2)$
, we have
$\hat {u}_{\textit{DF}} = k_xN^2/(2\omega k_z^2)\mathcal{A}(\lambda )\delta (\kappa )\delta (\mu )$
, where
$\mathcal{A}(\lambda ) = ({1}/{2\pi })\int _{-\infty }^{\infty } A^2(y)\exp (-i\lambda y)\ {\rm d}y$
is the Fourier transform of
$A^2(y)$
in the y direction. Substituting this into (C3), it is straightforward to inverse Fourier transform and find that
\begin{eqnarray} \boldsymbol{u}_{BF}= \frac {k_x}{\omega }\frac {N^2}{2k_z^2}s^2\begin{bmatrix} A^2 \\ 0 \\ 0 \end{bmatrix}={u}_{\textit{DF}}\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \end{eqnarray}
that is, the Bretherton flow is a uniform flow in the
$x$
direction with the amplitude given by the pseudo-momentum per unit mass.
