A new approach to understanding fluid mixing in process-study models of stratified fluids

. While well established energy-based methods of quantifying diapycnal mixing in process-study numerical models are often used to provide information about when mixing occurs, and how much much mixing has occurred, describing how and where this mixing has taken place remains a challenge. Moreover, methods based on sorting the density field struggle with under resolution and uncertainty as to the definition of the reference density when bathymetry is present. Here, an alternative method of understanding mixing is proposed. Paired histograms of user selected variables (which we abbreviate USP) are 5 employed to identify mixing fluid, and are then used to identify regions of fluid in physical space that are undergoing mixing. This paper presents two case studies showcasing this method: shoaling internal solitary waves and a shear instability in cold water influenced by the nolinearity of the equation of state. The USP method identifies differences in the mixing processes associated with different internal solitary wave breaking types, including differences in the horizontal extent and advection of mixed fluid. The method is also used to identify how density, and passive tracers are mixed within the core of the cold-water 10 Kelvin-Helmholtz instability.


Introduction
At the largest scales, the ocean is stably stratified, typically with warmer, less dense water overlaying cold, and more dense water (although in many locations, such as the polar regions the stratification is salinity dominated instead and warm water can underlay cold).Similar density stratifications also exist across other geophysical settings, such as lakes, estuaries, and the atmosphere.Thus the process by which the layers mix is of great interest to a range of disciplines, from biological scientists interested in the distribution of nutrients, plankton, and sediment, to physical oceanographers interested in the global (and local) distribution of heat and buoyancy.Mixing primarily occurs due to flows that induce small scale motions that stir the fluid and stretch density interfaces.Across the increased surface area, molecular diffusion is effectively sped up in the process of diapycnal mixing (Salmon, 1998).Since they have small scales, mixing processes occur below the grid scale of many oceanic and atmospheric models.Process-scale numerical modelling is thus a useful tool to better understand the routes to mixing, with a view to improving the way that climate-scale and regional models parameterise such processes through turbulent or eddy diffusivity.Large scale models may use sophisticated parameterisations of eddy diffusivity, dependent on local velocity shear and buoyancy, and examples include the well-known k − ϵ model for Reynolds-averaged Navier-Stokes (RANS) (Pope, 2000), and the Smagorinsky model for Large Eddy Simulations (LES) (Wyngaard, 2010).
The mechanical process of stirring is a geometric deformation of fluids, and does not in itself imply irreversible mixing.
However, the presence of a diffusion term in the equations of motion (equation 5) allows mixing (where the concentration of a tracer, or density, of a given fluid parcel is modified) to occur across density gradients that are stretched by stirring.Commonly, mixing in numerical simulations is quantified according to the framework set out by Winters et al. (1995).Under this framework, energy is partitioned into various components; Kinetic Energy (KE), Background Potential Energy (BPE), Available Potential Energy (APE) and Internal Energy (IE).APE is the potential energy available to be released to kinetic energy, and is computed by first adiabatically redistributing the density field to find the lowest possible energy state (this state in turn defines the BPE).The conversion between APE and BPE is considered the energy used to conduct irreversible diapycnal mixing, and is therefore an important concept.Crucially, calculating APE involves sorting the density field in order to identify the lowest energy state achievable by adiabatic redistribution (or under certain formulations, a far-field reference density profile (Lamb, 2008)).Other formulations for measuring mixing have been suggested, such as computing the Thorpe length scales (Thorpe, 1977), which have been applied to both laboratory experiments and numerical simulations (e.g.Carr et al., 2017).
The most widely used APE framework is valuable due to its uses in comparisons to the field, and parameterisation in larger scale models.However, the Winters et al. (1995) sorting process does little to inform us of interesting questions around where and how this mixing is taking place, or what happens to the mixed fluid (Moum et al., 2003;Carr et al., 2017).Exploration of the spatial distribution BPE density, and APE density, or local APE can give us some further information about the flow (e.g.Lamb, 2008;Scotti and White, 2014), but are rarely considered in comparison with the domain integrated, or bulk, values.This is in part due to the complexity of calculating these local values.
Graphical methods of understanding stratified turbulent systems have long been used, for example in the study of fossil turbulence (Caldwell, 1983;Gibson, 1999).Penney et al. (2020) considered the diapycnal mixing of passive tracers from a graphical point of view by using simple two-variable probability density histograms showing the statistical relationship between a tracer and density.Penney et al referred to their primary tool as "weighted density-tracer scatter plots", since we propose a methodology that is more general, we will adopt the acronym USP, for user-controlled scatter plot.The user specifies the two variables chosen, and the manner in which their ranges are set in order to focus on dynamical phenomena of interest.
The pair of fields studied in detail in Penney et al. (2020) is a subset of our more general methodology.
To demonstrate the efficacy of our methodology we choose two application areas, one in which the geometry is complex (shoaling internal waves) and one in which the equation of state is nonlinear (in water below the 4 °C temperature of maximum density).Grace et al. (2021) simulated the evolution of gravity currents in water below the 4 °C temperature of maximum density.They demonstrated profound asymmetries between cold gravity currents intruding into warm water, and warm gravity currents intruding into cold water.They labelled this temperature regime the "weak cabbeling" regime, because two water parcels mixed together yield a different density than the average density of the individual parcels.Grace et al found histograms of a single fluid quantity (e.g.temperature) to be a useful analysis tool in characterising the gravity currents.parcels based on the combination of two quantities, for example asking where the density is within one range, and the kinetic energy (KE) exceeds a threshold value simultaneously.An interactive Matlab tool is developed, and applied to two example flows which have previously been studied extensively.Firstly shoaling nonlinear Internal Solitary Waves (ISWs), that result from such waves propagating up-slope are studied, a process investigated in numerous process-studies for its role in transport and mixing of heat, nutrients, and sediment (e.g.Michallet and Ivey, 1999;Aghsaee et al., 2010;Sutherland et al., 2013;Arthur et al., 2017;Hartharn-Evans et al., 2022a).Secondly, the stratified shear flow and its Kelvin Helmholtz instability, a feature explored in many studies relating to mixing due to its fundamental importance (e.g.Winters et al., 1995;Caulfield and Peltier, 2000;Peltier and Caulfield, 2003;Caulfield, 2021), and here applied in three dimensions to a cold water setting, where the nonlinear equation of state alters the dynamics.Such features have been observed in a range of environmental settings, from enhancing mixing by ISWs on the Oregon Shelf (Moum et al., 2003), deepening the Arctic Ocean surface mixed layer (Lincoln et al., 2016), and affecting the highly stratified Connecticut River estuary (Geyer et al., 2010).
The joint probability histograms, and specifically the change in their form over time, allows identification of fluid quantities that are interesting.Additionally, the method allows us to determine the physical region of fluid where these quantitatively unambiguously identified "interesting things" are happening.Instead of the understanding of how much bulk mixing is taking place (as provided by the sorting algorithm), we can understand where, when and perhaps how mixing is taking place.Using this method of selecting regions based on user-selected paired histograms (USP) also has advantages over the Winters et al. (1995) sorting algorithm in that it is unaffected by small perturbations inherent to numerical methods (e.g.due to uncertainties in bathymetry, or general under resolution), and as there is no assumption of mass (or energy) conservation inside the region to be analysed, the spatial domain for analysis can be restricted to just the physical area we're interested in.This paper is organised as follows.In §2.1, the USP methodology is introduced.In §2.2 the numerical model SPINS is introduced, with specific model setups for the fissioning ISW and cold shear flows are introduced in § §2.2.1, 2.2.2 respectively.
The method is then applied to identify the causes of different mixing regimes according to wave breaking types in §3.1, followed by application of the method to understand diapycnal mixing in the cold shear instability in §3.2.The results conclude with the introduction of a passive, rather than active tracer to the cold shear instability §3.3, identifying the differences between active and passive tracers with differing diffusivity.

