1. Introduction
Stratified shear flows are commonplace in the ocean and the atmosphere. Understanding the mechanisms by which such flows become turbulent and lead to irreversible mixing due to the ultimate break down of different types of primary instabilities is vital in understanding diapycnal fluxes of heat and other important scalars such as salt and carbon (Caulfield Reference Caulfield2021). Various types of primary instabilities are possible depending on the particular structure and orientation of the background velocity and density profile, as well as whether (or not) there is external large-scale forcing (Linden Reference Linden1979; Caulfield Reference Caulfield2021). Here, we will focus on mixing in a class of flows associated with the now classical initial value problem of parallel stratified shear flow where the initial velocity and density profiles are chosen to have a smooth variation between two layers with uniform density, travelling at different horizontal velocities in the far field. In this class of flows, the initial background velocity and density distributions are often assumed to take the (symmetrical) form
\begin{eqnarray} \boldsymbol{U}^*_{0}\left (t^* = 0\right ) &=& \mathit{U^*_{0}} \tanh {\left (\frac {z^*}{d^*}\right )} \boldsymbol{e}_x, \quad \rho ^*\left (t^* = 0\right ) = -\rho ^*_{0} \tanh \left (\frac {z^*}{\delta ^*}\right ), \end{eqnarray}
such that the total velocity and density jumps are
$2U^*_{0}$
and
$2\rho ^*_{0}$
, respectively, superscript asterisks denote dimensional quantities and
$\boldsymbol{e}_x$
is the unit vector in the
$x$
-direction. Note that the frame of reference has been chosen so that the midpoint (
$z^*=0$
) of the ‘shear layer’ (i.e. the region over which the background velocity varies significantly) is stationary, and the density
$\rho ^*$
is the departure from the base hydrostatic state with characteristic density
$\rho ^*_r$
, with zero departure also chosen to be at the midpoint of the shear layer
$z^*=0$
. We will also assume that the (Oberbeck–) Boussinesq approximation (Oberbeck Reference Oberbeck1879; Boussinesq Reference Boussinesq1903) applies, with a linear equation of state relating temperature to density in a viscous incompressible fluid with constant kinematic viscosity
$\nu ^*$
and thermal diffusivity
$\kappa ^*$
with
$\rho ^*_0 \ll \rho ^*_{r}$
. Using
$\rho ^*_{0}$
,
$U^*_{0}$
and
$d^*$
as characteristic density, velocity and length scales, we obtain the non-dimensional governing (incompressible, yet non-constant density) Navier–Stokes equations:
\begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{u} &= 0, \end{align}
\begin{align} \frac {\partial \boldsymbol{u} }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla } \boldsymbol{u} &= -\boldsymbol{\nabla } p - Ri_{b}\rho \boldsymbol{e}_z+ \frac {1}{{\textit{Re}}} {\nabla} ^2 \boldsymbol{u}, \end{align}
\begin{align} \frac {\partial \rho }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla } \rho &= \frac {1}{\textit{RePr}} {\nabla} ^2 \rho , \end{align}
where
$p$
is the departure from the hydrostatic pressure associated with the base hydrostatic density profile and
$\boldsymbol{u}$
is the (total) non-dimensional velocity. The flow evolution thus depends on four non-dimensional parameters: the (bulk) Richardson number
$Ri_b$
; the Reynolds number
${\textit{Re}}$
; the Prandtl number
${\textit{Pr}}$
; and the length scale ratio
$R$
, defined as
\begin{equation} Ri_{b} := \frac {g^* \rho ^*_{0}d^*}{\rho ^*_{r}U_{0}^{*2}}, \quad {\textit{Re}} := \frac {U^*_{0}d^*}{\nu ^*}, \quad Pr := \frac {\nu ^*}{\kappa ^*}, \quad R := \frac {d^*}{\delta ^*} , \end{equation}
where
$g^*$
is the (assumed constant) acceleration due to gravity.
Motivated by oceanographic applications, we restrict attention to fluids with
$Pr \sim O(10)$
and flows with
${\textit{Re}} \gtrsim O(1000)$
. Such flows are well known to be prone to two qualitatively different kinds of primary instability for different values of
$Ri_b$
and
$R$
. If
$R \sim 1$
and
$Ri_b$
is sufficiently small, the flow is prone to the classical Kelvin–Helmholtz instability (KHI) (Helmholtz Reference Helmholtz1868; Kelvin Reference Kelvin1871), with the shear layer ‘rolling up’ into elliptical vortices or ‘billows’ that break down (at sufficiently high
${\textit{Re}}$
) to a ‘zoo’ of secondary instabilities, leading to vigorous ‘overturning’ mixing of the density distribution (Mashayek & Peltier Reference Mashayek and Peltier2012). Here, ‘overturning’ refers to fluid motion that advects dense fluid over light fluid near the centre of the interface where the density gradient is maximum. This concept has significant points of relation to the concept of ‘internal mixing’ as discussed by Turner (Reference Turner1973). In the simplest symmetrical case considered here, the phase speed of the primary KHI is zero, with associated critical layer located at the shear layer’s midpoint, and so the overturning mixing inevitably leads to an at least superficially diffusive-like broadening of the density interface.
For higher values of
$R$
(and hence a ‘sharp’ density interface embedded within a significantly thicker shear layer), the situation is somewhat more complicated (e.g. Caulfield (Reference Caulfield1994)). For sufficiently small
$Ri_b$
, such flows are still susceptible to KHI. However, for higher
$Ri_b$
, a bifurcation occurs (Smyth & Peltier Reference Smyth and Peltier1989, Reference Smyth and Peltier1991; Parker, Caulfield & Kerswell Reference Parker, Caulfield and Kerswell2019), with two instabilities appearing with non-zero phase speed with critical layers located (symmetrically in this simplest case) above and below the mid-plane
$z=0$
. These so-called ‘Holmboe wave instabilities’ (HWI) are associated at finite amplitude with counter-propagating vortices above and below the density interface, which induce counter-propagating cusped waves on the density interface, that ‘scour’ wisps of fluid away from the interface (Holmboe Reference Holmboe1962; Browand & Winant Reference Browand and Winant1973; Koop & Browand Reference Koop and Browand1979; Smyth, Klaassen & Peltier Reference Smyth, Klaassen and Peltier1988; Lawrence, Browand & Redekopp Reference Lawrence, Browand and Redekopp1991; Caulfield Reference Caulfield1994; Baines & Mitsudera Reference Baines and Mitsudera1994; Carpenter et al. Reference Carpenter, Tedford, Heifetz and Lawrence2011; Lefauve et al. Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018). Here, ‘scouring’ refers to fluid motion that advects dense fluid over light fluid at the edge of the interface, once again with some points of significant relation to the concept of ‘external mixing’ as discussed by Turner (Reference Turner1973). At sufficiently high
${\textit{Re}}$
, these instabilities are also known to be prone to secondary instabilities and lead to turbulent break down (Smyth & Peltier Reference Smyth and Peltier1991), although the ensuing mixing can still be appropriately characterised as ‘scouring’, leaving the density interface largely intact (Salehipour et al. Reference Salehipour, Caulfield and Peltier2016, Reference Salehipour, Peltier and Caulfield2018), and thus the mixing is not really possible to be described and parametrised in terms of a turbulent diffusivity.
As there is emerging evidence that density staircases (i.e. relatively deep, relatively well-mixed ‘layers’ separated by relatively thin ‘interfaces’ of enhanced density gradient) are commonplace in stratified turbulent fluids (Phillips Reference Phillips1972; Crapper & Linden Reference Crapper and Linden1974; Linden Reference Linden1979; Ivey & Corcos Reference Ivey and Corcos1982; Balmforth, Smith & Young Reference Balmforth, Smith and Young1998; Holford & Linden Reference Holford and Linden1999; Waite & Bartello Reference Waite and Bartello2004; Taylor & Zhou Reference Taylor and Zhou2017), it is of great interest to investigate in detail the different characteristics of these two broad types of shear-driven mixing events (see Woods et al. (Reference Woods, Caulfield, Landel and Kuesters2010), Caulfield (Reference Caulfield2021) for further discussion of these concepts), and it seems natural to analyse in detail the qualitative and quantitative differences between KHI-like ‘overturning’ mixing and HWI-like ‘scouring’ mixing.
Heretofore, although there have been many studies of the mixing induced by such stratified shear flows (e.g. Smyth & Peltier (Reference Smyth and Peltier1989), Smyth, Carpenter & Lawrence (Reference Smyth, Carpenter and Lawrence2007), Alexakis (Reference Alexakis2009), Salehipour, Caulfield & Peltier (Reference Salehipour, Caulfield and Peltier2016), Maffioli, Brethouwer & Lindborg (Reference Maffioli, Brethouwer and Lindborg2016)), they have typically been inherently and completely Eulerian. This has been the case particularly in numerical studies, where great advances have been made in considering the energetic pathways (and hence the ‘efficiency’) (Peltier & Caulfield Reference Peltier and Caulfield2003) of irreversible mixing. In the Oberbeck–Boussinesq, linear equation of state, initially statically stable situation considered here, irreversible mixing can be directly related to increases in the ‘background potential energy’ (Lorenz Reference Lorenz1955; Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995; Winters & D’Asaro Reference Winters and D’Asaro1996), i.e. the minimum possible potential energy of the system after volume-preserving (and hence adiabatic) rearrangement of fluid parcels (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995; Nakamura Reference Nakamura1996; Caulfield & Peltier Reference Caulfield and Peltier2000). Therefore, it is natural to consider the ‘efficiency’ of mixing (see the reviews of Peltier & Caulfield (Reference Peltier and Caulfield2003) and Caulfield (Reference Caulfield2021)), i.e. the proportion of the total amount of converted background kinetic energy that is converted irreversibly into potential energy (as opposed to being irreversibly lost to the internal energy reservoir via viscous dissipation).
There appears to be a clear signal of the difference between HWI-induced mixing and KHI-induced mixing in terms of this energetic measure (e.g. Salehipour et al. (Reference Salehipour, Caulfield and Peltier2016, Reference Salehipour, Peltier and Caulfield2018)), but the specific physical and mechanistic reasons for this observed difference are still somewhat unclear. A particular area of uncertainty is disentangling the significance of larger-scale advective (and at least notionally reversible) ‘stirring’ from the inherently smaller-scale, diffusive and irreversible ‘mixing’, not least because it is becoming increasingly clear that ‘history matters’ in understanding irreversible mixing (Tseng & Ferziger Reference Tseng and Ferziger2001; Villermaux Reference Villermaux2019; Caulfield Reference Caulfield2020, Reference Caulfield2021), in that it is necessary to understand in detail the spatio-temporal evolution of an individual fluid parcel (Sreenivasan Reference Sreenivasan2019).
Here, we take some first steps towards understanding the turbulent mixing ‘history’ of individual fluid parcels by taking a particular Lagrangian view of the mixing, enabled by having access to numerical simulations. We sample the spatio-temporal changes in fluid density sampled along Lagrangian (virtual) particle paths, advected by the spatio-temporally varying fluid velocity field as determined by solving the governing equations. This approach allows us to quantify irreversible mixing since the fluid density sampled along Lagrangian virtual particle paths only changes in response to molecular diffusion. This particular Lagrangian point of view allows us to understand how mixing occurs across a density interface by allowing us to observe and quantify the transport and associated irreversible density changes of fluid parcels transported across the interface induced by the (turbulent) motions in both KHI- and HWI-driven flows. We are particularly interested in understanding vertical transport and mixing, and so we follow the trajectories of particles located in streamwise-aligned lines located across the entire spanwise extent of the computational domain at different vertical locations. As we combine the data from the lines located at a given vertical location, we are not in a position to follow the inherently three-dimensional mixing mechanisms of (just to take one example) the streamwise-aligned horseshoe vortices discussed by Pham, Sarkar & Winters (Reference Pham, Sarkar and Winters2012). Furthermore, it is important to stress that our approach does not track the variation in density of specifically labelled Lagrangian finite volume parcels of fluid, but rather how the density changes at a particular virtual point particle whose location is advected by the evolving velocity field.
As we discuss in more detail later, we observe an aggregation phenomenon when there is sustained overturning motion inside the interface. We also examine the time evolution and final value of quantities related to these virtual particles’ buoyancy and find that for the same bulk Richardson number, the root mean square (r.m.s.) buoyancy (variation) of a set of Lagrangian particles starting with the same initial buoyancy is actually larger for HWI than KHI. This perhaps somewhat surprising observation suggests, if only in this highly specific Lagrangian sense, that HWI is a stronger ‘mixer’ even though such flows ‘scour’, rather than ‘overturn’ the density interface. Finally, we analyse the phenomenon that approximately
$65\,\%$
of particles starting near the mid-plane of the interface end up with a buoyancy with opposite sign of their initial buoyancy in the highly stratified case which we consider.
In summary, our objectives in this paper are hence two-fold. First, we wish to develop a physical, mechanistic picture of the time-evolving ways in which fluid parcels are both irreversibly mixed and advected (or stirred) in stratified shear flows in the vicinity of an interface. Second, we wish to investigate a quantitative understanding of the differences between ‘overturning’ and ‘scouring’ mixing events, as modelled by the turbulent break down of flows prone to primary KHI and HWI. To fulfil these objectives, the rest of the paper is organised as follows. In § 2, we present the details of the numerical simulations we consider. We then present our key results in §§ 3 and 4. We present a physical description of the time evolution of the simulations in § 3, as well as a qualitative consideration of the aggregation of particles in buoyancy space. In § 4, we focus on discussion of quantitative measures such as the time series and the final change in mean and r.m.s. of particle density, comparing flows associated with turbulent break down of KHI and HWI. We analyse the effect of stratification, and draw distinctions between overturning mixing and scouring mixing. We also examine the proportion of particles crossing the mid-plane of the interface, focusing on particles with starting position close to the mid-plane. A key result is that approximately
$65\,\%$
of particles starting close to the mid-plane cross the mid-plane for HWI at sufficiently high
$Ri_b$
and mixing occurs away from the Holmboe wave crests. The size of regions that either remain at the same side or cross the interface mid-plane is also analysed. Finally, we present our conclusions in § 5.
2. Numerical details
We consider stably stratified parallel shear flows with non-dimensional initial conditions given by
\begin{eqnarray} \boldsymbol{U}_{0}\left (t = 0\right ) &=& \tanh {\left (z\right )} \boldsymbol{e}_x ; \ \rho \left (t = 0\right ) = -\tanh \left (R z\right ), \end{eqnarray}
and solve (1.2)–(1.4) with a random perturbation to the velocity field with non-dimensional amplitude
$0.005$
to trigger instability, so that the initial perturbation kinetic energy is much smaller than the energy associated with the base velocity field. The boundary conditions are periodic in horizontal directions with rigid walls at the top and bottom where no tangential stress and no density flux boundary conditions are applied.
In all flows, we set
${\textit{Re}} = 1200$
. We also set
$Pr = 8$
, which is a characteristic value for the diffusion of heat in water. We set
$R = 1$
in flows which exhibit KHI and
$R = \sqrt {\textit{Pr}} = \sqrt {8}$
in flows which exhibit HWI (Smyth et al. Reference Smyth, Klaassen and Peltier1988). Note that the initial (in time) gradient Richardson number at the midpoint of the shear layer
$z=0$
is related to the bulk Richardson number and the length scale ratio:
\begin{equation} Ri_{0} := Ri_{g}\left (0\right ) = Ri_{b} R. \end{equation}
We investigate flows prone to primary KHI and HWI at three different values of bulk Richardson numbers
$Ri_{b} = ( 0.10, \ 0.15, \ 0.48 )$
. For
$Ri_{b} = 0.10$
and
$Ri_{b} = 0.15$
, we can construct initial base flows that are prone to either primary KHI and HWI simply by varying
$R$
. Conversely, for the highly stratified flow
$Ri_{b} = 0.48$
, the flow is only unstable to HWI, and so this flow is particularly well suited to investigate mixing due to ‘scouring’ alone.
The simulations use the open source code Oceananigans (Ramadhan et al. Reference Ramadhan, Wagner, Hill, Jean-Michel, Churavy, Souza, Edelman, Ferrari and Marshall2020). For simplicity, we set the domain size in the streamwise direction to be one wavelength of the most unstable linear mode, obtained by solving the viscous Taylor–Goldstein equation (Taylor Reference Taylor1931; Goldstein Reference Goldstein1931) numerically. Table 1 provides further information about the parameters used in the simulations. For all simulations, the resolution is sufficient to ensure that the grid spacing is always no more than
$3.3$
times the Batchelor length scale (Batchelor Reference Batchelor1959), defined as
\begin{equation} L^*_{B} = \frac {L^*_{K}}{\sqrt {\textit{Pr}}} = \frac {\left (\dfrac {\nu ^{*3}}{\epsilon ^*}\right )^\frac {1}{4}}{\sqrt {\textit{Pr}}}, \end{equation}
where
$L^*_{K}$
is the Kolmogorov length scale (Kolmogorov Reference Kolmogorov1962) and
$\epsilon ^*$
is the (dimensional) kinetic energy dissipation rate. A detailed discussion of the resolution of the simulations is given in Appendix A.
Details of the three-dimensional DNS performed. The domain size is
$(L_{x}, L_{y}, L_{z})$
,
$N_{x}$
,
$N_{y}$
,
$N_{z}$
are number of grid points in
$x$
,
$y$
,
$z$
direction, respectively, and
$t$
represents the non-dimensional duration of the simulation.
We run three simulations for each set of parameter values, with each simulation using a different random seed for the perturbation velocity field. This increases the sample size of turbulent events and improves the reliability of our particle statistics. In each simulation, to investigate the effect of the turbulent mixing on the particles, we track 13 planes of particles initially situated with different initial buoyancy. Within each plane, there are 10 lines equally spaced in the
$y$
-direction and each line parallel to the
$x$
-axis consists of 100 particles equally spaced in the
$x$
-direction. We have used a fifth-order weighted essentially non-oscillatory (WENO) method and third-order Runge–Kutta time stepper to evolve the flow. The particles are also advected using a third-order Runge–Kutta time stepper, and the velocities and buoyancy of particles are interpolated using a fifth-order accurate Lagrange polynomial. The 13 planes are initially positioned at heights associated with buoyancies:
\begin{equation} b_i=\left (\pm 0.95, \hspace {1pt} \pm 0.6, \hspace {1pt} \pm 0.4, \hspace {1pt} \pm 0.2, \hspace {1pt} \pm 0.1, \hspace {1pt} \pm 0.05, \hspace {1pt}\textbf {} 0\right ) \times Ri_b , \end{equation}
where non dimensional buoyancy is defined as
$b=-Ri_{b}\rho$
and hence
$2 Ri_b$
is the total jump in the buoyancy across the density interface. The associated initial
$z$
-coordinate of these lines of particles are then determined by
\begin{equation} z_{i} = \frac {1}{R} \text{arctanh}\left (\frac {b_i}{Ri_{b}}\right ) \end{equation}
for
$i=1,\ldots,13$
. This selection of buoyancy levels enables us to examine the effect of turbulent mixing on particles both within the density interface and outside the density interface, and provides us information about particles with initial position close to the mid-plane (
$z=0$
). To enhance clarity, we label the planes as plane
$-6, -5,\ldots ,+6$
in order of increasing initial buoyancy.
3. Flow evolution
The evolution of the buoyancy field and sample particles in all five flows we have considered are illustrated in figures 1–5, each showing four chosen snapshots to demonstrate the evolution of turbulence and mixing. The particles are coloured by their initial buoyancy. Figures 1–5 use the same scale bar and colour bar as shown below figure 5 and have a 1 : 1 aspect ratio, showing the difference in the wavelength of the most unstable mode (and hence horizontal extent) for each flow. Videos of the evolution of the five flows are provided as supplementary material.
Visualisation of the H10 flow are projected onto the mid
$xz$
-plane and are coloured by their initial buoyancy. The streamwise extent is the corresponding
$L_{x}$
specified in table 1, and the vertical extent is chosen to be 7 non-dimensional units centred at the mid-plane (
$z=0$
). The scale and colour bar are the same as shown for figure 5.
Visualisation of the K10 flow. The layout is the same as in figure 1.
Visualisation of the H15 flow. The layout is the same as in figure 1.
Visualisation of the K15 flow. The layout is the same as in figure 1.
There is a very rich collection of coherent structures in each flow. In the first stage, primary instabilities occur. For flows prone to primary KHI, there is a large primary billow at the centre of the interface that induces mixing by overturning, while for flows prone to primary HWI, the primary instability is due to nonlinear Holmboe waves on both sides of the interface which induce mixing by scouring. In the next stage, secondary instabilities occur. This includes braid instabilities that overturn the interface between neighbouring billows and secondary instabilities inside the primary billow induced by KHI (see Mashayek & Peltier (Reference Mashayek and Peltier2012) for a detailed discussion of the ‘zoo’ of secondary instabilities of KHI at sufficiently high
${\textit{Re}}$
). For flows susceptible to HWI, secondary shear instabilities develop inside the interface and mixing induced by the breakdown of Holmboe waves occurs on both sides of the interface (Smyth Reference Smyth2006; Salehipour et al. Reference Salehipour, Caulfield and Peltier2016). After the primary and secondary instabilities break down, small-scale turbulence develops across the shear layer in all flows except for the H48 flow, where small-scale turbulence mostly occurs outside the interface and is much weaker compared with that in the other four flows. Finally, the small-scale turbulence eventually decays and the flow becomes laminar again. As expected, increasing
$Ri_b$
reduces the turbulence intensity.
Although the time evolution of such stratified shear instabilities is now relatively well known, the way in which the turbulent breakdown advects virtual particles is not. The spatial distribution of the Lagrangian particles clearly also exhibits different patterns during the flow evolution in each flow. In the first stage when the primary instabilities are developing, it is apparent that Lagrangian particles cluster at certain locations in the streamwise direction. For the flows susceptible to primary KHI, particles inside the density interface are advected into the primary billow. However, for flows susceptible to primary HWI, particles on both sides of the interface cluster at the wisps induced by the nonlinear Holmboe waves. After small-scale turbulence develops in both canonical classes of flow, the particles become more uniformly distributed in the streamwise direction.
From the snapshot at
$t=155$
in figure 2, we see that for the K10 flow, there is an instability in the braid region (near the mid-point in the
$x$
-direction) where dense fluid sits above light fluid before the billow breaks due to the extremely strong density gradient generated away from the primary billow as it draws in a large amount of fluid (Mashayek & Peltier Reference Mashayek and Peltier2012). The subsequent mixing in the braid region does not affect many particles inside the density interface as almost all particles inside the interface have been advected into the large primary billow. The primary billow in the K15 flow is considerably smaller and weaker than that in the K10 flow (compare figures 2 and 4) but it is sustained for a much longer period of time. Unlike for K10, secondary instabilities do not develop in the braid region in K15 before the collapse of the primary billow.
From the second time snapshots at
$t=480$
and
$t=565$
respectively in figures 1 and 3, we see that there are two positions in
$x$
where particles inside the interface aggregate and overturning motion occurs in the H10 flow and the H15 flow. These two coherent overturns are due to the secondary shear instabilities inside the interface as the primary HWI only induces mixing on both sides of the interface by scouring. The coherent overturns inside the interface in the H15 flow are smaller in size and slower than in the H10 flow. At later times, the H10 flow becomes more turbulent with smaller overturns compared with the H15 flow.
For the H48 flow shown in figure 5, there are no secondary shear instabilities inside the interface and mixing is purely induced by the scouring of Holmboe waves on both sides of the interface. The asymmetry above and below the interface induced by the randomness of the perturbation velocity field is apparent from the profile at later times.
Plots showing aggregation of particles in buoyancy space. From top to bottom, the flows are H10, K10, H15, K15, H48, respectively. Panels (a,c,e,g,i) show the heatmap of
$\text{log}_{10}(\text{number of particles})$
in each bin over the duration of the simulation. The particle buoyancy is normalised by the corresponding
$Ri_{b}$
in each case. The limit of the colour bar is to aid comparison. Panels (b,d,f,h,j) show the buoyancy field at the time labelled
$t_{p}$
in the left column when a peak occurs in the heatmap. The flow visualisations have the same layout as in figure 1, except the particles are now chosen to be lines within planes
$\pm 3$
and plane
$0$
.
The time series of the buoyancy sampled at the particle positions also show a distinct signature when there are sustained overturning motions inside the interface. Figure 6 shows the heatmap of the logarithm (base 10) of the number of particles in each bin (rectangular regions in
$t{-}b$
space of the same size) over the duration of the simulation on the left and the corresponding flow visualisation on the right for each flow. These data are constructed from all the three simulations with given initial conditions, combining all the data from the 10 lines of particles located at different spanwise
$y\hbox{-}$
locations across the flow domain. Although the mixing is naturally inherently three-dimensional, we are interested here in understanding the extent of vertical transport induced by those three-dimensional (and often vortical) motions. The visualisations are shown at the times labelled
$t_{p}$
, which corresponds to the time when the largest number of particles have
$b \simeq 0$
. There are very clear signals of an increase in the number of particles with
$b \simeq 0$
for the flows H15 and K15 as can be observed in figure 6(e) and 6(g). The H10 and K10 flows show a much weaker peak and there is no such aggregation phenomenon for the H48 flow.
By examining the corresponding flow visualisation shown in figure 6(b,d,f,h,j), we see that the increase in particles with
$b \simeq 0$
corresponds to the formation and break up of coherent overturns inside the density interface. From figure 6(f), we see that for the H15 flow, particles near the centre of the interface accumulate in two small coherent overturns slightly offset from the midpoint of the interface. These coherent overturns are induced by secondary shear instabilities inside the interface. Similarly to the H15 flow, there are also two small coherent overturns inside the interface in the H10 flow as can be seen from figure 6(b), but these coherent overturns break down much faster due to secondary instabilities, and hence we do not see a significant increase in the number of particles with
$b \simeq 0$
. From figure 6(h), we see that for the K15 flow, almost all particles near the centre of the interface are advected into the primary billow centred on the interface associated with KHI. In this case, the primary billow is sustained for a long time as can be seen from figure 4. Therefore, for the K15 flow, the fluid inside the primary billow mixes thoroughly and we see a much more sustained increase in the number of particles with
$b \simeq 0$
.
Perhaps unsurprisingly, the buoyancy field in the K10 flow is much more mixed than in the other flows. As already observed in figure 2, density overturns in the braid region develop before the primary billow breaks down. Therefore, when the increase in particles with
$b \simeq 0$
occurs as the primary billow breaks down, the flow is already quite well mixed by the secondary instabilities in the vicinity of the braid. Hence, the primary billow breaks down quickly due to turbulence generated from the braid region, and this occurs before the fluid inside the billow becomes fully mixed. The peak in the number of particles with
$b \simeq 0$
is higher for
$Ri_{b} = 0.15$
than for
$Ri_{b} = 0.10$
, suggesting that the number of particles with
$b \simeq 0$
is a non-monotonic function of
$Ri_{b}$
. This is primarily associated with the more rapid break down of the coherent overturns when the flow is more turbulent.
4. Quantitative analysis
4.1. Comparison of mean and root mean square (r.m.s.) of particle buoyancy
For a purely laminar, diffusive flow, all particles initialised on a horizontal plane will share the same buoyancy for all times. In that case, the r.m.s. (variation) of the buoyancy for particles initialised on a horizontal plane will be zero. Therefore, any change in the r.m.s. buoyancy (variation) is a signature of flow-enhanced mixing (as opposed to laminar diffusion). In this subsection, we focus on the change in mean and r.m.s. buoyancy for the horizontal planes of particles. The particle buoyancy is normalised by half the buoyancy jump across the interface (i.e.
$Ri_b$
) to characterise the significance of the change within the density interface. Specifically, we examine the time evolution of r.m.s. buoyancy, and the final change in mean and r.m.s. buoyancy at the end of the simulations after the turbulence has decayed. The final values have been calculated by averaging over the last ten percent of the duration of the simulation to provide more robust statistics. We classify the particles into three classes: planes
$\pm 6$
are classified as being ‘outside’ the interface, planes
$\pm 5$
and planes
$\pm 4$
are classified as being on the ‘edge’ of the interface, and the remaining seven planes are classified as being at the ‘centre’ of the interface. This classification gives us insight into the behaviour of the flow in different vertical, horizontally averaged regions.
Comparison of time series of r.m.s. (variation) of normalised buoyancy of particles for the various flows. (a) Comparison between H10 and K10 (b) Comparison between H15 and K15 (c) Comparison between H10, H15 and H48 (d) Comparison between K10 and K15.
In figures 7(a) and 7(b), we first compare the time evolution of the r.m.s. particle buoyancy for two different values of
$Ri_b$
. Mixing induced by HWI starts at a later time than KHI, but the r.m.s. buoyancy is, in general, larger in the flows induced by primary HWI, except for the class ‘outside the interface’ for the H10 and K10 flows, where the braid instability in the K10 flow results in a larger r.m.s. buoyancy. For the more weakly stratified flows H10 and K10, mixing occurs almost simultaneously at all levels. However, for the H15 and K15 flows, particles with different initial buoyancy experience mixing at different times. In the H15 flow, particles outside the density interface (red lines) experience mixing slightly before particles near the centre of or on the edge of the density interface (plotted with blue and gold lines). Here, mixing is first induced by breaking of the primary Holmboe waves that develop on both sides of the interface in the H15 flow. The associated vortices that have developed above and below the interface break and mix before the overturning mixing induced by secondary shear instabilities inside the interface occurs. In the K15 flow, particles outside the density interface are mixed after the particles inside the density interface. Here, particles inside the interface experience mixing following the break down of the large primary billow centred at the interface, and particles outside the interface are only mixed when the turbulence in the vicinity of the interface spreads outward, as can be seen from the third panel at
$t=410$
of figure 4.
For flows susceptible either to primary KHI or primary HWI, as
$Ri_{b}$
increases, the r.m.s. buoyancy decreases (figure 7
c,d). For the H10, H15 and K15 flows, particles inside the interface experience more mixing than particles outside the interface as expected, but this is not true for the K10 and H48 flows. For the K10 flow, the braid instability induces significant mixing for particles outside the interface. For the H48 flow, the r.m.s. buoyancy is smaller than in the other flows at the end of the simulation, but the r.m.s. buoyancy for particles outside the interface is slightly larger than particles centred on the interface and particles at the edge of the interface are almost unaffected by the mixing. This is because mixing is mainly induced by the breaking of Holmboe waves on both sides of the interface and the mixing inside the interface is relatively weak as there are no coherent overturns inside the interface.
By comparing the length of time for which the r.m.s. buoyancy is increasing in figures 7(a) and 7(b), we see that turbulence in the flows associated with primary HWI lasts longer than in the KHI-associated flows, consistent with the findings of Salehipour et al. (Reference Salehipour, Caulfield and Peltier2016). This late change is nevertheless much weaker than the changes during the main turbulent mixing phase, although it still extends throughout our simulation. Furthermore, we can deduce that the early-time migration of lines towards the centre in the heatmaps in figure 6(a,c,e,g,i) is purely due to laminar diffusion as there is no change in r.m.s. buoyancy at early times, as shown in figure 7. The migration is strongest in the H48 flow and weakest in the K10 flow, as laminar diffusion is naturally determined by the spatial gradient of density, which is controlled by the total buoyancy jump across the interface as well as the density interface thickness.
Figure 8 shows the change in the mean of normalised buoyancy between the start and end of each flow, and the final r.m.s. of normalised buoyancy for each set of particles. Since the simulations are statistically symmetric about the mid-plane, the r.m.s. buoyancy for particles starting above and below the interface mid-plane with the same magnitude of initial buoyancy are averaged to give the averaged r.m.s. buoyancy in figure 8(b). The final r.m.s. buoyancy is larger for the HWI-associated flows than for the KHI-associated flows for all
$Ri_{b}$
at all levels except where there are braid instabilities for the K10 case. Since purely laminar diffusion will not affect the r.m.s. buoyancy, the final r.m.s. buoyancy is a particular measure of turbulent mixing. Since large values of the particle r.m.s. buoyancy imply that particles starting with the same buoyancy experience different mixing, the r.m.s. buoyancy can be viewed as a measure of the heterogeneity of mixing when viewed from a Lagrangian perspective. In this particular Lagrangian point of view considered in this paper, specifically quantifying flow-enhanced mixing through r.m.s. buoyancy, flows prone to primary HWI may be appropriately described as a stronger mixer than KHI for these two particular
$Ri_{b}$
. We stress that this classification is being made in the very specific sense that particles on average become more spread out in buoyancy space in flows prone to primary HWI than in flows prone to primary KHI. This behaviour is possibly unexpected since turbulent mixing is generally thought to be stronger in KHI compared with HWI using traditional Eulerian measures of mixing (e.g. Salehipour et al. (Reference Salehipour, Caulfield and Peltier2016)) and the overturning/scouring classification. It is important to note that, due to the sharper initial density interface for flows prone to primary HWI, enhanced spread in buoyancy space does not correspond to enhanced spread in physical space, further reinforcing the fact that the ‘stronger mixer’ property of HWI (compared with KHI) presented here is defined in a very specific way.
Comparison of the final change in mean and r.m.s. of normalised buoyancy of particles. The (red dashed) line labelled as ‘PM’ represents perfect mixing and hence has gradient
$-1$
. (a) Final change in the mean of normalised buoyancy (b) Final r.m.s of normalised buoyancy.
From figure 8(a), we see that for the K10 flow, the gradient is very close to
$-1$
(as shown in the figure with a red dashed line labelled ‘PM’ for perfect mixing) for the range of initial normalised buoyancy with magnitude no more than
$0.6$
, which suggests that the large primary billow which entrains almost all fluid inside the interface breaks symmetrically as expected. For all cases, the gradient is close to
$-1$
for the range of initial normalised buoyancy with magnitude no more than
$0.05$
. Despite this, there is a significant difference between the final r.m.s. buoyancy among the simulations for particles initialised close to the interface. The r.m.s. buoyancy behaves very differently, suggesting that near the centre of the interface, flow-enhanced mixing has a significantly greater influence on the spread of particle buoyancy rather than on its mean.
However, the difference in the final r.m.s. of normalised buoyancy is much more significant, both in terms of magnitude and structure. The r.m.s. of normalised buoyancy for HWI-associated flows is always larger than the r.m.s. for KHI-associated flows at almost all buoyancy levels for both values of
$Ri_{b}$
except where there is braid instability for the K10 flow. The final r.m.s. of normalised buoyancy for the H48 flow is much smaller than for the other four flows except outside the density interface as expected as the H48 flow induces much weaker mixing. Outside the density interface, the r.m.s. buoyancy for the K15 flow is smaller than for the H48 flow, which confirms our previous observations that flows prone to primary KHI are not as effective at mixing particles outside the density interface compared with flows prone to primary HWI, always remembering that the initial density interface in the flows prone to a primary HWI is much thinner in physical space than in the flows prone to primary KHI.
The high r.m.s. buoyancy of particles outside the density interface for the K10 flow is apparently due to secondary braid instabilities in the braid region as can be observed from the
$t=155$
panel in figure 2. Comparing HWI at three different
$Ri_{b}$
, we see that particles outside the interface experience a similar change in r.m.s. buoyancy as particles at the centre of the interface for the H48 flow, but particles outside the density interface are less mixed than particles inside the interface for the H10 and H15 flows. By observing figures 1 and 3, this would appear to be due to the effects of the overturns induced by the secondary shear instabilities inside the interface as noted in § 3. The overturns give rise to strong mixing inside and across the interface, and hence increase the r.m.s. of normalised buoyancy of particles inside the density interface for the H10 and H15 flows. The r.m.s. buoyancy for the planes
$\pm 5$
with initial buoyancy
$\pm 0.6 \times Ri_{b}$
is larger than the r.m.s. buoyancy of other planes for the H15 flow as it is at the edge of the density interface. Some particles are ejected out of the interface by the wisps induced by the Holmboe waves and some of them remain inside the interface, and hence create a (relatively) large r.m.s. buoyancy variation.
4.2. Comparison of the final distribution of particles in buoyancy space
In this section, we examine the final distribution of particle buoyancy for the three classes of particles identified previously, i.e. those near the ‘centre’ of the interface (shown in blue), on the ‘edge’ of the interface (shown in gold) and ‘outside’ the interface (shown in red). Figure 9 shows histograms of normalised buoyancy of particles for all five flows, combining all three simulations for each flow.
Histograms of final normalised buoyancy of particles located: near the centre of the density interface (plotted in blue); on the edge of the density interface (plotted in gold); and outside the interface (plotted in red) for the five flows. (a) H10 (b) H15 (c) H48 (d) K10 (e) K15.
The particles are the most uniformly spread out in buoyancy space in the H10 simulation, and the peak of the distribution in the H10 case is smaller than the peaks in the other four flows. A considerable amount of particles that start at the edge of the interface (shown in gold) have been ejected out of the interface by the breaking of the Holmboe waves. In the H15 flow, particles that are initially near the centre of the interface become more spread out in buoyancy space, but have distributions that remain centred at
$b=0$
. Particles on the edge of the interface have been spread out across the buoyancy space. There is some mixing for particles outside the interface, but most of the particles remain outside the interface.
In the H48 flow, the particles exhibit distinct peaks close to the initial buoyancy, reflecting weak mixing in this flow. For the particles on the edge of the interface, there is a shift in the buoyancy peaks towards
$b=0$
due to laminar diffusion. The buoyancy changes from
$ (\pm 0.6, \; \pm 0.4 ) \times Ri_b$
to approximately
$ (\pm 0.3, \; \pm 0.2 ) \times Ri_b$
. Combining information from figures 7(c) and 8(b), we conclude that this shift is due to laminar diffusion associated with a large density jump across a relatively thin density interface as there is no increase in r.m.s. buoyancy. However, the spread of particles outside the density interface in the H48 flow is comparable to the H15 flow and they have also shifted towards
$b=0$
due to laminar diffusion.
Turning our attention to the flows prone to primary KHI, in the K10 flow, the buoyancy distribution for particles initially inside the density interface has a distinct peak at
$b=0$
. Particles outside the interface are more uniformly distributed in buoyancy, but a small peak still occurs in the distribution near
$b=0$
. Figure 2 shows that at
$t=200$
, several particles that were initially outside of the interface (white and black dots) are found inside the turbulent patches left behind from the collapsing primary billows. This suggests that the peak in the particle buoyancy distribution near
$b=0$
is due to mixing that occurs within the primary KH billows.
In the K15 flow, the buoyancy distribution for particles near the centre of the interface has a peak near
$b=0$
as in the K10 flow, while the buoyancy distribution for particles at the edge of the interface is more uniform. Most particles outside the density interface retain their initial buoyancy, implying low levels of mixing. This suggests that the mixing induced by the roll-up and break down of the primary KH billow in the K15 flow is more confined vertically compared with the K10 flow, perhaps unsurprising due to the increased stratification.
By examining flows prone to HWI and KHI across our range of
$Ri_{b}$
, we see that at
$Ri_{b} = 0.10$
, the flow is highly turbulent, and the shape of the particle buoyancy distribution is not highly sensitive to the initial position of the particles for the K10 flow. However, for
$Ri_{b} = 0.15$
and
$Ri_{b} = 0.48$
, the initial location of particles plays a more important role in setting the final particle buoyancy distribution. Particles also behave differently depending on the flow structure at their initial position, for example, as wisps ejecting away from the interface outside the interface in the H15 and H48 flows compared with the primary billow centred on the density interface in the K10 and K15 flows. Comparing flows prone to primary HWI and prone to primary KHI at the same
$Ri_{b}$
, we again see that for both
$Ri_{b}=0.10$
and
$Ri_{b}=0.15$
, particles have a larger range of buoyancy values in flows prone to primary HWI compared with flows prone to primary KHI. Again, in this specific sense, particles appear to experience more mixing in flows prone to primary HWI than in flows prone to primary KHI.
4.3. Proportion of particles crossing the density interface mid-plane
In this section, we consider the fraction of particles which experience a change in sign in buoyancy between the start and end of the simulation. Since both the density interface and the velocity interface are centred at
$z=0$
, when a particle ends with a different sign in buoyancy, its final velocity after the turbulent mixing event also changes sign, meaning that the particle will travel in the opposite streamwise direction, a qualitatively significant change. A particle is identified as having crossed the interface mid-plane if its buoyancy averaged over the last ten percent of the duration of the simulation has a different sign compared with its initial buoyancy. Using buoyancy instead of vertical position to determine whether particles cross the mid-plane of the interface helps to eliminate the influence of long-lived Holmboe or internal waves. The number of particles crossing the interface is undefined for particles with
$b=0$
at
$t=0$
, and these particles are thus excluded from the analysis in this section. Statistics for particles starting the same distance above and below the mid-plane have been combined, exploiting the statistical symmetry of the flow.
Figure 10 shows the averaged proportion of particles crossing the density interface mid-plane for the five cases considered here. The maximum proportion of particles crossing the mid-plane start near the mid-plane as expected. With the exception of the H48 flow (which will be discussed separately later), approximately half of the particles that start with a buoyancy magnitude less than
$0.2$
cross the interface. This is consistent with the observation that particles inside the density interface are entrained into relatively coherent overturns, which then act as vigorous and thorough mixers. Compared with the flows with
$Ri_b=0.10$
, fewer particles with initial normalised buoyancy greater than
$0.2$
cross the mid-plane in flows with
$Ri_b=0.15$
and
$Ri_b=0.48$
.
Comparison of the proportion of particles crossing the density interface mid-plane in all five flows.
The H48 flow is anomalous. Only
$0.3\,\%$
of particles with initial normalised buoyancy
$b = 0.2$
cross the interface by the end of the simulation. However, interestingly, more particles with relatively small initial values of buoyancy cross the interface for the H48 flow than for any of the other flows, with approximately
$65\,\%$
of the particles with an initial buoyancy magnitude of 0.05 crossing the interface. This result is remarkable as it implies that particles that start above the interface are more likely to end up below the interface, and vice versa.
To analyse the mechanism responsible for the large fraction of particles crossing the interface in the H48 flow, figure 11 shows the
$x$
and
$b$
values for two lines of particles with an initial buoyancy magnitude of 0.05 above and below the interface mid-plane. Several times are shown and the particles are labelled according to whether their buoyancy changes sign. The class ‘stay above’ refers to particles that have positive initial and final buoyancy, and the class ‘cross from above’ refers to particles with positive initial buoyancy and negative final buoyancy (analogous definitions apply for particles with negative initial buoyancy). In the snapshot at time
$t=565$
, we see that the particles tend to aggregate at two
$x$
locations (approximately at
$x\approx 2.5$
and
$x\approx 7.0$
). This is also observed in the H10 and H15 flows (not shown here), suggesting that this is common for flows prone to primary HWI. This is a result of different flow velocities near and away from the crests of the Holmboe waves on both sides of the interface as discussed by Smyth & Winters (Reference Smyth and Winters2003), Hogg & Ivey (Reference Hogg and Ivey2003), Smyth (Reference Smyth2006). If the flow is sufficiently turbulent (as in the H10 or H15 flows), coherent overturns as discussed in § 3 will be generated at these
$x$
locations.
Buoyancy space evolution of particles with initial buoyancy
$\pm 0.05 Ri_b$
for the H48 flow.
Another surprising phenomenon can be observed from the snapshot at time
$t=605$
in figure 11. If we focus on the particles marked with squares, which indicates that both their initial and final buoyancy is negative (i.e. they are in the ‘stay below’ class), we see that a considerable proportion of square particles have positive buoyancy at this specific intermediate time. This shows that some square particles cross the mid-plane and have positive buoyancy during their evolution, and then they cross the mid-plane again and return to negative buoyancy. The return of the square-marked particles to negative buoyancy and the change of the upward-pointing-triangle-marked particles to positive buoyancy happens at a very similar time. Therefore, some particles that ‘stay below’ the interface cross the interface twice, while particles that ‘cross from below’ only cross the interface once. A similar phenomenon occurs for the circle-marked ‘stay above’ particles that have positive initial buoyancy and positive final buoyancy. Some of these particles develop negative buoyancy at some point during the simulation before the buoyancy becomes positive again.
Chain length of particles crossing the mid-plane and staying on the same side for particles with initial normalised buoyancy of magnitude
$0.05$
for the five different flows.
Snapshot at
$t=496$
for particles with initial buoyancy
$\pm 0.05 Ri_{b}$
in the H48 flow. The horizontal extent is the corresponding
$L_{x}$
and the vertical extent has been chosen to be two non-dimensional units centred at the interface. The colour of the particle indicates its instantaneous rate of change of buoyancy, as given by the vertical colour bar.
From the initial snapshot at
$t=0$
in figure 11, we see that there are coherent regions in which particles either cross the mid-plane (denoted ‘cross’) or stay on the same side (denoted ‘stay’). We refer to a series of particles which share the same interface crossing behaviour as a ‘chain’. We can then quantify the size of the region by chain length by counting the number of contiguous particles (i.e. with consecutive initial
$x\hbox{-}$
locations) with the same crossing behaviour. This gives us a way to examine an aspect of the coherence of the motion of fluid particles near the interface mid-plane.
Table 2 shows the average chain length across the three different simulations for all five flows for particles with initial normalised buoyancy of magnitude
$0.05$
. Particles which start above and below the mid-plane have been included in the average. For the K10 and K15 flows, the chain length for the ‘cross’ and ‘stay’ classes of particles are similar and close to
$2$
, which is the expectation value of chain length for a randomly arranged chain with probability
$0.5$
of crossing for either class. This suggests that the breakup of the primary billow in KHI-associated flows distributes particles that start close to the mid-plane randomly to either side of the interface. For both HWI-associated flows and KHI-associated flows, as
$Ri_{b}$
decreases, the flow becomes more turbulent and the average chain length decreases towards the random limit
$2$
as expected. For the H48 flow with high
$Ri_{b} = 0.48$
, the chain lengths of ‘cross’ and ‘stay’ are markedly larger, so there is clearly a highly coherent structure of mixing for the highly stratified flow. To identify the origin of such long chains of ‘cross’ particles, we turn our attention to the buoyancy change of particles in physical space in an effort to identify precisely where mixing takes place.
Figure 12 shows a snapshot of particles in physical space with particles coloured by their instantaneous rate of change of buoyancy at a time when the opposing Holmboe waves are passing each other. This specific time is chosen so that the contrast of rate of change of buoyancy at different locations with respect to the Holmboe waves on both sides of the interface is most apparent. The symbol type is the same as in figure 11, except the colour of the particle now indicates its instantaneous rate of change of buoyancy. Particles with an
$x$
location near the sharp crest of the Holmboe waves have a small rate of change in buoyancy, while particles further from the Holmboe waves have a larger rate of change of buoyancy. At the time shown in figure 12, some particles at the left end of the domain with similar
$x$
and
$z$
values have different instantaneous rates of change of buoyancy, suggesting variability in the spanwise (
$y$
) direction and an effect of three-dimensional mixing. We deduce that mixing occurs away from the sharp crest of the Holmboe waves. For the H48 flow, there are no vigorous secondary coherent overturns inside the interface that can overturn the interface and destroy the particle aggregation through turbulent mixing, so larger coherent regions survive. In summary, formation and breaking of billows or other spanwise vortices tend to make the particles more randomly scattered, while coherent Holmboe waves will make the particles aggregate in the streamwise direction due to different flow velocity near and away from the crest and, importantly, irreversible mixing occurs away from the sharp crest of the Holmboe waves. This latter effect appears to be associated with the compression of buoyancy gradient there, leading to an enhanced diffusive buoyancy flux.
5. Conclusion
In this study, we have used direct numerical simulations to study mixing in a stratified shear layer from a (very specific) Lagrangian point of view. We have considered initial conditions with co-located midpoints of the velocity and density distributions that are unstable either to primary Kelvin–Helmholtz instability (KHI) or to primary Holmboe wave instability (HWI) with different bulk Richardson numbers. The resulting flows exhibit mixing by ‘overturning’ of the density interface, and ‘scouring’ due to turbulent motions above and below the density interface. We have tracked the spatio-temporally evolving fluid buoyancy at the location of virtual point particles seeded at various vertical locations relative to the density interface, whose Lagrangian trajectory we track through use of the evolving velocity field.
For flows prone to primary KHI, the Lagrangian particles aggregate in buoyancy space. Specifically, mixing within the large primary KH billow centred on the density interface increases the number of particles with buoyancy
$b \simeq 0$
. However, this effect is transient. As the billow breaks down, fluid that was formerly inside the billow mixes with fluid outside the billow, decreasing the number of Lagrangian particles with
$b\simeq 0$
. For the KHI-associated flow with
$Ri_b=0.10$
(K10), the billow is quickly destroyed by turbulent motion generated by the braid instability, and hence the increase in the number of particles with
$b\simeq 0$
is much less significant. For flows prone to primary HWI, the situation is more subtle. For the HWI-associated flow with
$Ri_b=0.15$
(H15), two sustained coherent overturns inside the interface lead to a transient increase in the number of particles with buoyancy
$b \simeq 0$
, similar to the KHI-associated flows. For the HWI-associated flow with
$Ri_b=0.10$
(H10), there are also two coherent overturns induced by secondary shear instabilities inside the interface, but the coherent overturns are quickly destroyed by secondary instabilities as the flow is more turbulent, and there is no significant increase in the number of particles with buoyancy
$b \simeq 0$
. For the HWI-associated flow with
$Ri_b=0.48$
(H48), there are no such coherent overturns formed inside the interface and the mixing is purely scouring, so there is also no aggregation in buoyancy space.
We have then examined the time evolution and final value of the r.m.s. (variation) and mean of normalised buoyancy of particles starting with the same initial buoyancy, using r.m.s. buoyancy as an indicator for turbulent mixing. We have analysed the response to mixing in three different regions of the flow: near the centre of the density interface; at the edge of the density interface; and far above (and below) the density interface. We have found distinct features for all five flows. The onset time of mixing identified as a change in the r.m.s. or mean particle buoyancy depends on the particle’s initial vertical position. Consistent with the previous work of Salehipour et al. (Reference Salehipour, Caulfield and Peltier2016), the HWI-associated flows exhibit long-lived mixing. The fact that the final change in r.m.s. buoyancy at various buoyancy levels is always larger for flows prone to primary HWI than for flows prone to primary KHI suggests that flows prone to primary HWI are stronger ‘mixers’ than flows prone to primary KHI, at least in this specific Lagrangian point of view in buoyancy space, showing that ‘scouring’ (and the associated Eulerian persistence of a relatively sharp density interface) does not imply that Lagrangian transport of fluid parcels is suppressed.
The HWI-associated and KHI-associated flows with
$Ri_b=0.10$
(the smallest value considered here) exhibit very different mixing mechanisms. As noted previously, mixing within the primary KH billows increases the number of particles with
$b\simeq 0$
. This is true for particles that start inside and outside the interface. However, breaking Holmboe waves in the
$Ri_b=0.10$
flow lead to the ejection of fluid from the interface. In the H10 flow, particles are spread out across the buoyancy space, while in the K10 flow, the particles end up with the same distribution in buoyancy space, regardless of their starting position class. For the higher bulk Richardson number flows K15, H15 and H48, the final distribution in buoyancy space depends largely on the starting position.
An especially instructive measure to understand mixing from our specific Lagrangian view has proved to be the proportion of the particles crossing the mid-plane of the interface. Approximately
$50\,\%$
of the particles that start inside the interface in the K10 flow end on the opposite side of the mid-plane (
$z=0$
). In this flow, the large primary billow centred on the interface mixes the fluid inside the interface regardless of its initial position. Remarkably, for the H48 flow, approximately
$65\, \%$
of the particles starting close to the mid-plane were transported across the mid-plane.
Particles that end up on the same side of the density interface as their original vertical location actually tend to cross the mid-plane twice in the H48 flow. In this flow, despite the fact that there are no coherent overturns induced by the secondary shear instabilities inside the interface, particles still aggregate at two
$x$
-coordinates during the evolution due to the variation of flow speed, especially close to and far from the propagating Holmboe wave crests. Particles with similar interface crossing behaviour occur in coherent regions at
$t=0$
. Analysis of the size of coherent regions by a measure of chain length shows that the distribution of coherent regions becomes more random as the flow gets more turbulent.
Mixing occurs away from the sharp crests of Holmboe waves on both sides in the H48 flow. The H48 flow presents a very large coherent region of particles crossing the mid-plane, which is distinctive from the other flows. The coherent regions are associated with the particle aggregation due to different flow speed near and away from the Holmboe waves on both sides of the interface and the fact that mixing happens away from the Holmboe wave crests. In the H48 flow, the coherent regions survive as there are no billows or coherent overturns inside the interface that can break and overturn the interface to break up the particle chains.
This study has multiple interesting aspects for further investigation. The number of Lagrangian particles that can be placed is restricted by current computing capacity. Tracking more Lagrangian particles would provide more information about the nature of the instabilities and their mixing properties. In particular, initialising with more particles close to the mid-plane of the interface would help us to understand the skewness of mixing towards the interface for flows prone to primary HWI and could also conceivably help to explain the relatively high proportion of particles crossing the mid-plane when there is primary HWI in highly stratified flows. Ideally, by employing enough particles, it should be possible to construct a probability density function relating the initial and final particle buoyancy, which would provide a full measure of the cumulative mixing from our particular Lagrangian point of view. Just to take one example, it would be particularly interesting to consider the spatial dependence of mixing, particularly in the spanwise direction associated with secondary vortical motions. Clearly, it would be useful to study the effect of other parameters such as
${\textit{Pr}}$
and
${\textit{Re}}$
, and examine more cases with different
$Ri_{b}$
. This would be an important step for extending our results to the wider range of parameters that occur in environmental flows. Finally, the technique we have used here to trace buoyancy evolution along Lagrangian trajectories could be generalised to consider other flow properties, such as turbulent kinetic energy, available potential energy and their associated dissipation rates, perhaps illuminating the key question as to if and when mixing can be viewed as a diffusive process.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10626.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Justification of simulation resolution
To justify that the resolutions of the simulations are high enough so that we have a DNS for all five flows, we compare the largest grid spacing with the Batchelor scale over the flow evolution. In the H10 flow, the largest grid spacing is
$\Delta x = 2.49 \times 10^{-2}$
. In the other four flows, the largest grid spacing is
$\Delta y = 2.34 \times 10^{-2}$
. The Batchelor scale has been calculated by
\begin{equation} L_{B} = \left (\frac {\nu ^3}{\max _{z} \langle \epsilon \rangle _{xy}}\right ) ^ {\frac {1}{4}} {\textit{Pr}} ^ {-\frac {1}{2}}, \end{equation}
where
$\langle \epsilon \rangle _{xy}$
denotes the horizontally averaged kinetic energy dissipation rate. For the H10 flow, the maximum value of
$({\Delta x}/{L_{B}})$
is
$3.13$
. The maximum values of
$({\Delta y}/{L_{B}})$
are
$3.27, 2.90, 2.44, 2.45$
for the K10, H15, K15 and H48 flows, respectively. Previous studies have suggested that a maximum grid spacing of
$3{-}6$
times of the Batchelor scale (Smyth & Moum Reference Smyth and Moum2000; Salehipour et al. Reference Salehipour, Peltier and Mashayek2015, Reference Salehipour, Caulfield and Peltier2016) is sufficient to resolve mixing in DNS of stratified shear layers.
Finally, we calculate the averaged perturbation potential energy (PPE) at the mid
$xz$
-plane (
$y=0$
), defined as
$\overline {\langle { \hat {b'} \hat {b'}^{*}}/ {2 N^2} \rangle }$
, where
$\hat {b'}$
is the Fourier transform of the perturbation buoyancy in the
$x$
-direction,
$^*$
indicates complex conjugation,
$N$
is the buoyancy frequency,
$\langle \boldsymbol{\cdot }\rangle$
represents averaging in the
$z$
-direction within the turbulent region (chosen to be 5.23 non-dimensional units centred at the mid-plane
$z=0$
) and
$\overline {\boldsymbol{\cdot }}$
represents averaging over the main phase of the turbulent mixing. The spectra of the averaged PPE are plotted in figure 13 as a function of streamwise wavenumber,
$k_x$
. The spectra are shown for the flows H10 and K10 with the lowest spatial resolution in terms of the Batchelor scales. The averaging time window for the H10 and K10 flows are chosen to be
$250 \leqslant t \leqslant 750$
and
$100 \leqslant t \leqslant 300$
, respectively. The spectra exhibit roll-offs at high wavenumbers for both cases, and the buoyancy variance at the grid scale is small. The blips at the tail of the spectrum are associated with the cut-off Nyquist frequency. These spectra further support that the resolution of the simulations are sufficient for our purposes.
Streamwise spectrum of perturbation potential energy for the H10 and K10 flows. (a) Streamwise spectrum of PPE for the H10 flow (b) Streamwise spectrum of PPE for the K10 flow.
Open access
