1. Introduction
Over the past three decades, Antarctic ice loss has contributed
$7.4 \pm 1.47$
mm to global mean sea level (Shepherd et al. Reference Shepherd2018; Otosaka et al. Reference Otosaka2023), and represents a major uncertainty in future sea level rise (Edwards et al. Reference Edwards2021; Hill et al. Reference Hill, Rosier, Gudmundsson and Collins2021; Fricker et al. Reference Fricker, Galton-Fenzi, Walker, Freer, Padman and DeConto2025). Antarctic outflow glaciers, grounded on land, flow into the ocean, where they become floating ice shelves, and once floating, their mass contributes to the sea level. Ice shelves are eroded by a combination of thinning by ocean-driven melting from below and calving at the edges (Depoorter et al. Reference Depoorter, Bamber, Griggs, Lenaerts, Ligtenberg, van den Broeke and Moholdt2013; Greene et al. Reference Greene, Gardner, Schlegel and Fraser2022), processes that reduce their ability to restrain (buttress) the flow of ice sheets into the ocean (Fürst et al. Reference Fürst, Durand, Gillet-Chaulet, Tavard, Rankl, Braun and Gagliardini2016; Reese et al. Reference Reese, Gudmundsson, Levermann and Winkelmann2018) and that are consequently the main cause of recent grounded ice losses over Antarctica (Pritchard et al. Reference Pritchard, Ligtenberg, Fricker, Vaughan, van den Broeke and Padman2012; Rignot et al. Reference Rignot, Jacobs, Mouginot and Scheuchl2013; Reese et al. Reference Reese, Gudmundsson, Levermann and Winkelmann2018). Quantifying the rate at which ice shelves melt into the underlying ocean, now and in the future, is therefore an urgent challenge that must be tackled if we are to improve sea level projections (Edwards et al. Reference Edwards2021; Hill et al. Reference Hill, Rosier, Gudmundsson and Collins2021). The rate at which ice melts from below is critically dependent on mixing processes that bring deeper warm waters up towards the ice. Yet, the technical and logistical challenges of direct observations limit our current understanding of this process, leaving a critical knowledge gap.
Theoretical models (MacAyeal Reference MacAyeal1985; Jenkins Reference Jenkins1991, Reference Jenkins2021), large eddy simulations (Vreugdenhil & Taylor Reference Vreugdenhil and Taylor2019; Begeman et al. Reference Begeman, Asay-Davis and Van Roekel2022; Vreugdenhil et al. Reference Vreugdenhil, Taylor, Davis, Nicholls, Holland and Jenkins2022; Anselin et al. Reference Anselin, Holland, Jenkins and Taylor2024) and limited observations (Davis & Nicholls Reference Davis and Nicholls2019; Rosevear, Galton-Fenzi & Stevens Reference Rosevear, Galton-Fenzi and Stevens2022; Davis et al. Reference Davis2023; Schmidt et al. Reference Schmidt2023) have indicated that, at the interface between the ocean and ice shelf, a stratified shear flow can form, regulating the vertical transfer of heat towards the ice. Freshwater released into the ocean by ice ablation is more buoyant than the ambient, salty (but relatively warm) ocean, and flows along the gently tilting underside of the ice shelf (figure 1 a). This buoyant freshwater forms a stratification which means that fluid parcels perturbed vertically feel a restoring (buoyant) force, and the turbulent transport of scalars across stratified interfaces is suppressed. At the same time, the shear generated by this buoyant flow acts to produce turbulence and destabilise the flow, causing warm, salty ambient ocean water to entrain (mix) into the ice-shelf-adjacent meltwater-modified flow. In order to model the effect of this flow on ice melting in the real world, a fundamental understanding of the flow under idealised conditions is needed; however, there is currently no such fundamental understanding of this situation of a stratified shear instability developing in proximity to a constant buoyancy interface.
Schematic of (a) the ice-shelf–ocean boundary layer, where a fresh, cold meltwater-modified layer is formed from mixing of the warm, salty ocean with the fresher, colder products of ice-shelf melting. The upper layer propagates along the ice underside upslope. (b) Indicative density (
$\rho$
) profiles (background and blue/white line) and horizontal velocity (
$u$
) profiles (black/white line) from the growth phase of the buoyant boundary shear current. Direction of the gravity vector,
$g$
is shown.
Figure 1. Long description
The figure consists of two schematic panels, labelled (a) and (b), illustrating the physical and numerical setup of a stratified flow problem. Panel (a) depicts a tilted rectangular domain representing a cross-section beneath an ice shelf. Three distinct fluid regions are identified by text labels: ‘Fresh, cold ice shelf’ at the top, ‘Fresh, cold meltwater’ in an intermediate layer, and ‘Warm, salty ocean’ below. Panel (b) illustrates the boundary conditions applied to the domain. The upper boundary has a fixed density condition (ρ at z = L_z = constant), the lower boundary has a fixed density condition (ρ at z = 0 = constant), and the lateral (left and right) boundaries are periodic. The domain is inclined to the horizontal, and a gravity vector g is decomposed into two components: g cos θ acting perpendicular to the domain floor, and g sin θ acting along it (parallel to the slope). Coordinate axes are shown, with z pointing perpendicular to the inclined base and x along it.
Motivated by such a geophysical setting, this paper uses numerical simulations at the laboratory scale to study a stratified shear flow that develops a unique instability under the following conditions:
-
(i) The density at the upper and lower boundaries is fixed at a constant value (Dirichlet boundary conditions), representing the boundary formed by the melting interface between the ice shelf (upper boundary) and far-field ocean (lower boundary).
-
(ii) The velocity at the upper and lower boundaries is 0 (no-slip boundary conditions), representing a solid upper interface (the ice shelf at the upper boundary). Matching lower interface boundary conditions are required by the model set-up, but this condition has little impact as it is far from any dynamic processes.
-
(iii) The domain is tilted, in effect, providing a gently sloping upper interface in a computationally efficient manner.
-
(iv) The domain is periodic in the along-slope (
$x$
) dimension, representing conditions far from a lateral boundary.
Studying flows at this scale allows for the smallest scale of the dynamics that would be missed or parametrised in ocean models to be fully resolved. In turn, these simulations can provide insight into the calibration of sub-grid-scale models used in large eddy simulations (LES) and parameterisations in larger-scale ocean models.
In general, it is the balance between shear and stratification that determines where flow instability and associated mixing may occur. The relationship between shear and stratification can be characterised by the gradient Richardson number
\begin{align} \textit{Ri}_g(z) \equiv \dfrac {-g}{\rho _0}\dfrac {\mathrm{d}\rho /\mathrm{d}z}{(\mathrm{d} u/\mathrm{d}z)^2} , \end{align}
where
$g$
is the gravitational acceleration,
$\rho _0$
a reference density,
$\mathrm{d}\rho /\mathrm{d}z$
the vertical density gradient and
$\mathrm{d}u/\mathrm{d}z$
the vertical gradient in the streamwise (along-slope) velocity, and where the condition
${Ri}_g\lt 1/4$
is a necessary, but not sufficient criterion for shear instability to occur (Miles Reference Miles1961). Under idealised conditions, stratified shear flows are categorised into three canonical unstable modes. Firstly, the Kelvin–Helmholtz instability (KHI) (Helmholtz 1868; Kelvin Reference Kelvin1871) describes flow where the interface rolls up into spiralling billows, where density interfaces overturn and mixing occurs by convective and secondary instabilities. Secondly, the Holmboe wave instability (HWI) describes a flow where the stratified interface is continually scoured away (Holmboe Reference Holmboe1962). The final of these modes, the Taylor–Caulfield instability (Taylor Reference Taylor1931; Caulfield et al. Reference Caulfield, Peltier, Yoshida and Ohtani1995) is not relevant for this study, but is fully described in Eaves & Balmforth (Reference Eaves and Balmforth2019). In each case, mixing characteristics of the flow are uniquely determined by overturning, scouring and the development of secondary and/or convective instabilities that can transition the flow to a turbulent state.
Our focus on under-ice meltwater flows is motivated by the occurrence of an unusual combination of characteristics that shape the flow and control the development of instability into a mixed-mode instability. Because of the gentle boundary slope, the stratification that develops adjacent to it is a source of both static stability and dynamic instability because of the sheared flow that develops (figure 1). Forcing on the flow arises from the buoyancy generated at the boundary itself, so is strongest at the boundary and decays away from it. Meanwhile, velocity shear is generated through both the sloping isopycnals and viscous drag, since the meltwater flows along a solid boundary, and the far-field velocity is zero in quiescent regions with no background or tidal currents. This differs from canonical boundary layer structures where a far-field velocity decays monotonically to zero at the boundary. A complex relationship exists between shear and buoyancy profiles that evolves with distance from the interface, while the development of instability is influenced by the proximity of the boundary. In the present study we leave aside other features of real-world, under-ice flows. These complexities could include implementing two buoyancy-determining properties, temperature and salinity, which have widely differing diffusivities and are both influenced by melting; and the rapid planetary rotation of the Earth. Instead, we focus on the fundamentals of the interplay between buoyancy and shear profiles and how the resulting instabilities develop in proximity to the solid boundary.
This work explores the fundamental dynamics of an unusual but nevertheless important class of flows, which develop in specific geophysical settings such as under melting ice shelves or where the atmospheric boundary layer is cooled over a sloping terrain (the process that gives rise to katabatic winds). We make use of idealised laboratory-scale numerical simulations (described in
$\S$
2.1) of a tilted rectangular domain (described in section
$\S$
2.2), alongside linear stability analysis (described in
$\S$
2.3) of profiles from the numerical simulations. In
$\S$
3, the dynamics of the idealised ice–ocean flow are presented from growth of the buoyant boundary shear current (
$\S$
3.1), through primary (
$\S$
3.2) and secondary instability phases (
$\S$
3.3) to the long-term evolution of the flow. A mechanistic explanation of the flow (
$\S$
3.2.3) is discussed. Finally, the implication of these results is discussed (
$\S$
4), specifically what these results indicate about marginal instability in meltwater flows beneath ice shelves.
2. Methods
2.1. Numerical model (SPINS)
Numerical simulations were carried out with the pseudospectral code SPINS described by Subich, Lamb & Stastna (Reference Subich, Lamb and Stastna2013), designed to simulate stratified fluids with high (spectral) accuracy and efficiency when parallelised. The model solves the stratified incompressible Navier–Stokes equations subject to the Boussinesq approximation
\begin{align} \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} = -\frac {1}{\rho _0} \boldsymbol{\nabla } P + \nu \boldsymbol{\nabla} ^2 \boldsymbol{u} -\frac {\rho g}{\rho _0} \hat {e}_g , \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \rho }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla } \rho = \kappa \boldsymbol{\nabla} ^2 \rho , \\[-9pt] \nonumber \end{align}
where
$\boldsymbol{u}$
is the velocity,
$t$
is time,
$P$
is the pressure,
$\rho$
is the density and
$\rho _0$
is the reference density of the fluid
$1000\,\textrm {kg m}^{-3}$
. The physical parameters are gravity
$g$
(set at
$9.81$
ms
$^{-2}$
), the shear viscosity
$\nu$
(set at
$10^{-6}$
m
$^2$
s
$^{-1}$
, unless otherwise stated, chosen to be consistent with the physical value) and the scalar diffusivity
$\kappa$
set so that the Schmidt number,
$ \textit{Sc} = \nu /\kappa = 7$
for all cases (equivalent to the molecular value of
$\textit{Pr} = \nu /\kappa$
for temperature (Nayar et al. Reference Nayar, Sharqawy, Banchik, Lienhard and John2016)). The unit vector in the vertical direction is denoted by
$\hat {e}_g$
. The model solves the equations using a spectral collocation method, using a third-order multi-step method developed to allow adaptive time steps. An exponential filter is used to ensure numerical stability near the grid scale according to default settings (see Subich et al. Reference Subich, Lamb and Stastna2013,for details).
The code has been thoroughly validated using theoretical results and physical laboratory experiments in a number of different configurations including shear instabilities (Subich et al. Reference Subich, Lamb and Stastna2013), boundary layer instabilities (e.g. Harnanan, Stastna & Soontiens Reference Harnanan, Stastna and Soontiens2017), interaction with topography (e.g. Deepwell et al. Reference Deepwell, Stastna, Carr and Davies2017) and shoaling internal solitary waves (e.g. Hartharn-Evans et al. Reference Hartharn-Evans, Carr, Stastna and Davies2022). It is available for download through its online manual: https://wiki.math.uwaterloo.ca/fluidswiki/index.php?title=SPINS_User_Guide.
2.2. Case set-up and experimental design
The numerical tank is an initially two-dimensional (2-D) rectangular (unmapped) domain with length
$L = 1$
m, and depth
$H = 0.3$
m. One simulation is extended for 3-D simulation in § 3.3. Sloping upper and lower boundaries were implemented by ‘tilting’ the domain, adding a term to
$\hat {e}$
in the buoyancy forcing term in (2.1) of the form
$\hat {e} = (-\sin \theta , \cos \theta )$
. Since the grids remain oriented with the boundaries of the tank, this method remains computationally efficient by avoiding the need for a mapped grid. Throughout this paper, references to the streamwise and vertical directions strictly refer to the
$x$
-axis (along the tank) and
$z$
-axis (vertical in the frame of the tank) (figure 1
a). Periodic boundary conditions were applied in the streamwise (
$x$
) direction. Grid resolution in the
$x$
dimension was a uniformly spaced
$N_x = 1024$
points, giving
$\mathrm{d}x = 0.98$
mm. Both 2-D and 3-D simulations were carried out without Coriolis rotation.
Numerical simulations were carried out with no-slip boundary conditions (see table 1) at the upper and lower boundaries and employed a Chebyshev grid in the vertical, implying a clustering of points near both the upper and lower boundaries that scales with the number of points in the vertical squared. For these simulations, grid resolution in the vertical was
$N_z = 512$
, giving variable
$\mathrm{d}z$
between
$2.8\, {\unicode{x03BC}} \textrm {m}$
at the upper and lower boundaries and
$920 \, {\unicode{x03BC}} \textrm {m}$
near mid-depth, such that the region of the tank that becomes turbulent has high vertical resolution. We demonstrated mesh independence in Appendix A.
Parameters for simulations used in this paper; density difference
$\Delta \rho$
, slope
$\theta$
and viscosity
$\nu$
. Parameters that differ from those in ISO_Base are highlighted in bold.
Table 1. Long description
7 simulations are presented with the density difference (kg m^-3), slope (degrees) and viscosity (m^2/s)ISO_Base simulation has density difference of 10, slope of 2 and viscosity 10^{-6}ISO_0.5x_rho simulation has density difference of 5, slope of 2 and viscosity 10^{-6}ISO_2x_rho simulation has density difference of 20, slope of 2 and viscosity 10^{-6}ISO_5_slope simulation has density difference of 10, slope of 5 and viscosity 10^{-6}ISO_10_slope simulation has density difference of 10, slope of 10 and viscosity 10^{-6}ISO_2x_nu simulation has density difference of 10, slope of 2 and viscosity 2x10^{-6}ISO_0.5x_nu simulation has density difference of 10, slope of 2 and viscosity 5x10^{-7}.
Upper boundary conditions representing the ice–ocean boundary were implemented via Dirichlet (fixed value) vertical boundary conditions for density
\begin{align} \rho (z = 0) = \rho _0, \\[-28pt] \nonumber \end{align}
\begin{align} \rho (z = H) = \rho _1. \\[0pt] \nonumber \end{align}
An initial density profile was set using a very thin tanh profile to meet the smoothness condition for the spectral model at initialisation
\begin{align} \rho (z) = \rho _0 + \frac {\Delta \rho }{2}\left (1-\tanh \left (\frac {2\left (z-(H - z_{\textit{pyc}})\right )}{h_{\textit{pyc}}}\right )\right ) , \end{align}
where
$\Delta \rho = \rho _0 - \rho _1$
is the density difference,
$z_{\textit{pyc}}$
the pycnocline location (distance beneath the upper interface) and
$h_{\textit{pyc}}$
is the pycnocline thickness. Both
$h_{\textit{pyc}}$
and
$z_{\textit{pyc}}$
are set throughout these simulations at
$0.001$
m to be as small as possible whilst maintaining the smoothness criteria in initial conditions (identified via test simulations). In figure 5
a, we show these profiles evolve in line with the expected analytic profiles, showing the simulations are insensitive to the initial condition.
The velocity profiles
$u(z), w(z)$
were initialised at rest with small white noise perturbations (with amplitude of
$10^{-3}$
) added at the initialisation of all simulations to trigger instability.
Simulations carried out are summarised in table 1, where prefixes of ISO (and throughout this paper) refer to ice-shelf–ocean simulations. Unless otherwise stated, results in the ensuing sections are based on ISO_Base, with additional sensitivity simulations presented later.
2.3. Linear stability analysis
By applying theoretical approaches to the outputs of time-resolved numerical simulations, fundamental processes driving the dynamics can be understood and identified. In particular, we aim to understand the stability of the time-evolving flow through the use of linear stability analysis (LSA), linearising about snapshots of the slowly evolving 2-D simulation profiles, averaged in
$x$
. The 2-D simulation profiles are essentially one-dimensional until instabilities grow and dominate the flow; averaging in
$x$
is performed to mitigate these small-scale instabilities to enable LSA to provide reasonable predictions during the early stages of unstable growth.
The implicit assumption here is that LSA solutions are valid when the flow evolution is slow compared with the growth rate of instabilities. In other words, once unstable growth begins in the simulations, it dominates over other flow evolution driven by either gravity or diffusion. The validity of these assumptions will be shown in
$\S$
3.2.1.
The base flow profiles at time
$t$
are denoted
$\boldsymbol{U}(z) = (U(z),0,0)$
, for velocity where
$U(z) = \langle (u(x,z)\rangle _x$
,
$P(z) = \langle P \rangle _x$
for pressure and
$B(z) = - \langle \rho (x,z) g / \rho _0 \rangle _x$
for buoyancy, where
$\langle \rangle _x$
represents averaging in the streamwise direction. Perturbations to these base flow profiles are given by
$\boldsymbol{u}'$
,
$p'$
and
$b'$
. Following the procedure outlined in Drazin & Reid (Reference Drazin and Reid2004), Liu, Thorpe & Smyth (Reference Liu, Thorpe and Smyth2012) and Lian, Smyth & Liu (Reference Lian, Smyth and Liu2020), linearised equations can be formulated depending only on the vertical velocity and buoyancy perturbations. General normal-mode solutions are subsequently sought of the form
$w' = \hat {w}(z) \text{exp}(ik_x x + \lambda t)$
and
$b' = \hat {b}(z) \text{exp}(ik_x x + \lambda t)$
, allowing for a general vertical dependence of the complex amplitudes or eigenfunctions
$\hat {w}$
and
$\hat {b}$
. Here,
$k_x$
represents the streamwise wavenumber, and
$\lambda = \sigma - i\omega$
represents a complex growth rate or eigenvalue related to the real growth rate
$\sigma$
and the temporal frequency
$\omega$
. Consistent with the simulations, we restrict our analysis to 2-D modes with
$k_y = 0$
. Substitution into and manipulation of the governing equations leads to the viscous Taylor Goldstein (vTG) eigenvalue problem (Liu et al. Reference Liu, Thorpe and Smyth2012; Lloyd et al. Reference Lloyd, Dorrell and Caulfield2022; Lloyd & Dorrell Reference Lloyd and Dorrell2024)
\begin{align} \lambda \begin{bmatrix} \Delta & 0 \\ 0 & I \end{bmatrix} \begin{bmatrix} \hat {w} \\ \hat {b} \end{bmatrix} = \begin{bmatrix} - i k_x U \Delta + i k_x \frac {\text{d}^2 U}{\text{d} z^2} + \widehat {D}_w & - k_x^2 \\ - \frac {\text{d} B}{\text{d} z} & - i k_x U + \widehat {D}_b \end{bmatrix} \begin{bmatrix} \hat {w} \\ \hat {b} \end{bmatrix}, \end{align}
where
$I$
represents the identity matrix,
$\Delta = ({\text{d}^2 }/{\text{d} z^2}) - k_x^2$
and the diffusive operators are defined as
\begin{align} \widehat {D}_w = \nu \left ( k_x^4 - k_x^2 \frac {\text{d}^2 }{\text{d} z^2} + \frac {\text{d}^4 }{\text{d} z^4} \right )\ \ \text{and} \ \ \widehat {D}_b = \kappa \left (-k_x^2 + \frac {\text{d}^2 }{\text{d} z^2} \right ). \end{align}
Note that we have neglected the streamwise component of the gravity unit vector
$\boldsymbol{k}$
due to the low slope angles (
$\theta \lt \sin \theta \ll \cos \theta$
), where the destabilising effect of
$\sin \theta$
on the streamwise component is negligible (Atoufi et al. Reference Atoufi, Zhu, Lefauve, Taylor, Kerswell, Dalziel, Lawrence and Linden2023).
For a given wavenumber
$k_x$
, and specified profiles of
$U$
and
$B$
, the differential eigenvalue problem (2.7) can be (numerically) solved to obtain the complex growth rates (eigenvalues)
$\lambda$
and associated structure functions (eigenvectors)
$\hat {w}$
and
$\hat {b}$
. The eigenvalue problem (2.7) is numerically solved, with discrete operators derived using Chebyshev polynomials with a resolution of
$N_z=128$
grid points (note that solutions are insensitive to further increase in resolution). Six boundary conditions are required to close the system. Consistent with the nonlinear simulation boundary conditions we specify
$\hat {w} = 0$
,
$\hat {w}_z = 0$
and
$\hat {b} = 0$
at the lower and upper boundaries, implemented using a ‘give-back’ matrix, following the procedure outlined by Lian et al. (Reference Lian, Smyth and Liu2020).
Every 10 s of compute time we take the planar-averaged profiles from simulations and compute vTG solutions, using a robust peak-finding algorithm to maximise the growth rate
$\sigma$
and its corresponding wavenumber
$k_x$
. After some a priori testing we find that the maximum growth rates occur with
$k_x \approx 18$
. To find the precise
$k_x$
that maximises
$\sigma$
, to within two decimal places, we first adopt a coarse
$k_x$
grid between 0 and 30, then iteratively refine using a bisection method.
As we will show in § 3.2.1, the assumptions made for LSA produce modal structures, wavelengths and growth rates in excellent agreement with full numerical simulations.
3. Flow characteristics
3.1. Time-evolving base flows
After initialisation at rest, buoyancy (low density) diffuses from the upper boundary, which forms a continually stratified layer, the depth of which increases entirely by diffusion of density down the density gradient (figure 2 a). The density evolves in a way that is insensitive to the initial condition (beyond a requirement for a smooth, thin transition), quickly deviating in form from the tanh profile ((2.6), figure 2 a) to an error function. The density interface consists of a quasi-linearly stratified layer at the upper surface, overlaying a well-mixed lower layer (figures 2 a and 2 e). This density gradient weakens over time as the depth of the stratified layer increases (without a corresponding increase in magnitude of the density difference).
With this source of upslope buoyancy forcing, the flow begins to accelerate, forming a scythe shaped velocity profile (figure 2
b). At the upper boundary, the velocity is 0 due to the no-slip condition, and similarly the ambient fluid is at rest owing to a lack of forcing at depth, in between the two is the velocity maximum. Between the upper boundary and the peak in velocity at
$z\approx 0.29$
m, a strong shear layer (
$\mathrm{d}u/\mathrm{d}z \gt 0$
) forms, and between the ambient fluid and the velocity maximum, a further shear layer of the opposite polarity (
$\mathrm{d}u/\mathrm{d}z \lt 0$
) forms (figures 2
b, 3
a– and 3
b). The flow is both laminar and parallel to the slope (and therefore also to isopycnals), and so drives no additional mixing (as confirmed by comparison of density profiles in a simulation with
$\theta = 0$
). During this phase (until
$t \approx 570$
s) the flow is entirely one-dimensional (uniform in the streamwise direction) (figure 3
a).
The evolution of flow in simulation ISO_Base. In (a–c) the evolution of profiles of (a) density, (b) streamwise velocity at
$t =$
5225, 450 and 570 s and (c) profile of Richardson number at
$t = 570$
s prior to instability forming, where shaded regions indicate depths where
${Ri}_g (z)\lt 0.25$
. Panels (d–m) show the evolution of density fields at times
$t =$
400, 570, 586, 594 and 600 s (d–h) prior to secondary instability, and
$t =$
610, 630, 660, 690 and 770 s (i–m). In panels (d–g) the white line is provided as the isopycnal at the velocity maximum to aid visualisation of the internal wave. Supplementary material and movies are available at 1.
Figure 2. Long description
The figure has two sections: a top row of three line-graph panels (a–c) and a lower grid of ten pseudocolour panels (d–m). Panel (a) – Density profile (ρ, kg m-3) vs. height z (m, 0.26–0.30). Three curves in black, dark grey, and light grey show a sharp pycnocline near z = 0.285 m, where density decreases from ∼1000 kg m-3 in the lower layer to ∼990 kg m-3 in the upper layer. In successive curves, the curve is deeper and les sharp. Panel (b) - Horizontal velocity profile (u, m s-1) vs. z. All three curves are negative (−0.12 to 0), indicating flow in the negative x-direction throughout the shown depth range, with maximum speed near the pycnocline and decreasing toward the upper boundary. Panel (c) - Gradient Richardson number (Ri_g) vs. z, plotted from 0 to 1. A vertical dashed line marks Ri_g = 0.25 (the canonical shear-instability threshold). Values rise steeply from near zero at the surface to well above 0.25 in the upper layer, before dropping back to near-zero at greater depths below the pycnocline. Panels (d–m) – Pseudocolour images of density (ρ, kg m-3) in the x–z plane (x: 0.1–1.0 m; z: 0.26–0.30 m), arranged in two columns of five rows. A shared colourbar ranges from 990 (light blue) to 1000 (dark purple) kg m-3. The left column (d–h) and right column (i–m) represent two different experimental conditions or time sequences. Panel (d): Nearly flat, undisturbed interface visible as a thin white contour line near z = 0.285 m; dense purple fluid below, light blue above. Panel (e): Interface remains largely flat with a slight waviness beginning to develop. Panel (f): Interface shows gentle, long-wavelength undulations; otherwise laminar. Panel (g): Visible wave-like distortions along the interface; amplitude increasing. Panel (h): Clear Kelvin–Helmholtz-like billows forming; interface broken into rolling structures. Panel (i): Interface already disrupted; irregular density structures visible throughout the domain. Panel (j): Strongly undulating interface with overturning features across the full width. Panel (k): Intense mixing; complex interleaving of dense and light fluid, with coherent vortical structures visible. Panel (l): Advanced turbulent mixing; density field highly disordered with filamentary structures. Panel (m): Partially restratified after mixing event; large-scale structure reduced but density field remains irregular.
Time sequence of vorticity (
$\zeta = \boldsymbol{\nabla }\times \boldsymbol{u}$
) outputs from simulation ISO_Base, where panels (a–j) match figure 2 panels (d–m) respectively.
Figure 3. Long description
Ten pseudocolour panels (a–j) showing the spanwise vorticity field ζ (s^-1) in the x–z plane (x: 0.1–1.0 m; z: 0.26–0.30 m), arranged in two columns of five rows. A shared diverging colourbar runs from –10 s^-1 (green) through 0 (white) to +30 s^-1 (dark purple), with positive (cyclonic) vorticity in red-purple and negative (anticyclonic) in green. Panel (a): Nearly quiescent flow. A thin, intense red-pink band of positive vorticity sits just below z = 0.30 m (the upper boundary shear layer). The remainder of the domain is uniformly near-zero (white-pale green). A black isopycnal contour sits at approximately z = 0.285 m, nearly flat. Panel (b): Similar to (a) but with a slightly thicker upper shear layer and a faint red tinge developing near the interface. The isopycnal remains nearly horizontal with very slight waviness. Panel (c): The upper shear layer remains prominent. The isopycnal shows gentle, long-wavelength undulations across the domain. A weak, diffuse region of positive vorticity has spread slightly below the interface. Panel (d): Growing interfacial waves are visible in the isopycnal, now with moderate amplitude. Positive vorticity extends further into the interior below the upper boundary layer, with a banded structure suggestive of wave-induced shear. Panel (e): The upper boundary vorticity layer remains, and the interior shows broad, patchy regions of positive vorticity extending across the full horizontal extent. The flow is approaching a critical state for instability. Panel (f): Strong Kelvin-Helmholtz billows are evident. Concentrated red/purple patches of high positive vorticity alternate with green (negative) regions across the domain, indicating overturning structures. The upper boundary layer is disrupted and irregular. Panel (g): Fully developed Kelvin-Helmholtz-like billows dominate. Alternating red and green vorticity cores are arrayed quasi-periodically in x, with clear cat’s-eye structures indicating coherent overturning. The flow is strongly nonlinear. Panel (h): Advanced turbulent breakdown. The organised billow structure is fragmenting; intense red vorticity cores persist but are surrounded by fine-scale green filaments, indicating secondary instabilities and the onset of three-dimensional turbulence. Panel (i): Turbulent mixing phase. Vorticity is distributed across a wide range of scales, with interspersed red and green filamentary structures throughout the domain. No coherent large-scale billow organisation remains. Panel (j): Post-mixing, partially restratified state. Vorticity magnitudes are reduced overall; broad regions of weak positive (red-pink) vorticity dominate, with residual green patches. The flow is decaying toward a more quiescent state.
Following Carpenter et al. (Reference Carpenter, Tedford, Rahmani and Lawrence2010b ), it is useful to identify interfaces on which vorticity and gravity waves can propagate. The density profile indicates there is a single density interface that supports gravity waves that extends across the entire depth, however, gravity wave propagation will be concentrated over the region of stratification near the surface. The velocity profile indicates the presence of two vorticity interfaces (figure 2 b), one approximately coincident with the density interface, and one just below it, with opposing signs.
Using the classical critical gradient Richardson number
${Ri}_g \lt 1/4$
, we identify two regions of flow susceptible to shear instability from as early as
$t = 70$
s. A near-surface region with
${Ri}_g \lt 1/4$
is driven by high shear at the no-slip boundary and a further deeper region where
${Ri}_g \lt 1/4$
from
$z \approx 0.2725{-}0.285$
m occurs where the density gradient is gentle, and shear is moderate (figure 2
a–c). These two regions are separated vertically by a stable region around the velocity maximum (where shear is therefore lowest and, at the velocity maximum itself, identically zero). However, despite the necessary, but not sufficient, critical Richardson number for shear instability being reached at this point, the flow remains stable. One dimensional mean flow profiles continue to dominate. Once the velocity structure has developed (by around
$t = 70$
s) the flow always has regions where
${Ri}_g \ll 1$
and regions where
${Ri}_g \to \infty$
, and hence
${Ri}_g \lt 0.25$
(or any other critical value) is a poor parameter for assessing stability of this flow.
3.2. Primary instability of the base flow
In numerical simulations, small perturbations on these interfaces grow into waves of finite amplitude on the interface (figures 2 f and 3
b) propagating in the same direction as the main flow (right to left) at
$-0.0859\,\textrm {m s}^{-1}$
(slower than
$U_{\textit{max}} = -0.1263\,\textrm {m s}^{-1}$
) from
$t \approx 580$
s. From the cusps of the waves, long wisping structures emerge in a paired system centred on the velocity maximum. On each resulting vorticity interface these wisping structures emerge, with smaller-amplitude wisps (geometrically constrained by the upper boundary) forming upward and clockwise, and larger-amplitude wisps that develop into overturning billows forming downward and anti-clockwise (figures 2
g, 2
h, 3 d, 3
e).
3.2.1. Linear stability of the time-evolving base flow
To explain the instability, and the conditions required to go unstable, we first explore the time dependence of the dominant modes arising from solutions to the vTG system, subject to temporally evolving buoyancy and velocity background states. The time dependences of the dominant mode eigenvalue with
$\text{max}(\sigma )$
and its corresponding wavenumber
$k_x(\sigma = \text{max}(\sigma ))$
and frequency
$\omega (\sigma = \text{max}(\sigma ))$
are shown in figure 4. First, note that the LSA predicts a stable flow until
$t \approx 340$
s, at which point a single dominant mode emerges with
$\sigma \gt 0$
, with
$\sigma$
growing approximately linearly with time beyond this point. Note that we do not perform computations after
$t = 560$
s, at which point planar-averaged simulation profiles are dominated by the growing instability and LSA is no longer useful (see e.g. figures 2
e and 2
f). After
$t = 340$
s,
$\sigma$
and
$\omega$
change fairly linearly with time (
$\sigma$
increasing,
$\omega$
decaying), while the dominant
$k_x$
slowly decreases with time, dropping from approximately 18 to 17 over 200 s. The inset of figure 4 (a) shows the single dominant mode with
$\sigma \gt 0$
, computed with profiles at
$t=400$
s with
$k_x = 17.9$
. This is indicative of all solutions for
$340 \lesssim t \lesssim 560$
s. Note that the frequency
$\omega$
of the dominant instability is negative, such that instabilities propagate in the direction of the mean flow.
Time series of the wavenumber
$k_x$
(a), maximum growth rate
$\text{max}(\sigma )$
(b) and frequency
$\omega$
(c) associated with the dominant mode emerging from vTG solutions (mode with the maximum growth rate) subject to simulation flow profiles. Panel (a) inset shows frequency growth rates associated with the
$k_x$
marked by the cross, demonstrating a single mode with
$\sigma \gt 0$
. Lines in panel (b) represent estimated growth rates of turbulent kinetic energy (TKE) for the simulation ISO_Base (solid), and the same simulation but initialised at
$t=300$
with the LSA dominant mode eigenfunctions (dashed).
To compare LSA solutions against simulations we have also estimated TKE growth rates following Howland, Taylor & Caulfield (Reference Howland, Taylor and Caulfield2018). Instantaneous growth rates from simulation data are calculated as
$\sigma _{\textit{sim}} = ({1}/({2E}))({\mathrm{d} E}/{\mathrm{d} t})$
, where
$E$
is the perturbation kinetic energy calculated
$({1}/{2})\langle (\boldsymbol{u}')^2\rangle$
. Growth rates of ISO_Base are reported in figure 4(b) (solid red line). We see that unstable growth does not occur until
$t \approx 520$
s for ISO_Base, but when the flow does go unstable the growth rates compare well against the LSA predictions, with
$\sigma \approx 0.1$
. The simulation TKE growth rate peaks at
$t \approx 590$
s, at which point strong overturning dominates the dynamics (figures 2
f and 2
g). However, when the simulations are seeded with the appropriate perturbations, unstable growth occurs at the time predicted by the LSA (see figure 4
b). Here, the simulation was restarted at
$t=300$
s with low-amplitude perturbations added to the flow profiles, taken from the eigenfunctions associated with the most unstable modes predicted by LSA (see Appendix B for details). In this case, the growth of
$E$
is consistent with LSA predictions until nonlinearity dominates (
$t \gtrsim 400$
s), indicating that, while growth is possible during these earlier stages of flow development, it requires a particular form of perturbation to experience instability.
To develop our understanding of why the slowly evolving profiles go unstable at
$t \approx 340$
s, we first revisit the planar-averaged simulation profiles, rescaling them by the time-varying maximum absolute velocity,
$U_{\textit{m}}$
and the boundary current thickness,
$\delta _m$
, which we define as the distance from the upper boundary down to the point beneath the velocity maximum where the absolute flow velocity reduces below 1 % of its maximum value. The planar-averaged time-varying simulation profiles of streamwise velocity, buoyancy and associated spatial derivatives, are reported in figure 5, where velocities are scaled by
$U_{\textit{m}}$
and spatial coordinates by
$\delta _m$
. When the vertical coordinate is taken as
$(z-H)/\delta _m$
, all appropriately scaled profiles collapse for
$ t \lesssim 560$
s, after which instabilities dominate the simulation (see e.g. the
$t=600$
s profile in figure 5). This clearly demonstrates that the flow is self-similar during its slow evolution phase, until instability dominates for
$t \gtrsim 560$
s. Interestingly,
$\delta _m$
appears to be an appropriate length scale to collapse the buoyancy fields, despite the comparatively thin density interface when compared with shear. This indicates that the ratio
$\delta _b/\delta _m$
is approximately constant over this time range (figure 7
d), where
$\delta _b$
represents the thickness of the buoyant jet, which we define as the distance from the upper boundary to the location where
$B$
is 1 % of its maximum value.
Time dependence of planar-averaged velocity, buoyancy and respective spatial derivatives, from simulation data, scaled by the absolute maximum velocity
$U_{\textit{m}}$
and the boundary current thickness
$\delta _m$
. Note that self-similar profiles from
$t = 50$
to
$560$
s overlap and are hidden by the
$t = 560$
s profile. The markers of panels (a) and (d) represent the analytical solutions for a convective flow developing on a vertically oriented heated plate (Ke et al. Reference Ke, Williamson, Armfield, McBain and Norris2019).
Figure 5. Long description
Five vertical profile panels (a–e) sharing a common y-axis labelled (z – H)/δ_m, ranging from 0 at the top to –1.0 at the bottom, representing depth below the boundary normalised by the boundary current thickness δ_m. Each panel shows planar-averaged profiles at eight time steps from t = 50 s to t = 600 s, indicated by a shared legend at the top. Line colours progress from warm (orange/yellow at early times) through pink/magenta to black at later times; the t = 600 s profile is dashed black. A note in the caption states that profiles from t = 50–560 s are self-similar and overlap, so earlier profiles are largely hidden by the t = 560 s solid black line. Panel (a): x-axis is U/U_m (normalised horizontal velocity), ranging from –1 to 0. Profiles show a boundary current structure: velocity magnitude increases from zero at the surface (y = 0) to a maximum near (z – H)/δ_m ≈ –0.2, then decays back toward zero at (z – H)/δ_m = –1.0. Red circular markers overlie the profile, representing the analytical solution for a convective boundary current on a vertically heated plate (Ke 2019). All time steps collapse well onto a single curve, confirming self-similarity. Panel (b): x-axis is U_z δ_m / U_m (vertical shear of velocity, normalised), ranging from approximately –1 to 10+. The shear profile peaks sharply near the surface and decays with depth, again showing good collapse across time steps. Panel (c): x-axis is U_zz δ_m-2 / U_m (second vertical derivative of velocity, normalised), ranging from approximately –10 to 100+. This profile exhibits a sharp positive peak very close to (z – H)/δ_m = 0 and a negative lobe beneath, indicating the curvature structure of the boundary current velocity profile. Panel (d): x-axis is B (buoyancy), ranging from 0 to 1 (normalised). Profiles show buoyancy decreasing monotonically from a maximum at the upper boundary to near zero at (z – H)/δ_m = –1.0. Red circular markers again show the Ke (2019) analytical solution, with good agreement. Profiles from all time steps collapse tightly. Panel (e): x-axis is B_z δ_m (vertical buoyancy gradient, normalised), ranging from 0 to approximately 5+. The buoyancy gradient peaks near the surface and decays with depth. Unlike the other panels, the t = 600 s dashed black profile shows a modest but visible departure from the earlier self-similar collapse, suggesting the buoyancy gradient structure begins to evolve at late times.
Figure 5 also includes the analytical solution for a temporally evolving convective boundary layer that develops on an isothermal vertically aligned plate (Illingworth Reference Illingworth1950; Schetz & Eichhorn Reference Schetz and Eichhorn1962; Ke et al. Reference Ke, Williamson, Armfield, McBain and Norris2019). The present ice-sheet problem is equivalent to the convective boundary layer case if the wall were aligned with the gravity vector (
$\theta = 90 ^\circ$
). Clearly, during initial flow evolution, flow development is identical between the two cases. The flows evolve due to diffusion from the wall, and through acceleration due to the streamwise buoyancy force (in the present case,
$g_x = - g \sin \theta$
). Until nonlinearity develops, it is unsurprising that both flows collapse to the same velocity and buoyancy profiles, subject to appropriate scaling.
It is informative to compare our LSA with that of Ke et al. (Reference Ke, Williamson, Armfield, McBain and Norris2019), who performed equivalent analysis on the vertically aligned heated plate. In figure 6 we present our solutions using the same scaling as figure 2 of Ke et al. (Reference Ke, Williamson, Armfield, McBain and Norris2019), which is based on the diffusive length and time scales
\begin{align} \ell ^\nu = \big (\kappa ^2 / g_r \big )^{1/3}, \ \ \text{and}, \ \ \tau ^\nu = \big (\kappa / g_r^2 \big )^{1/3}, \end{align}
where
$g_r$
represents the reduced streamwise gravitational acceleration:
$g_r = (\Delta \rho / \rho _0) g \sin \theta$
. Definitions of the diffusive time
$t^\nu = t / \tau ^\nu$
, wavenumber
$\alpha = 2 \sqrt {t^\nu } k_x \ell ^\nu$
and growth rate
$\sigma \tau ^\nu$
, are consistent with the definitions of Ke et al. (Reference Ke, Williamson, Armfield, McBain and Norris2019), where the wavenumber time dependence arises due to the time dependence of the base flow. Instability is also presented as a function of the time-dependent Grashoff number
$Gr_\delta = g_r \delta _u^3 /\nu ^2$
, where
$\delta _u = \int _{-H}^0 (U/U_m) \ \mathrm{d} z$
is the integral boundary layer thickness. As noted by Ke et al. (Reference Ke, Williamson, Armfield, McBain and Norris2019), the flow can be assumed to vary sufficiently slowly for LSA when the Grashoff number is large.
Linear stability of the temporally evolving stratified shear flow, presented with a diffusive scaling. Here, the wavenumber and growth rates are scaled by the diffusive length and time scales
$l^\nu$
and
$\tau ^\nu$
, and instability is presented as a function of dimensional time
$t$
, diffusive time
$t^\nu$
and the Grashoff number
$Gr_{\delta }$
.
Figure 6. Long description
A single scatter/bubble plot showing the linear stability properties of a temporally evolving stratified shear flow as a function of time, with three aligned x-axes and one y-axis. Y-axis: α = 2(t^v) k_x l^ν, the diffusively scaled wavenumber, ranging from 0.0 to approximately 1.25. This is the normalised streamwise wavenumber of the most unstable mode. X-axes (three, vertically stacked below and above the plot): Top axis: Gr_δ (Grashof number based on boundary current thickness), ranging from approximately 20,000 to 90,000, increasing left to right. Middle axis (primary): t [s] (dimensional time), ranging from approximately 200 to 560 s. Bottom axis: t^ν (diffusive time), ranging from approximately 1000 to 2250, increasing left to right. All three x-axes are linearly related and represent the same temporal progression of the flow. Colour: Each point is coloured by the normalised growth rate στ^ν, with a colourbar on the right ranging from 0.00 (light yellow) through green to approximately 0.02 (dark teal/black). Higher values indicate faster-growing instabilities. Data description: No unstable modes are present at early times (t ≲ 300 s; Gr_δ ≲ 30,000). Unstable modes first appear around t ≈ 300–320 s with small wavenumbers (α ≈ 0.2–0.4) and low growth rates (pale yellow). As time progresses, the band of unstable wavenumbers broadens and shifts to higher α, while growth rates increase. By t ≈ 400–560 s the unstable region spans α ≈ 0.1–1.2, with the highest growth rates (darkest points, σr^ν ≈ 0.02) concentrated at intermediate wavenumbers around α ≈ 0.3–0.7. The overall envelope of unstable modes grows and fills out progressively, consistent with a flow that becomes increasingly susceptible to shear instability as the boundary current strengthens and the Grashof number rises.
The present flow is unstable over a range of wavenumbers
$\alpha \lesssim 0.6$
, and for
$t^\nu \gtrsim 1400$
or
$Gr_\delta \gtrsim 40\,000$
. This is considerably more stable than the vertically oriented flow of Ke et al. (Reference Ke, Williamson, Armfield, McBain and Norris2019), who observed instability for
$\alpha \lesssim 2$
,
$t^\nu \gtrsim 50$
and
$Gr_\delta \gtrsim 500$
. Nevertheless, there are similarities between the linear stability of figure 6 and that of Ke et al. (Reference Ke, Williamson, Armfield, McBain and Norris2019), particularly the magnitudes of the growth rates and unstable wavenumbers, and the qualitative shape of the marginal stability curve. It is expected that, as the slope angle
$\theta$
increases, the stability diagram would transition towards that of the vertically aligned problem of Ke et al. (Reference Ke, Williamson, Armfield, McBain and Norris2019).
This disagreement is of course due to the role of the vertical buoyancy force. We argue that the balance between vertical buoyancy forces and vertical shear (i.e. a Richardson number) controls the transition to instability. To this end, we reformulate the eigenvalue problem based upon time-invariant base profiles, with appropriate scaling by
$U_m$
and
$\delta _m$
(figure 5). Time invariance of the base profiles further justifies our use of LSA to describe the dynamics for this temporally and spatially evolving flow. Rescaling the governing equations leads to the following system:
\begin{align} \frac {\partial \boldsymbol{u}^*}{\partial t^*} + \left (\boldsymbol{u}^* \boldsymbol{\cdot }{\boldsymbol{\nabla }}^* \right ) \boldsymbol{u}^* = - {\boldsymbol{\nabla }}^* p^* + b^* \boldsymbol{e}_g + \frac {1}{Re_u}\Delta ^* \boldsymbol{u}^* , \\[-28pt] \nonumber \end{align}
\begin{align} {\boldsymbol{\nabla }}^* \boldsymbol{\cdot }\boldsymbol{u}^* = 0, \\[0pt] \nonumber \end{align}
and
\begin{align} \frac {\partial b^*}{\partial t^*} + \left (\boldsymbol{u}^* \boldsymbol{\cdot }{\boldsymbol{\nabla }}^* \right ) b^* = \frac {1}{Re_u Sc}\Delta ^* b^* , \end{align}
where the
$*$
superscript represents a dimensionless variable. Crucially, the buoyancy is now redefined as
$b^* = {Ri}_u b$
; the governing dimensionless parameters are therefore a Schmidt number, and time-varying Reynolds and Richardson numbers
\begin{align} \textit{Sc} = \frac {\nu }{\kappa }, \ \ \textit{Re}_u = \frac {U_m \delta _m}{\nu }\ \ \text{and} \ \ {Ri}_u = \frac {\Delta \rho g \delta _m}{U_m^2} = \frac {\delta _m}{H} \frac {u_b^2}{U_m^2} , \end{align}
where
$u_b = \sqrt {\Delta \rho gH/\rho _0}$
and
$H$
is the domain depth. Flow evolution, and therefore transition to an unstable dynamics, is therefore controlled by the time dependence of
$ \textit{Re}_u$
and
${Ri}_u$
rather than evolution of base flow profiles; likely either through a reduced ability of viscosity to mitigate shear instability through an increased Reynolds number, or through a reduced ability of buoyancy forces to suppress shearing, through a reduced Richardson number.
Time dependence of key dimensionless parameters obtained from planar-averaged simulation profiles.
Figure 7. Long description
Four panels (a-d) arranged in a 2 × 2 grid, each with x-axis t ranging from 0 to 600 s. Panel (a): y-axis is Re_u (velocity-based Reynolds number) on a logarithmic scale, ranging from approximately 10^1 to 10^4. The single black curve grows rapidly from near zero at t = 0 following an approximately power-law increase, and appears to be approaching but not yet reaching a plateau by t = 600 s. Panel (b): Dual y-axes. The left y-axis (black, solid curve) shows U_m/u_b (peak velocity normalised by a buoyancy velocity scale), ranging from 0.00 to 0.75, increasing approximately linearly with time. The right y-axis (blue, dashed curve) shows (δ_m/H)^2 (squared ratio of boundary current thickness to domain height), ranging from 0.00 to 0.09, also increasing roughly linearly. Panel (c): y-axis is Re_b (buoyancy-based Reynolds number) on a logarithmic scale, ranging from approximately 10^0 to 10^2. The single black curve decreases monotonically from a high value near t = 0, decaying smoothly over the full simulation period. Panel (d): y-axis shows two thickness ratios, ranging from approximately 0.35 to 0.45. The solid black curve shows δ_b/δ_m (ratio of buoyancy layer thickness to boundary current thickness), which fluctuates around approximately 0.42-0.43 with a slight upward trend and a small step-like increase near t ≈ 500 s. The dashed blue curve shows δ_u/δ_m (ratio of velocity layer thickness to boundary current thickness), which is nearly flat around 0.38-0.39 throughout, with minor fluctuations. Both ratios remain close and approximately constant.
The time dependence of these parameters is reported in figure 7. First, note that the ratios
$\delta _b / \delta _m$
and
$\delta _u/\delta _m$
are approximately constant over the region of the time series where profiles are self-similar, fluctuating about values of
$\delta _b / \delta _m = 0.41$
and
$\delta _u/\delta _m =0.375$
. Secondly, note that both
$U_m / u_b$
and
$\delta _m / H$
have a clear and trivial time dependence, with the velocity scale growing linearly with time and the boundary current width growing like the square root of time. As a result, we see that
$ \textit{Re}_u \sim t^{3/2}$
and
${Ri}_u \sim t^{-3/2}$
. The result of this is that both the Reynolds number and the Richardson number change by more than two orders of magnitude in the range
$0 \lt t \lt 600$
s.
Figure 7 shows that, when instability develops at
$t \approx 340$
s,
$ \textit{Re}_u$
has reached approximately
$5000$
. More interestingly,
${Ri}_u \approx 1$
at
$t \approx 340$
s, indicating that transition to instability occurs due to increased shear relative to suppressive buoyancy forces. To separate the influence of these two parameters we reformulate the vTG eigenvalue problem, instead assessing the linear stability of the rescaled time-invariant background state,
$U^* = U/U_m$
and
$B^* = {Ri}_u B$
. While the values of
$ \textit{Re}_u$
and
${Ri}_u$
arise naturally through flow evolution, here, we explore the wider parameter space to understand their role in instability. The resultant dimensionless eigenvalues and eigenvectors are denoted
$\lambda ^* = \sigma ^* - i \omega ^*$
,
$\hat {v}^*$
and
$\hat {b}^*$
respectively. The growth-rate dependence of this rescaled system on
$ \textit{Re}_u$
and
${Ri}_u$
is reported in figure 8(a), where the peak-finding algorithm has been adopted to find the maximum growth rate (
$\text{max}(\sigma ^*)$
) and associated scaled wavenumber
$k^*_x$
for each of the marked points. The dashed line shows the parameter values that arise through simulations, combining data of figure 7. The unstable region of the phase space is bounded by
$ \textit{Re}_u \gtrsim 250$
and
${Ri}_u \lesssim 1$
. For
$ \textit{Re}_u\gtrsim 1000$
the Reynolds number dependence is reasonably weak, although note that the flow is most unstable for
$ \textit{Re}_u \approx 1000$
, indicating that the mode is destabilised by viscosity, a behaviour similar to that of plane Poiseuille flow (Smyth & Carpenter Reference Smyth and Carpenter2019). There is a clear boundary between stable and unstable flow at
${Ri}_u \approx 1$
, indicating that the balance between shearing and buoyancy forces plays a key role in the dynamics, with the maximum growth rate increasing with decreasing
${Ri}_u$
. The region where the simulation crosses into the unstable regime is fairly insensitive to the Reynolds number, but highly sensitive to
${Ri}_u$
. For most of the flow development the Reynolds number is high enough to enable unstable growth (figures 7 and 8), indicating that it is the crossing of the
${Ri}_u$
threshold that finally enables flow transition.
(a) Reynolds and Richardson number dependence of the maximum growth rate arising through the appropriately scaled vTG solutions, with parameters defined in (3.5). The dashed line represents parameter values arising through 2-D simulations. The maximum growth rates of panel (a) are obtained using the peak-finding algorithm. Panel (b) shows the wavenumber and Richardson number dependence of the maximum growth rate for
$ \textit{Re}_u = 1000$
.
Figure 8. Long description
Two panels (a) and (b) sharing a common y-axis and colourbar, showing the maximum normalised growth rate max(σ*) from viscous Taylor-Goldstein (vTG) stability solutions. Shared colourbar: max(σ*) ranging from 0.00 (pale green, stable) through yellow and pink to approximately 0.15 (dark teal/black, most unstable). Shared y-axis: Ri_u (Richardson number) on a logarithmic scale, ranging from 10^-1 to 10-2, with instability concentrated at low Ri_u and stability (pale green, no unstable modes) at high Ri_u. Panel (a): x-axis is Re_u (Reynolds number) on a logarithmic scale from 10-2 to 10-4. Each point represents the maximum growth rate over all wavenumbers for a given (Re_u, Ri_u) pair. Unstable modes (warm-coloured points) are confined to low Ri_u (roughly Ri_u ≲ 1) and are absent at high Ri_u regardless of Re_u. At fixed low Ri_u, growth rates increase with Re_u, with the highest values (darkest points, max(σ*) ≈ 0.15) appearing at Re_u ≳ 10^3 and Ri_u ≲ 0.2. A dashed black line runs diagonally from upper-left to lower-right across the panel, representing the locus of (Re_u, Ri_u) values realised in the 2-D simulations; this line passes through the unstable region at low Ri_u and high Re_u and into the stable region at lower Re_u and higher Ri_u. Panel (b): x-axis is k_x* (normalised streamwise wavenumber), ranging linearly from 0 to 4, at fixed Re_u = 1000. Each point shows the maximum growth rate for a given (k_x*, Ri_u) pair. Unstable modes are again confined to Ri_u ≲ 1. At low Ri_u, instability spans a broad range of wavenumbers (roughly k_x* ≈ 0.5–3), with the highest growth rates centred around k_x* ≈ 1–2. Growth rates decrease toward both low and high wavenumber limits, consistent with a band-limited instability. At Ri_u ≳ 1 the flow is stable across all wavenumbers shown
For clarity the
$k_x^*$
dependence of the maximum growth rate is also reported in figure 8(b), where the Reynolds number is fixed to
$ \textit{Re}_u = 1000$
. This clearly shows a single peak in growth rate at
$k_x^* \approx 1.8$
.
The structure of the dominant mode is reasonably insensitive to Reynolds and Richardson numbers, so long as their values fall in the unstable regime. The mode is visualised in figure 9 for parameter values of
$ \textit{Re}_u = 10\,000$
,
${Ri}_u = 0.5$
and
$k^*_x = 1.36$
. Figure 9 also directly compares LSA eigenfunctions with simulation perturbations for case ISO_Base during the early stages of unstable growth at
$t=570$
s. There is clear strong agreement between the two datasets, providing further evidence that LSA is a valuable tool for assessing the primary instability of this flow. It is important to stress that (as with other figures) the vertical axes of figure 9 have been stretched to better visualise the strongly sheared modal structure near the upper surface. As a result, the lower vortical structures beneath the lower critical layer are more strongly affected by shear than presented in the figure. This is important when interpreting the nature of this instability.
Visualisation of the dominant unstable mode, obtained through LSA, for
$ \textit{Re}_u = 10\,000$
,
${Ri}_u = 0.5$
and
$k^*_x = 1.36$
. Vorticity contours overlayed by arrows indicating flow direction, scaled and coloured by the velocity magnitude are shown in panel (a), while the streamwise and vertical velocity components are shown in panels (b) and (c). Panel (d) shows the buoyancy structure. Note that all quantities are scaled by respective maxima. The horizontal dashed lines represent critical levels where the wavespeed
$c^* = \omega ^* / k_x^* = U/U_m$
. Modal structure is reported over one wavelength; note the stretched scaling of the
$y$
-axes to better visualise the strongly sheared near-wall structures. Numerical simulation data perturbations for case ISO_Base at
$t=570$
s, taken as instantaneous values with planar means subtracted, are reported in panels (e–h): spanwise vorticity (e), streamwise velocity (f), vertical velocity (g) and buoyancy (h). For reference, critical levels of LSA solutions overlay all simulation perturbations.
Figure 9. Long description
Eight pseudocolour/contour panels arranged in a 2×4 grid. The top row (a-d) shows the dominant unstable mode. The bottom row (e-h) shows the corresponding perturbation fields from DNS (case ISO_Base, t = 570 s, instantaneous with planar mean subtracted). Columns show the same physical quantity in LSA and DNS respectively. All quantities are normalised by their respective maxima. Both rows share a common y-axis: (z – H)/δ_m, ranging from 0.00 at the top (wall) to –1.00 at the bottom, with a stretched scale emphasising the near-wall region. The x-axis in all panels is x/δ_m*, ranging from 0.0 to approximately 3.0 (one wavelength). A horizontal dashed line in each panel marks the critical level where the wave phase speed c* = ω*/k_x* = U/U_m. Four colourbars (top of figure): Ω’_y(x,z): diverging red-white-blue, ±0.7 (spanwise vorticity perturbation) u’(x,z): diverging yellow-white-teal, ±0.9 (streamwise velocity perturbation) w’(x,z): diverging teal-white-purple, ±0.9 (vertical velocity perturbation) b’(x,z): diverging blue-white-red, ±0.9 (buoyancy perturbation) Panel (a and e), Ω’_y: Spanwise vorticity contours overlaid with velocity direction arrows scaled and coloured by velocity magnitude. Strong red (positive) vorticity is concentrated in a thin layer at the wall (z – H ≈ 0), with alternating red and blue patches just below, centred near the critical level. The arrows reveal a recirculating cat’s-eye flow structure straddling the critical level. Panel (b and f), u’: Streamwise velocity perturbation shown as filled contours. A pair of lobes of opposing sign are arranged vertically, where the maximum u’ displacement is situated at the lower critical level, and the lobe is skewed towards positive x with increasing height, so the lobes appear vertically stacked with a lateral tilt. A further but less intense sign reversal of the structure at greater depths below the critical level is seen. Panel (c), w’: Vertical velocity perturbation. A single dominant lobe of large magnitude (deep purple/teal) is centred near, but below the critical level and the mid-domain, flanked by lobes of opposite sign. Above the lower critical level the structure is sheared strongly towards positive x giving the contours a pronounced lateral tilt relative to the more symmetric structure below. Panel (d) – LSA, b’: Buoyancy perturbation. Similar lobe structure to u’, with opposing-sign regions above and below the critical level. The buoyancy perturbation is concentrated near each of the critical levels, with rapidly decaying amplitude away from both. Top and bottom panels are visually indistinguishable, showing the excellent agreement between DNS and LSA.
These parameter values lead to an unstable mode with
$\sigma ^* = 0.074$
and
$\omega ^* = -0.94$
, although we stress that the dominant modal structure is reasonably insensitive to parameter values so long as the mode is unstable. Critical levels are also shown, where the wavespeed of the mode matches the base flow speed, and indicate regions we may expect wave breaking. We see that these critical levels approximately denote the regions of strongest buoyancy and streamwise velocity perturbations. Vortical structures are most clear beneath the critical levels which agree well with the 2-D simulation solutions where overturning is present (figure 2). However, the highest magnitude of velocity perturbations occurs in the region of highest mean velocity, bounded between the two critical levels with a peak at
$(z-H)/\delta _m \approx -0.2$
. The influence of shear is imprinted on the modes above the lowest critical level. Note also that the buoyancy and vertical velocity are out of phase.
3.2.2. Sensitivity of instability to governing parameters
To identify the ubiquity of these flow features to varying real-world conditions, a number of sensitivity simulations were carried out where the governing parameters (the slope,
$\theta$
, and Reynolds number,
$ \textit{Re} = u_{b} H/\nu$
) were adjusted. Two methods of altering the
$ \textit{Re}$
were applied. Firstly,
$\Delta \rho$
was doubled and halved, which in turn increases
$ \textit{Re}$
(via
$u_b$
) by factors
$\sqrt {2}$
and
$1/\sqrt {2}$
, respectively. Secondly,
$\nu$
was doubled and halved, with
$\kappa$
halved and doubled, respectively, such that the final governing parameter
$ \textit{Sc}$
is kept constant. Whilst the overall effect should be governed only by the overall result on
$ \textit{Re}$
, the two methods may give the reader better insight into the effect of differing boundary/far-field conditions (expressed via
$\Delta \rho$
), and scaling (expressed via
$\nu$
).
Snapshots of the overturning billow structures are presented in figure 10. Overall, the flow characteristics are robust across the range of forcing slopes and density differences tested (figure 10). As density difference changes, in each case as the flow becomes unstable, the reflected wisping structure is reproduced. Reducing the density difference (and
$ \textit{Re}$
) (figure 10
a–c), has the effect of allowing the density interface to develop over a longer time before going unstable, and therefore the density interface becomes thicker before instability begins. Overall, the effect on billow size (both in amplitude and wavelength) is to decrease as the density difference increases. Similarly, the effect of increasing viscosity (and therefore
$ \textit{Re}$
) also results in smaller (amplitude and wavelength) billows forming at earlier times at higher Reynolds numbers. A rough estimate of the relevant field-scale Reynolds number, assuming that flow speeds remain at approximately the same scale (Jenkins Reference Jenkins1991), with
$H = 100{-}1000$
m, and
$\nu$
the same as in ISO_Base, giving
$ \textit{Re} = 10^7{-}10^8$
. Although it is impossible to extrapolate the relatively small changes in
$ \textit{Re}$
explored here to these much higher values, it is reasonable to assume the vertical extent of the instabilities (as a proportion of total depth) would continue to decrease, with more rapid cycles of instability and re-stability.
Single time outputs of density from simulations ISO_0.5x_rho at
$t = 1030$
s (a, Supplementary material and movies are available at 2), ISO_Base at
$t = 1030$
s (b,d,h), ISO_2x_rho at
$t = 592$
s (c, Supplementary material and movies are available at 3), ISO_5_slope at
$t = 370$
s (e, Supplementary material and movies are available at 4), ISO_10_slope at
$t = 220s$
s (f, Supplementary material and movies are available at 5), ISO_2x_
$\nu$
at
$t = 835$
s (g, Supplementary material and movies are available at 6) and ISO_0.5x_
$\nu$
at
$t = 470$
s (i, Supplementary material and movies are available at 7).
Figure 10. Long description
Nine pseudocolour panels of density ρ (kg m^-3) in the x-z plane, arranged in three columns representing three sets of sensitivity experiments. All panels share a common colourbar ranging from light blue (low density, ∼990 kg m^-3) to dark purple (high density, ∼1000 kg m^-3), with x-axis spanning 0-1 m and z-axis spanning 0.25-0.30 m. Each column contains a sensitivity case above and below a Base simulation, allowing direct visual comparison. All panels show the flow at a time selected to capture active billowing at the density interface near z ≈ 0.285-0.295 m. Left column - stratification sensitivity (a-c): • Panel (a), 0.5ρ (Re/2): Halved stratification case. Two large, well-developed billows are visible, with broad overturning structures and significant lateral extent. The interface is strongly deformed. • Panel (b), Base: Base case. Higher billow count and smaller scale to (a), serving as the reference for both stratification and viscosity comparisons. • Panel (c), 2ρ (2 Re): Doubled stratification case. More numerous, smaller-amplitude billows compared to (a) and (b), consistent with a shorter instability wavelength at higher stratification. Centre column - slope angle sensitivity (d-f): • Panel (d), Base: Identical to panel (b); reproduced as the reference for slope comparisons. • Panel (e), 5°: Five-degree slope case. Billows are present but appear somewhat more irregular and compressed compared to the Base. • Panel (f), 10°: Ten-degree slope case. Billows are again clearly visible, with an even earlier onset time. The structures appear more tightly spaced and the interface more sharply defined, consistent with more vigorous shear at this slope. Right column - viscosity sensitivity (g-i): • Panel (g), 2ν (0.5 Re): Doubled viscosity (halved Re) case. Billows are present but appear somewhat smoother and larger in scale relative to the Base. • Panel (h), Base: Identical to panels (b) and (d); reproduced as the reference for viscosity comparisons. • Panel (i), 0.5ν (2 Re): Halved viscosity (doubled Re) case. Billows are visible and appear more sharply defined with finer internal structure compared to the Base.
Changing slope steepness has a similar effect, although the relative contributions of changing gravity forcing on
$u$
and
$w$
becomes evident by the 10
$^\circ$
simulation. Fundamentally the same instability is produced for the steeper slopes of 5
$^\circ$
and 10
$^\circ$
(figure 10
d–f), and as with increasing the density difference, increasing the slope reduces billow dimensions, but additionally also produces more deformed billows (figure 10
f). As the gravitational vector strays considerably from vertical, the system switches from one where density is able to stabilise the flow at small slopes, to one where density is destabilising and tending to overturning at steep slopes. Since these regimes with steeper slopes have been studied elsewhere (e.g. Carey & Gebhart Reference Carey and Gebhart1981; Cenedese & Gatto Reference Cenedese and Gatto2016; McConnochie & Kerr Reference McConnochie and Kerr2018), here, we focus only on shallower slopes, and note that slopes beyond
$5 ^\circ$
show a transition towards the convection-driven instability regime (Ke et al. Reference Ke, Williamson, Armfield, McBain and Norris2019).
Schematic representation of profiles producing the canonical Holmboe (a), KHI (b) and the paired instability presented here (c), with horizontal velocity profile in blue, density profile in black, and a sense of the vorticity in the background colours (marked with clockwise, CW, and counter-clockwise, CCW). Circular arrows indicate the effect of the vorticity on shear interfaces, and the relative depths of shear and density interfaces are marked via the vertical arrows in blue and black for velocity and density respectively for each panel.
Figure 11. Long description
Three side-by-side schematic panels labelled Holmboe (left), Kelvin-Helmholtz (centre), and Mixed-Mode (right), each depicting the vertical profiles and vorticity structure associated with the named instability type. All panels share a common y-axis labelled z, ranging from 0 at the bottom to H at the top, with H/2 marked on the Holmboe panel. A shared colourbar on the right indicates background vorticity: clockwise (CW, red/pink) at the top and counter-clockwise (CCW, teal/dark) at the bottom, with white/neutral in between. In each panel: The blue curve shows the horizontal velocity profile (S-shaped, representing a shear layer). The black curve shows the density profile (also S-shaped, representing a pycnocline). Grey circular arrows indicate the sense of rotation induced by the vorticity at the shear interfaces. Vertical double-headed arrows mark the depth of the shear interface (blue arrow) and the density interface (black arrow), illustrating their relative vertical positions. Holmboe (left): The density interface (black profile) is sharp and centred at the centre of the shear layer (blue profile), which is broader. The black vertical arrow is short, indicating a thin density interface, while the blue arrow is longer, indicating a thicker shear layer. The two interfaces are vertically offset, with the density interface sitting near z = H/2 and the shear layer centred higher. The background vorticity is weak and diffuse. Circular arrows appear on both sides of the density interface, indicating the counter-propagating wave interaction mechanism characteristic of Holmboe instability. The vorticity background is pale, reflecting weak shear-induced rotation. Kelvin-Helmholtz (centre): The density and velocity profiles are co-located, both centred at the same intermediate depth (z ≈ H/2). The blue and black vertical arrows are equal in length and aligned, indicating that the shear and density interfaces coincide. The background vorticity is strongly concentrated in a horizontal band at mid-depth, with intense red (CW) above and teal (CCW) below the interface, and fading to neutral away from it. Circular arrows are arranged symmetrically on either side of the interface, consistent with the single-interface KH roll-up mechanism. Mixed-Mode (right): The shear layer has a turning point at around the mid-point of the profile, and is asymmetric around it, with a broader profile beneath and a sharp one above it. The density profile (black) is sharp and at the centre of the profile. The blue vertical arrow is longer than the black, again indicating the shear layer is thicker than the density interface. The background vorticity shows a strong CW (red) band associated with the shear layer above, and a weaker CCW (teal) region below, mostly below the density interface, with the intense clockwise vorticity region offset upward from the density interface. Circular arrows appear both near the shear layer above and near the density interface below, indicating that both KH-type and Holmboe-type interaction mechanisms are active simultaneously.
3.2.3. Mechanisms of the shear instability
The stratified shear instability described above is unusual compared with the canonical instabilities described throughout most of the literature (figures 11
a and 11
b). Key to understanding the instability is the opposing directions of shear above and below
$z(u = U_{\textit{max}} )$
, with a stabilising density gradient throughout that velocity structure. This results in the sense of rotation of vortices each side of the velocity maximum being in opposing directions. Two interacting modes of instability are observed each side of the velocity maximum, and unlike canonical versions of mode interactions, these modes are both propagating in the same direction, due to the opposing directions of shear (figure 11
c). The wave resonance interpretation of instability precludes the instability of certain jet-like channel flows such as these (Smyth & Carpenter Reference Smyth and Carpenter2019), but as shown in Appendix C, the flows are able to go unstable where the jet is bounded by fluid with constant (including 0) velocity.
A further complicating factor is the presence of a solid upper boundary which directly intersects both shear and density interfaces. Given these complications, identifying the forms of these modes against descriptions and definitions based on idealised flows is difficult. Whilst a full diagnosis is beyond the scope of this study (and may not yield additional insight into the useful dynamics of the flow), some considerations to these definitions are considered here.
Invoking descriptions based on the form of the instability itself, rather than based on profile-based measures, Parker, Caulfield & Kerswell (Reference Parker, Caulfield and Kerswell2020) define KHI as shear instability involving overturning of the shear layer, whereas HWI involves propagating vortices either side of the shear layer. For the instability below the velocity maximum, vortices indeed form below the shear layer, indicating that instability is a HWI-type instability. For the instability above the velocity maximum, the shear continues throughout the layer, perhaps suggesting that the only available mode in this region is KHI indicating overturning of the shear layer itself.
Other work has used the ratio,
$R_i$
, of velocity variation length scale (
$h_{u,i}$
) to density variation length scale (
$h_{\rho , i}$
), as a diagnostic for Holmboe instability, specifically where
$R \gt 2$
indicates HWI is present (Alexakis Reference Alexakis2007; Carpenter et al. Reference Carpenter, Tedford, Rahmani and Lawrence2010b
). Split into each instability structure by the velocity maximum, the upper shear layer (
$i = 1$
)
$h_{u, 1} = h_{\rho , 1} = H - z(u = U_{\textit{max}} )$
, by definition returns
$R_1 = 1$
, whilst for the lower shear interface (
$i = 2$
),
$h_{u, 2} = \delta _u - h_{u, 1}$
and
$h_{\rho , 2} = \delta _\rho - h_{\rho , 1}$
, giving
$R_2 = 3.27$
. Again, this criterion indicates mixed instability types, with KHI at the upper interface, and HWI on the lower layer.
Finally, we consider a descriptive diagnostic based on the eigenfunctions of dominant modes (see figure 9, and Zhu et al. (Reference Zhu, Atoufi, Lefauve, Kerswell and Linden2024, figure 5). In figure 9, the modal structures beneath the velocity maximum show characteristics of the HWI, specifically two pairs of counter-rotating roll cells, centred at the lower critical level, whilst above the velocity maximum, alternating (strongly sheared) bands dominate the structure in both vorticity and density (figures 9 a, 9 e, 9 d and 9 h).
Application of the wave interaction-based diagnostic from Carpenter et al. (Reference Carpenter, Balmforth and Lawrence2010a ) may provide further details as to the modes of instability. However, despite remaining uncertainty as to the mode of the upper instability, this work identifies complex systems of interacting modes of primary stratified shear instabilities. Paired interactions affect the leading-order behaviours of the secondary instability, and subsequent turbulent mixing with the dynamics playing out differently to two unpaired stratified Couette flows at the point of secondary instability due to this pairing between vorticity interfaces.
3.3. Secondary instability of the flow and the transition to turbulence
3.3.1. Secondary instability and long-term evolution in two dimensions
Once billows begin to overturn, the dynamics can only be represented by time-evolving nonlinear simulations, so we return to qualitative descriptions of 2-D SPINS simulations in order to test the mechanisms by which the flow may undergo cyclic or marginal instability equivalent to those proposed by Smyth, Nash & Moum (Reference Smyth, Nash and Moum2019). Here, we define the system as undergoing cyclic instability if: (a) there is a cycle of instability, turbulent growth, turbulent decay and forcing towards instability, (b) multiple cycles of instability that occur in the same manner, indicated here by self-similarity of profiles during the forcing stage, meeting
$ \textit{Re}_u$
and
${Ri}_u$
thresholds at instability, and similarity of the unstable mode.
As the overturning pattern develops, further shorter waves develop on the interface between the counter-rotating wisps/billow (figures 2 h, 3 i and 3 e, 3 f). Specifically, as these secondary waves form, the anticlockwise vorticity of the lower layer drives shedding of parcels of upper layer fluid downwards (across the leading edge of the lower billow) (figures 2 j–l, 3 g–i). Secondary instability appears to form from these secondary waves – smaller-scale billowing structures, that drive the true degeneration into a state where turbulent mixing is active (figures 2 k–l and 3 h–i).
The long-term evolution of ISO_Base, showing Hovmöller plots of density (a), and vorticity (b) as well as evolution of
$ \textit{Re}_u$
and
${Ri}_u$
(c). Hovmöller plots are each taken for a vertical profile at the mid-point of the domain. Lines in b indicate timings of panels shown in figures 2 and 3, the reference line in (c) shows
${Ri}_u = 1$
, with two dimensions (solid line) and three dimensions (dashed line) shown. Panel (d) shows evolution of TKE for ISO_Base simulation over the initial cycle of instability (solid line) and subsequent evolution (dotted line).
Figure 12. Long description
Four stacked panels sharing a common x-axis of dimensional time t [s] from 0 to 1900 s. Panel (a) – Hovmöller plot of density ρ (kg m^3): z-axis spans 0.24-0.30 m. Colourbar ranges from ∼990 (light blue) to 1000 (dark purple) kg m^3. From t = 0 to approximately t = 450 s the density field is quasi-steady: a sharp horizontal pycnocline near z ≈ 0.285-0.290 m separates light fluid above from dense fluid below, with little temporal variation. Around t ≈ 500-600 s a marked transition occurs: the pycnocline broadens and lightens visibly (lighter colours appear at mid-depth), indicating vertical mixing and erosion of the density interface. After t ≈ 700 s the density field partially restratifies, with the interface reforming but at a lower density contrast, and remains broadly steady through to t = 1900 s with only slow evolution. Panel (b) - Hovmöller plot of vorticity ζ (s^1): z-axis spans 0.24-0.30 m. Colourbar ranges from –5 (teal, CCW) to +25 (dark purple, CW) s^1. Prior to t ≈ 450 s the vorticity is concentrated in a thin band near the upper boundary (z ≈ 0.29-0.30 m), with near-zero values elsewhere, consistent with a laminar boundary current. Short vertical tick marks along the top of the panel indicate the timings of snapshots shown in companion figures. Around t ≈ 500-650 s a dramatic broadening of the vorticity field occurs: intense vorticity spreads across a wide vertical range, with strong positive values filling much of the water column, indicating the onset and development of instability and turbulent mixing. After t ≈ 700 s the vorticity field contracts back toward the upper boundary but remains broader and more structured than the pre-instability state, slowly reorganising through t = 1900 s. Panel (c) - Re_u and Ri_u time series: Dual y-axes on a logarithmic scale. The black curve (left axis) shows Re_u growing monotonically from ∼10^2 at t = 0 to ∼10^5 by t ≈ 500 s, before the instability onset causes a slight inflection; it then continues to evolve at high values. The orange curve (right axis) shows Ri_u decreasing from high values (well above 10^0) at early times, crossing the horizontal reference line at Ri_u = 1 (marked by a thin horizontal line) at approximately t ≈ 400-450 s, then dropping further before the instability event causes it to spike sharply and irregularly around t ≈ 500-650 s. After the mixing event Ri_u recovers and stabilises above 1. Solid lines represent the 2-D simulation and dashed lines the 3-D simulation; the 3-D simulation peaks quicker around the instability event. Panel (d) – Turbulent kinetic energy (TKE, J m^-3): y-axis on a logarithmic scale from ∼10^-10 to ∼10^-5 J m^-3. TKE remains at very low (near-zero) levels from t = 0 until approximately t ≈ 450-500 s, then rises sharply by several orders of magnitude over a short interval, peaking near t ≈ 600-650 s at ∼10^-5 J m^-3. This is shown as a solid line, representing the initial instability cycle. After the peak, TKE decays irregularly (shown as a dotted line for the subsequent evolution), fluctuating around intermediate values and slowly declining through t = 1900 s, consistent with intermittent turbulence following the primary mixing event.
The first criterion appears to be met for at least these two cycles (figures 12
a and 12
b). Once turbulence has been produced from the secondary instability, turbulent diffusion (mixing at the small scales) smooths shear and buoyancy gradients, so that the production of turbulence by the transfer of kinetic energy from the mean flow decreases and dissipation dominates (figures 2
m, 3
j and 12
b, 12
d). This decaying of turbulence over around
$150$
–
$200$
s allows the buoyancy forcing at the upper boundary to re-stratify the fluid, producing renewed forcing for the mean flow to grow. Over time, the re-stratifying density boundary condition maintains a strong top to bottom density difference, and the fluid re-forms to a stratified boundary current state resembling earlier times (figure 13
a–e). The flow develops to instability once again (figure 13 f-i) in a very similar manner to that described in
$\S$
§ 3.2 and 3.3. Throughout this second cycle of instability, profiles are broadly self-similar and meet similar criteria of
$ \textit{Re}_u$
and
${Ri}_u$
(figure 12
c), and meet the second criteria for cyclic instability.
Over time, this 2-D mixing process is an efficient means of vertical mixing, and the pycnocline thickens substantially, exerting an upslope buoyant forcing throughout the water column. A more diffuse lower shear region forms after the second cycle of this instability after
$t \approx 1200$
s (figure 12
a). Once this dynamics begins to be affected by the lower boundary, the stability criteria
${Ri}_u$
and
$ \textit{Re}_u$
based on
$\delta _u$
is no longer applicable, since
$\delta _u$
is now bounded by the lower boundary. It is possible that cyclic behaviour could continue beyond this time if the vertical extent of the domain were longer, but in these simulations, the instabilities stop due to constraint of the flow by the lower boundary, and the dynamics deviate from the geophysically relevant problem.
As in figure 2, but for the second instance of the instability forming. (a–d) show scaled forms of the respective profiles at times
$t = 570$
s (black) and
$t = 894$
s (blue) for assessment of self-similarity in later iterations. Panels e-n show times
$t = 850$
,
$894$
,
$914$
,
$934$
,
$954$
,
$966$
,
$994$
,
$1014$
and
$1054$
s respectively.
Figure 13. Long description
The figure mirrors the structure of the first instability figure, with four profile panels (a-d) at the top and ten density pseudocolour panels (e-n) below, arranged in two columns of five. The density colourbar is shared across all pseudocolour panels, ranging from ∼990 (light blue) to 1000 (dark purple) kg m^3, with x spanning 0.1-0.9 m and z spanning 0.25-0.30 m. Panels (a-d) – vertical profiles: The common y-axis is (z – H)/δ_m, ranging from 0 to approximately –2.5, noting the extended range relative to the first instability cycle, reflecting a thicker boundary current at this later time. Two curves are shown in each panel: black at t = 570 s (from the first instability cycle, for self-similarity comparison) and blue at t = 894 s (pre-instability state of the second cycle). Panel (a), U/U_m: Both profiles show the same S-shaped boundary current structure. The blue curve is broader than the black, extending to greater normalised depth, but the overall shape is similar, suggesting approximate self-similarity between the two cycles. Panel (b), U_z δ_m/U_m: Vertical shear profiles. Both curves peak near the surface and decay with depth; the blue curve is slightly broader, consistent with the thicker boundary layer. Panel (c), U_zz δ_m^2/U_m: Second derivative profiles, showing the curvature structure of the velocity profile. Both curves exhibit a sharp positive peak near the wall and a negative lobe below, again broadly similar between the two times. Panel (d), B: Buoyancy profiles. Both curves show a monotonic decrease from maximum at the wall to near-zero at depth. The blue curve decays more gradually, indicating a somewhat more diffuse buoyancy interface at the later time. Panels (e-n) – density snapshots: Ten pseudocolour panels showing the evolution of the density field through the second instability cycle, at times t = 850, 894, 914, 934, 954, 966, 994, 1014, and 1054 s (panels e-n respectively, with one time shared across both columns). Panel (e), t = 850 s: Near-laminar state. A thin white isopycnal contour near z ≈ 0.285 m is nearly flat, with only very slight undulations. The density field is largely undisturbed. Panel (f), t = 894 s: Long-wavelength interfacial waves are developing. The isopycnal shows gentle, coherent undulations across the domain. Panel (g), t = 914 s: Wave amplitude has grown; the interface shows clear sinusoidal deformation with a well-defined wavelength. Panel (h), t = 934 s: Billows are beginning to form. The interface is strongly deformed and the first signs of overturning appear. Panel (i), t = 954 s: Active KH billowing. Large overturning structures are visible in the left portion of the domain, with a rolled-up billow clearly evident. Panel (j), t = 966 s: Billows are well developed across the domain, with multiple overturning structures of comparable scale visible. Panel (k), t = 994 s: Billows are breaking down. The organised structures are fragmenting and the density field shows increased small-scale complexity and mixing. Panel (l), t = 1014 s: Advanced turbulent mixing. The interface is highly disrupted, with complex interleaving of dense and light fluid and filamentary structures throughout. Panel (m), t = 1054 s: Post-mixing, partial restratification. Large-scale billow structures have collapsed; the density field is disordered but beginning to reorganise, with a diffuse interface reforming near the upper boundary.
Early 3-D evolution of the secondary instability (centre, right) compared with the 2-D simulation ISO_Base at the same time (a,d,g). Shown for a
$x{-}z$
slice for density (b,e,h) where contours representing the density isosurfaces shown in c,f and i are added in white. Top panels show
$t = 604$
s middle panels show
$t = 610$
s and lower panels show
$t = 618$
s. Note
$z$
axis exaggerated from previous figures.
Figure 14. Long description
Nine panels arranged in a 3×3 grid. Rows correspond to times t = 604 s (top), t = 610 s (middle), and t = 618 s (bottom). Columns show: the 2-D simulation ISO_Base x-z slice (left), the 3-D simulation x-z slice with isosurface contour lines overlaid in white (centre), and a 3-D perspective volume rendering of density isosurfaces (right). The z-axis spans 0.24-0.30 m, exaggerated relative to earlier figures, and x spans 0-0.3 m in the 2-D/slice panels. The shared colourbar ranges from ∼990 (light blue) to 1000 (dark purple) kg m^3. The 3-D panels use two isosurface colours: teal/blue for the light (low density) isosurface and dark red/maroon for the dense isosurface, with semi-transparency revealing internal structure. Top row – t = 604 s: Panel (a), 2-D x-z slice: A large, well-developed billow dominates, with a tightly wound spiral of white density contours in the centre-left of the domain. Light blue fluid (low density) occupies the upper portion, with the billow core showing intermediate densities. The structure is smooth and coherent, characteristic of the 2-D roll-up. Panel (b), 3-D x-z slice with white isosurface contours: Very similar billow structure to (a), with the spiral contours closely matching the 2-D case. The density field in the slice is nearly identical, indicating that at this early time the 3-D simulation has not yet diverged significantly from 2-D behaviour in this plane. Panel (c), 3-D isosurface rendering: Two smooth, coherent isosurfaces are visible. The teal surface (light fluid) forms a broad, gently curved sheet near the top of the domain. The dark red surface (dense fluid) wraps into a clean, elongated spiral roll – the billow core – that is uniform across the spanwise (y) direction, confirming the essentially 2-D character of the instability at this time. Middle row – t = 610 s: Panel (d), 2-D x-z slice: The billow has tightened and the spiral is more compact. The light blue region at the top has been partially entrained into the roll, and the density gradients within the billow core are steeper. The overall structure remains a coherent, smooth roll. Panel (e), 3-D x-z slice with white isosurface contours: The billow structure in the slice is similar to (d), but subtle differences in the contour positions are beginning to emerge, hinting at the onset of 3-D effects. Panel (f), 3-D isosurface rendering: The isosurfaces show the beginning of spanwise (y-direction) deformation. The teal surface remains broadly sheet-like but shows slight corrugation. The dark red spiral isosurface, while still recognisably a roll, exhibits visible waviness along the y-axis, indicating the growth of a spanwise secondary instability on the primary billow. Bottom row – t = 618 s: Panel (g), 2-D x-z slice: The 2-D billow remains coherent, with a tightly wound spiral core and smooth density gradients. The 2-D simulation shows no sign of breakdown. Panel (h), 3-D x-z slice with white isosurface contours: In strong contrast to (g), the 3-D slice shows dramatic disruption. The lower portion of the domain is filled with fine vertical streaks and irregular density structures, indicating the rapid onset of 3-D turbulent breakdown. The upper billow region retains some coherence but the interior is highly disordered. Panel (i), 3-D isosurface rendering: The 3-D structure has broken down dramatically. The teal isosurface at the top remains partially intact as a fragmented sheet. The dark red isosurface has completely lost its coherent spiral form, fragmenting into a dense, irregular tangle of filamentary structures filling much of the lower domain – indicating the transition to three-dimensional turbulence through the rapid amplification of the spanwise secondary instability seen nascently in panel (f).
3.3.2. Secondary instability and long-term evolution in three dimensions
To test the effects of spanwise variability on the instability, a 3-D extension simulation of ISO_Base was conducted. The outputs of ISO_Base at
$t=570$
s, just prior to the primary instability forming, were extended into the spanwise direction, where
$L_y = 0.128$
m, and
$N_y = 128$
, with white noise added in the spanwise direction to velocity fields to trigger any instabilities (in the same manner as at initialisation for 2-D simulations). In addition, the
$x$
domain was reduced by one billow wavelength, so that
$L = 0.66$
m and
$N_x = 684$
, to optimise the computational efficiency. The evolution of the primary instability was confirmed to occur without significant spanwise variability or 3-D effects (figure 14
a–c), so 2-D simulations are sufficient until this stage. Spanwise variability only emerges once the upper and lower billows are significantly overturning and the flow is undergoing secondary instability, at which point the 2-D and 3-D simulations rapidly diverge (figures 14
d and 14
e). In contrast to the continued overturning and formation of waves on the interface between the billows observed in 2-D simulations (figure 14
d), in 3-D simulations, the billows themselves undergo a transition towards a turbulent state (figures 14
e, 14
f, 14
h and 14
i), in line with theory on 3-D overturning. Small scale structures emerge, but unlike in the 2-D case, these are not cast downwards into the ambient fluid, instead the turbulence decays over a relatively rapid time scale (figure 15
e).
Late 3-D evolution of the secondary instability. Shown for a
$x{-}z$
slice for density (c–f) and horizontal velocity (g–j). Corresponding profiles of horizontally averaged density (a) and along-slope velocity (b) are shown. Panels (c–f) and (g–j) and the profiles in increasing darkness are at times
$t =$
625, 645, 655 and 688 s respectively. Dashed lines in (a, b) show reference profiles from
$t = 560$
s (as in figure 5).
Figure 15. Long description
Ten panels showing the late-stage 3-D evolution of the secondary instability, combining profile plots (top) with pseudocolour x-z slices (bottom two columns). All slice panels share x spanning 0.1-0.7 m and z spanning 0.25-0.30 m. Four times are shown – t = 625, 645, 655, and 688 s – represented by increasing line darkness in the profiles and panels (c)/(g) through (f)/(j) respectively. Panel (a) – Horizontally averaged density profile: z-axis spans 0.24-0.30 m; x-axis shows ρ from 990-1000 kg m^3. Four solid curves of increasing darkness show the evolving mean density profile; a dashed curve shows the reference profile at t = 560 s (pre-instability, from figure 5). At t = 625 s (lightest solid) the profile still shows a relatively sharp pycnocline near z ≈ 0.285 m, though already broader than the dashed reference. With increasing time the pycnocline progressively broadens and the density gradient weakens across a growing depth range, indicating irreversible diapycnal mixing. By t = 688 s (darkest solid) the profile is considerably more diffuse, with the density transition spread over nearly the full z range shown. Panel (b) – Horizontally averaged along-slope velocity profile: z-axis spans 0.24-0.30 m; x-axis shows u from –0.12 to 0 m s^1. The dashed curve (t = 560 s reference) shows a well-defined boundary current with peak velocity near z ≈ 0.290 m. The solid curves at the four later times show progressive broadening and weakening of the velocity maximum, with the profile becoming more diffuse at greater depths. The velocity peak shifts slightly and the shear layer thickens, consistent with turbulent momentum redistribution during the mixing event. Left column – Density x-z slices, panels (c-f): Colourbar ranges from 990 (white/light) to 1000 (dark red) kg m^3. Panel (c), t = 625 s: The density field shows strong turbulent disruption throughout the domain. Vertically oriented streaks and filaments of mixed-density fluid fill the interior, with no coherent large-scale billow structure remaining. The upper boundary retains a thin layer of light fluid. Panel (d), t = 645 s: The turbulent structure is somewhat less intense. Broad, irregular density variations persist across the domain but the fine vertical streaking has reduced, suggesting the turbulence is beginning to decay. Panel (e), t = 655 s: Further decay of turbulent structure. The density field is becoming more horizontally layered, with lighter fluid above and denser fluid below, though considerable disorder remains. Panel (f), t = 688 s: The density field has largely restratified. A diffuse but recognisable horizontal interface has re-established near z ≈ 0.28-0.29 m, with light fluid above and dense fluid below. Small-scale structure has largely dissipated. Right column – Horizontal velocity x-z slices, panels (g-j): Colourbar ranges from –0.08 (dark blue, strong negative flow) to +0.01 (yellow, weak positive) m s^1. Panel (g), t = 625 s: The velocity field shows significant spatial variability. The boundary current structure near the top is disrupted, with irregular patches of varying velocity magnitude across the domain. Moderate negative velocities dominate but with substantial horizontal inhomogeneity. Panel (h), t = 645 s: The velocity field is becoming more organised. A clearer horizontal stratification of velocity is re-emerging, with stronger negative velocities near the upper boundary and weaker velocities below, though patchiness persists. Panel (i), t = 655 s: Further reorganisation. The velocity field is broadly horizontally layered, with the boundary current signature strengthening near the top of the domain and near-zero velocities below. Panel (j), t = 688 s: The boundary current has largely re-established. Strong negative velocities (dark blue) are concentrated near z ≈ 0.285-0.30 m in a relatively coherent layer, with weaker velocities below, consistent with a reformed but somewhat thicker boundary current following the mixing event.
It is important to note that these 3-D simulations do not include a Coriolis acceleration term which would deflect the base flow (to the form shown by e.g. Jenkins Reference Jenkins2021) and likely result in more 3-D instabilities than those shown here. The more realistic mixing processes in 3-D simulations over a single cycle of instability result in a broadly similar process to their 2-D equivalents studied here. However, the improved representation of turbulent dissipation and mixing means the pycnocline does not thicken in the same manner, and shear remains constrained to the upper part of the water column (figure 15 b). Whilst it remains unclear how well the similarity of the instability is maintained through subsequent cycles of instability from just a single cycle of instability, this indication of suppressed long-term deepening of the pycnocline indicates cyclic instability is more likely to occur when the 3-D dynamics is considered.
4. Discussion
This flow with a re-stabilising buoyancy flux across a tilted interface, which produces a unique dynamics, is likely to be a key feature of geophysical settings where a stable boundary layer forms over surface topography, but possibly unique in the ocean to the sub-ice-shelf environment and highly influential in driving mixing of temperature and salinity. Similar flows exist in the idealised form of stratified plane Poiseuille flows (e.g. Lloyd et al. Reference Lloyd, Dorrell and Caulfield2022; Lloyd & Dorrell Reference Lloyd and Dorrell2024), and turbulent gravity currents in a channel (e.g. Cantero et al. Reference Cantero, Balachandar, Cantelli, Pirmez and Parker2009). Consideration of gravity currents by Wells, Cenedese & Caulfield (Reference Wells, Cenedese and Caulfield2010) discusses the two sources of turbulence in density currents as interfacial entrainment at the boundary between the different density layers (expressed in plume-type models as the entrainment ratio,
$E$
) and drag at the solid lower boundary (expressed through the drag coefficient,
$C_d$
). In the flows considered in this paper, it is clear that the hard cutoff with the no-slip upper boundary generates a shear layer that is interacting with the generation of turbulence and shear at the diffuse boundary to the ambient fluid, whilst the continuous gradient of density throughout a significant portion of these two shear layers produces different dynamics to that observed in the earlier studies. The combination of continuous and multiple interfaces along with interactions of the flow with solid boundaries (Holt Reference Holt1998; Baglaenko Reference Baglaenko2016; Liu, Kaminski & Smyth Reference Liu, Kaminski and Smyth2023) impact these instabilities, and as is shown in this paper, drive a new dynamics.
Whilst laboratory experiments and direct numerical simulations have been carried out for similar ice–ocean boundaries with steeply tilting ice (e.g. Kerr & McConnochie Reference Kerr and McConnochie2015; Gayen, Griffiths & Kerr Reference Gayen, Griffiths and Kerr2016; McConnochie & Kerr Reference McConnochie and Kerr2018; Mondal et al. Reference Mondal, Gayen, Griffiths and Kerr2019) (see recent review from McCutchan & Johnson (Reference McCutchan and Johnson2022)), studies for gently sloping ice surfaces are restricted to LES (e.g. Vreugdenhil & Taylor Reference Vreugdenhil and Taylor2019; Begeman et al. Reference Begeman, Asay-Davis and Van Roekel2022; Anselin et al. Reference Anselin, Holland, Jenkins and Taylor2024) which, due to underlying assumptions in the sub-grid-scale model, may not resolve this new shear dynamics. Across various approaches, the general form of vertical profiles of density and streamwise velocity are consistent (Jenkins Reference Jenkins2016; Mondal et al. Reference Mondal, Gayen, Griffiths and Kerr2019; Begeman et al. Reference Begeman, Asay-Davis and Van Roekel2022; Patmore et al. Reference Patmore, Holland, Vreugdenhil, Jenkins and Taylor2023), indicating the key research question is what the routes to mixing and turbulence are in such flows.
In this case, evolution of these mixing events is considerably different to those observed in the canonical shear instabilities of Holmboe or Kelvin–Helmholtz. The initial instability is Holmboe-like, with a continual scouring of the pycnocline, rather than complete overturning of the pycnocline (due to separation of buoyancy and vorticity interfaces). This supports the theory of continuous marginal instability, which cannot be sustained for vertical wall systems (Ke et al. Reference Ke, Williamson, Armfield and Komiya2023). However, in contrast with those earlier results, later development of the flow becomes more vigorous, and almost analogous to later stages of KHI-like instability, with a complete decay to a turbulent state (although via a different mechanism).
Previous work using 1-D analytic models incorporating rotation has identified a marginally stable state of the pycnocline (Jenkins Reference Jenkins2021; Anselin et al. Reference Anselin, Holland, Jenkins and Taylor2024) under-ice shelves, a feature that may be ubiquitous in geophysical settings (Salehipour, Peltier & Caulfield Reference Salehipour, Peltier and Caulfield2018; Smyth & Carpenter Reference Smyth and Carpenter2019). The criterion for this marginal stability for the ice-shelf–ocean boundary remains an open question, with deviation around a critical
$Ri_g$
gradient through the pycnocline suggested by Jenkins (Reference Jenkins2021), whilst Anselin et al. (Reference Anselin, Holland, Jenkins and Taylor2024) instead found the boundary current was governed by a critical shear (itself dependent on basal slope), without the same dependence on buoyancy gradients. Here, our analysis of the linear stability of profiles identifies the instability of these profiles is dependent on both
$ \textit{Re}_u$
and
$Ri_u$
. If the system over longer time scales is assumed to deviate around the point of instability, this analysis may point towards the interplay of a combination of the shear and buoyancy gradients.
The instability of the flow is found to be poorly predicted by
${Ri}_g \lt 1/4$
, instead LSA of the relevant profiles reveal patterns of instability based on
$ \textit{Re}_u$
and
${Ri}_u$
. In particular, a region with
${Ri}_g \ll 1/4$
will exist from very early in the simulation (once shear begins to develop) due to the thick shear layer compared with the buoyancy (controlled through
$ \textit{Sc}$
), whilst a region with
${Ri}_g \gg 1/4$
will also exist at the velocity maximum (where shear is zero). The use of LSA enables the identification of such regimes, whilst its use is validated both by excellent agreement between vertical structures (figure 9) and growth rates (figure 4). The validity of LSA is confirmed by the self-similarity of profiles during the flow evolution, and the relative insensitivity of the structure of the dominant mode to Reynolds and Richardson numbers, meaning the flow profiles are effectively steady state, with a time dependence felt through
$ \textit{Re}_u$
and
${Ri}_u$
.
Salehipour et al. (Reference Salehipour, Peltier and Caulfield2018) argue that marginal instability in strongly stratified flows is more likely linked to the continuous, ‘slow burn’ features of HWIs, than the episodic ‘flaring’ KHIs. Subsequent work calls into question the exact mechanisms by which marginal instability can be achieved for stratified shear (Zhou Reference Zhou2022; Stastna et al. Reference Stastna, Bhavsar, Hartharn-Evans and Castro-Folker2025). These simulations identify further mechanisms by which a marginally stable state might be achieved, via a combination of slower, more continuous scouring in the primary instability, followed by episodic overturning mixing in the secondary instability. We show evidence of cyclic instability where energetic turbulence diffuses away the sharp velocity gradients thus increasing
${Ri}_u$
, and re-starting the cycle. It is important to note that periodic simulations may be more susceptible to cyclic behaviour than large domains relevant to geophysical settings (Smith, Caulfield & Taylor Reference Smith, Caulfield and Taylor2021; Vieweg & Caulfield Reference Vieweg and Caulfield2026), and further work is required to understand the sensitivity of mechanisms for cyclic behaviour presented here.
Now the fundamental dynamics is understood in the simplest case, future work is needed to identify the role of various processes present in the real world ice shelves. Firstly to consider the sensitivity of simulations to varying
$ \textit{Sc}$
from the
$ \textit{Sc} = 7$
used here for computational efficiency to higher
$ \textit{Sc}$
representative of salt, the main stratifying agent under-ice shelves. Only when two properties (temperature and salinity) with widely differing diffusivities are considered can simulations include the effects of double diffusion, which can be a significant factor in ice-shelf melting regimes under particularly quiescent conditions (Middleton et al. Reference Middleton, Vreugdenhil, Holland and Taylor2021, Reference Middleton, Davis, Taylor and Nicholls2022; Rosevear et al. Reference Rosevear, Gayen and Galton-Fenzi2021, Reference Rosevear, Galton-Fenzi and Stevens2022), in turn defining the profiles of buoyancy forcing. Furthermore, planetary rotation exerts a strong control over the development of the shear profile. Rotation of the velocity vector leads to non-zero shear at the speed maximum, and is therefore likely to produce truly 3-D instabilities that are not represented by simulations herein. Meanwhile viscous effects are limited to the Ekman layer, beyond which the flow is close to thermal wind balance. Incorporating such effects would be critical to compare simulations analogous to these with past LES studies.
5. Conclusions
In this paper, we identify a complex paired instability that emerges and drives the transition to turbulent mixing under-ice shelves using high-fidelity numerical simulations. This flow exhibits features of both the slowly evolving HWI, with scouring-based mixing typically associated with marginally stable flows; and of the flaring KHI associated with short lived, but vigorous overturning-based mixing. Such features emerge after a long phase of slowly evolving, diffusion-based flow, and are well predicted by LSA of profiles, with the onset of instability predicted by a combination of
${Ri}_u$
and
$ \textit{Re}_u$
, rather than
${Ri}_g\lt 0.25$
used for many other stratified shear flows. Cyclic behaviour in the instabilities forming, mixing by turbulence, decay of turbulence, and then flow building back to instability is evidenced in 2-D simulations by changes in
${Ri}_u$
(figure 12); instabilities develop for
${Ri}_u\lt 1$
which reduce shear and increase
${Ri}_u$
above 1, where turbulence decays and flow begins to sharpen again.
As the flow evolves, increasing levels of complexity are required in the modelling approach. The early evolution of the flow is entirely one-dimensional, and the initial instability is well predicted by LSA. That analysis provides valuable insight into the regimes of instability, and the mechanisms driving instability, but the flow rapidly evolves outside the range of validity for small-amplitude perturbations, at which point time-resolved numerical simulations in two dimensions are required. Such simulations are effective at simulating the initial transition from a 1-D to 2-D state under the primary instability, but once overturning has occurred, 3-D secondary instabilities form, requiring 3-D simulations. Such a combined approach of LSA of 1-D profiles, 2-D simulations and extensions of these into the third dimension provide a valuable combination of insights into the behaviour of these flows. Building on the approach in this manuscript, further valuable insights into these systems could be identified via a LSA generalised to any slope angle (from the near-horizontal slopes in this study, to the vertical walls of Ke et al. (Reference Ke, Williamson, Armfield, McBain and Norris2019)) following the approach of Atoufi et al. (Reference Atoufi, Zhu, Lefauve, Taylor, Kerswell, Dalziel, Lawrence and Linden2023). Additionally, work to implement a sorting-based method similar to Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995) for these periodic, tilted domains with surface forcing could yield interesting insights.
The key research question underpinning any investigation into the ice-shelf–ocean boundary is to understand better the rates at which oceanic heat can be mixed towards the ice and lead to melting, a pre-cursor to which includes improving the resolution towards full turbulence-resolving simulations. Further work is required to identify the applicability of existing parameterisations to flows with these paired instability features, but given the difficulties observing such features in situ, this work presents an important step towards understanding the mechanisms and regimes and of shear-driven mixing at this interface.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.11840.
Acknowledgements
We thank C. Subich, A. Grace and M. Stastna for their contributions to developing this case in SPINS. This work used Northumbria University’s Oswald High-Performance Computer.
Funding
CJL was supported by an Early Career Fellowship funded by The Leverhulme Trust.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study is openly available at Northumbria University’s data repository, accessible at https://doi.org/10.25398/rd.northumbria.c.7941254.
Author contributions
All authors contributed to development of the project, review and editing the manuscript. AJ and SGHE were responsible for project conceptualisation, SGHE developed, carried out, analysed and visualised the numerical simulations and wrote the original draft of the manuscript. CJL conducted the LSA and prepared the figures and initial draft for this section.
Appendix A. Grid sensitivity
Grid sensitivity simulations with
$Nz = 256$
,
$Nz = 512$
and
$Nz = 1024$
were carried out between t = 290 and 650 s. Simulations with
$Nz = 256$
rapidly diverge from the base simulation (
$Nz = 512$
) due to numerical instability that arise from insufficient diffusion at low grid resolutions, indicating lower grid resolutions are unsuitable. Simulations with
$Nz = 512$
and
$Nz = 1024$
are compared with quantitative time series of domain-integrated density measures the density variance (
$\langle \rho '\rho '\rangle$
) and density variance dissipation rate (
$\chi = \kappa \langle (\partial \rho '/\partial {x_i} ) (\partial \rho '/\partial {x_i} )\rangle$
) (figures 16
a and 16
b), along with kinetic measures of kinetic energy (KE) and dissipation (
$\epsilon$
) (figures 16
c, and 16
d). The two simulations show near identical-integrated KE (within 6 %) and dissipation rates (within 20 % prior to
$t = 650$
s). This error (
$KE_{\textit{sim}}-KE_{N_z=512})/KE_{N_z=512}$
) is minimal during early periods and only peaks during active small-scale turbulence after
$t\approx 640$
s. To confirm grid sensitivity for the evolution of density (since
$\textrm {Sc} \lt 1$
), the simulations show good agreement in domain-integrated
$\langle \rho '\rho '\rangle$
and
$\langle \chi \rangle$
. Figure 17 qualitatively compares the two resolutions as instability is developing, and show excellent agreement in the dynamics between each simulation.
Small-scale differences between solutions emerge due to sub-grid-scale motions, where the smallest Batchelor scale (
$\lambda _B = (({\nu \kappa ^2})/{\epsilon } )^{1/4}$
, where
$\epsilon$
is the rate of dissipation of TKE) decreases to approximately 20 % of the vertical grid size during the most energetic periods of flow evolution. Despite this, our simulations demonstrate only a small degree of mesh dependence. Given our primary focus is the development of instability (figure 16
c–f), we deem these grid resolutions appropriate.
Grid sensitivity of simulation ISO_Base shown with the domain-integrated density variance (a), density variance dissipation rate (b), KE (c) and dissipation rate (d) for the
$400$
s sensitivity simulation for sensitivity simulations with
$Nz = 512$
and
$Nz = 1024$
. Note that density variance and
$\chi$
are output every
$10$
s for the high resolution case, and every
$2$
s for the base case whilst
$\epsilon$
and
$KE$
are output every simulation time step. Inset to panel c zooms into the region marked by the box.
Figure 16. Long description
Four panels comparing two grid resolutions – Nz = 512 (blue/grey) and Nz = 1024 (orange, marked with × symbols) – over t ≈ 300-680 s. The legend in panel (c) applies to all panels. Panel (a) – Domain-integrated density variance 〈ρ’ρ’〉: y-axis ranges from 0 to ∼3.5 × 10^-4. Both curves remain near zero and closely overlapping from t = 300 s until approximately t = 500 s, indicating a laminar, low-variance state. Around t ≈ 520-530 s both curves rise sharply, peaking near t ≈ 610-620 s at ∼3 × 10^-4, then declining. The two resolutions agree closely throughout, with only minor divergence near the peak, indicating the density variance evolution is well-resolved at both grid spacings. Panel (b) – Density variance dissipation rate log10χ: y-axis spans approximately –17 to –8. Both curves begin at very low values (∼–16 to –17) and rise slowly and step-wise through t = 300-520 s, then increase more rapidly during the instability, reaching a peak near log10χ ≈ –8 at t ≈ 620-640 s before declining. The × markers for Nz = 1024 are sparser (output every 10 s) compared to the more densely sampled Nz = 512 curve. An inset zooms into the boxed region t ≈ 575-690 s, showing that the two resolutions track each other closely through the peak, with small differences in the decay phase, confirming adequate resolution of the dissipation rate. Panel (c) – Kinetic energy KE: y-axis ranges from 0 to ∼0.20. Both curves are nearly indistinguishable throughout the full time range, rising smoothly from near zero at t = 300 s to a broad peak of ∼0.19 near t ≈ 600-610 s, followed by a sharp drop and a secondary smaller peak near t ≈ 640 s, then declining. The excellent agreement between Nz = 512 and Nz = 1024 confirms that the kinetic energy evolution is fully converged with respect to vertical resolution. Panel (d) – Dissipation rate log10ε: y-axis spans approximately –3.2 to –2.5. Both curves rise from ∼–3.1 at t = 300 s, accelerating around t ≈ 500 s and peaking sharply near log10ε ≈ –2.6 at t ≈ 600-610 s. A sharp drop follows, with a secondary peak near t ≈ 640 s. The two resolution curves agree closely up to and through the primary peak, with a small but visible divergence in the secondary peak and decay phase at t ≳ 630 s, where the higher-resolution Nz = 1024 case shows slightly higher dissipation, as expected given its ability to resolve finer scales.
Qualitative comparisons of the instability dynamics for grid sensitivity of simulation ISO_Base showing
$\rho$
(a,b) and
$\zeta$
(c,d) for the base case (left) and double resolution case (right) for
$t = 610$
s.
Figure 17. Long description
Four pseudocolour panels arranged in a 2×2 grid, comparing base resolution (Nz = 512, left column) and double resolution (Nz = 1024, right column) at t = 610 s. All panels share x spanning 0.1-1.0 m and z spanning 0.25-0.30 m. Panels (a) and (b) – Density ρ (kg m^3): Colourbar ranges from ∼990 (light blue) to 1000 (dark purple) kg m^3. Both panels show four to five large KH billows arranged quasi-periodically in x, with dark purple dense fluid forming deep, elongated downwelling lobes and light blue low-density fluid occupying the upper portion of the domain. The billow cores are clearly defined in both cases, with similar spatial positions, scales, and overall morphology. The structures in (b) appear marginally sharper and more clearly delineated than in (a), consistent with the higher resolution capturing finer interface detail, but the large-scale billow geometry is effectively identical between the two resolutions. Panels (c) and (d) – Vorticity ζ (s^1): Colourbar ranges from –10 (green, CCW) to +30 (dark red, CW) s^1. Both panels show a broad band of strong positive (red) vorticity concentrated near the upper boundary (z ≈ 0.29-0.30 m), with the vorticity field modulated in x by the billow structures below. Green (negative) vorticity regions are visible within and between the billow cores at mid-depth. The spatial organisation and magnitude of the vorticity field are closely matched between (c) and (d), with the double-resolution case again showing slightly crisper gradients at the vorticity interfaces. The overall agreement between both resolution pairs confirms that the instability dynamics at this time are well-resolved by the base grid.
Appendix B. Perturbed simulations
In order to assess the ability of the flow to transition to an unstable state at the points identified by the LSA, two versions of ISO_Base simulation were restarted with added perturbations representing the dominant unstable mode identified via LSA, as discussed in
$\S$
§ 2.3 and 3.2.1. Each field was re-initialised to one dimensional mean profiles with 2-D structure arising entirely from the added perturbations
\begin{align} u(x, z) = \langle \left (u(x, z)\right )\rangle _x + u' , \\[-28pt] \nonumber \end{align}
\begin{align} w(x, z) = w' , \\[-28pt] \nonumber \end{align}
\begin{align} \rho (x, z) = \langle \left (\rho (x, z)\right )\rangle _x + \rho '. \\[0pt] \nonumber \end{align}
The amplitude of the perturbations for each field (
$B, \boldsymbol{u}$
) was set to 1 % of the magnitude of that field in the original profiles. One simulation (ISO_Base+LSA_t340) was initialised from profiles at
$t= 340$
s, using the dominant mode from
$t = 340$
s, to represent the scenario predicted from figure 4, where the flow becomes unstable at this time. A further simulation (ISO_Base_LSA_t300) was initialised from simulation profiles at
$t = 300$
s, using the mode from this same time with the wavenumber matching the dominant mode at
$t = 340$
s (since there is no dominant or unstable mode at this earlier time).
In ISO_Base+LSA_340, the flow immediately becomes unstable as a result of the perturbation (figure 4
b), evolving in a near-identical manner to ISO_Base at the later time steps (figure 18
j–l). It is notable that the evolution of the flow is slower, as predicted by the slower
$\sigma$
at earlier times (figure 4). In ISO_Base+LSA_300 there is an initial decay of the perturbation between
$t = $
300 and 340 s due to dissipation. When the flow becomes unstable, the instability is slower to grow owing to a smaller initial state at the point of instability, but nonetheless evolves in the same manner as ISO_Base+LSA_340 and ISO_Base (figure 18
e–h).
The evolution of perturbations in the density field for simulations ISO_Base (left), ISO_Base+LSA_t300 (centre) and ISO_Base+LSA_t340 (right). Top panels indicate the initial state with seeded perturbations visible in (e) and (i), the second row shows the emergence of the unstable mode (as in figure 9d) with the lower two panels showing the emergence and evolution of the overturning structures. Note the time intervals are not necessarily equal between simulations, and the
$x$
-axis is cropped for clarity.
Figure 18. Long description
Twelve pseudocolour panels of density perturbation ρ’ arranged in a 4×3 grid. Columns correspond to three simulations: ISO_Base (left), ISO_Base+LSA_t300 (centre), ISO_Base+LSA_t340 (right). Rows show four evolutionary stages, with times differing between simulations reflecting their different perturbation seeding times. All panels share x spanning 0–0.5 m and z spanning 0.25–0.30 m, with a common diverging colourbar from –1×10^3 (blue) through 0 (white) to +1×10^3 (dark red). A black isopycnal contour is overlaid in each panel. Row 1 – Initial state: Panel (a), Base, t = 340 s: The density perturbation field is essentially featureless – uniform white/near-zero throughout – indicating no organised perturbation structure at this early time. The isopycnal is flat near z ≈ 0.285 m. Panel (e), Base+LSA_t300, t = 300 s: Small but visible perturbations are present, seeded from the LSA mode at t = 300 s. Alternating weak red and blue patches are visible near the interface, with a spatial structure consistent with the dominant unstable wavenumber. The isopycnal shows very slight undulation. Panel (i), Base+LSA_t340, t = 340 s: Similar to (e) but seeded at t = 340 s. Weak alternating perturbation patches are visible near the interface, slightly larger in amplitude than in (e) reflecting the stronger base flow at the later seeding time. Row 2 – Emergence of the unstable mode: Panel (b), Base, t = 570 s: The perturbation field remains weak and relatively unstructured, with only faint, small-amplitude variations near the interface. The isopycnal is nearly flat, consistent with the flow still being in a pre-instability laminar state. Panel (f), Base+LSA_t300, t = 370 s: A clear, organised perturbation pattern has emerged. Alternating positive (red) and negative (blue) lobes are arranged periodically in x near the interface, matching the buoyancy structure of the LSA mode (as in figure showing LSA mode visualisation, panel d). The isopycnal shows coherent wave-like undulations. Panel (j), Base+LSA_t340, t = 350 s: Similar organised modal structure to (f), with alternating red/blue lobes near the interface. The pattern is slightly more developed than (f) at the equivalent stage, consistent with the stronger shear at the later seeding time accelerating growth. Row 3 – Onset of overturning: Panel (c), Base, t = 590 s: The perturbation field now shows strong, large-amplitude structures. Broad positive (red) and negative (blue) regions dominate, with the isopycnal beginning to overturn. The spatial scale is larger than the seeded cases, reflecting natural (noise-seeded) growth to a later nonlinear stage. Panel (g), Base+LSA_t300, t = 440 s: Clear overturning is developing. The perturbation field shows incipient cat’s-eye structures with the characteristic positive/negative lobe arrangement beginning to roll up. The isopycnal is strongly deformed. Panel (k), Base+LSA_t340, t = 420 s: Similar overturning onset to (g), with comparable perturbation amplitude and spatial organisation, occurring slightly earlier in time than the t300 case. Row 4 – Advanced overturning: Panel (d), Base, t = 600 s: Fully developed overturning. Large-amplitude positive and negative perturbation regions are arranged in a clear billow pattern, with the isopycnal rolled into a tight spiral. The structures span the full vertical extent of the panels. Panel (h), Base+LSA_t300, t = 450 s: Advanced billow structure closely resembling (d), with large red and blue lobes and a strongly overturning isopycnal. The morphology is qualitatively identical to the Base case despite occurring ∼150 s earlier. Panel (l), Base+LSA_t340, t = 430 s: Similar advanced overturning to (h), occurring ∼20 s earlier than the t300 case and ∼170 s earlier than the unforced Base simulation, confirming that LSA-mode seeding at either time substantially accelerates the onset of nonlinear instability.
Appendix C. Wave interaction interpretation of instability applied to the compound instability
Using the interpretation of stratified shear instability as interacting vorticity and gravity waves, the fundamental tendency of this flow to go unstable is analysed, and contrasted with the inherently stable jet flow discussed in Smyth & Carpenter (Reference Smyth and Carpenter2019). To simplify this analysis for a fundamental overview of the instability, the flow is reduced to a jet flow, symmetric about the velocity maximum and far from boundaries (figure 19 a). Each of those assumptions is likely to impact the exact dynamics observed.
The wave field diagram (figure 19) shows potential interactions between the four interfaces (three vorticity interfaces, one buoyancy interface) in a piecewise profile representative of these flows. Following Lloyd & Dorrell (Reference Lloyd and Dorrell2024), the coincident density and buoyancy interface at the mid-depth forms a combined vorticity-gravity wave, and hereafter is considered a single interface. The analysis shows two pairs of interactions between waves, firstly the lower vorticity interface forms a resonant interaction with the central vorticity-gravity wave, where the nodes of maximum vertical velocity perturbation are able to coincide with sympathetic peaks and troughs in the wave, acting to amplify them (i.e. the Rayleigh condition that there exists an inflection point somewhere in the flow (Smyth & Carpenter Reference Smyth and Carpenter2019) is met for this interaction). The intrinsic direction of propagation for the lower vorticity wave is left, whilst the intrinsic direction of propagation for the middle vorticity wave is right, the resonant internal gravity wave is assumed to take this same propagation direction (and is the direction used to identify resonance), and so this pair of interfaces matches the Fjørtoft condition that, somewhere in the flow, the second derivative of velocity is of the opposite sign to the velocity (Smyth & Carpenter Reference Smyth and Carpenter2019). In much the same manner, the upper vorticity interface and mid vorticity-gravity interface meet both the Rayleigh condition and Fjørtoft condition, overall indicating the susceptibility in qualitative terms for this flow to become unstable with the vorticity-gravity wave forming interactions with each of the outer vorticity interfaces.
Wavefield diagram for a stratified jet system. Left panel is streamwise velocity profile, centre is vorticity, right is a representation of the waves at each interface, marked on with intrinsic propagation directions, a general sense of the anomaly in flow rotation (vorticity, reduced to
$\text{d}u/\text{d}z$
for a 1-D flow), and vertical velocity perturbation at the nodes.
Figure 19. Long description
Three panels sharing a common vertical z-axis, illustrating the wave interaction mechanism for instability in a stratified jet profile. Left panel – Streamwise velocity profile u: A piecewise-linear profile with u ranging from 0 to 1. The profile increases linearly from zero at the bottom to a maximum at mid-height, then decreases symmetrically back to zero at the top, forming a triangular (tent) shape. This represents a jet-like velocity profile with a single maximum at mid-depth and zero velocity at both boundaries. Centre panel – Vorticity profile ω: A piecewise-constant (step) profile consistent with the piecewise-linear velocity above. Three distinct levels are shown: a negative vorticity value (left step) in the lower portion, zero in the mid-section, and a positive vorticity value (right step) in the upper portion – or equivalently, the profile shows two vorticity jumps corresponding to the two kinks in the velocity profile (at the lower and upper inflection points of the jet). This is the du/dz representation, with the sign of the vorticity reversing between the lower and upper shear layers. Right panel – Wave representations at three interfaces: Three horizontal rows of wave diagrams are shown, each representing the wave activity at one of the three vorticity interfaces (corresponding to the top boundary, the jet maximum, and the bottom boundary from top to bottom). Each row contains: A grey sinusoidal wave shape (the interface displacement), with a wavelength spanning roughly half the x-domain (0–1). Two circular arrows indicating the sense of the vorticity anomaly induced by the wave displacement: both arrows in each row rotate in the same sense, consistent with the local vorticity sign at that interface. Black horizontal arrows at the left edge of each wave row indicating the intrinsic phase propagation direction of the wave along the interface (leftward for the upper and lower interfaces, rightward for the middle interface, reflecting counter-propagating waves). Blue vertical arrows at the wave nodes (zero-crossing points) indicating the sign and relative magnitude of the vertical velocity perturbation w’ at those locations: upward (blue arrow pointing up) at one node and downward (blue arrow pointing down) at the adjacent node, consistent with the phase relationship between interface displacement and vertical velocity in a propagating wave. The diagram illustrates the resonance mechanism for shear instability: the upper and lower interface waves propagate in opposite directions but can phase-lock and mutually amplify through their induced velocity fields, while the middle interface wave interacts with both. The relative phasing of the vertical velocity perturbations (blue arrows) at each interface determines whether the wave interactions are constructive (leading to growth) or destructive (neutral/stable).
The nature of the paired interactions, each being between an internal gravity interface and a vorticity interface, indicates that the flow falls within the broad categorisation of HWI, however, such a categorisation is only relevant for this symmetric form.
ρ
u
g
Δρ
θ
ν
t=
t=570
Rig(z)<0.25
t=
t=
ζ=∇×u
kx
max(σ)
ω
kx
σ>0
t=300
Um
δm
t=50
560
t=560
lν
τν
t
tν
Grδ
Reu=1000
Reu=10000
Riu=0.5
kx∗=1.36
c∗=ω∗/kx∗=U/Um
y
t=570
t=1030
t=1030
t=592
t=370
t=220s
ν
t=835
ν
t=470
Reu
Riu
Riu=1
t=570
t=894
t=850
894
914
934
954
966
994
1014
1054
x−z
t=604
t=610
t=618
z
x−z
t=
t=560
400
Nz=512
Nz=1024
χ
10
2
ϵ
KE
ρ
ζ
t=610
x
du/dz
Open access