Paired Histograms
This study employs pseudocolour plots showing the characteristics of fluid parcels in a two-dimensional (2-D) space where the two dimensions, instead of being spatial dimensions, are chosen to be fluid characteristics, similar to the weighted densitytracer scatter plots in Penney et al. (2020).Each USP plot shows the characteristics of fluid parcels in terms of this redefined 2-D space, with two fluid properties as their coordinates, and the colour showing the proportion of the fluid within each discretized coordinate bin.The change of a scatter plot, reminiscent of temperature-salinity plots widely used to trace water masses in oceanography (Helland-Hansen, 1916;Imasato et al., 1993), to the paired histogram by discretising into bins and applying weightings is shown in figure 1, and described in the following.
Following the method of Penney et al. (2020), the algorithm employed is applied as follows.For two variables ϕ and θ (which could for example represent density and kinetic energy), the domains are subdivided into N ϕ and N θ bins, with sizes: For each grid cell, the nearest bin centre is identified: For each grid cell, the I ij (ϕ, θ) is calculated: and the total weight for each given bin with centre (ϕ i , θ j ) is given by: Thus for a given bin the weight W ij is the total volume of grid cells that had (ϕ i − 1 2 δϕ ≤ ϕ < ϕ i + 1 2 δϕ) and (θ j − 1 2 δθ ≤ θ < θ j + 1 2 δθ).The end result of the process is a bivariate weighted histogram.The weighting ensures that the Probability Density Function (PDF), reflects the probability a volume of fluid has given properties in mapped domain cases.A code for this applied to the numerical model SPINS (introduced in §2.2) can be found at: https://github.com/HartharnSam/SPINS_USP.This code generalises the formulation in Penney et al. (2020) to any two fluid variables, and to a non-uniform grid (including mapped grids).
The second new tool introduced herein is the technique of Region of Interest (ROI) plots.Here, the variable is plotted in physical space, but only where conditions identified on the paired histogram are met, for example that 1040 < ρ < 1050 kgm −3 and 0.1 < KE < 0.2 Jkg −1 .In relation to mixing in the shoaling wave simulations, the kinematic quantity investigated is enstrophy density, Ω = 1 2 |ω| 2 (where ω = ∇ × u) which describes the rotational energy linked to dissipation in 2-D flows, and so is a quantity associated with (but not a measure of) mixing.The Kinetic Energy density, KE = 1 2 u 2 is used as a measure of energy of translation in the flow.
In order to unambiguously identify time periods of interest in the various case studies reported on below we have utilised the methodology of Shaw and Stastna (2019).Briefly, this technique computes the Empirical Orthogonal Functions (EOFs) for a particular physical quantity using the built in svds function in Matlab.The EOFs are used to create reconstructions of the original physical field, and the infinity norm of the difference between the reconstruction and original field is plotted as a function of time and the number of modes used in the reconstruction.Time periods for which a jump in the number of EOFs needed to obtain an approximation below a given error tolerance are taken as "periods of interest" and hence selected for detailed examination.This technique is attractive because it builds on a standard technique (EOFs), while addressing its well-known shortcoming (the fact that individual EOFs are not physical quantities).

