1. Introduction
The transition to turbulence from Kelvin–Helmholtz (KH) instability in stably stratified, parallel shear flows is a long-standing open problem in fluid dynamics with broad geophysical and astrophysical relevance. Classically, this transition is assumed to occur after normal mode perturbations grow to finite amplitude, producing quasi two-dimensional vortex structures (referred to as billows) that subsequently collapse, resulting in the rapid onset of three-dimensional motions. Importantly, the properties of subsequent turbulence and mixing are sensitive to the details of the transition (Mashayek & Peltier Reference Mashayek and Peltier2013; Caulfield Reference Caulfield2021).
At high Reynolds number, the finite-amplitude billow structure is subject to a wide range of secondary instabilities (Mashayek & Peltier Reference Mashayek and Peltier2012a , Reference Mashayek and Peltierb ) that may be grouped broadly into three classes. The first class consists of vortex pairing instabilities (e.g. Klaassen & Peltier Reference Klaassen and Peltier1989) occurring at subharmonics of the most unstable primary streamwise wavenumber. The second class extracts their energy from unstable density gradients in the billow core and are referred to as secondary convective instabilities (e.g. Klaassen & Peltier Reference Klaassen and Peltier1985). The third class instead extracts their energy from the narrow, inclined strips of vorticity connecting adjacent billows, referred to as braids (Corcos & Sherman Reference Corcos and Sherman1976; Staquet Reference Staquet1995; Smyth Reference Smyth2003). Though a range of braid instabilities are possible (Mashayek & Peltier Reference Mashayek and Peltier2012a , Reference Mashayek and Peltierb ), this study focuses specifically on secondary shear instability (SSI), which produces trains of vortices that are advected along the braid and towards the cores. Turbulent transition is governed by a competition between the growth rates of these instabilities, which change substantially as the primary billow grows and saturates in amplitude.
In the ocean, KH instability is considered a primary pathway between submesoscales and stratified turbulence (Smyth & Moum Reference Smyth and Moum2012), with coherent billow structures observed from the deep ocean to the thermocline, continental shelves and estuaries. In some cases, detailed acoustic, velocity and density observations (Geyer et al. Reference Geyer, Lavery, Scully and Trowbridge2010; Vladoiu et al. Reference Vladoiu2025; Lefauve et al. Reference Lefauve, Bassett, Plotnick, Lavery and Geyer2026) show intense turbulent mixing concentrated in the braids, with large overturns rarely present. A compilation of observations from the past four decades (Appendix A of Lefauve et al. Reference Lefauve, Bassett, Plotnick, Lavery and Geyer2026) suggests that such braid-dominated mixing is widespread in continuously forced shear flows across stratification types (temperature, salinity and turbidity interfaces) and shear-layer thicknesses spanning
${\approx} 0.5$
–
$200$
m, corresponding to Reynolds numbers
$\approx 10^5$
–
$10^8$
. Providing a mechanistic and predictive basis for the hypothesis of Geyer et al. (Reference Geyer, Lavery, Scully and Trowbridge2010) that mixing is braid-dominated at these Reynolds numbers, the present study clarifies why and when SSI is expected to control the transition to three-dimensional turbulence, thereby shaping subsequent mixing across a range of ocean flows.
Seminal early work on SSI by Corcos & Sherman (Reference Corcos and Sherman1976) analysed the local flow in a tilted, braid-aligned coordinate system, assuming the resulting instabilities to be of the classical inflectional shear type. Subsequently, Staquet (Reference Staquet1995) drew on the stability of a strip of vorticity under strain (Dritschel et al. Reference Dritschel, Haynes, Juckes and Shepherd1991) to propose a modified braid stability criterion based on the ratio of strain to shear, later adopted by Smyth (Reference Smyth2003) and Mashayek & Peltier (Reference Mashayek and Peltier2012a , Reference Mashayek and Peltierb ). Importantly, these later studies assumed either a steady or frozen-in-time background flow following saturation of the primary billow. Here, we relax this assumption and explicitly model the time-evolving steepening of the braid. This distinction is essential because, as we will show, at geophysically relevant parameters SSI can develop well before billow saturation and must therefore be predicted from the evolving braid.
To this end, we introduce an inviscid model for the time evolution and stability of the braid in § 2, followed by two-dimensional numerical simulations in § 3. In § 4, we present results for the onset of SSI across a range of Richardson and Reynolds numbers and examine the agreement with the inviscid model. Finally, we conclude in § 5.
2. Theory
2.1. Governing equations and parameters
We consider a two-dimensional, stably stratified, vertically sheared layer evolving in time
$t$
with streamwise coordinate
$x$
and vertical coordinate
$z$
. The flow is incompressible, with small density variations under the Boussinesq approximation. The dimensionless governing equations and initial conditions are
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol u =0,\qquad \frac {{\rm D} \boldsymbol u}{{\rm D}t} = -\boldsymbol{\nabla }p+ {Ri} \,b\, \boldsymbol{ \hat {z}} + \frac {1}{\textit{Re}} {\nabla} ^2 \boldsymbol u, \qquad \frac {{\rm D}b}{{\rm D}t} = \frac {1}{ \textit{Re}\textit{Pr}} {\nabla} ^2 b, \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol u(z,\; t=0)= \tanh (z)\,\boldsymbol{\hat x}, \qquad b(z,\;t=0) = \tanh (z) .\\[0pt] \nonumber \end{align}
All lengths are scaled by the shear layer half-depth
$d$
, whilst velocities and buoyancy are scaled with half the dimensional velocity and buoyancy differences across the shear layer
$u_0$
and
$b_0$
. Equations (2.1) and (2.2) include three dimensionless numbers, the Reynolds number
$ \textit{Re}$
, the Richardson number
$ \textit{Ri}$
and the Prandtl number
$ \textit{Pr}$
:
\begin{align} \textit{Re} = \frac {u_0 \,d}{\nu }, \qquad \textit{Ri} =\frac { b_0\,d}{u_0^2}, \qquad \textit{Pr}=\frac {\nu }{\kappa }=1, \end{align}
where
$\nu$
and
$\kappa$
are the kinematic viscosity and buoyancy diffusivity, respectively.
2.2. Inviscid model of braid evolution
Following Corcos & Sherman (Reference Corcos and Sherman1976) (hereafter CS76), we begin by defining a braid-aligned coordinate system tilted at an average angle
$\phi (t)$
to the horizontal. The along-braid variation of the angle is insignificant (as we will see). The along-braid coordinate
$x'$
and cross-braid coordinate
$z'$
are defined by
\begin{align} x' = x\cos (\phi )+z\sin (\phi ), \qquad z' = -x\sin (\phi )+z\cos (\phi ), \end{align}
where the inflection point of the braid is located at the origin
$(x,z)=(x',z')=(0,0)$
. Neglecting the time-dependent rotation of the coordinate system, the governing equations (2.1) in braid coordinates can be written as
\begin{align} \frac {\partial b}{\partial t} &= -u'\frac {\partial b}{\partial x'}-w'\frac {\partial b}{\partial z'} + \frac {1}{\textit{Re}\textit{Pr}}{\nabla} ^{\prime 2}b, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \varOmega }{\partial t}&=-u'\frac {\partial \varOmega }{\partial x'}-w'\frac {\partial \varOmega }{\partial z'}-{Ri} \cos \phi \frac {\partial b}{\partial x'}+{Ri} \sin \phi \frac {\partial b}{\partial z'}+\frac {1}{{Re}}{\nabla} ^{\prime 2}\varOmega , \\[0pt] \nonumber \end{align}
where
$\varOmega =\partial u'/\partial z'-\partial w'/\partial x'$
is the vorticity, and
$u'$
and
$w'$
refer to the braid-aligned velocities. In this long and thin braid flow, we assume gradients in
$z'$
are greater than those in
$x'$
, thus neglecting the
$\partial /\partial x'$
terms, and identify vorticity as nearly equal to the cross-braid shear. The baroclinic amplification of shear is due to the term
$ \textit{Ri} \sin \phi \, (\partial b/\partial z')$
in (2.6). We further assume a uniform compressive strain field
$w' = -\gamma z'$
local to the braid (to be validated later by simulations), where
$\gamma (t)$
is the local strain rate. These assumptions are particularly accurate far away from the cores and reduce (2.5) and (2.6) to
\begin{align} \frac {\partial b}{\partial t} \approx \gamma z'\frac {\partial b}{\partial z'} + \frac {1}{\textit{Re}\textit{Pr}}\frac {\partial ^2 b}{\partial z^{\prime 2}}, \qquad \frac {\partial \varOmega }{\partial t} \approx \gamma z'\frac {\partial \varOmega }{\partial z'} + {Ri} \sin \phi \frac {\partial b}{\partial z'} + \frac {1}{{Re}}\frac {\partial ^2 \varOmega }{\partial z^{\prime 2}}. \end{align}
We seek similarity solutions
\begin{align} b(z^{\prime},t) = F(\xi _b),\quad \varOmega (z^{\prime},t) = S(t) \, G(\xi _\varOmega ), \end{align}
where
$\xi _b = z'/\delta _b(t)$
and
$\xi _\varOmega =z'/\delta _\varOmega (t)$
for time-dependent buoyancy and velocity layer thicknesses
$\delta _b$
and
$\delta _\varOmega$
. The homogeneous vertical shear
$S(t)$
in the expression for
$\varOmega$
accounts for the baroclinic acceleration on the right-hand side of (2.7b
). We further consider the inviscid limit
$ \textit{Re} \to \infty$
in (2.7a–b
), expected to hold as long as
$\delta _b$
and
$\delta _\varOmega$
are sufficiently large compared to the length scales at which diffusion of buoyancy and velocity are
$O(1)$
. The above similarity solutions satisfy (2.7a–b
) provided that (i)
$ \textit{Re} \to \infty$
; (ii)
$G=\partial F/\partial \xi$
, where
$\xi \equiv \xi _b = \xi _\varOmega$
or equivalently
$\delta (t)\equiv \delta _b(t)=\delta _\varOmega (t)$
, valid with initial conditions (2.2); and (iii)
$\delta (t)$
and
$S(t)$
satisfy
\begin{align} \begin{aligned} \frac {{\rm d} \delta }{{\rm d} t}=-\gamma \delta , \qquad \frac {{\rm d} S}{{\rm d} t}= \frac {{Ri} \sin \phi }{\delta }. \end{aligned} \end{align}
Letting
$t=0$
in (2.8), the function
$F$
is given by
$b(z^{\prime},0) = F(z^{\prime}/\delta (0))$
, which corresponds to
$\delta (0) = 1/\partial b/\partial z^{\prime}|_{z^{\prime}=0}$
and
$\partial F/\partial \xi (0) = 1$
. In practice,
$b(z^{\prime},0)$
is the sum of the initial condition and the normal mode perturbation. We will show that (2.9a–b
) accurately predicts the diagnostic quantities governing braid stability. Further details on the validity of the similarity-solution assumptions are given in supplementary material §§ 1 and 2.
For simplicity, here we assume
$\gamma (t)=\gamma _1t+\gamma _0$
and
$\phi (t)=\phi _1t+\phi _0$
during the initial period of billow roll-up (to be validated later by simulations), where
$\gamma _0$
,
$\gamma _1$
,
$\phi _0$
and
$\phi _1$
are anticipated to be empirical functions of
$ \textit{Ri}$
and of the initial conditions. This approach may be compared with CS76, who derived an inviscid model for the advection of vorticity from the braid into the cores that was then solved numerically to compute
$\gamma (t)$
and
$\phi (t)$
. Note that, for comparison with our simulations, the initial angle
$\phi _0$
and strain
$\gamma _0$
are both small but non-zero. The evolution of braid thickness
$\delta (t)$
is then given by
\begin{align} \delta (t) = \delta (0)\exp \left [-\int _0^t \gamma (\tau ) \;{\rm d}\tau \right ] = \delta (0)\exp \left [- \frac 1 2 \gamma _1 t^2 -\gamma _0 t\right ]\!, \end{align}
whilst the evolution of braid shear
$S(t)$
is given by
\begin{align} S(t) = S(0) + \int _0^{t}{Ri} \,\sin \left (\phi _1 \tau + \phi _0 \right ) \,\delta (0)^{-1}\exp \left [\frac 1 2 \gamma _1 \tau ^2 +\gamma _0 \tau \right ] \! \; {\rm d} \tau . \end{align}
2.3. Braid stability criterion
CS76 viewed the braid as a steady parallel stratified shear layer, where a necessary criterion for instability is that the minimum gradient Richardson number at the centre is
\begin{align} {Ri}_{g}(0,t) = {Ri} \frac { (\partial b/\partial z') \cos \phi }{\varOmega ^2} \Big |_{z'=0} = \frac {{Ri}\cos \phi }{\delta (t) S(t)^2}\lt 1/4. \end{align}
They argued that, under this condition, SSI may develop, assuming that the timescale of the local braid flow (which dictates the roll-up time of subsequent shear instabilities) is much shorter than that of the background shear flow. In our notation, this corresponds to
$S\gg 1$
. However, identifying a threshold value is difficult as the assumptions underlying standard linear stability analysis break down in the braid. In particular, the strain field acts to damp growing perturbations.
To address this limitation, Staquet (Reference Staquet1995) used results from the study of a strip of vorticity under uniform straining by Dritschel et al. (Reference Dritschel, Haynes, Juckes and Shepherd1991) to argue that a stricter criterion for SSI onset should involve the ratio of strain to vorticity (shear):
\begin{align} \frac {\gamma }{\varOmega }\Big |_{z'=0} = \frac {\gamma _1t +\gamma _0}{S(t)} \lt \varGamma ^*. \end{align}
In this picture, SSI is in fact the product of rapid transient (non-normal) growth that can reach significant amplitude before the strain is able to damp perturbations. Critical values of
$\varGamma ^*$
between
$0.02$
and
$0.1$
have been suggested (Staquet Reference Staquet1995; Smyth Reference Smyth2003; Mashayek & Peltier Reference Mashayek and Peltier2012a
). Importantly, however, these studies neglected the time-evolving braid and considered the stability of the flow at a later time when braid thickness
$\delta$
is maintained in an equilibrium in which straining is balanced by viscous diffusion. The validity of this approach was verified at the time for simulations with Reynolds numbers up to
$O(10^4)$
. The viscous braid equilibrium occurs approximately at time
$t_{2D}$
, defined here as the time at which the perturbation kinetic energy first reaches a local maximum in the absence of SSI. This time also corresponds approximately to saturation of the primary billow and the onset of a plateau in
$\gamma (t)$
. We will show that the threshold (2.13), marking the dominance of baroclinic shear over strain, can be reached on a timescale considerably shorter than
$t_{2D}$
, thereby allowing early SSI onset.
3. Numerical simulations
3.1. Setup
We run direct numerical simulations (DNS) using Oceananigans.jl, a GPU-accelerated finite-volume solver (Wagner et al. Reference Wagner2025) with a centred second-order advection scheme and constant isotropic diffusivities. We present a set of 16 simulations with geophysically relevant
$ \textit{Ri} = [0.05, 0.10, 0.15, 0.20]$
and
$ \textit{Re} = [10^4,10^5, 10^6, 10^7]$
. Simulations are limited to
$ \textit{Re}\ge 10^4$
as lower values do not exhibit early SSI onset.
The flow is horizontally periodic in the
$x$
domain
$[-L_x/2,\; L_x/2]$
and is subject to free-slip (
$\partial u/\partial z=0$
) and no-flux (
$w=0$
and
$\partial b/\partial z=0$
) boundary conditions at
$z =-L_z/2,\;L_z/2$
. The length
$L_x$
is chosen to be the wavelength of the fastest-growing linear mode for each
$({Ri}, {Re})$
– typically
$L_x\approx 14$
– thus restricting the domain to one billow length and excluding the pairing instability (Smyth Reference Smyth2003). The domain height
$L_z=20$
is chosen sufficiently large to minimise boundary effects on instability growth. A fixed stretched grid was used in
$z$
to increase resolution in the central region where gradients are sharpest, with grid spacing at
$z=0$
five times finer than the average. To initiate KH instability growth, the flow is perturbed with the fastest-growing linear mode corresponding to each
$({Ri}, {Re})$
with energy
$10^{-4}$
of the base flow. All cases were run on a grid of
$8000\times 8000$
points, which was sufficient to ensure convergence (verified in the supplementary material, § 3).
Figure 1 shows the evolution of a sample simulation, illustrating the initial condition (figure 1
a), primary billow growth (figure 1
b) and subsequent onset of SSI on the braid (figure 1
c). Figure 1(d) shows the subsequent advection of secondary instabilities towards the cores and the ultimate breakdown of the primary and secondary billows. These dynamics are beyond the scope of this study as they will clearly be dominated by three-dimensional instabilities (Mashayek & Peltier Reference Mashayek and Peltier2012a
,
Reference Mashayek and Peltierb
). Therefore, simulations are considered up to the time of primary billow saturation
$t_{2D}$
or the time of SSI onset
$t_S$
, whichever occurs first – the valid 2-D regime. Note that
$t_S$
is diagnosed as the time at which
$ \textit{Ri}_{g}$
reaches a minimum, as SSI growth causes a subsequent increase.
An example of the numerical braid analysis for a simulation with
$({Ri},{Re})=(0.15,10^6)$
, showing (a–d) the transient evolution of the buoyancy field
$b$
. The initial profiles for
$b$
,
$u$
and
$w$
are shown in (a). The white dashed lines in (b) and (c) illustrate the location of the braid, and a white box shows the local braid field. The local buoyancy field
$b(x',z')$
is shown in (e) and (f). (e) also overlays the
$x'$
-averaged profiles
$\overline b(z')$
,
$\overline u'(z')$
and
$\overline w'(z')$
. The black dashed line in (e) shows the linear approximation for
$\overline w'$
used to diagnose the strain
$\gamma (t)$
.
Figure 1 Long description
Panel A: A contour plot shows the initial profiles for buoyancy (b), horizontal velocity (u), and vertical velocity (w) at time t = 0. The x-axis ranges from -5 to 5, and the z-axis ranges from -5 to 5. Panel B: A contour plot displays the buoyancy field at time t = 30. The x-axis ranges from -5 to 5, and the z-axis ranges from -5 to 5. White dashed lines illustrate the location of the braid, and a white box shows the local braid field. Panel C: A contour plot shows the buoyancy field at time t = 51. The x-axis ranges from -5 to 5, and the z-axis ranges from -5 to 5. White dashed lines illustrate the location of the braid, and a white box shows the local braid field. Panel D: A contour plot displays the buoyancy field at time t = 60. The x-axis ranges from -5 to 5, and the z-axis ranges from -5 to 5. Panel E: A contour plot shows the local buoyancy field b(x’, z’) with overlaid x’-averaged profiles for buoyancy (b̄), horizontal velocity (ū’), and vertical velocity (w̄’). The x’-axis ranges from -1 to 1, and the z’-axis ranges from -1 to 1. A black dashed line shows the linear approximation for w̄’ used to diagnose the strain γ(t). Panel F: A contour plot displays the local buoyancy field b(x’, z’) at a later time. The x’-axis ranges from -1 to 1, and the z’-axis ranges from -1 to 1.
3.2. Numerical braid diagnostics
To allow for comparison between the predictions of the inviscid braid evolution model, the simulation outputs were analysed as follows. First, the initial linear mode perturbation was positioned so that the braid inflection point is always located at the origin
$(x,z)=(0,0)$
. Second, we defined the braid following the fitting procedure outlined in Smyth (Reference Smyth2003), fitting a cubic polynomial to the points of strongest shear, as seen by white dashed lines in figure 1(b–c). Third, we calculated the average braid angle
$\phi$
from the gradient of the braid locus. Fourth, transects taken along the braid inclined at angle
$\phi$
reconstructed the buoyancy
$b$
, along-braid shearing velocity
$u'$
and cross-braid straining velocity
$w'$
. Figures 1(e) and 1(f) illustrate
$b$
where the braid is now centred around
$z'=0$
. Before the growth of secondary instabilities, the flow is approximately parallel along
$x'$
, and so we average along
$x'$
to construct averaged profiles of buoyancy
$\overline b(z')$
, and velocities
$\overline u'(z')$
and
$\overline w'(z')$
; such profiles are overlaid in figure 1(e). The braid is similar to the initial base flow in figure 1(a) but is subject to additional straining. These averaged profiles are used to calculate
$ \textit{Ri}_{g}$
and
$\gamma /\varOmega$
, as defined in (2.12) and (2.13), respectively. Finally, a linear fit
$\overline w'=-\gamma z'$
(figure 1
e) validates our earlier assumption of uniform strain and provides the time-evolving strain
$\gamma (t)$
to be used in the inviscid model.
4. Results
4.1. Overview and SSI onset in
$({Ri},{Re})$
space
Figure 2 presents an overview of the full parameter range using snapshots of the buoyancy field at
$t_{2D}$
, where
$t_{2D}$
is calculated for each
$ \textit{Ri}$
using the corresponding value for
$ \textit{Re}=10^4$
(for which SSI does not develop). Qualitatively, figure 2 confirms the well-known result that as
$ \textit{Ri}$
increases (from left to right), the amplitude of the primary billows decreases (Thorpe Reference Thorpe1973). As stratification increases with
$ \textit{Ri}$
, greater work is required to raise the denser fluid in the billow, resulting in flatter and weaker cores. The role of increasing
$ \textit{Re}$
(from top to bottom) is less apparent but clearly
$ \textit{Ri}$
-dependent; at
$ \textit{Ri}=0.05$
, the billows appear identical across all
$ \textit{Re}$
since there is no onset of SSI. However, above this
$ \textit{Ri}$
, increasing
$ \textit{Re}$
favours SSI. At higher
$ \textit{Ri}$
, SSI occurs at a lower
$ \textit{Re}$
. Figure 2 shows that high
$ \textit{Re}$
and
$ \textit{Ri}$
are generally favourable to the onset of SSI before primary billow saturation (
$t_S\lt t_{2D}$
), motivating the use of a predictive instability criterion.
Snapshots of the buoyancy field
$b$
at the time of primary billow saturation,
$t_{2D}$
, over a range of
$ \textit{Re}$
(rows) and
$ \textit{Ri}$
(columns). The values of
$t_{2D}$
are
$[33.8, 49.6, 61.3, 106.8]$
and the horizontal domain
$L_x$
is
$[14.3,14.3,14.0,13.4]$
for
$ \textit{Ri}=$
$[0.05, 0.10, 0.15, 0.20 ]$
, respectively. Symbols in the upper-right corners represent the stability outcome of the braid before
$t_{2D}$
; stars denote SSI, crosses denote marginal SSI (instigated by additional noise, not shown here) and circles denote the absence of SSI.
The billow growth is quantified in figure 3, where the growth rates and maximal values of both
$\phi$
(a proxy for core amplitude) and
$\gamma$
(a proxy for core rotation rate) decrease with
$ \textit{Ri}$
. This billow growth is
$ \textit{Ri}$
-dependent and essentially
$ \textit{Re}$
-independent at
$ \textit{Re}\ge 10^4$
. When SSI occurs, higher
$ \textit{Re}$
results in an earlier SSI onset (smaller
$t_S$
). The red dashed lines in figure 3 represent the linear fits
$\phi = \phi _1 t+ \phi _0$
and
$\gamma =\gamma _1 t +\gamma _0$
used as inputs in the inviscid model (exact values in supplementary material, § 4). The slopes
$\phi _1$
and
$\gamma _1$
(given in red) are found to decrease approximately linearly with
$ \textit{Ri}$
.
In the following sections, we assess the stability of the braid by considering the evolution of
$ \textit{Ri}_{g}$
and
$\gamma /\varOmega$
in each simulation (§ 4.2). We consider three key timescales: the primary billow saturation time
$t_{2D}$
, the SSI onset time
$t_S$
and the time
$t^*$
at which the instability criterion (2.13) is first satisfied. The prediction of
$t^*$
from the inviscid model is considered in § 4.3, where we will establish the relevant range of
$ \textit{Re}$
for which the inviscid model and its constituent assumptions hold.
The temporal evolution of braid angle
$\phi$
and strain
$\gamma$
at
$ \textit{Re}=[10^4,10^5,10^6,10^7]$
, as given by the legend, and
$ \textit{Ri}=[0.05,0.10,0.15,0.20]$
, as given by column titles. The primary billow saturation is marked at
$t_{2D}$
by a circle (see corresponding snapshots in figure 2), and SSI onset is marked at
$t_S$
by a star. Red dashed lines and slopes
$\phi _1,\gamma _1$
show the linear fits to be used in the inviscid model (see figure 4).
Figure 3 Long description
Eight line graphs depict the temporal evolution of braid angle (phi) and strain (gamma) at different Reynolds numbers (Re) and Richardson numbers (Ri). The graphs are organized into two rows and four columns. The top row (a-d) shows the evolution of braid angle (phi) over time (t) for different Ri values: 0.05, 0.10, 0.15, and 0.20. The bottom row (e-h) shows the evolution of strain (gamma) over time (t) for the same Ri values. Each graph includes data for Re values of 10^4, 10^5, 10^6, and 10^7, as indicated by the legend. The primary billow saturation is marked by a circle at t_2D, and the onset of secondary shear instability (SSI) is marked by a star at t_S. Red dashed lines and slopes (phi_1, gamma_1) represent linear fits used in the inviscid model.
4.2. Braid evolution and stability
The evolution of the minimum gradient Richardson number
$ \textit{Ri}_{g}(0,t)$
calculated using the middle expression in (2.12) is shown in figure 4 (top row). Early on, the braid is not dynamically distinct from the primary billow and
$ \textit{Ri}_{g}$
is effectively independent of
$ \textit{Re}$
. During this stage, straining thins the braid, increasing local stratification and raising
$ \textit{Ri}_{g}$
, effectively stabilising the flow. However, once shear begins to intensify, this trend reverses and
$ \textit{Ri}_{g}$
decreases rapidly, at which point differences between
$ \textit{Re}$
cases become apparent.
The temporal evolution of the minimum gradient Richardson number
$ \textit{Ri}_{g}$
and ratio of strain to shear
$\gamma /\varOmega$
. Symbols indicate the braid stability outcome as in figure 2 (black crosses indicate SSI onset with the addition of white noise). The red dashed curves are the inviscid model predictions for
$ \textit{Ri}_{g}$
and
$\gamma /\varOmega$
, respectively. In the lower row, horizontal red dashed lines show the threshold
$\varGamma ^* =0.02$
, while vertical red dashed lines indicate the time
$t^*$
at which the inviscid model attains this threshold value. Vertical black dotted lines indicate the time
$t_{2D}$
. Note that the separation
$t_{2D}-t^*$
increases with
$ \textit{Ri}$
for
$ \textit{Ri}\ge 0.10$
.
While
$ \textit{Ri}_{g}$
is informative, it is insufficient to diagnose SSI onset as it remains below
$1/4$
throughout, and there is no other apparent useful critical value. All simulations attain similarly low values of
$ \textit{Ri}_{g}\lesssim 0.05$
before
$t_{2D}$
, emphasising the weakly stratified nature of the braid. This motivates the comparison to the stability of a strained strip of vorticity and the choice of a stability criterion based on
$\gamma /\varOmega$
, whose evolution is shown in the bottom row of figure 4. Again, all growth is initially
$ \textit{Re}$
-independent; the initial increase of
$\gamma /\varOmega$
is approximately linear and governed solely by the increasing strain from the primary billow. As baroclinic acceleration then dominates,
$\gamma /\varOmega$
decreases. The onset of SSI occurs below the threshold
$\varGamma ^*=0.02$
across the parameter space studied here.
In our stability model, early SSI onset at
$t_S$
requires more than satisfying the threshold
$\varGamma ^*$
: the time interval between threshold crossing and primary-billow saturation must also be long enough for transiently growing perturbations to reach finite amplitude. Thus, early SSI requires a sufficiently large window
$t_{2D}-t^*$
. This window increases with
$ \textit{Ri}$
, because
$t_{2D}$
grows with
$ \textit{Ri}$
while
$t^*$
remains nearly unchanged (see the lower row of figure 4). The relative independence of
$t^*$
with
$ \textit{Ri}$
in the range
$ \textit{Ri}=0.1{-}0.2$
can be understood from (2.9a–b
). Although the baroclinic shear-production term contains an explicit factor
$ \textit{Ri}$
, higher
$ \textit{Ri}$
also implicitly slows the angle growth
$\phi (t)$
and the strain-driven thinning of the braid
$\delta (t)$
(through
$\gamma (t)$
), both of which are needed for shear to intensify. As a result,
$\varOmega (t)$
evolves, counterintuitively, slightly more slowly with increasing
$ \textit{Ri}$
, approximately offsetting the simultaneous slowing of
$\gamma$
, such that
$\gamma /\varOmega$
reaches
$\varGamma ^*$
at comparable
$t^*$
across
$ \textit{Ri}=0.1$
–
$0.2$
.
To further confirm the transient growth interpretation in the
$t^*\lt t\lt t_{2D}$
window, all simulations that reached
$\varGamma ^*$
but did not exhibit SSI were restarted
$\approx 10$
time units before
$t^*$
with the addition of divergence-free white noise with modest energy
$10^{-4}$
of the initial flow. This was sufficient to trigger SSI in two marginal cases, marked with crosses on figures 2 and 4:
$({Ri},{Re})=(0.10,10^5)$
and
$(0.20,10^4)$
. By contrast, this additional noise did not trigger visible SSI prior to
$t_{2D}$
at the lower values of
$({Ri}, {Re}) = (0.10,10^4)$
and
$(0.15,10^4)$
, where transient growth did not have sufficient time to reach finite amplitude.
4.3. Agreement with the inviscid model at finite
$ \textit{Re}$
Finally, figure 4 compares the simulations at increasing
$ \textit{Re}$
with the inviscid model predictions (red dashed curves) for
$ \textit{Ri}_{g}$
(the right-hand side of (2.12)) and
$\gamma /\varOmega$
(from (2.13), calculated using the linear fits for
$\phi ,\gamma$
). As expected, the agreement improves with increasing
$ \textit{Re}$
. The inviscid predictions for
$t^*$
(vertical red dashed line) are accurate for
$ \textit{Re}\ge 10^5$
. The lower
$ \textit{Re}=10^4$
cases deviate from the inviscid evolution only after the braid has thinned enough for viscous diffusion to act and limit further thinning, resulting in weaker baroclinic shear amplification than predicted. In cases with intermediate
$ \textit{Ri}$
, for which
$t^*$
is only shortly before
$ t_{2D}$
already in the inviscid limit, even relatively modest viscous effects are sufficient to prevent the early onset of SSI, as is observed for
$ \textit{Ri} = 0.10$
and
$0.15$
at
$ \textit{Re}=10^4$
.
5. Conclusions
We studied the mechanism leading to the early onset of SSI on the braids between KH billows over an extended range of Richardson and Reynolds numbers overlapping with ocean observations that frequently show braid-dominated turbulent mixing. To do so, we constructed a semianalytical inviscid model (see (2.10), (2.11)) for the transient evolution of the braid, comparing the results to high-resolution two-dimensional DNS.
As the braid steepens and thins under the influence of the growing primary billow, its stability reflects a competition between the stabilising strain and the destabilising baroclinic amplification of braid-aligned shear (see (2.9a–b
)). This balance is captured by the ratio of strain to braid-aligned shear,
$\gamma /\varOmega$
(see (2.13)). The simulations confirm that crossing the threshold
$\gamma /\varOmega \lt \varGamma ^*=0.02$
is necessary for SSI onset. The inviscid model further shows that the associated threshold-crossing time,
$t^*$
, is largely independent of
$ \textit{Ri}$
. In the inviscid limit, early SSI nevertheless requires
$ \textit{Ri} \gtrsim 0.10$
, for which sufficient separation develops between
$t^*$
and the primary-billow saturation time
$t_{2D}$
. Since
$t_{2D}$
increases steeply with
$ \textit{Ri}$
, this provides a longer window for SSI perturbations to grow to finite amplitude.
As summarised in figure 2, simulations confirm early SSI onset for
$ \textit{Ri} \gtrsim 0.10$
at
$ \textit{Re}=10^6$
and
$10^7$
, where the braid dynamics remain essentially inviscid and well predicted by our model. At
$ \textit{Re}=10^5$
, viscosity slightly slows braid thinning and thus shear amplification, delaying
$t^*$
and narrowing the window
$t_{2D}-t^*$
available for SSI growth (
$t_{2D}$
being independent of
$ \textit{Re}$
). This renders
$ \textit{Ri} = 0.10$
marginal, with onset requiring additional noise, likely present in natural flows. At
$ \textit{Re}=10^4$
, viscosity further slows and nearly halts braid thinning, thus preventing SSI at
$ \textit{Ri} = 0.10$
and
$0.15$
, and leaving only a marginal case at
$ \textit{Ri}=0.20$
.
We note that our simulations are two-dimensional and permit only a single primary billow wavelength, thus precluding vortex pairing and convective instabilities. Nonetheless, when SSI outpaces primary billow saturation, it should precede both types of instabilities. The growth rate of vortex pairing decreases rapidly with
$ \textit{Ri}$
(Klaassen & Peltier Reference Klaassen and Peltier1989) and is sufficiently small at
$ \textit{Ri} = 0.10$
that pairing generally occurs well after
$t_{2D}$
. Similarly, Mashayek & Peltier (Reference Mashayek and Peltier2012a
,
Reference Mashayek and Peltierb
) show that convective instabilities become significant only near or after
$t_{2D}$
. As noted in the introduction, other types of braid instabilities are possible beyond SSI. Although SSI is the first to emerge in all of our simulations, dynamics resembling that of the braid stagnation point instability described by Mashayek & Peltier (Reference Mashayek and Peltier2012a
) are observed by
$t_{2D}$
in some cases (e.g. panels k and o of figure 2).
An avenue for future work is to investigate the turbulence generated by the early onset of SSI, and its impact on the subsequent evolution of KH billows in three-dimensional DNS. As first noted by CS76, rapid mixing in the braid may strongly affect billow dynamics, and field observations at
$ \textit{Re} = O(10^5{-}10^7)$
indeed show intense braid mixing and weak mixing in billow cores, suggesting that braid turbulence may arrest primary billow growth. However, even in the absence of fully developed turbulence, rapid braid thinning during primary KH billow growth poses a significant computational challenge at high
$ \textit{Re}$
. Poorly resolved interfaces can trigger spurious instabilities that vanish with increasing resolution (McNally et al. Reference McNally, Lyra and Passy2012; Lecoanet et al. Reference Lecoanet, McCourt, Quataert, Burns, Vasil, Oishi, Brown, Stone and O’Leary2016; Lewin & Caulfield Reference Lewin and Caulfield2021; Briard et al. Reference Briard2024). The simple framework introduced here to estimate the interface thickness
$\delta$
provides a criterion for the minimum grid spacing required to capture SSI onset. In our simulations,
$\delta \gt 10\Delta z$
at the time of SSI onset. In three-dimensional DNS, this constraint should be considered alongside, and may be more restrictive than, resolution criteria based on turbulent Kolmogorov or Batchelor scales.
Similarly to the baroclinically accelerated KH braid, shear instabilities are also commonly observed as a route to turbulence in tilted interfaces present in a variety of other geophysical flow configurations, such as standing and progressive interfacial waves (Fringer & Street Reference Fringer and Street2003; Kalinichenko Reference Kalinichenko2005; Cavelier et al. Reference Cavelier, Gréa, Briard and Gostiaux2022), solitary waves (Barad & Fringer Reference Barad and Fringer2010; Xu et al. Reference Xu, Stastna and Deepwell2019) and gravity currents (Scotti Reference Scotti2008). As illustrated in these studies, local stability analyses have served as a useful guide for the onset of instability, though critical values of the Richardson number much less than
$1/4$
have often been reported as necessary for instability. The inclusion of a local ratio of strain to shear criterion in addition to the usual Richardson number may also prove useful in these contexts for making more accurate predictions regarding the onset of shear instability.
Finally, oceanic flows have
$ \textit{Pr} = O(10{-}10^3)$
– well above the
$ \textit{Pr} = 1$
considered here. At geophysically relevant
$ \textit{Re} \ge 10^6$
, we showed that the braid dynamics are effectively inviscid, so our results for SSI onset across
$({Ri},{Re})$
space should be independent of
$ \textit{Pr}$
. At lower
$ \textit{Re} = O(10^4{-}10^5)$
, however,
$ \textit{Pr} \gg 1$
can influence braid evolution and SSI onset, particularly near marginal
$({Ri},{Re})$
values. A higher
$ \textit{Pr}$
will reduce the interface thickness at which diffusion hinders strain-induced thinning, thus promoting closer agreement with the inviscid model and favouring SSI. However, supplementary simulations (supplementary material, § 5) show that it also induces non-trivial changes in the temporal evolution of
$ \textit{Ri}_{g}$
and
$\gamma /\varOmega$
. These high-
$ \textit{Pr}$
effects likely arise from the decoupling of buoyancy and shear interface thicknesses, which is not captured by our model and remains to be investigated.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2026.12012.
Acknowledgements
We thank A. Kaminski for helpful discussions and A. Atoufi and J. R. Taylor for sharing code that was adapted for the simulations.
Funding
The authors acknowledge the support of the Geophysical Fluid Dynamics Program at the Woods Hole Oceanographic Institution, supported by the NSF and ONR. E.R.B. was supported by Saint-Gobain. A.L. was supported by a NERC Independent Research Fellowship NE/W008971/1. This work used the Delta advanced computing and data resource through allocation PHY250188 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by NSF grants 2138259, 2138286, 2138307, 2137603 and 2138296. Delta is a joint effort of the University of Illinois Urbana–Champaign and its National Center for Supercomputing Applications, supported by the NSF (award OAC 2005572) and the State of Illinois.
Declaration of interests
The authors report no conflict of interest.
(Ri,Re)=(0.15,106)
b
b
u
w
b(x′,z′)
x′
b¯(z′)
u¯′(z′)
w¯′(z′)
w¯′
γ(t)
b
t2D
Re
Ri
t2D
[33.8,49.6,61.3,106.8]
Lx
[14.3,14.3,14.0,13.4]
Ri=
[0.05,0.10,0.15,0.20]
t2D
ϕ
γ
Re=[104,105,106,107]
Ri=[0.05,0.10,0.15,0.20]
t2D
tS
ϕ1,γ1
Rig
γ/Ω
Rig
γ/Ω
Γ∗=0.02
t∗
t2D
t2D−t∗
Ri
Ri≥0.10
Open access