Numerical model
Simulations were carried out with the pseudospectral code SPINS described in Subich et al. (2013).The code has been thoroughly validated using theoretical result and physical laboratory experiments in a number of different configurations including, shear instabilities, boundary layer instabilities (e.g.Harnanan et al., 2017), interaction with topography (e.g.Deepwell et al., 2017), and shoaling ISWs (e.g.Hartharn-Evans et al., 2022a).It is available for download through its online manual: https://wiki.math.uwaterloo.ca/fluidswiki/index.php?title=SPINS_User_Guide The model solves the stratified Navier-Stokes equations subject to the Boussinesq approximation: where u is the velocity, t is time, P is the pressure, Ξ i is the passive tracer field, ρ is the density and ρ 0 is some reference density of the fluid, and throughout this paper density is reported using the density perturbation ρ ′ = ρ − ρ 0 .The physical parameters are gravity g (set at 9.81ms −2 ), the shear viscosity ν (set at 10 −6 m 2 s −1 , chosen to be consistent with the physical value) and scalar diffusivity κ ρ (or κ Ξi for tracer).The unit vector in the vertical direction is denoted by k.The boundary conditions, domain dimensions and other varying parameters for each simulation are described for each case in §2.2.1 and 2.2.2.

Shoaling Internal Solitary Waves
To illustrate the USP method, simulations from previously published papers on ISW shoaling will be utilised, namely simulations 27_111120 (here the surging case), 26_091120 (here the collapsing case), 24_071020 (here the plunging case) from No slip boundary conditions were applied at the flat upper, and mapped lower boundaries to satisfy model requirements.A mapped Chebyshev grid is employed in the vertical, implying a clustering of points near both the upper and lower boundary that scales with the number of points in the vertical squared, and that vertical resolution improves over the slope.Free-slip boundary conditions were applied at the vertically oriented left and right ends of the computational domain, the grid spacing of which was regularly spaced.Grid resolution was 4096 points in the x and 256 grid points in the z coordinate, giving dx = 1.7 mm in all cases except the fissioning case, where dx = 3.5 mm.Away from the slope, dz varies between 0.124 mm at the upper and lower boundaries and 1.8 mm near mid-depth, simulations were two-dimensional and κ ρ = 10 −7 m 2 s −1 .
The vertical density profile in the main tank was set according to the smoothed two-layer, or hyperbolic tangent profile: set for a system like that studied previously in the literature consisting of a thin, linearly stratified pycnocline (h pyc = 0.015 m) sandwiched between homogeneous layers, here referred to as "thin tanh stratification".The density change, ∆ρ = 20 kgm −3 , the reference density ρ 0 = 1026 kgm −3 and the depth of the pycnocline, z pyc = 0.07 m.Simulations are carried out for a slope of s = 0.2 in all cases except the fission case at s = 0.033, and in all cases the slope height is 0.3 m, the full depth of the water column.The bottom boundary follows the form of Lamb and Nguyen (2009): and d (= 0.03) represents a characteristic distance for the transition from 0 to a constant slope of 1, and L s the length of the slope (1.5 or 9m for the fission case).The function smooths the transition from the flat bed to the slope, which is necessary for the spectral code.The surging (a = 0.009 m), collapsing (a = 0.048m) and plunging (a = 0.076 m) cases represent increasing wave amplitudes (a), and the fissioning (a = 0.063m) case has an intermediate wave amplitude.

Cold shear instability
A Kelvin-Helmholtz instability was produced for a three-dimensional (3-D) simulation in a domain 0.512 m in the x direction, and 0.128 m in the y and z directions by initialising a temperature stratified shear flow perturbed with white noise.Free slip boundary conditions were applied at the upper and lower boundaries, with regular grid spacing in all dimensions, and periodic boundary conditions in the x and y domains.To produce a stratified shear flow, the temperature field, T (z) and horizontal velocity field, u(z) were initialised as follows: A nonlinear equation of state suitable for cold, fresh water was chosen.The equation of state to calculate densities from temperature and salinity is the polynomial fit from Brydon et al. (1999), assuming salinity and excess pressure are 0 everywhere.
The background temperature and velocity gradient regions are thus in a similar location, but the density stratification is asymmetric (figure 2 b).Salinity was set to 0 everywhere, and v, w were set at rest, with small white noise perturbations to initialise the instability.The interface depth, z mix and thickness, h mix , were 0.064m and 0.01m respectively.The minimum initial Richardson number, Ri 0 is defined as: The configuration was set such that Ri 0 ≪ 0.25 (= 0.087).The low value was chosen because Ri = 0.25 is the necessary, but not sufficient, criterion for shear instability.For this temperature range, the density range is very small, so for the cold shear instability cases, the density perturbation, ρ ′ is a useful measure.
Such flows characterise situations where the temperature and momentum are linked, rather than density and momentum.For example, a warm river flow entering a cold lake.
To establish the interplay between passive and active tracers, a passive tracer (i.e.dye in a laboratory setting) was additionally initialised in two further simulations with the same initial conditions as temperature (equation 12), but with each case performed at different tracer diffusivity to the temperature tracer, as described in §3.3.The evolution of shoaling ISWs in this stratification takes one of four forms: surging, collapsing, plunging or fissioning (for more details see e.g.Sutherland et al., 2013;Hartharn-Evans et al., 2022a).Each of these manifests different mixing and dissipation characteristics, as identified via bulk mixing measures (Arthur and Fringer, 2014).Figure 3 shows the evolution of one of these examples, collapsing, which occurs with moderate steepness waves over moderate slopes, will be explored here in more detail (each of the other examples is presented in Supplementary Figures 1-3).A more complete description of these processes can be found in Hartharn-Evans et al. (2022a) §4.3, but here the focus is on the evolution of USP during this process.
First, the reference USP case for an ISW prior to interacting with the slope, as shown in figure 1, is explored.At this point, most of the fluid in a three-layer stratification is at either extreme of density (ρ = 1015 or 1035kgm −3 ) (figure 1 b, d), with only a small amount of fluid in intermediate densities across a thin pycnocline.Meanwhile, the KE is strong in both the upper and lower layers around the wave itself (figure 1 a), and is small at the centre of the pycnocline due to shear here.The overall result is an arcing pattern in USP space.This pattern is asymmetric, due to the asymmetric stratification (in the vertical), which concentrates the KE in the upper layer.In order to investigate mixing due to the shoaling of this wave, we will now explore how this pattern changes during the shoaling process, and identify the fluid regions which these new fluid properties apply to, instead using enstrophy, Ω as a measure of energy in the flow available for mixing.Almost inverse to the KE USP in figure 1, Ω is highest at the highest shear region across the pycnocline in the wave, and at the upper and lower boundaries around the wave (figure 3 i).
As the wave begins to interact with the slope and steepen (figure 3 a, b, e, f), a separation bubble forms on the slope due to boundary layer separation Carr et al. (2008); Boegman and Ivey (2009), meanwhile enstrophy remains constrained to the pycnocline (now extended), and the separation bubble, reflected in much higher enstrophy across mid-densities (figure 3 f, j).
As the wave continues to propagate upslope and concentrate into a shallower region, the separation bubble strengthens and produces what is sometimes termed a "global" instability (figure 3 c, g).USP reveals higher enstrophy primarily associated with higher densities, with remnants of the high enstrophy at the upper and lower boundaries still visible in USP as the high enstrophy at high and low density limits (figure 3  A region of USP space can be identified at early output times, which includes only Ω higher than that associated with a stable ISW, and intermediate densities, and therefore is here taken to represent active mixing regions (white box figure 3 i).
The fluid that has the properties contained within such a region can be identified, and is plotted in figure 4 e-h, for each type of wave shoaling.For a collapsing wave at a late time step, this shows a large continuous patch of actively mixing fluid, extending over around 0.7m (figure 4 f), concentrated at intermediate densities (the pycnocline), although slightly skewed towards higher densities (figure 4 j).
Following this same process for four example waves across the four breaking regimes reveals their mixing properties.Surging is observed for small amplitude waves, and the resulting process is a non-turbulent surge of dense fluid propagating upslope (figure 4a).A single region of high enstrophy forms as the combination of that at the pycnocline and the lower boundary (figure 4 e), and this is advected upslope.Little mixing occurs during this process (figure 4 i), and the small active mixing region is advected upslope in the pulse of dense fluid.Plunging, which takes place at high slope and wave steepness values, results in an anvil-shaped instability plunging forward from the rear steepening face of the wave.As a result, overturning occurs, and this process has been associated with high levels of mixing.Although this process produces high enstrophy across the density As an overview of this pattern, as the wave steepness increases, and therefore the wave breaking regime moves from surging through collapsing to plunging, the region of actively mixing fluid increases.Fissioning takes this a step further, with a similar input wave to the collapsing case, the gentle slope allows the wave to mix over a longer distance.

Cold shear instability
The unstable shear flow undergoes the typical evolution through the formation of Kelvin-Helmholtz (K-H) instability, as follows.Initially, the flow is a shear flow in stably stratified fluid (figure 5 a, e), where the Kinetic Energy (KE) is high in both the upper and lower layers, but at the interface of the layers is zero, at the inflection point of the horizontal velocity (figure 2 b).As the fluid is a two layered shear flow, one would expect the ρ ′ /KE USP to be a symmetrical curve, where the highest kinetic energy is associated with the highest and lowest density anomalies away from the interface (which is also where most of the fluid resides), whilst intermediate densities represent the shear layer where KE is lower.Here, the ρ ′ /KE USP show a "crooked smiley face", with most the fluid concentrated at the extremes of density, but due to the temperature stratification and non-linear Equation of State, the curve is shifted towards higher densities (figure 5 i).Despite kinetic energy being high at this stage in the other directions, there is very little kinetic energy in the y direction (KE v = v 2 2 ), indicating the flow features at this stage are not dependent upon 3-D elements of the flow (figure 5 i, m).
It is worth noting that while the evolution of the shear instability is generic, the nonlinearity of the equation of state means that the symmetry across the centre of the co-located shear and density layers is broken.The visual manifestation of the billows is effectively slightly biased toward one side.
Small fluctuations in the flow result in stationary waves growing on the interface, which due to the shear, grow, and roll up forming braids and billows (figure 5 b, f).As the billows form and stretch, an increasing amount of the flow's KE builds around the billows, and therefore associated with the intermediate densities (figure 5 f, j).Until now, the flow is two-dimensional, with negligible kinetic energy in the transverse direction (figure 5 m, n), and little transverse variability.During the braiding and billow growth, the pycnocline is stretched thin (figure 5 b), enhancing the potential for mixing by stretching density gradients.
The roll up of these billows is such that alternating layers of high and low density form, with the density now statically unstable, leading to turbulence (and in turn spanwise instability) and mixing.These billows pair and coalesce by t = 225 s (figure 5 c).
The kinetic energy associated with intermediate densities continues to build as the billows coalesce, but by this stage, an increased volume of fluid that represents fluid in the ρ ′ = −0.03 to −0.01kgm −3 range is evident (figure 5k), indicating that 265 mixing has occurred in the fluid.When the billows coalesce, KE v is high, and almost all of this is associated with the same density range as that of the newly mixed fluid layer at late stages (figure 5 o).This indicates the importance of these transverse flow features in the mixing produced by billows formed from a K-H instability.The skew in the density of newly-mixed fluid, as a result of the cold temperature nonlinearity in density, is to some extent visible (but easily missed) in the ρ ′ plots (figure 5 c), but is only truly apparent in USP space (figure 5 k, l, o, p).The USP region in which active mixing related to 3-D processes is identified in figure 5 o.This is where KE v is higher than any of the flow in early time steps, and with intermediate density.This region of the flow for t = 225s is isolated in figure 6.
The transverse average KE v (figure 6 a) shows that this active 3-D mixing is occurring right at the core of the K-H billow, whilst there is considerable variability throughout the y direction (figure 6 c, d).Peaks in KE v are particularly localised (figure 6 b), indicating the small scale processes at play responsible for the mixing.Such results are indicative of the importance of 3-D simulations for understanding these mixing processes.By isolating the fluid undergoing these processes, we can identify regions of the flow worthy of future study.

Cold shear instability with passive tracers
Two additional simulations were carried out with a passive tracer, Ξ i , with different prescribed diffusivities.A less diffusive case (case K-H Dye 1) had κ Ξ1 = 10 −8 m 2 s −1 (case K-H Dye 1) and a more diffusive case (case K-H Dye 2) had κ Ξ2 = 10 −6 m 2 s −1 .These represent diffusivity approximately one order of magnitude lower and greater than the active temperature field, κ T = 1.4 × 10 −7 m 2 s −1 .A snapshot taken of the time step at which billows coalesce (t = 225s), reveals the emergence of fine structures in the K-H Dye 1 case within the billow core, and with sharp interfaces (figure 7 a), and a much more diffuse   8d) shows which gradient is diffused faster, here indicating the density gradient is diffusing faster than the tracer, that is that the extremes of salinity, associated with the upper and lower layer fluids initially, now map onto a range of densities.Meanwhile, the skew of this figure to the K-H Dye 2 case (figure 7 f) shows the tracer gradient is diffusing faster than the density, that is that the extremes of density, associated with the upper and lower layer fluids, now map onto a range of densities, and indicating that if the fluids were sorted again by density, each layer would contain varying levels of the tracer (that is to say, the tracer has mixed).
Applying the same region of active 3-D mixing (based on the density of that fluid at the initiation of the simulation), as in figure 6 (see region in figure 7 h), produces figures 8 and 9.The mixing of the less diffusive passive dye indicates a pattern of highly localised, filamentous structures that are regions of active dye mixing (figure 8 b-d).In contrast, the diffusive passive dye is mixing in broader-scale structures (figure 9 c).However, the large-scale regions in which this passive dye is mixing is similar across both simulations (figure 8 a, c 9 a, c), indicating that there is an important interplay between the large scale flow structures (the K-H billow, and overturn of the isopycnals shown as the black isosurface) and filamentous fine-scale structures, across which the diffusive dye is able to diffuse.This result highlights the small scales at which diffusion plays an important

Discussion
By pairing a kinematic measure (KE, KE v , or Ω) with a conserved tracer (density, temperature or passive dye), the role of both turbulence in stirring the fluid, and diffusive processes in irreversibly blurring density gradients are shown for several example stratified flows.This process has also been used to identify which fluid constituents are being mixed more effectively by a given process, whether the low density upper layer fluid is being mixed into the pycnocline, or higher density, lower layer fluid mixed upwards, and the eventual fate of these fluids.In the case of the shoaling ISWs, informing whether mixed fluids are retained locally to the site of mixing, advected upslope, or away from the slope as an intermediate (or nepheloid) layer.Each of these questions is difficult to answer using existing bulk energetics models.The USP method has an additional advantage over existing methods for understanding in its ease of implementation, and its reliability, in cases where the model is under-resolved, or has non-trivial bathymetry.
The results from shoaling ISWs reveal considerably different mixing regimes between breaking types, for example the horizontal extents of mixing between collapsing and plunging, where past studies have identified similar results between breaking types based on the energetics-based mixing measures alone.However, whilst it may be the bulk measure that is most useful in parameterising global-scale models, the locality of mixing, and transport of this fluid, which may carry with it nutrients, heat (i.e.thermal refuge for biota), and sediment (depending on the local geomorphological setting), is crucial to understand the effects of such processes on ecosystems.The features outlined by ROI in some shoaling ISW examples (e.g.figure 3 h) have striking similarity to the intermediate nepheloid layers identified by Bourgault et al. (2014), formed due to the shoaling of ISWs, and the interaction of sediment resuspension and mixing.Similar features have been observed in various locations, but due to their episodic nature, their formation mechanisms can be unclear (Schulz et al., 2021).Understanding how these features form and evolve may help understand how such features may contribute to larger-scale sediment fluxes, and as a result carbon exports (Schulz et al., 2021).To investigate these processes, the USP method could be applied with an active sediment tracer.
For the cold shear instability, these results also reveal new insights into the dynamics of shear flows at low temperatures (below the density maximum).The shear flow presented could represent a cold river flow into a lake, and as such the position of the velocity gradient and temperature gradient are linked, whilst the density gradient position is offset.The result of this in the shear flow is to produce asymmetric K-H features (e.g.figure 5b), and mixing bias towards one density (figure 5o).These asymmetric stratified shear flows have attracted recent attention following field observations (e.g Tu et al., 2020;Olsthoorn et al., 2023).Previous work on K-H billows has highlighted the role of 3-D instabilities beyond the point at which billows coalesce.Here, by isolating the fluid where these 3-D (v) flows are playing an active part in mixing, it shows that not only are 3-D processes important at this stage, but they remain within the coalesced billow region, and small scale (figure 6).Beyond that stage, this v component of the KE remains high around the active mixing region, but also increasingly across the entire flow.Furthermore, investigating the mixing of passive tracers is also not possible using energetics models, and has been presented here, identifying which densities the mixing of passive tracers is associated with, and the role of 3-D instabilities in the passive tracer mixing processes.

Conclusions
The USP method identifies how the different regimes of mixing associated with different breaking types (and therefore with differing wave and slope characteristics) manifest themselves.ISWs breaking with surging behaviour produce a small mixing region, which is advected upslope, indicating transport of the mixed fluids, potentially to regions of high biological activity.
In contrast, collapsing ISWs produce a large patch of actively mixing fluid over a long horizontal area, particularly with lower layer fluid mixed up into the pycnocline.Plunging ISWs, which initially appear to have similarities to collapsing ISWs, are instead dominated by upper layers mixed down into the pycnocline, and are more horizontally constrained.Finally, ISWs undergoing fissioning produce active mixing over a long horizontal region, advected upslope in the form of quasi-turbulent dense pulses (i.e.boluses).
Temperature-stratified shear flow simulations in the cold, nonlinear equation of state regime reveal important differences with a density-stratified shear flow.USP at early time steps immediately reveals the asymmetry of the distribution of density and kinetic energy, an asymmetry which goes on to play an important role in the mixing induced by the shear instability.Once density gradients have been stretched, and transverse flow at medium-high densities forms, mixing occurs in an asymmetric manner, producing a late stage state with the layer of newly mixed fluid skewed towards medium-high densities.Simulations of the shear flow with passive tracers of varying diffusivity highlight the nature of flow structures associated with mixing fluid.At low diffusivities, filamentous structures in the tracer are observed, heavily associated with 3-D elements of the flow, whilst at high diffusivity, tracer gradients are blurred.Both scenarios show interplay between the dynamic drivers in large-scale flow features and diffusive effects at filamentous fine-scales once large gradients are produced by stretching.Such fine-scale structures are revealed by the new ROI plots based on the USP double histograms.Overall, the new USP method provides valuable insights into the evolution and dynamics of sheared flows, and the role of temperature-based stratification within the nonlinear EOS regime, and diffusivity on mixing processes in particular.
The USP method provides a method of interpreting and investigating mixing that is not blind, i.e. it invites the researcher to understand the ongoing processes in a way that is ignored by the relative black box of typical (bulk) mixing methods.The tagging of particles based on variables provides a methodology that lies somewhere between a strictly Eulerian and strictly Lagrangian approach to understanding mixing.In the case of the passive tracer the method has been used to trace the motion of specific parcels of fluid.This could be further refined by introducing tracer at different times in the simulation into localized regions of interest.Future work could also set up higher dimensional histograms (although these become increasingly difficult to visualise), and/or more complex selection criteria based on multiple quantities; effectively a language to classify aspects of the fluid motion.

Figure 1 .
Figure 1.Schematic of the creation of USP diagrams.The combination of two fluid properties (left) are represented in variable-variable space (centre), and then summed into discrete bins (right) Hartharn-Evans et al. (2022a) and the Thin-20L case (here the fissioning case) fromHartharn-Evans et al. (2022b).Full details of the simulations can be found in the respective papers.ISWs were simulated in a rectangular tank with waves initiated using the lock gate technique, and a slope at one end of the tank as shown in the schematic of the model set-up in figure2 (a).The length of the tank L x = 7 m (extended to L x = 14.5 m for the fissioning case), and the depth of the tank L z = 0.3 m.A hyperbolic smoothing function produces the gate region by a numerical step in density 0.3 m from the left end of the tank, from which an ISW is produced and propagates from left to right.

Figure 2 .
Figure 2. (a) Shoaling ISW schematic diagram showing the numerical domain used for the shoaling internal solitary waves simulations in this study.(b) Cold shear instability schematic indicating the initial vertical profiles of the density, temperature, and velocity fields for the cold shear instability simulations (horizontal scales adjusted), the dashed horizontal line indicates the mid-tank, and centre of the interface for T and u.

Figure 3 .
Figure 3.Time sequence of collapsing wave shoaling process shown as ρ ′ (top), Ω (middle), and USP for the same variables (bottom).t = [55, 60, 70, 75]s g, k).This evolves into a bolus of fluid moving upslope (figure3 d, h).Although plots of enstrophy and density indicate mixing is occurring around the pycnocline at the head of the bolus (figure3 d, h), USP is informative of where mixing is occurring.As the bolus propagates upslope, the skew of high enstrophy fluid towards higher densities (lower layer fluid) becomes clearer in the USP diagram (figure 3 l).

Figure 4 .
Figure 4. Late time step density (top) region of interest (ROI, middle) where the fluid meeting the criteria is marked in dark red, and USP (lower) plot for four example waves, representing different breaking regimes.The white box in the lower panels indicate the ROI selected, based on USP at earlier time steps.

Figure 5 .
Figure 5.Time sequence of the cold shear instability, shown as density (a-d) and Kinetic Energy density (e-h) fields, and the USP for density vs KE density (i-l) and for density vs KEv (the KE density for just the v component) (m-p).t =[5, 170, 225, 275]s.

270
Beyond this time, mixing across the braids continues, and leads to a new state (figure 5 d) where a quasi-turbulent intermediate density layer forms, which represents a new steady state, where the Ri is reduced by a smoothed density gradient.As the flow begins to settle back towards the new steady state, kinetic energy at the intermediate densities falls, with the presence of the layer of intermediate density evident at ρ ′ ≈ −0.02kgm −3 .

Figure 6 .
Figure 6.Regions of interest of the cold shear instability case at t = 225s, based on the region indicated in the white box of figure 5 o.Shown as the ROI in the transverse mean KEv field (a), the KEv x, z slice most representative of the mean (b), the proportion of the y domain within the ROI for the (x, z) plane (c), and the boundaries of the ROI (red) with the ρ ′ = −0.0454kgm −3 isopycnal for reference (black) (d).

Figure 7 .
Figure 7. Late time step (t = 225s) plots for K-H Dye 1 case (left), the base case (centre) and K-H Dye 2 case (right).Top shows the tracer field (tracer for K-H Dye Case 1 and 2, temperature for base case), centre shows the density vs tracer USP, and lower row shows tracer vs KEv.

Figure 8 .
Figure 8. ROI plots as in figure 6 but for K-H Dye 1 case

Figure 9 .
Figure 9. ROI plots as in figure 6 but for K-H Dye 2 case, for the ROI indicated by the white box in figure 7 h