Jared Penney

Yves Morel

Peter Haynes

Francis Auclair

Cyril Nguyen

Diapycnal mixing of passive tracers by Kelvin-Helmholtz instabilities

This paper considers the diapycnal transport of passive tracers during a Kelvin-Helmholtz mixing event. Numerical simulations of a traditional Kelvin-Helmholtz (KH) configuration of a stratified shear flow are extended to include layers of passive tracer at different locations relative to the shear layer. The evolution of the tracers during the simulation is followed and is analysed using different theoretical approaches. One is to consider the evolution via the distribution in density-tracer space which clearly reveals how the tracers are redistributed across isopycnals by the mixing driven by the growing and saturating KH billow. The shape of the distribution places constraints on the redistribution of the tracer and, for this problem of symmetrically stratified shear, it is shown that the distribution typically tends to a compact, piecewise-linear form. The redistribution across isopycnals is also considered via an "effective diffusion equation" for the tracer relative to density-based coordinates. The equation is a generalization of an equation previously derived for transport of density in these coordinates and includes an extra "eddy term" that depends on the relative geometry of the density and tracer surfaces. Under certain circumstances and at later stages of the flow, the eddy term can be neglected, and the evolution of the mean tracer profile can be adequately represented using a simple diffusion equation where diffusivity is defined as the effective diffusivity of density, scaled by the molecular diffusion of the tracer.

Introduction

On the large scale, both the atmosphere and ocean have stable density stratifications, and processes by which fluid properties are mixed in the vertical direction are crucial to both the circulation and to the distribution of chemical and biological species. Shear between interacting masses of fluid is an integral component of the transfer of energy from the largest scales of flow down to the smallest dissipation and mixing scales through the formation of instabilities and turbulence. Away from boundaries vertical mixing is likely to occur through intermittent events triggered by dynamical shear instabilities, such as Kelvin-Helmholtz (KH) (e.g., [START_REF] Smyth | Ocean Mixing by Kelvin-Helmholtz Instability[END_REF]. The evolution of the flow in simple configurations for KH instability has been been studied in detail, both in the laboratory [START_REF] Thorpe | Experiments on instability and turbulence in a stratified shear flow[END_REF][START_REF] Caulfield | Secondary instability and three-dimensionalization in a laboratory accelerating shear layer with varying density differences[END_REF][START_REF] Patterson | Timedependent mixing in stratified Kelvin-Helmholtz billows: Experimental observations[END_REF]) and through numerical simulation [START_REF] Klaassen | The onset of turbulence in finite-amplitude KelvinHelmholtz billows[END_REF][START_REF] Klaassen | The role of transverse secondary instabilities in the evolution of free shear layers[END_REF][START_REF] Scinocca | The Mixing of Mass and Momentum by Kelvin-Helmboltz Billows[END_REF][START_REF] Alexakis | Stratified shear flow instabilities at large Richardson numbers[END_REF][START_REF] Carpenter | Identifying unstable modes in stratified shear layers[END_REF]Mashayek & Peltier 2012a,b, 2013). Shear instabilities have also been observed and recorded in nature. For example, in the atmosphere, clear air turbulence (CAT) represents a specific hazard for aircraft, and is thought to be largely driven by shear instabilities [START_REF] Browning | Observations of Clear Air Turbulence by High Power Radar[END_REF][START_REF] Fritts | Convective and dynamical instabilities due to gravity wave motions in the lower and middle atmosphere: Theory and observations[END_REF], while "roll-up patterns" in clouds are known to be caused by KH instabilities [START_REF] Fritts | Convective and dynamical instabilities due to gravity wave motions in the lower and middle atmosphere: Theory and observations[END_REF] and the references therein). In the ocean, studies of the effect of KH instabilities have been motivated by the likely importance of mixing processes on the overall density stratification and hence the large-scale circulation (e.g., [START_REF] Wunsch | Vertical Mixing, Energy, and the General Circulation of the Oceans[END_REF]. Measurements indicative of KH instabilities have been recorded near regions of shear along oceanic thermoclines [START_REF] Woods | Wave-induced shear instability in the summer thermocline[END_REF][START_REF] Marmorino | Observations of Small-Scale Mixing Processes in the Seasonal Thermocline, Part II: Wave Breaking[END_REF], during the downwelling tidal phase near seamounts [START_REF] Van Haren | A deepocean KelvinHelmholtz billow train[END_REF], and the surface mixed layer has been observed to deepen due to shear instabilities along its base [START_REF] Lincoln | Surface mixed layer deepening through wind shear alignment in a seasonally stratified shallow sea[END_REF]. KH instabilities have also been observed along interfaces in estuaries [START_REF] Geyer | Shear Instability in a Highly Stratified Estuary[END_REF][START_REF] Geyer | Mixing by shear instability at high Reynolds number[END_REF], which can be chemically and nutritionally rich due to material in the river outflows.

Density is a dynamically active tracer in the sense it plays an active role in driving the flow. However, processes such as Kelvin-Helmholtz instabilities are important in vertical mixing of passive tracers, which have no direct effect on the flow, as is the case with certain low concentration or neutrally buoyant species [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF][START_REF] Canuto | Vertical diffusivities of active and passive tracers[END_REF]), but which are important for other reasons. These include, in the ocean, nutrients and microscopic biological species (e.g., [START_REF] Vaquer-Sunyer | Thresholds of hypoxia for marine biodiversity[END_REF][START_REF] Brierley | Impacts of Climate Change on Marine Organisms and Ecosystems[END_REF] and, in the atmosphere, chemical species that are radiatively active or affect human health (e.g., [START_REF] Seinfeld | Atmospheric chemistry and physics: From air pollution to climate change[END_REF][START_REF] Yang | Aerosol transport and wet scavenging in deep convective clouds: A case study and model evaluation using a multiple passive tracer analysis approach[END_REF], though in addition to dynamical mixing events, cross-isentropic tracer transport in the atmosphere may also be driven by radiative heating and cooling [START_REF] Sparling | Diabatic cross-isentropic dispersion in the lower stratosphere[END_REF]. In stably-stratified flows, passive tracers are mostly advected adiabatically along isopycnal surfaces. However, there is potentially diapycnal transport of tracers associated with the intermittent mixing events, and understanding the details of these events and their effect on tracers, and representing that effect through parameterization in atmospheric and oceanic models, is thus of major importance. In the present paper, numerical simulations are used to investigate this problem.

In most atmospheric and oceanic models, the details of intermittent diapycnal mixing events are not simulated directly. Such mixing is typically regarded as a sub grid-scale process that must be parameterized, often through a turbulent or eddy diffusivity (see, for example, Seinfeld & Pandis (1998, Chapter 18) for a summary of functional forms of eddy diffusivities in atmospheric models, or Griffies (2004, Chapter 7) for a discussion of parameterized phenomena in ocean models). The magnitude of the required eddy diffusivity has been estimated in various ways, for example, with ocean tracer release experiments to estimate diffusivity based on spreading across isopycnals (e.g., [START_REF] Ledwell | Evidence for slow mixing across the pycnocline from an open-ocean tracer-release experiment[END_REF], and extracting estimates from atmospheric radar data (e.g., [START_REF] Fukao | Seasonal variability of vertical eddy diffusivity in the middle atmosphere: 1. Three-year observations by the middle and upper atmosphere radar[END_REF]. Parameterizations have been drastically improved and are able to reproduce the dynamics of the mixed layer from high frequency to seasonal time scales. However, in most developments and existing parameterizations, turbulent diffusivities of passive and active tracers are assumed for simplicity to be equal (or proportional). Whether or not this assumption is justified is not yet clear. One aspect of this uncertainty is the effect of different molecular diffusivities between passive and active species or indeed between different passive species. This certainly needs to be taken account when considering turbulent mixing in the ocean, where diffusivities of heat and salt, both of which may contribute to density, differ by a factor of about 100. The effect of differing molecular diffusivities on vertical transport in KH mixing events has been considered by [START_REF] Smyth | Differential diffusion in breaking KelvinHelmholtz billows[END_REF] (for practical reasons the ratio of diffusivities between a passive scalar and density was taken to be to 10). But in addition the assumption of equal turbulent diffusivities may be an oversimplification if the passive tracer and the density have different large-scale sources and sinks. The geometry of the two may then be different at small scales and the resulting differences in molecular diffusive transport may potentially lead to differences in turbulent diapycnal transport at macroscales. There has been some investigation of this topic using numerical simulation [START_REF] Nagata | The difference in turbulent diffusion between active and passive scalars in stable thermal stratification[END_REF] and some discussion of the potential importance in the atmospheric boundary layer [START_REF] Li | Monin-Obukhov Similarity Functions for the Structure Parameters of Temperature and Humidity[END_REF]. Further work is needed to evaluate under which circumstances turbulent diffusivities are equal or similar for all tracers and under which they are not and to assess the implications for representation of diapycnal fluxes of different tracers in oceanic and atmospheric models.

The research reported in this paper examines the vertical mixing of passive tracers in KH instability, focusing in particular on whether the extent to which what has been learned already regarding mixing of density can also be applied to passive tracers, and what other factors need to be taken into account. Numerical simulations of KH instability in a standard flow configuration, including a set of passive tracers, are presented (sections 2 and 3), and analyzed using various techniques (sections 4 and 5). One is to use a density-tracer scatter plots to examine the relative distribution of density and tracer and to relate this to the mixing. The shape of the scatter plot places constraints on the redistribution of the tracer and in particular it is shown that for the flow configuration considered the relation between the density and the tracer is often piecewise linear or close to piecewise linear in the end state. Another technique used in previous studies is to exploit a tracer-based coordinate system in which the effect of mixing can be represented completely by an effective diffusivity. The question addressed here is whether the effective diffusivity for density also usefully represents the mixing of other tracers, which is tested by varying the initial distribution of tracer (section 7).

Numerical model and flow configuration

Governing equations and numerical model

The governing equations are non-dimensionalized using scales typical for studies of Kelvin-Helmholtz instabilities, with tildes denoting the dimensional forms of distance x = (x, y, z), time, t, velocity, u = (u, v, w), density, ρ, tracer concentration, φ, and pressure, p, x = x/h, y = ỹ/h, z = (z -z 0 ) /h, t = t U 0 /h, u = ũ/U 0 , ρ = (ρ -ρ 0 ) /∆ρ, φ = φ/∆φ, p = p/ρ 0 U 2 0 .

The dimensional parameters used in this non-dimensionalization are the value defining the width of the pycnocline and the shear layer, h, the midpoint of the stratification and shear layer, z 0 , the maximum magnitude of the initial velocity field, U 0 the density at the midpoint of the density distribution, ρ 0 , the half of the change in density across the two layers, ∆ρ, and the maximum initial tracer concentration, ∆φ. Using these dimensionless variables, the continuity equation and equations of conservation of momentum, density and passive tracer can be written in their incompressible and Boussinesq forms as

U0 1 m/s ρ0 1017 kg/m 3 ∆ρ 1.2 kg/m 3 h 10 m z0 -143 m ν 5 × 10 -3 m 2 /s κρ 5 × 10 -3 m 2 /s κ φ 5 × 10 -3 m 2 /s g 9.81 m/s 2
∇ • u = 0, (2.1a) ∂u ∂t + u • ∇u = -∇p -Ri 0 ρ k + Re -1 0 ∇ 2 u, (2.1b) ∂ρ ∂t + u • ∇ρ = (Re 0 Pr) -1 ∇ 2 ρ, (2.1c) ∂φ ∂t + u • ∇φ = (Re 0 Sc) -1 ∇ 2 φ, (2.1d )
with k denoting the unit vector of the vertical axis. The dimensionless parameters included in these equations are the initial Reynolds number, Re 0 = U 0 h/ν, the Prandtl number, Pr = ν/κ ρ , and the Schmidt number, Sc = ν/κ φ , where ν is the kinematic viscosity, and κ ρ and κ φ are the prescribed molecular diffusivity constants for density and the passive tracer, respectively. Ri 0 is the minimum initial Richardson number, which follows the standard definition for the gradient Richardson number,

Ri = - g ρ 0 ∂ ρ ∂ z ∂ ũ ∂ z 2 , (2.2)
and can be expressed as

Ri 0 ≈ g∆ρ h ρ 0 U 2 0 , (2.3)
based on the initial conditions presented in the following section. The configuration parameters are chosen so that the necessary criterion for stratified shear instability, Ri 0 < 1/4, is satisfied. For the sake of completeness, the dimensional parameters prescribed in the simulations presented in this paper are given in table 1.

The numerical simulations presented in this paper were performed using the nonhydrostatic, non-Boussinesq version of the Coastal and Regional Ocean Community model (CROCO). This model was adapted from the Regional Ocean Modeling System (ROMS, [START_REF] Shchepetkin | The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model[END_REF] to include non-hydrostatic and compressible effects [START_REF] Auclair | A non-hydrostatic non-Boussinesq algorithm for free-surface ocean modelling[END_REF]. While numerical simulations of Kelvin-Helmholtz instabilities are often considered in a periodic domain with free-slip rigid lid conditions for the upper and lower boundaries (e.g., Mashayek & Peltier 2012b, 2013), the implementation presented here utilizes a free-surface upper boundary, and a flat, solid, bottom boundary, with periodic lateral boundary conditions in the x-and y-directions (streamwise and spanwise directions, respectively). The existence of a free-surface and compressibility adds two dynamical processes (surface and acoustic waves) compared to more traditional studies of Kelvin-Helmholtz instabilities in incompressible, unbounded or rigid lid flows. It has been verified that, with the chosen configurations where the instability develops far from the vertical boundaries, the impact of these additional processes on the results is negligible when compared to the traditional configuration (see Appendix A for further details regarding the implementation of CROCO and the discussion in the concluding section).

Initial conditions

The initial conditions (ICs) for the simulations presented in this paper follow those from previous numerical studies of Kelvin-Helmholtz instabilities (e.g., Salehipour et al. 2015).

In their dimensionless forms, the extents of the domain in the streamwise, spanwise, and vertical directions are given by L x , L y , and L z , respectively (their dimensional forms can be obtained by multiplying each parameter by h). Both 2D and 3D configurations are examined, with periodic boundary conditions in the horizontal directions, a flat, freeslip, rigid bottom and a free surface in the vertical direction. The ICs are presented in their dimensionless forms. The initial density distribution is defined as a two-layer stratification with a weak linear background state, ρ (x, 0) = -βz -tanh (z) .

(2.4)

The linear background term βz (where β = h/10L z ) is a minor modification to the standard configuration for the purpose of defining the scatter plots in density-tracer space discussed in section 4.1 over the entire vertical domain in physical space (i.e., each density value initially corresponds to a unique value of height). Note that the definition of the initial Richardson number (2.3) ignores the weak linear background term of the initial density profile, since it was confirmed to have a negligible influence on the development of the instability when compared to the initial stratification without the weak linear background (β ≈ 3.5 10 -3 1 in this study). The initial velocity field is given by,

u (x, t = 0) = U (z) + u (x), v (x, t = 0) = v (x), w (x, t = 0) = w (x).
(2.5) U (z) is the initial background flow providing the shear and is defined by a hyperbolic tangent profile, with the upper layer moving leftward, and the lower layer moving rightward,

U (z) = tanh (z) . (2.6)
u , v , and w are small amplitude perturbations, required to kickstart the instability, defined as

u (x, y, z) = f (z) sin 2πn x L x x 1 + 3D sin 2πn y L y y , v (x, y, z) = 0, w (x, y, z) = -f (z) 2πn x L x cos 2πn x L x x 1 + 3D sin 2πn y L y y . (2.7)
This choice of functional form for the perturbation ensures that the initial velocity field is non-divergent, which is important in the context of this numerical implementation. The function

f (z) = 1 -tanh 2 z α ,
(2.8) ensures that the initial perturbation is localized within the region of the shear, with α = 4.29. = 0.01 is a dimensionless parameter that sets the magnitude of the perturbation in the streamwise direction, while n x and n y set the wavelengths of the initial perturbation in the x-and y-directions, respectively. 3D sets the magnitude of the spanwise perturbations, and for 3D configurations, 3D = 0.2. For 2D configurations, both n y = 0 and 3D = 0, so that the spanwise perturbation remains null.

In addition to the dynamical fields, passive tracer fields are defined by an initial profile of the form

φ (x, 0) = a sech 2 (az -b) , (2.9)
where a is a dimensionless parameter that determines the width of the tracer distribution, while ensuring the total amount of tracer remains constant between simulations, and b indicates the offset of the tracer from the pycnocline. The tracer is thus initially distributed in a layer, with its maximum initial value located at z = b. Both the vertical position and the width of the tracer layer will be varied in simulations presented later in this paper, but for the different dynamical configurations discussed first in section 2.3, the tracer layer is collocated with the midpoints of the stratification and the shear (b = 0), and the maximum initial tracer value is 1 (a = 1). Sketches of the initial vertical profiles of the shear, density, and passive tracer fields are depicted in figure 1, with associated dimensional parameters.

Description of dynamical configurations

The primary experiment presented in this section compares simulations of three different dynamical configurations, all with identical initial background fields and prescribed parameters, as listed in tables 1 and 2. The nondimensional wavelength of the fastest growing mode as predicted by the Taylor-Goldstein equations for these given initial background velocity and density profiles is λ KH ≈ 14.3. As such, the length of the domain is the streamwise direction was chosen as L x = λ KH , so that a single KH billow develops, or L x = 2 λ KH , so that two KH billows develop, possibly leading to pairing in 2D configurations. N x , N y , and N z are the number of grid points in the streamwise, spanwise, and vertical directions, respectively. Combined with the length of the domain in each direction, this leads to a consistent dimensionless grid spacing in the x-, y-, and z-directions of ∆x = ∆y = ∆z = 0.0559. ∆t is the barotropic time-step. Table 3 presents the parameters that vary between each dynamical configuration. Additional parameters relevant to the non-Boussinesq components of the model are listed in Appendix A. The three dynamical configurations are defined as follows: (i) Configuration I is two-dimensional, with L x = λ KH and n x = 1, n y = 0 (so that v = 0 at all times). The initial perturbation corresponds to the most unstable wavelength and only one billow develops in this simulation.

(ii) Configuration II is also two-dimensional, except that L x = 2 λ KH = L z and n x = 2. The length of the domain and wavelength of the perturbation allows for the formation of two billows, which eventually gives way to pairing and additional mixing.

(iii) Configuration III is three-dimensional, with spanwise parameters L y = 0.375 λ KH , n y = 4, while the streamwise and vertical parameters are the same as configuration I (L x = λ KH , n x = 1). This allows for the generation of another type of mixing, with the development of spanwise secondary instabilities that lead to the breakdown of the primary KH billow.

Due to the inclusion of the third dimension and resulting spanwise instabilities, configuration III provides a more realistic simulation than the other two cases. Configurations I and II, while less physically realistic, provide additional pathways to turbulence and mixing, despite showing nearly identical initial behaviour. Comparing the three configurations therefore permits the evaluation of the impact of the route to turbulence.

Diapycnal mixing of a passive tracer by Kelvin-Helmholtz billows: comparison of dynamical configurations

This section presents two-and three-dimensional numerical simulations of stratified turbulence developing from shear-induced Kelvin-Helmholtz billows. Each simulation undergoes the same initial two-dimensional growth as predicted by linear theory, with later stages leading to different forms of turbulence and mixing.

Evolution of the different dynamical configurations

In order to follow the time evolution of the three dynamical configurations, the mean background and perturbation kinetic energies are defined in the usual way (e.g., [START_REF] Mashayek | Shear-induced mixing in geophysical flows: does the route to turbulence matter to its efficiency[END_REF]. This requires definitions for the mean background and perturbation velocity fields. The mean background velocity field is a function of z only, and is given as

u(z) = 1 L x L y u dx dy. (3.1)
The 2D velocity perturbation is a function of x and z, and is defined as

u 2D = 1 L y (u -u (z)) dy. (3.2)
Finally, the 3D velocity perturbation can be expressed using the background velocity and the 2D velocity perturbation,

u 3D = u -u 2D -u (z) . (3.3)
Following these definitions, the background kinetic energy, kinetic energy due to 2D perturbations, and kinetic energy due to the 3D perturbations are defined as

K = 1 L z u 2 2 dz, (3.4) K 2D = 1 L x L z u 2 2D 2
dx dz, (3.5)

K 3D = 1 L x L y L z u 2 3D 2
dx dy dz, (3.6) respectively, noting that the mean kinetic energy can be expressed as the sum of these three values,

K = 1 V u 2 2 dV = K + K 2D + K 3D . (3.7)
The general evolution of the dynamical configurations is tracked in figure 2, which plots their mean kinetic energy as a function of time. Note that although each of the simulations are performed on domains with different volumes, the initial mean kinetic energy is identical for all three cases. This provides a straightforward depiction of the divergent behaviour between the three cases. The 2D and 3D perturbation kinetic energies will be examined later in section 4.2 when discussing tracer mixing.

Figures 3 and4 present the evolution of the density and passive tracer, respectively, for the three dynamical configurations. Each row corresponds to a specific time depicted on the KE time series of figure 2, which indicates an important stage in KH development, or points at which the development of the configurations diverge. The distinct phases enumerated in figure 2 are as follows:

(i) The initial period of linear growth, as predicted by inviscid Taylor-Goldstein theory. This phase is nearly identical for all three simulations, though some weak spanwise effects are visible in the 3D simulation. Each simulation exhibits stirring of the low and high density regions by 2D KH billows with a wavelength close to that of the fastest growing mode predicted by Taylor-Goldstein theory.

(ii) The branching point between 2D and 3D simulations due to the onset of secondary instabilities. All three configurations exhibit secondary shear instabilities (SSI) along the tilted pycnocline, and secondary convective instabilities (SCIs) within the vortices, while spanwise instabilities begin to emerge in the 3D configuration. The density field develops alternating layers of high and low density, while the passive tracer takes the shape of ellipses with maximum values centred within the billows, gradually decreasing to no tracer at the billow exterior. A thin strand of low tracer concentration occurs along the braids, while even lower concentrations of tracer are visible between the horizontal edges of the billows.

(iii) First onset of small scale features, resulting to greater mixing in all simulations. In the 2D configurations, this leads to the density within the billows becoming nearly uniform. In the 3D configuration, this leads to a shallower pycnocline that is no longer unstable to shear instabilities. The 2D cases keep their tracer maxima focused at the centre of the billows, though some tracer is pulled from the high concentration regions towards the outsides of the billows. Meanwhile, the density field experiences greater mixing due to the appearance of the small scale features, and becomes nearly homogeneous within the billows. The secondary instabilities in configuration III destroy the maximum at the centre of the billow, redistributing the tracer from the high concentration region throughout the rest of the pycnocline.

(iv) The onset of pairing in configuration II, which results in much greater mixing further into the low and high density layers, while mixing the concentrated regions of tracer within the two billows with the regions external to the billows without tracer. The single billow of configuration I maintains a nearly constant shape and size.

(v) Second onset of small scales due to enhanced stretching from the pairing in configuration II. Configuration III is in an essentially steady-state.

(vi) The nearly-steady final phase of the 2D configurations. The single billow of configuration I continuously rolls and experiences slow diffusion of density at the edges of the vortex. The pairing of configuration II induces stretching, injecting small scale instabilities, enhancing homogenization of tracer fields within the billows. This leads to formation of a single larger billow that eventually remains mostly steady.

Density-tracer scatter plots

This section presents the first diagnostic technique used to describe and quantify the diapycnal fluxes of a passive tracer. It is synoptic and is based on density-tracer scatter plots. 2.4) and (2.9) where the pycnocline and tracer layer are collocated (solid line). Hypothetical evolution of the scatter plot after a mixing event showing diapycnal fluxes of tracer (dashed line). In the new profile, mixing has led to a flux of tracer toward lower densities.

Weighted density-tracer scatter plots

For convenience, mixing is here defined as an irreversible process during which the density or tracer concentration of a given fluid parcel is modified. Stirring is defined as a process by which a fluid parcel is displaced in space without altering its material properties. Stirring is a reversible process. In this manuscript, a process is denoted adiabatic or diabatic if it is reversible or irreversible, respectively. To characterize the diapycnal mixing of tracers during dynamical event, the scatter plots of tracer concentration against density are employed. In such diagrams, a fluid parcel will remain at the same position in density-tracer (ρφ) space whatever its displacement in physical space, provided it is not affected by diabatic processes such as diffusion. Any modification of the scatter plot is therefore indicative of irreversible mixing. Similar diagrams are often used to characterize mixing in geophysical flows, such as temperature-salinity diagrams used to qualify the mixing of large water masses in the ocean (e.g., [START_REF] Tomczak | A multi-parameter extension of temperature/salinity diagram techniques for the analysis of non-isopycnal mixing[END_REF][START_REF] Teramoto | Deep Ocean Circulation: Physical and Chemical Aspects[END_REF]; tracersalinity estuarine mixing curves used to identify whether an estuary may act as a source or sink of a given tracer (e.g., [START_REF] Loder | The dynamics of conservative mixing in estuaries[END_REF][START_REF] Officer | Dynamics of mixing in estuaries[END_REF]; or tracertracer diagrams used to examine compact relationships between different atmospheric tracers (e.g., [START_REF] Tilmes | Development of tracer relations and chemical ozone loss during the setup phase of the polar vortex[END_REF][START_REF] Plumb | Tracer interrelationships in the stratosphere[END_REF]. Scatter plots are thus convenient for the purpose of the analysis presented here as they permit the straightforward identification of diapycnal fluxes of tracer, as they are simply indicated by the generation of new points in density-tracer space (see figure 5).

Additionally, density-tracer scatter plots do not correspond to unique density and tracer profiles or flow distributions in physical space. If every individual point is plotted in a tracer-tracer scatter plot, it is possible for information to be lost when specific points have the same value or overlap due to graphical limitations. This can provide a misleading representation of the relative amounts of tracer in the domain. A value similar to the formulation of the density-tracer probability function (PDF) presented in Appendix D of [START_REF] Plumb | Tracer interrelationships in the stratosphere[END_REF], here called the "weight", is defined to quantify the amount of fluid in a given density-tracer bin, and plotted at discrete points in the scatter plot. The procedure to calculate the weight is outlined as follows. The density and tracer domains are subdivided into N ρ and N φ individual bins, respectively, with sizes

δρ = ρ Max -ρ Min N ρ , (4.1) δφ = φ Max -φ Min N φ , (4.2)
where the subscripts Max and Min indicate the maximum and minimum values of the density and the tracer. The centre of a given bin is defined by the point (ρ i , φ j ), where

ρ i = ρ Min + 2i -1 2 δρ, i = 1, 2, . . . , N ρ , (4.3) φ j = φ Min + 2j -1 2 δφ, j = 1, 2, . . . , N φ . (4.4)
The weight corresponding to a given bin with centre (ρ i , φ j ), W ij (t), is calculated as

W ij (t) = W (ρ i , φ j , t) = 1 V I ij (ρ, φ, t)∆V, (4.5)
where ∆V = ∆x × ∆y × ∆z, and

I ij (ρ, φ, t) = 1, (ρ (x, t) -ρ i , φ (x, t) -φ j ) ∈ -1 2 δρ, 1 2 δρ × -1 2 δφ, 1 2 δφ 0, otherwise, (4.6)
and its value is represented in colour on the scatter plot. Note that the total weight is conserved (i.e., the total is always 1), and the integral of the weight multiplied by the tracer concentration or density is also conserved in the absence of sources and sinks. Since graphical limitations can make it difficult to discern when a scatter plot has converged to a compact relationship, it is useful to define a diagnostic that acts as a measure of the scatter,

R (t) = (φ(x, t) -φ * (z, t)) 2 dV φ(x, t) -φ 2 dV . (4.7)
This value will be referred to as the scatter variance. It relates the total variance of the tracer from the isopycnal mean φ * (as defined by equation (C 11) in Appendix C) to the total variance of tracer from the mean over the whole domain φ. As such, it presents a time evolution of the scatter plots, with larger values of R indicating greater relative variance over a given density bin, and smaller values indicating a convergence towards functional, compact relationships between ρ and φ.

Scatter plot evolution for the different dynamical configurations

Figure 6 provides the evolution of the weighted density-tracer scatter plots for each of the dynamical configurations, with each of the rows corresponding to those in figures 3 and 4. The weight of each density-tracer bin is represented in colour, with white indicating no fluid occupies that region of density-tracer space, dark blue indicating a small amount, and yellow to red indicating the regions most occupied by the fluid. Prior to t = 25, all scatter plots remain close to their initial shape, reflecting that the evolution is mostly adiabatic and advective. Starting around t = 25 (figures 6(a,g,m)), the scatter plots begin to spread just below the top of the initial curve, corresponding to the initial roll-up of the pycnocline and slow diffusion of ρ and φ near the pycnocline. Most of the fluid remains located at the edges of the scatter plots, which is indicative of the near-constant high and low density layers below and above the pycnocline, while the rest of the curve indicates the pycnocline itself, which occupies relatively small region of the domain. The greater spreading of the scatter plots at t = 56 relates to the irreversible mixing of both density and passive tracer which starts when the roll-up of the billow has significantly stretched the interface between the regions of high and low density. While the scatter plots of the 2D configurations appear identical, the slight deviation in shape of the 3D configuration scatter plot is due to the introduction of the spanwise instabilities. The development of these secondary instabilities is the point of divergence between the evolution of the 2D and 3D configurations. The spanwise instabilities lead to the rapid breakdown of the 3D billows, and thus rapid homogenization of the density and passive tracer. This leads to a tent-like scatter plot with a rounded top around ρ = 0, and nearly linear branches on either side, as shown in figure 6(q). The convergence of the scatter plots towards more compact curves indicates that the range of tracer concentrations for a given density value has decreased. The 2D billows continue to rotate uninhibited, which mixes and homogenizes the density within the billow, while trapping local maxima of the passive tracer within the core of the billows. These tracer maxima localized in regions of relatively constant density are visible in the scatter plot as narrow vertical protrusions centred at ρ = 0, as in figures 6(c) and (i). As the density field homogenizes, these vertical protrusions converge to more compact vertical lines, as in figures 6(d) and (j), indicating a wide range of tracer values in a region of nearly constant density. However, while configuration I equilibrates, configuration II experiences new instabilities and undergoes a new turbulent phase associated with billow pairing. This generates large meanders protruding deep into the upper and lower layers of the fluid, which involves mixing over a wide density range. A new single billow is formed that eventually stabilizes. The scatter plot evolves to a new shape where the central branch has vanished.

The scatter variance is plotted with different components of the perturbation kinetic energy, in order to compare the scatter to the dynamical evolution of the instabilities. Note that until t ≈ 50, the evolution of the 2D KE is essentially identical for all three configurations, while it remains the same for the 2D configurations until t ≈ 200. Each configuration shows an increase in scatter variance that occurs just after the initial increase in 2D KE, with both quantities depicting similar rates of increase. The first peak in the plots of scatter variance correspond to the period during which secondary instabilities have started to develop. For configuration I, the scatter variance undergoes a rapid decrease as the scatter plot begins to converge to its three branch shape, with a slight increase around t = 150. After this, the decrease in scatter variance is quite slow, as is the decrease in 2D KE, as the flow has reached a stable state. Configuration II observes similar behaviour in the evolution of its scatter variance, but see a small rapid increase after the onset of pairing (depicted by the rapid increase in 2D KE around t = 300), corresponding to the spreading of the scatter plot into away from its three branch structure to the filled triangle structure, as depicted in figures 6(j) and (k), respectively. The scatter variance then undergoes a relatively rapid decrease to near zero as the scatter plot reaches its stable compact shape. The scatter variance of configuration III experiences the same rapid increase as the other cases after the initial billow roll-up, but rapidly decreases after the organized three-dimensional secondary instabilities develop (the 2D KE sees a rapid decrease during the development of the spanwise instabilities, which appears to relate to a brief breakdown of the billow). The rate of decrease of the scatter variance slows from around t = 90 to 130, following a rapid increase in the 2D KE (related to a brief reformation of the coherent billow which occurs prior to the complete breakdown of the primary instability and further development of turbulence). It quickly reaches zero around t = 150 as the scatter plot begins to reach its compact tent-like shape.

Convergence principle for scatter plots

Because density and the passive tracer are both governed by advection-diffusion equations without sources or sinks, there exists an important constraint on the evolution of the density-tracer scatter plot (see [START_REF] Plumb | Tracer interrelationships in the stratosphere[END_REF] and references therein. [START_REF] Lauritzen | Evaluating advection/transport schemes using interrelated tracers, scatter plots and numerical mixing diagnostics[END_REF] used this constraint to determine if the mixing in numerical models is physical). In a general sense, scatter plot evolution can be understood as follows. As previously mentioned, stirring will not modify the scatter plot, but bring fluid parcels from different regions closer together (e.g., the green points in figure 8). Provided the molecular diffusivities of the density and tracer are equal, mixing will homogenize the density-tracer characteristics of parcels within a cell whose size is determined by the diffusivity coefficient. The resulting density-tracer characteristics of this cell are then the averaged values of the initial parcels, weighted by their volumetric ratio (e.g., the red point in figure 8). This implies that the scatter plot after stirring and mixing will contained within the convex envelope of the initial distribution (region within the red dashed curve in 8). Assuming the mixing does behave realistically, certain properties can be inferred from this constraint:

• Whatever the route to turbulence and mixing, the tracer concentration in a given density range remains within the interval determined by the initial convex envelope. This limits the possible diapycnal fluxes of the tracer. • During mixing, the scatter plot evolves continuously. Thus, at any given time, a scatter plot must lie within the convex envelope of every preceding scatter plot. In addition, extreme values of tracer or density are eroded by mixing. As a result, the convex envelope reduces with time and may eventually converge to a more compact scatter plot.

• The convex envelope of a straight line is simply the same straight line. If fluid parcels along a given straight line mix, they will remain confined to that line.

Note that these properties are only valid if both the tracer and the density diffuse at the same rate, or equivalently that the homogenization cell is the same for density and tracer. The convex envelope constraints can be broken if the density and the tracer have different diffusivities, though the importance of this effect warrants further investigation.

The scatter plot evolution depicted in figure 6 follows each of the constraints listed above. The large scale dynamical mixing due to the KH billows (both the initial billows and pairing) reduces the overall size of the convex envelopes and thus the scatter plots, and the maximum possible value of the tracer concentration in specific density ranges. After a certain amount of time, the maximum tracer value slowly decreases, showing that mixing at this point is no longer dynamically active between the layers, but acts at smaller, localized scales between adjacent points in density-tracer space. This localized mixing appears responsible for the formation of the nearly-linear regions of the converged scatter plots on either side of ρ = 0 because the previous arguments can be applied to restricted portions of the scatter plots provided that mixing acts locally in ρφ space or at reduced scales.

Background density and tracer profile evolution

This section presents the second approach used to describe and quantify the evolution of the tracer field and contrast it with the evolution of the density. It is more quantitative and based on the evolution of mean density and tracer profiles obtained by an adiabatic rearrangement of the density field.

Background density profiles and effective diffusivity

The traditional approach to representing transport and mixing of tracers in turbulent flows is via a diffusive formulation, i.e., a turbulent diffusivity is sought that represents the effect of the turbulent flow on the tracer. One of the fundamental limitations of this approach is that the diffusive representation of random walks, and by extension, of the effect on tracers of quasi-random flows, describes evolution that is the result of many small random steps. This condition cannot be justified for flows that are spatially inhomogeneous, such as the KH instability considered here, where the typical fluid particle displacement by a single eddy is comparable to the length scale on which the properties of the flow change.

One way of overcoming this is to a move to a tracer-based coordinate system. This is the approach in the effective-diffusivity formalism devised by [START_REF] Nakamura | Two-Dimensional Mixing, Edge Formation, and Permeability Diagnosed in an Area Coordinate[END_REF] and [START_REF] Winters | Diascalar flux and the rate of fluid mixing[END_REF]. This formalism applies to systems where tracers are advected and diffused. Contours (in 2D) or surfaces (in 3D) of tracer concentration are used to define coordinate surfaces or contours, but the latter are labelled not by the value of the tracer concentration but the area (in 2D) or volume (in 3D) enclosed by the surface. If the tracer has some kind of geometric organization, then the coordinate system and the variation of tracer concentration within that system represent that organization. For the KH instability and for other flows in density-stratified fluids, the natural tracer is the density and the corresponding tracer-based coordinate, z * , represents a vertical coordinate. By construction, the density is a function of one space variable, z * , and time t alone, written as ρ * (z * , t), which may be shown to satisfy the diffusion equation

∂ρ * ∂t = ∂ ∂z * K ρ ∂ρ * ∂z * , (5.1)
where K ρ is the effective diapycnal diffusivity of density, defined as

K ρ (z * , t) = κ ρ |∇ρ| 2 z * ∂ρ * ∂z * -2 , (5.2)
where κ ρ is the molecular diffusivity of density, and • z * denotes an appropriately defined average over a z * surface (as defined by equation (C 11) in Appendix C). Note that the rearranged (or background) density field ρ * (z * , t) may be calculated from the threedimensional simulation by rearrangement to construct a monotonic profile, i.e. fluid elements, each with a specified infinitesimal volume, are ordered by their density, giving density as a function of cumulative volume (i.e., the volume of fluid with density less than a given value), and then the volume is converted to a vertical coordinate z * by dividing by the horizontal area of the fluid domain. The coordinate z * is then a decreasing function of ρ * . The value of z * for a given ρ * is therefore proportional to the integral, down to ρ * of the weight defined previously in section 4.1. † The left-hand column of figure 9 shows the time evolution of the density profiles sorted from the simulations of the different dynamical configurations (solid lines), and the profiles calculated from equation (5.1), using the effective diffusivity calculated from those three simulations. The right-hand column plots the final profiles based on the rearrangement of the simulation results and the diffusion equation. The area around the pycnocline is magnified to show the area of interest in better detail. As suggested by the evolution of the contour plots in figure 3, the profiles of all three simulations undergo similar evolution until the roll-up of the primary KH billows, at which point the pycnocline of configuration III undergoes greater spreading than the 2D configurations due to the mixing from secondary instabilities. The profiles for the 2D configurations continue to evolve identically until the onset of pairing, at which point the pycnocline of the configuration II profile widens significantly. After mixing, for both 2D configurations, the edges of the pycnocline widen slightly until the end of the simulation. The final state of configuration I is an approximately three-layer profile, with a centre layer near ρ * = 0, and rapid changes to the upper and lower layers. Configuration II has a near constant ρ * = 0 layer of similar width to configuration I, with less steep changes in profile toward the upper and lower layers. The final profile of configuration III differs by not having a constant density middle layer, instead showing the density continuously decrease with height.

Equation (5.1) describes the transport which occurs solely through molecular diffusion of ρ, and hence ρ * , across z * surfaces (there is no advective component). As discussed by [START_REF] Nakamura | Two-Dimensional Mixing, Edge Formation, and Permeability Diagnosed in an Area Coordinate[END_REF] and [START_REF] Winters | Diascalar flux and the rate of fluid mixing[END_REF] (and in subsequent papers that exploit this formalism), the effective diffusivity K ρ is determined by the geometry of the full three-dimensional density field. The dimensionless factor |∇ρ| 2 z * (∂ρ * /∂z * ) -2 has minimum value 1 when the ρ surfaces are planes, and increases as the ρ surfaces become more complex. Thus whilst there is no advective transport required in equation (5.1), the indirect effect of advection is to deform the surfaces of ρ and hence to increase the effective diffusivity.

The approach followed here is to investigate whether equation ( 5.1) provides a quantitatively useful expression of the effect of transport and mixing on density and whether it can be extended to other tracers. Whilst equation (5.1) follows exactly from the partial differential equation for advection-diffusion, with a specified value of molecular diffusivity κ ρ without any approximation, the KH instability simulations rely on a numerical implementation of this partial differential equation, and it is first important to establish whether or not this implementation (5.1) remains quantitatively accurate. It is straightforward to solve the one-dimensional diffusion equation (5.1) numerically by providing the history of K ρ (z * , t) calculated using the density field from the numerical simulations, and to compare the predicted evolution of ρ * (z * , t) with that predicted by the numerical simulation itself. In Appendix B it is shown that applying this approach using the value of κ ρ specified for the numerical simulation under-predicts the mixing of density in the evolution of the KH instability. The explanation is that the numerical approximations to each of the derivative terms of the governing equations, which have been designed with certain properties (for example, controlling spurious oscillations near discontinuities [START_REF] Shu | High Order ENO and WENO Schemes for Computational Fluid Dynamics[END_REF])), in effect augment the specified molecular diffusivity. An estimate of the extra diffusivity provided by the numerical schemes is provided by considering the ratio between the time rate of change of the density variance and the variance of the density gradient (see equation (B 3)). It may be concluded that by this measure, in the simulation depicted (configuration III), the "numerical" diffusivity is up to 100% of the specified molecular diffusivity (see figure 22). When the numerical diffusivity deduced on this basis is added to κ ρ to give a time-varying "net diffusivity" κ Net ρ (t) the agreement between ρ * (z * , t) as evaluated directly from the numerical simulation, and the ρ * (z * , t) obtained by solving (5.1) with κ ρ replaced with κ Net ρ (t) is significantly improved (as depicted in figure 23). Note that this agreement is obtained by adjusting a single quantity, κ ρ , but that the improved agreement is seen in the predicted ρ * as a function of z * . Figure 9 indicates that both the computed and rearranged profiles are in very good agreement for the duration of the simulations, especially towards the edges of the pycnoclines. The greatest disparity appears for the ρ * = 0 isopycnal (green lines) in the configuration I (figure 9(a)), but remains modest.

For reference, the corresponding time evolution of K ρ (divided by κ ρ ) for configuration III is presented in figure 10. As has been demonstrated in previous papers, and in applications to different flows (e.g., [START_REF] Nakamura | Two-Dimensional Mixing, Edge Formation, and Permeability Diagnosed in an Area Coordinate[END_REF][START_REF] Winters | Diascalar flux and the rate of fluid mixing[END_REF][START_REF] Shuckburgh | Diagnosing transport and mixing using a tracerbased coordinate system[END_REF], the distortion of ρ-surfaces leads to a substantial enhancement of K ρ relative to κ ρ , and in the case shown, a factor of several hundred at certain stages of the flow evolution. Note also that the difference between κ Net ρ (t) and κ ρ is not very significant in this enhancement but, again, that it is significant in giving precise quantitative agreement between the evolution of ρ * (z * , t) predicted by equation (5.1) and the evolution according to the full numerical simulation.

Effective tracer diffusivity

Once ρ * is determined, a tracer in the flow can be sorted relative to this density profile by calculating tracer means over isopycnal surfaces,

φ * (z * , t) = φ (x, t) z * ,
(5.3) as defined by equation (C 11). Following the process described in Appendix C, a complete governing equation for this background (or rearranged) tracer profile can be derived. A preliminary step during the derivation provides a formulation that holds for any φ and ρ,

( φ z * ) t = φ t z * + ∂ ∂z * ρ t z * φ z * -φρ t z * ∂ ρ z * ∂z * -1
, (5.4) provided the time derivatives φ t and ρ t are supplied via specific evolution equations. The final step in this derivation that is relevant for the purposes presented in this article supplies these time derivatives via the advection-diffusion equations (2.1c) and (2.1d),

∂φ * ∂t = ∂ ∂z * κ φ κ ρ K ρ ∂φ * ∂z * + ∂ ∂z * ∂ρ * ∂z * -1 κ φ ∇φ • ∇ρ -κ ρ φ ∇ 2 ρ z * , (5.5)
where φ and ρ have diffusivities κ φ and κ ρ , respectively. Here, φ = φ -φ * represents the perturbation from the mean value of tracer for a given density. † An intermediate formulation could be derived by supplying ρ t via an advection-diffusion equation, but allowing for a more general form of the equation for φ. Focusing on equation (5.5), the first term on the right-hand side is straightforwardit simply captures the molecular diffusion of φ * taking account of the geometry of the ρ-surfaces. The second term on the right-hand side takes account of the variation of φ over ρ-surfaces, represented by the quantity φ . Estimating this "eddy term" therefore presents a closure problem, since information about φ is not available from φ * . There are two distinct contributions to this second term, both requiring non-zero diffusivity. The first, proportional to κ φ ∇φ is associated with a diffusive flux of φ across ρ-surfaces. The second, proportional to -κ ρ φ ∇ 2 ρ, is associated with the motion of ρ-surfaces relative to the fluid. Rather than considering the closure problem further, the importance of the tracer eddy term will be investigated by considering a "virtual tracer" φ v (z * , t) which satisfies (5.5) with the tracer eddy term neglected,

∂φ v ∂t = ∂ ∂z * κ φ κ ρ K ρ ∂φ v ∂z * , (5.6)
where initially φ v (z * , t = 0) = φ * (z * , t = 0). Figure 11 shows the time evolution of φ * (solid lines) as calculated from the different dynamical configurations using equa- tion (5.3), and φ v as calculated from equation (5.6) (dashed lines), for the tracer layer collocated with the pycnocline. During the evolution, there exists strong differences between both profiles, starting immediately after the first mixing event. However, both configurations II and III show an increase in agreement between the background profiles and the virtual profiles, with quite good agreement at the end of the simulation (Figure 11(e-f)), despite the virtual profiles slightly underestimating the tracer concentration toward the centre of the domain, and slightly overestimating the amount of tracer at the edges of the profiles. In contrast, configuration I shows marked disagreement between the background and virtual profiles, with φ v being wider than φ * , and significantly underestimating the amount of tracer at the centre of the domain. This does indicate that, provided there are adequate mixing events, φ v and φ * can eventually match, suggesting that the mean tracer profile resulting from mixing can be predicted using a simple diffusive equation with diffusivity given by the density effective diffusivity. This is not a trivial result as it requires that, during the evolution, the cumulative effect of the eddy term vanishes. The fact that the eddy term vanishes at the end of the simulation is not enough.

Relative importance of the eddy term

In this section, the eddy term is examined in greater detail to assess its importance in tracer profile evolution. First note that practically calculating φ can be complicated. For discrete profiles from numerical simulations or observational data, a mean discrete value of φ * can been ascribed to a corresponding ρ * value, with density bins evenly divided. This mean φ * is assumed to be the mean value of each density cell of the physical domain that is associated with the density bin in the background profile. In reality, φ * may still fluctuate over different cells in a given bin. This can pollute the computation of φ , and introduce error into the overall calculation of the eddy term. To avoid this problem, equation (5.5) can rewritten in a way that omits the calculation of φ or its gradient, facilitating the numerical computation of the eddy term,

∂φ * ∂t = ∂ ∂z * κ φ K ρ ∂φ * ∂z * + ∂ ∂z * ∂ρ * ∂z -1 κ φ ∇φ • ∇ρ z * -κ φ ∂φ * ∂ρ * ∇ρ • ∇ρ z * -κ ρ φ∇ 2 ρ z * + κ ρ φ * ∇ 2 ρ z * ,
(5.7)

For the simulations presented here, it has been verified that calculating the evolution of φ * using equation (5.7) yields good results. The profiles obtained from the equation closely match the background profiles sorted from the 2D or 3D fields, even for configuration I, when the profile does not converge to the virtual, purely diffusive, profile. Figures 13(a) and (b) presents the time evolution of the profile of the diffusive term and eddy term, respectively, for the three-dimensional configuration III. The diffusive term first appears non-zero around t = 30, while a non-zero contribution by the eddy term is not visible until about t = 50. Following the initial roll-up of the billow and the during the development of secondary instabilities (from about t = 50 to 130), both terms are of the nearly the same magnitude and of opposite sign, indicating that the eddy term is not negligible. It redistributes tracer across density surfaces but the diffusive term acts so as to compensate this effect. After t = 130, the contribution of the eddy term is again reduced to zero, while there is a negative contribution from the diffusive term near z * = 0 until about t = 150.

Preliminary interpretations of results

For configurations II and III, the reasoning behind the relatively good agreement between the stable state profiles φ v and φ * is not entirely clear. This section presents possible behaviours and features that may be indicative of a tendency toward good agreement between φ v and φ * , as well as situations that may yield discrepancies between the two.

In section 4.2, it is observed that the scatter plots for the fully mixed configurations II and III (figure 6(l,q)) tend to compact forms where there are two outer regions in which the density-tracer relationship appears linear (as suggested in section 4.3), and a central region where the relation lies on a curve. The highest values of tracer concentration are located in the curved central region, with tracer values decreasing along the piecewise linear regions towards the edges. Consideration of the weights in density-tracer space shows that each of these regions correspond to a non-negligible volume of the fluid. Note that there is a strong relationship between the eddy term of equation (5.7) and the scatter plots, in that perfect compactness ensures that φ = 0, and thus the eddy term is null. The behaviour depicted in figure 13 can thus be interpreted as a tendency for the eddy term to generate small scale perturbations on the mean profile φ * and for diffusive term to compensate this effect, as long as the eddy term remains active. After the strong mixing phase, the eddy term diminishes, reflecting the convergence to a compact relationship.

The density-tracer scatter plot is now compared to the virtual tracer plotted as a function of density, φ v (ρ * ) for configuration III. Figure 14 plots the end time (t = 329) density-tracer scatter plot (indicated by primarily blue dots), which has converged to a compact relationship that recreates φ * (ρ * ) (solid grey line overlapping the points of the scatter plot). This convergence is to be expected, since φ * represents the mean value of the tracer over a given density surface. φ v (ρ * ) is plotted at different times as dashed grey lines, showing that the shape of the virtual tracer curve converges to a shape similar to the compact scatter plot relationship and φ * (ρ * ). Both figures 12(o) and 14 indicate that although φ v (ρ * ) underestimates the amount of tracer at the central region around ρ * = 0, and overestimates the amount of tracer along the two outer regions, φ * and φ v are in relatively good agreement.

As is the case for the scatter plots (and thus φ * (ρ * )), φ v (ρ * ) has tendency toward nearly piecewise linear relationships. Indeed, since φ v is a 1D profile, it can be rewritten as

φ v (z * , t) = F (ρ * , t).
(6.1)

Subbing equation ( 6.1) into equation (5.6), a new equation in terms of F can be written as

∂F ∂t - ∂ 2 F ∂ρ 2 * K ρ ∂ρ * ∂z 2 = 0. (6.2)
It is reasonable to assume that after enough time has elapsed (as demonstrated in figure 14), mixing will cause φ v to converge to some final steady shape (only slowly evolving because of the molecular diffusivity). Thus, if F tends toward a steady state, ∂F ∂t can be neglected from equation (6.2). This requires that, provided the effective diffusivity and gradient of the density profile are non-zero, ∂ 2 F ∂ρ 2 = 0, implying that F should be a linear function of ρ * . This hypothetical tendency towards linearity can be used as a framework in which to compare the final stages of φ v (ρ * ) and φ * (ρ * ).

The mixing event sets the distribution of density weights within a given density range, say [ρ L , ρ H ], which will be affected by turbulent mixing over the course of the entire event. Within this range, the total amount of passive tracer is conserved. Therefore, assuming symmetric mixing across the mid-density, a unique ideal piecewise linear tent shape can be entirely determined by the distribution of the density and the total amount of tracer within this range. It can be defined as

φ I (ρ * ) =      (ρ * -ρ L ) 2φMax ρ H -ρ L , ρ * ∈ [ρ L , ρ M ] , -(ρ * -ρ H ) 2φMax ρ H -ρ L , ρ * ∈ [ρ M , ρ H ] , (6.3)
where the maximum tracer concentration value φ Max occurs at the mid-density defined as ρ M = (ρ L + ρ H ) /2, and is given by solving

ρ H ρ L φ I (ρ)W (ρ)dρ = ρ H ρ L φ(ρ)W (ρ)dρ. (6.4)
Figure 14 depicts φ I (ρ * ) for configuration III as solid black lines. Though there are important differences that must be considered, the general shape of φ I (ρ * ) is close to that of φ v (ρ * ) and φ * (ρ * ). Primarily, the maximum of the tracer concentration is greater for the ideal tent shape. Given the conservation of total amount of tracer in the mixing region, this discrepancy requires that part of the external branches of φ v (ρ * ) and φ * (ρ * ) lie above the linear regions of the ideal tent shape. This indicates that φ v (ρ * ) and φ * (ρ * ) have slightly convex shapes. The general idea is that the convergence could be associated with a tendency for scatter plots and profiles to stay close to the ideal piecewise linear tent shape, which is uniquely determined by the initial condition and the mixing process. However, this ideal shape is difficult to achieve through physical mixing processes, because there is mixing across the central region, which leads to a smooth transition between branches. In addition, the symmetry assumed by piecewise linear shape is not exact in these simulations and in general. There thus exist differences between all profiles, and what happens in the central region is of particular importance for determining the shape of the late time profiles and their differences. For example, the lack of convergence between φ v (z) and φ * (z) for configuration I can be explained by the fact that, in the 2D simulation, the scatter plot does not tend to a compact relationship. There is a prominent vertical branch at the mid-density (in the central region) corresponding to incomplete homogenization of tracer along that density surface. Such a feature is impossible to achieve for the virtual profile. However, configuration I seems anomalous for the dynamical configurations tested here, as configurations II and III (the most physically relevant case) demonstrate convergence toward similar vertical profiles. An important consequence is that the effective diffusivity of density provides an excellent approximation for the tracer diffusivity in these flows.

To conclude, a no more precise interpretation behind the convergence to similar compact shapes and the physics involved can be offered. But it seems related to both the systematic compensation of the eddy term by the diffusive term in equation (5.7) (see figure 13) and the tendency to converge to compact scatter plots with nearly linear side branches and similar to the unique ideal tent shape. The latter being entirely determined by the total amount of tracer in the mixing region, the details of the initial tracer distribution should not strongly influence the final profile. This idea is tested in the following sections, which present additional simulations where the initial tracer structure is varied.

Sensitivity to variations in the initial tracer field

The comparison of the three different dynamical configurations only considers a tracer that is centred along and of similar width to the pycnocline. The tracer is fully distributed within the density range affected by turbulence and mixing, and its initial gradient is similar in magnitude and geometry to the density. This section examines the mixing of tracer layers not initially collocated with or of similar width to the pycnocline in the three-dimensional dynamical configuration III. The background shear and density profiles are maintained as before, but the width and vertical position of the initial tracer profiles are varied by modifying parameters a and b in equation (2.9), respectively.

Sensitivity to vertical extent of the tracer

As argued in section 6, the final profile should only depend on the initial values of the tracer at ρ = ρ H and ρ = ρ L , and the integral quantity of the tracer within the density range [ρ L , ρ H ] affected by the dynamics of the flow. To test this principle, the width and maximum of the initial tracer profile are varied following

φ (x, t = 0) = a sech 2 (az) , (7.1)
where a is a dimensionless parameter that determines the width of the tracer distribution, while ensuring the total amount of tracer remains constant when a is varied. In the reference configuration for the tracer (i.e., the configuration presented in previous sections), a = 1. When a > 1, the width of the tracer layer is narrower than the that of the reference configuration, while the initial maximum is greater. Therefore, the tracer in these simulations initially resides within the region subject to mixing. The left-hand column of figure 15 presents the evolution of the background tracer profiles (solid lines) and virtual profiles (dashed lines) for three passive tracers with initial fields given by equation (7.1)) with a = 1, 2, and 5 (top to bottom), all subject to the same density mixing obtained by configuration III. The right-hand column depicts the initial and final profiles for the rearranged (solid blue and black lines, respectively) and final virtual profiles (dashed red line), while profiles at specific times are depicted in figure 16.

During the development of the fastest growing instabilities, and prior to the formation of secondary instabilities around t = 56, the virtual and background profiles agree almost perfectly for all three tracers. After the onset of the secondary instabilities, the effective diffusivity overestimates the width of the tracer layers, and the virtual tracer layer is broader than the background tracer layer for all three tracers. At intermediate times (e.g., t = 86), the virtual profiles are smoother than the background profiles, which have slowly changing, low concentration values at the edges, with sharp increases leading towards flatter maximums towards the middle of the domain. This trend continues for the duration of the simulation, though the background and virtual profiles show much better agreement by the end. In each case, the virtual profile tends to underestimate the tracer concentration at the middle of the domain, and overestimate the amount at the edges of the layer. Despite the variation in initial profile widths and maximums, the final profiles for each of the tracers are all of similar sizes and shapes, as predicted by the arguments presented in section 6.

Figure 17 presents the evolution of density-tracer scatter plots for each tracer layer with different initial widths. The tracer axis has been normalized by the maximum value for the wide tracer layer (a = 1) at each time step to better compare scatter plots. Additionally, scatter plots are magnified in the vertical direction between the depicted times. During the initial roll-up of the KH billows, as depicted in figure 17(a,g,m), there is an increase in scattering along φ-space as a increases. By t = 56 at the onset of secondary instabilities, there is significant scattering induced by the turbulent mixing, filling the area below the initial curves. The maximum concentration of each tracer has been reduced by approximately 20 %, and there are strong discrepancies between all scatter plots. The outer edges of the scatter plots collapse to protrusions mostly centred at ρ = 0 (t = 85) while maintaining similar maximum tracer values to the previous times. By t = 140, the difference in maximum tracer values of each of the scatter plots has greatly reduced. The amount of scatter over a given density bin varies, however, Figure 17. Density-tracer scatter plots for tracer layers centred along the pycnocline with varying widths, but identical tracer totals (left to right: wide (a = 1), middle (a = 2), and narrow (a = 5)). Note that the tracer axis is magnified after each time step as the scatter plots collapse.

increasing as the initial profile profile becomes narrower (a increasing). Convergence of the scatter plots to similar shapes becomes apparent by the end of the simulation (t = 342), with each scatter displaying compact relationships in much closer agreement, with minor differences in tracer maximum and concavity. By this point any effects of turbulent mixing have finished (see figure 10), and only diffusive mixing remains to homogenize the tracers. The density-tracer scatter plots with the narrow and medium distributions show important cross-isopycnal fluxes of tracers, with a final distribution extending to a wider density range. The amount of tracer affecting new density ranges can be calculated from the knowledge of the density evolution and the total initial tracer amount. Finally, while all three scatter plots tend to nearly piecewise linear tent shapes, the wide (reference) case (a = 1) is slightly convex on either side of the tracer maximum, while the scatter plots for the other two tracers are concave on either side of the maximum, with the narrow case (a = 5) being more concave than the medium case (a = 3). Additionally, the tracer concentration maximum increases slightly as a increases. This relationship between scatter plot concavity and tracer concentration maximum reflects the observations made about deviations away from the ideal tent shape made in section 6.

Sensitivity to vertical position of the tracer

In the results presented in previous sections, all tracer profiles are initially concentrated within the region subject to turbulent mixing. As a result, the amount of tracer initially for densities less than ρ = ρ L and greater than ρ = ρ H is null. In this section, the validity of the convergence principle is tested for configurations where the tracer structure is offset above the main turbulent region, so that its value at ρ = ρ H remains nearly constant. These initial tracer profiles are of the form,

φ (x, t = 0) = sech 2 (z -b) , (7.2)
where b indicates the offset of the tracer from the pycnocline. For this simulation, b = 0 sets a layer centred along the pycnocline (i.e., the reference tracer configuration), b = 2.86 is offset above the pycnocline by 10% of the vertical extent of the domain, and b = 5.72 is offset above the pycnocline by 20% of the vertical extent of the domain. The evolution of profiles tracers subject to these initial conditions for the three-dimensional configuration III are presented as contour slices in figure 18. From left to right, the tracers are aligned with the midpoint of the pycnocline (b = 0), slightly offset from the pycnocline (b = 2.86), and completely offset from the pycnocline (b = 5.72). The corresponding density evolution is unchanged, and is depicted in the right-hand column of figure 3. The slightly offset tracer initially resides on the edge of the fastest-growing vortex, so that when the KH billow first develops, some of the tracer is entrained into the pycnocline by the vortex, and the tracer layer is significantly distorted. With the onset of small scale features from the secondary instabilities, the bottom half of this tracer layer is eroded and a thick layer of low tracer concentration develops, extending to the lower edge of the pycnocline. These secondary instabilities affect the interior and upper parts of the high concentration region, but do not significantly mix regions of low and high tracer concentration (e.g., figure 18(k)), and there remains a layer of relatively high tracer concentration mostly unaffected by mixing. Along the upper edge of this layer, molecular diffusion is the primary source of slow mixing, but this mixing is weak enough to be considered negligible when compared to the strong mixing happening below. The layer of tracer completely offset from the pycnocline is weakly distorted by the outer edge of the billow, but experiences almost no mechanical mixing (some strands of tracer pulled into the secondary instabilities are visible in figure 18(q)), mostly a spreading of tracer due to diffusion. The effect of the localized mechanical mixing away from the offset tracer layers is visible in the background tracer profiles. The complete evolution of the background and virtual tracer profiles is presented in the left-hand column of figure 19, with the initial and final profiles presented in the right-hand column. Profiles at specific times corresponding to the contour plots of figure 18 are presented in figure 20.

Compared to the tracers collocated with the pycnocline presented in previous sections, the offset background tracer profiles typically show greater agreement with the virtual profiles throughout the evolution of the flow. The slightly offset tracer shows some disparity at the lower edge of the virtual and background profiles after the onset of the secondary instabilities, and the virtual profile underestimates the maximum value of the profile towards the midpoint of the simulation (e.g., t = 140). For the completely offset tracer, the virtual and background profiles appear to agree perfectly for the duration of the simulation. This can be explained by considering that the tracer is located mostly outside the region of strong turbulent mixing and is essentially subject only to slow molecular diffusion. The tracer eddy term is almost null for this tracer distribution, since both φ and the gradient of density are nearly zero over the density levels in that region. The rearranged tracer profile equation therefore reduces to the virtual tracer equation in this case.

Figure 21 presents the scatter plots of each of the tracers in figure 18 relative to density. The first thing to note is the offset tracer scatter plots take a different form than the centred tracer scatter plot. They are asymmetric, and have different convex envelopes (essentially triangles with corners at the top left, bottom left, and bottom right corner of the graphs). The slightly offset tracer (b = 2.86) experiences most of its scattering above the concave region of the plot, and tends to a final shape that is nearly piecewise linear, with lines of two different slopes, and a change of slope at the density mean. The tracer maximum occurs along an isopycnal with a density lower than the mean. The cross-isopycnal flux of tracer is strong in this case, since at the end of the simulation, a significant quantity of the tracer is distributed over a much wider density range. For most of the flow, the completely offset tracer is modestly affected by mixing and does not show much scattering, though some redistribution of the scatter plot is visible at the lower edge of the tracer layer where it is more significantly eroded by mixing. As for the experiences with different tracer width, the present offset configurations show that mixing can strongly modify tracer concentration and entrain significant amounts of tracer into new density ranges.

Conclusions

This paper uses the classical problem of turbulence arising from Kelvin-Helmholtz instabilities to investigate the mixing of passive tracers. Three different stratified shear configurations providing different routes to turbulence are simulated, and the mixing of different passive tracers layers is examined. Simulations examining the mixing of tracer layers with different widths and vertical positions are performed. Weighted density-tracer scatter plots are proposed as a method through which to analyze diapycnal transport of tracer. When tracer and density diffusivity are equal, a convex envelope constraint can be placed on the evolving scatter plot, which limits on the maximum tracer values achievable in specific density ranges. When the fluid motion is adiabatic, there is no modification of the scatter plots. However, adiabatic rearrangement during the development of the billows leads to stretching of the density and tracer fields extends the surface area over which regions of different density and tracer concentrations can interact. Diabatic mixing then yields irreversible modifications of the density-tracer scatter plot over the density range affected by mixing. For the configurations depicted here, this first generates strong scattering, but there is a convergence towards tighter relationships following each scattering event. The small-scale stretching (e.g., due to secondary instabilities) appears to play an important role in the tightening of scatter plots, as this creates regions of tracer variation closer to the scale of the homogenization cell set by the diffusivity coefficients. The formulation for the scatter variance as introduced by equation (4.7) and presented in figure 7 provides an additional diagnostic to compare the evolution of the the tracer layer collocated with the pycnocline in each of the three dynamical configurations. Times at which the diagnostic is near zero can also be used to identify if the relationship between the density and tracer has become compact. This can supplement the information available in the weighted scatter plots, as it can be difficult to ascertain whether a compact density-tracer relationship has been reached based only on graphical analysis.

The use of effective diffusivity as a model for the mixing of the passive tracer is also examined. It appears that when the scatter plots display less variation in tracer-space along a given density bin (less scattering), the effective diffusivity approach gives a reliable estimate of the effect of mixing on the tracer after the mixing event as a whole, but there are more significant differences between the prediction and simulation at intermediate times.

The approach taken in this paper is to investigate, for φ representing the concentration of a passive tracer, to what extent most of the evolution of φ * is captured by (5.5), with the eddy term ignored. A question beyond the scope of this paper is whether, for cases where the eddy term cannot be ignored, can it be estimated using a closure assumption and whether this can be accomplished whilst retaining the advantages that have been demonstrated for the standard use of a tracer-based coordinate and the corresponding effective diffusivity. Note that [START_REF] Shuckburgh | Diagnosing transport and mixing using a tracerbased coordinate system[END_REF] consider whether particle transport can be approximated using an effective diffusivity calculated from a tracer field, i.e., to κ φ = 0, corresponding to Example 2 in Appendix C. Some success (at least semi-quantitative) was demonstrated. In terms of the new results given above, this would be equivalent to This manuscript examines a specific set of tracer transport experiments in which the diffusivity of the tracer is equal to the diffusivity of the density. In the case where κ ρ > κ φ (for example, in an oceanic environment where the density is controlled by temperature, and the passive tracer is some dissolved chemical species), preliminary simulations have shown a decent agreement between the virtual and background tracer profiles in the sense that it it still provides a valuable estimate, though there remains room for improvement. More work is needed to investigate the effect of differential molecular diffusivity. However, equation (5.5) shows that setting κ ρ = κ φ can change the dominant term in the equation, which can alter how well the virtual tracer represents the background tracer evolution. Additionally, when the density and tracer diffuse at different rates, this can drastically alter the diapycnal fluxes of tracers, and in particular breaks the convex envelope constraint placed on the scatter. In the hypothetical limit where only density diffuses and the passive tracer is non-diffusive, the displacement of the points on the scatter plot would be strictly horizontal, eventually converging to a vertical line located at the mean density. Similarly, if only the passive tracer were diffusive and the density is non-diffusive, the displacement of points would be entirely vertical, leading to a horizontal line at the tracer mean. It therefore remains to be determined what constraints can be placed on scatter plots for tracers with different diffusivities.

∂ ∂z * -φ κ ρ ∇ 2 ρ z * (∂z * /∂ρ) ∂ ∂z * K ρ ∂ φ z * ∂z * . ( 8
It also remains to be determined if these scatter plots provide a useful analysis in the event where the tracer is an active component of the density. For example, if ρ is a function of φ only, then a ρφ-scatter plot would simply depict the functional form of ρ = ρ (φ). However, if ρ is a function of two tracers, ρ = ρ(φ 1 , φ 2 ), it remains to be seen what information could be derived from ρφ 1 or ρφ 2 scatter plots. It is also worth considering a situation in which one is concerned with the mixing of a passive tracer when density is a function of two different active tracers, such as for double-diffusive flows (for example in the ocean, where salt and temperature are the active tracers). Because double-diffusion allows for the creation of new density anomalies, this immediately presents the possibility of breaking the convex envelope constraint presented here.

A complete understanding of convergence to compact density-tracer relationships has not been achieved in this study. The physics of tracer redistribution along streamlines have been examined from a theoretical perspective for steady and oscillatory 2D flows by [START_REF] Rhines | How rapidly is a passive scalar mixed within closed streamlines[END_REF]. They indicate that during a rapid mixing phase, passive tracers will tend to homogenize along streamlines. Thus the circulation patterns may play a significant role in the structure of tracers (including density) and the convergence of scatter plots to compact forms. However, these arguments cannot be easily extended to the configurations studied here, as density is not a passive tracer, there are strong 3D effects (for configuration III), and convergence is achieved during the chaotic mixing phase (see figure 7) when no clear streamlines can be identified.

As a final note, the presence of a free-surface and a compressible context are specific aspects of the numerical model used in the present study. It potentially adds two dynamical processes to the flow: surface waves and acoustic waves. The dynamical configurations have been selected (with the development of KH instability far from the top free-surface) so that there is no significant impact of these waves on these results (see Appendix A). However, surface and acoustic waves are realistic physical processes, which can play a role in modifying the turbulent cascade in certain circumstances (see, for example, [START_REF] Shete | Effect of free surface on submerged stratified shear instabilities[END_REF]). one can solve for the diffusivity at a given time by multiplying the equation by the tracer field, then integrating over the entire domain, giving

κ ρ = -V ∂ρ 2 ∂t dV 2 V |∇ρ| 2 dV . (B 2)
For an ideal numerical simulation with perfect numerical schemes, the resulting value of κ ρ would be constant and equal to the prescribed value at all times. In reality, numerical schemes are imperfect, and the resulting effects due to extra terms can act as extra diffusion inherent to the schemes. The resulting "net diffusivity" will vary with time, and depend on the schemes used, the grid size, the time-step, and the flow geometry,

κ Net ρ (t) = -V ∂ρ 2 ∂t dV 2 V |∇ρ| 2 dV . ( B 3)
As an example, the net diffusivities as calculated by equation (B 3) and normalized by the prescribed molecular value (5 × 10 -3 m 2 /s) are presented as a function of time in figure 22 for the density and passive tracer fields with initially different vertical positions (as described in section 7.2) of configuration III. The net estimates all begin at the prescribed molecular value, with the time-series for the density and the b = 0 tracer showing similar overall trends: an initial increase of approximately 100% corresponding to the stretching of the field interfaces by the initial billow formation before decreasing to near the prescribed value as the secondary instabilities collapse the billow. The evolution of the net diffusivity of the middle tracer (b = 2.86) is similar to that of the density and b = 0 tracer, though it only sees a maximum increase of approximately 70%. In contrast, the upper layer of tracer (b = 5.72) experiences a much smaller increase in net diffusivity of about 10% when the billow initially distorts the field. This increase away from the maximum is much less significant and shorter lived, as there is little meaningful interaction being the upper tracer and the dynamic mixing process. For this study, the net diffusivity was calculated in order to account for the effects of numerical processes when calculating profile evolution. Figure 23 depicts the complete time evolution (with figure 24 depicting profiles at individual times) of the background density profile as a calculated from the rearrangement of the density from simulations (solid lines) versus the profiles calculated from equation (5.1) using the constant prescribed diffusivity value κ ρ for the simulation (left panel), and the time-dependent net diffusivity κ Net ρ (t) computed from the simulation using equation (B 3). As shown in figure 23(a), the isopycnals calculated from the diffusion equation typically lie within the sorted isopycnals, indicating that the the constant diffusivity leads to the spreading of the pycnocline being underestimated. Using the time-dependent net diffusivity shows better agreement between the calculated and sorted isopycnals. This is also apparent in the later times shown in figure 24, which also demonstrates that the profiles calculated using the prescribed diffusivity underestimate the steepness of transition at ρ * = 0. with the velocity field u(x, t) being non-divergent. Note that for convenience, all variables presented in this appendix are given in their dimensional forms. Following [START_REF] Nakamura | Two-Dimensional Mixing, Edge Formation, and Permeability Diagnosed in an Area Coordinate[END_REF] and [START_REF] Winters | Diascalar flux and the rate of fluid mixing[END_REF], the density field may be used to define a new vertical coordinate z * , as follows. For a specified value ρ * of the density, define (t) (dashed red lines), and calculated from equation (5.1) using the prescribed, constant diffusivity κρ (dashed blue lines).

z * (ρ * , t) = 1 A d ρ(x,t) ρ * dV. (C 2)
Note that z * is a monotonically decreasing function of ρ * , so the inverse function ρ * (z * , t) is well-defined. Now for any field φ(x, t) define This is the case considered by [START_REF] Nakamura | Two-Dimensional Mixing, Edge Formation, and Permeability Diagnosed in an Area Coordinate[END_REF] and [START_REF] Winters | Diascalar flux and the rate of fluid mixing[END_REF]. Note that (C 11) κ ρ |∇ρ| dA , (C 16) where the last equality follows by first using the divergence theorem to transform the volume integral to a surface integral. Then, using the fact that for the term including u † Defining a mean of one quantity over isosurfaces of a different quantity, as suggested here, is a technique often employed in atmospheric sciences (e.g., [START_REF] Butchart | The Area of the Stratospheric Polar Vortex as a Diagnostic for Tracer Transport on an Isentropic Surface[END_REF][START_REF] Nakamura | Modified Lagrangian-Mean Diagnostics of the Stratospheric Polar Vortices. Part I. Formulation and Analysis of GFDL SKYHI GCM[END_REF].

Φ(ρ * , t) = ρ(x,t) ρ * φ(x, t)dV. ( C
the factor ρ is constant on the surface S(z * ) and can be taken outside the integral, the incompressibility condition can be applied so that the integral vanishes. The final step is to use (C 7) to introduce a factor ∂ρ * /∂z * , giving The question is how to organize (C 13) in this case. The key difference is that φ varies on the coordinate surfaces defined by z * . Therefore define φ = φ -φ z * to capture this variation. Assume that φ satisfies an advection-diffusion equation, i.e., φ t + u • ∇φ = κ φ ∇ 2 φ, allowing the diffusivity of φ to be different from κ ρ . Consider The term inside the brackets on the right-hand side of the first equality comes from writing the volume integral of κ φ ∇ 2 φ -u • ∇φ as a surface integral using the divergence theorem. In the second equality, the term containing u φ z * is zero by the same reasoning as applied in the previous section. In the second term in brackets on the right-hand side the quantity ∇ φ z * can be rewritten as (∂ φ z * /∂ρ * )∇ρ and then, again, using the same reasoning as in the previous section, the entire term can be written in terms of K ρ .

Inserting in (C 13) implies that This is the governing equation for φ z * . Some checks on the form of the above equation are as follows.

Example 1: When κ ρ = 0 and κ φ = 0 both ρ and φ are materially conserved and therefore, since ρ-surfaces are material surfaces, there will be no change in the total amount of φ within each ρ-surface. Furthermore the volume within each ρ-surface remains constant. Therefore for each z * , φ z * is constant in time. Equation (C 22) is consistent with this.

Example 2: When κ ρ = 0 and κ φ = 0 then φ is materially conserved. The rate of change of the total amount of φ within the volume z * is therefore given by

Figure 1 .

 1 Figure 1. (Left to right) Initial vertical profiles of the velocity, density, and passive tracer fields for the described simulations with relevant dimensional parameters. Ri0 Re0 Pr Sc 0.1158 2000 1 1

Figure 2 .

 2 Figure 2. Time evolution of the mean kinetic energy K for configuration I (dashed line), configuration II (solid line), and configuration III (dash-dotted line). The labelled times correspond to to the rows of contours plotted in figures 3 and 4.

Figure 3 .

 3 Figure 3. Evolution slices of the density field for configurations I (left, (a)-(e)), II (middle, (f)-(j)), and III (right, in the streamwise (k)-(n) and spanwise (o)-(r) directions) at the times indicated in figure 2.

Figure 4 .

 4 Figure 4. Evolution slices of the passive tracer field for configurations I (left, (a)-(e)), II (middle, (f)-(j)), and III (right, in the streamwise (k)-(n) and spanwise (o)-(r) directions) at the times indicated in figure 2.

Figure 5 .

 5 Figure5. Sketch of the initial scatter plot of tracer as a function of density associated with the profiles given by equations (2.4) and (2.9) where the pycnocline and tracer layer are collocated (solid line). Hypothetical evolution of the scatter plot after a mixing event showing diapycnal fluxes of tracer (dashed line). In the new profile, mixing has led to a flux of tracer toward lower densities.

Figure 6 .

 6 Figure 6. Weighted density-tracer scatter plots for configurations I, II and III (left to right). Blue points indicate the density-tracer pairs occupy a low percentage of the physical domain, yellow points indicate an intermediate percentage, and red points indicate a high percentage.

Figure 7 .

 7 Figure 7. Perturbation kinetic energies (black lines), and scatter variance (grey lines) for (a) configuration I, (b) configuration II, and (c) configuration III. Configurations I and II depict only the 2D kinetic energy, while configuration III depicts both the 2D kinetic energy (solid line) and 3D kinetic energy (dotted line). The values for the kinetic energy are depicted on the left-hand axes, while the values for the scatter variance are depicted on the right-hand axes.

Figure 8 .

 8 Figure 8. Diagram illustrating the convex envelope constraint on the evolution for a typical density-tracer scatter plot.

Figure 9 .

 9 Figure 9. Time evolution of the background density profile from the simulations (solid lines) and as calculated from the diffusion equation (5.1) (dashed lines) for the different dynamical configurations (a,b,c). The final density profiles are presented in the right-hand column (d,e,f), with the background profile given by the solid black line, and the profile from the diffusion equation given by the solid red line. The initial profile is depicted by the blue line. The solid black and dashed red lines overlap almost perfectly.

Figure 10 .

 10 Figure 10. Example of dimensionless effective diffusivity computed from the sorted density profile for configuration III. The colour indicates the dimensionless values of effective diffusivity, while the white lines indicate different values of the background density profile. The effective diffusivity presented here is not scaled by the net molecular diffusivity.

Figure 11 .

 11 Figure11. Time evolution of the virtual passive tracer profiles (dashed lines) and the background passive tracer field from simulations (solid lines) (a,b,c). The final passive tracers profiles are presented in the right-hand column (d,e,f), with the rearranged profile given by the solid black line, and the virtual profile given by the solid red line. The initial profile is given by the blue line.

Figure 12 .

 12 Figure 12. Background tracer profiles (solid black lines) and virtual tracer profiles (dashed red lines) for the three dynamical configurations (left to right: configuration I, configuration II, and configuration III).

Figure 13 .

 13 Figure 13. Time evolution of the profile of the diffusive term (top) and complete eddy term (bottom) for the 3D dynamical configuration (III).

Figure 14 .

 14 Figure 14. Evolution of the virtual tracer profiles (dashed lines) as a function of rearranged density for the default tracer field of configuration III. The blue and yellow points indicate the density-tracer scatter plot at t = 329, with the solid grey line indicating the rearranged tracer profile as a function of background density at t = 329. The solid black line indicates the ideal tent shape at t = 329.

Figure 15 .

 15 Figure15. (Left) evolution of the background (solid lines) and virtual (dashed lines) tracer profiles with time (top to bottom: wide (a = 1), medium (a = 2), and narrow (a = 5)). (Right) Initial tracer profiles (blue), final background profile (black), and final virtual profile (red).

Figure 16 .

 16 Figure16. Background tracer profiles (solid black lines) and virtual tracer profiles (dashed red lines) for the variable tracer layer width simulations (left to right: wide (a = 1), medium (a = 2), and narrow (a = 5)).

Figure 18 .

 18 Figure 18. Contour plots of the evolution of tracer layers centred at various depths (left to right: along the pycnocline (b = 0), slight offset (b = 2.86), and completely offset (b = 5.72)).

Figure 19 .

 19 Figure 19. (Left) evolution of tracer profiles with time (top to bottom: along the pycnocline (b = 0), slightly offset from the pycnocline (b = 2.86), and completely offset from the pycnocline (b = 5.72)), as sorted from the simulation (solid lines) and as predicted from the diffusion equation. (Right) Initial tracer profiles (blue), final background tracer profile (black), and final virtual profile (red).

Figure 20 .

 20 Figure 20. Background tracer profiles (solid black lines) and virtual tracer profiles (dashed red lines) for the simulation where the passive tracers are offset from the pycnocline (left to right: along the pycnocline (b = 0), slightly offset from the pycnocline (b = 2.86), and completely offset from the pycnocline (b = 5.72)).

Figure 21 .

 21 Figure 21. Density-tracer scatter plots for tracer layers centred at various depths (left to right: along the pycnocline (b = 0), slight offset (b = 2.86), and completely offset (b = 5.72)).

Figure 22 .

 22 Figure 22. Evolution of the estimate of global effective diffusivity for the density (solid line) and passive tracer (dashed lines) fields for each of the offset tracers described in section 7.2 normalized by the prescribed molecular diffusivity.

Figure 23 .

 23 Figure 23. Evolution the density profiles calculated from equation (5.1) using the constant prescribed value of diffusivity (left) and a time-dependent value of the diffusivity as calculated from equation (B 3).

Figure 24 .

 24 Figure 24. Background density profiles as calculated from simulation outputs (solid black line), calculated from equation (5.1) using the time-dependent, net diffusivity κ Net

  ρ

  implies the identityf (ρ)ψ z * = f ( ρ z * ) ψ z * = f (ρ) ψ z * (C 14)for any function f (•) and for any quantity ψ. In particular ρ t z * ρ z * = ρρ t z * and hence from equation (C 13)( ρ z * ) t = ρ t z * . (C 15)Because ρ z * and ρ are the same, and ρ * is the notation used for the value of ρ on the ρ-surface which is labelled by z * , ρ * will be used instead of ρ z * .The standard result on effective diffusivity follows if ρ is the solution of an advectiondiffusion equation, with incompressible flow field u(x, t) and diffusivity κ ρ . Then from the above results,

(

  diffusivity K ρ is defined by the second equality and hence, via (C 7) and (C 11), is given byK ρ = κ ρ |∇ρ| 2 z * (∂z * /∂ρ * ) 2 . (C 18)C.3. The non-standard case: φ = ρ.

  ∇φ -φ u] • dA + S(z * ) κ φ ∇ φ z * • dA (C 19)

(

  (z * ) [κ φ ∇φ -φ u] • dA .(C 20) Now consider the two quantities in the bracket in the second term on the right-hand side, which can be combined as a surface integral, using (C 7) and (C 11), plus the fact that dA = ∇ρ dA/|∇ρ|, to give-φ ρ t z * ∂z * ∂ρ * + 1 A d S(z * ) [κ φ ∇φ -φ u] • dA = 1 A d S(z * ) -κ ρ φ ∇ 2 ρ + φ u • ∇ρ + κ φ ∇φ • ∇ρ -φ u • ∇ρ |∇ρ| dA = ∂z * ∂ρ * (κ φ ∇φ • ∇ρ -κ ρ φ ∇ 2 ρ) z * (C 21)Hence, substituting into (C 19),( φ z * ) t = ∂ ∂z * κ φ κ ρ K ρ ∂ φ z * ∂z * + ∂ ∂z * ∂z * ∂ρ * (κ φ ∇φ • ∇ρ -κ ρ φ ∇ 2 ρ) z * . (C 22)

  S -u) • ∇ρ dA |∇ρ| = φ(u S -u) • ∇ρ z * (∂z * /∂ρ * ),(C 23) where u S is the local normal velocity of the surface S(z * ). Since S(z * ) is a surface of constant ρ (but the value of ρ on the surface varies in time), on S(z * ),ρ t + u S • ∇ρ = κ ρ ∇ 2 ρ + (u S -u) • ∇ρ = α(t), (C 24)where α(t) is to be determined by the condition that the volume within the surface S(z * ) stays constant, requiring thatS(z * ) u S • dA = S(z * ) u S • ∇ρ dA |∇ρ| = u S • ∇ρ z * (∂z * /∂ρ * ) = 0. (C 25)Incompressibility implies (e.g., replace u S by u in the preceding equation) that u • ∇ρ z * = 0. Therefore applying • z * to (C 24) implies thatα(t) = κ ρ ∇ 2 ρ z * (C 26)and it follows also from (C 24) that (u S -u)• ∇ρ = κ∇ 2 ρ z * -κ∇ 2 ρ. Additionally (u S -u) • ∇ρ z * = 0 implies that φ(u S -u) • ∇ρ z * = φ (u S -u) • ∇ρ z * = φ ( κ ρ ∇ 2 ρ z * -κ ρ ∇ 2 ρ) z * = -φ κ ρ ∇ 2 ρ z * . (C 27) It follows from (C 23) that ( φ z * ) t = ∂ ∂z * -φ κ ρ ∇ 2 ρ z * (∂z * /∂ρ * ) , (C 28) consistent with (C 22) when κ φ = 0. Example 3: When (formally) κ ρ = 0 but κ φ = 0, then as in Example 1 ρ-surfaces are material surfaces, therefore there is no advective flux of φ across ρ-surfaces. Therefore the rate of change of the total amount of φ within a ρ-surface is equal to the diffusive flux across the surface. φ = φ(ρ) z * + φ and hence ∇φ = (∂ φ z * )/(∂ρ)∇ρ + ∇φ . It follows that d dt z * φ dV = κ φ ∂ φ z * ∂ρ S(z * ) ∇ρ • dA + S(z * )

Table 1 .

 1 Prescribed dimensional parameters consistent across all simulations presented in this paper.

Table 2 .

 2 Prescribed

dimensionless parameters consistent across all simulations presented in this paper.

Table 3 .

 3 Relevant computational parameters of the three dynamical configurations presented in this paper.

  , as an independent variable, i.e. to regard ρ * as a function of z * , and t. Again, using (C 4) and (C 5), where the final equality follows from (C 7), noting that ∂z * /∂ρ * = {∂ρ * /∂z * } -1 , where both partial derivatives are evaluated for constant t. † To reflect the fact that z * is from now on going to be an independent variable, rather than ρ * , the notation S(ρ * ) is replaced by S(z * ) (but recall that ρ = ρ * on the surface S(z * )).Note that this implies thatA -1 d S(ρ * ) φ dA/|∇ρ| = φ z * (∂ ρ z * /∂z * ) -1 , so (C 10) is , t) = ρ t z * φ z * -φρ t z * (∂ ρ z * /∂z * ) -1 + 1 A d ρ(x,t) ρ * φ t dV. (C 12)Now differentiating with respect to z * gives, using (C 11) as an identity that holds for any field, ( φ z * ) t = φ t z * +

			∂ ∂t	Φ(z S(ρ * )	φ	dA |∇ρ|	S(ρ * ) ρ t dA/|∇ρ| S(ρ * ) dA/|∇ρ|	(C 10)
								+	ρ(x,t) ρ *	φ t dV +	S(ρ * )	φ ρ t	dA |∇ρ|	.
	A natural definition of the mean value of φ over a z * -co-ordinate surface is given by
	φ(z * , t) z * =	1 A d	∂ ∂z *	Φ(z * , t) =			S(z * ) φ dA/|∇ρ| S(z * ) dA/|∇ρ|	= -	1 A d	∂ρ * ∂z * S(z * )	φ dA/|∇ρ|, (C 11)
	equivalent to
	Simple vector calculus implies that ∂ ∂t Φ(ρ * , t) = ρ(x,t) ρ * 1 A d ∂ ∂t Φ(z ∂ ∂z * ρ t z * φ z * -φρ t z * φ t dV + S(ρ * )	φ ρ t ∂ ρ z * dA |∇ρ| ∂z *	,	-1	.	3) (C 4) (C 13)
	and
					∂ ∂ρ *	Φ(ρ S(ρ * )	φ	dA |∇ρ|	,	(C 5)
	where the surface S(ρ 1/A d
	that
					∂ ∂t	z * (ρ * , t) =	1 A d S(ρ * )	ρ t	dA |∇ρ|	,	(C 6)
	and
					∂ ∂ρ *	z * (ρ * , t) = -	1 A d S(ρ * )	dA |∇ρ|	.	(C 7)
	The next step is to consider z ∂ ∂t Φ(ρ(z * , t), t) = -S(ρ * ) φ |∇ρ| dA	∂ρ * ∂t	+	ρ(x,t) ρ *	φ t dV +	S(ρ * )	φ ρ t	dA |∇ρ|	. (C 8)
					S(ρ * )		ρ t	dA |∇ρ|	+	∂ρ * ∂t	S(ρ * )	dA |∇ρ|	= 0.	(C 9)

* , t) = -* ) is the boundary of the volume {x|ρ(x, t) ρ * } and dA is the area element on that surface. Note that it follows from (C 4) and (C 5) with φ = * (ρ * , t), rather than ρ * Now note that the condition on ∂ρ * /∂t for z * (ρ * , t) to remain constant in time is that -Hence using (C 9) to substitute for ∂ρ * /∂t in (C 8) it follows that, regarding Φ as function of the two variables z * and t, * , t) = -* C.2. The standard case: φ = ρ.

  Finally, differentiating with respect to z * , dividing by A d and re-arranging the right-hand side,

	( φ z * ) t = =	∂ ∂z * ∂ ∂z *	κ φ |∇ρ| 2 κ φ κ ρ K ρ ∂ φ z * z * ( ∂z *	∂z * ∂ρ * +	) 2 ∂ φ z * ∂z * ∂ ∂z * κ φ ∇φ • ∇ρ z * + ∇φ • ∇ρ z * ∂ρ ∂z * ∂z * ∂ρ * ,	(C 31)
	again consistent with (C 22) in the relevant limit.

∇φ • dA = A d κ φ ∂ φ z * ∂ρ |∇ρ| 2 z * + ∇φ • ∇ρ z * ∂z * ∂ρ .

(C 30)

† Email address for correspondence: jared.penney@legos.obs-mip.fr

† Note that typically this equation is written in terms of ρ and z * (for example,[START_REF] Smyth | Differential diffusion in breaking KelvinHelmholtz billows[END_REF], with ∂ρ ∂z * implicitly referring to the adiabatic rearrangement of the density field. The decision was made to explicitly write ρ * in parts of the equation that refer to the background density profile, while writing ρ refers to the density field in Cartesian coordinates, such as in the numerator of the right-hand side of equation (5.2).

† This is calculated in physical space by determining the mean value of tracer corresponding to the density at a given point, φ * = φ(z * , t) = φ(ρ * , t), then subtracting it from the value of the tracer at that point.

Appendix A. Non-Boussinesq, non-hydrostatic system of equations (CROCO model)

The numerical models used in the study of stratified shear generally consider that volume (rather than mass) is conserved, and that the flow satisfies the Boussinesq assumption. For non-hydrostatic, incompressible processes such as Kelvin-Helmholtz instabilities, the system of Boussinesq equations (2.1) accurately represents the flow. From a numerical point of view, the non-hydrostatic implementation of an incompressible flow degenerates to an elliptic system of equations (with respect to acoustic waves). A 3D Poisson equation must thus be solved to obtain a pressure "avatar" for the resulting flow and several types of algorithms (based on, for example, pressure projection (e.g., [START_REF] Subich | Simulation of the Navier-Stokes equations in three dimensions with a spectral collocation method[END_REF] or pressure correction (e.g., Jiang 2019)) can be implemented to simulate non-hydrostatic flows. When dealing with geophysical flows, whether in the atmosphere or in the ocean, the re-introduction of acoustic waves can lead to efficient numerical algorithms for non-hydrostatic flows (see for instance [START_REF] Janjic | An Alternative Approach to Nonhydrostatic Modeling[END_REF] for atmospheric flows and [START_REF] Auclair | A non-hydrostatic non-Boussinesq algorithm for free-surface ocean modelling[END_REF] for oceanic flows). Several reasons can explain this at-firstglance paradoxical conclusion: energy transfers are simulated in a presumably consistent way [START_REF] Tailleux | On the energetics of stratified turbulent mixing, irreversible thermodynamics, Boussinesq models and the ocean heat engine controversy[END_REF], and the resulting computations are "local" and they are thus well-adapted to massively-parallel implementations. In the ocean, the re-introduction of acoustic waves additionally provides a convenient way to simulate non-hydrostatic flows with time-splitting algorithms [START_REF] Auclair | A non-hydrostatic non-Boussinesq algorithm for free-surface ocean modelling[END_REF].

The numerical simulations presented in this article were performed using the CROCO ocean model. CROCO is a free-surface, non-hydrostatic, non-Boussinesq model. It inherited from ROMS and ROMS-AGRIF [START_REF] Shchepetkin | The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model[END_REF][START_REF] Debreu | Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation[END_REF]) their numerically-efficient two-mode leapfrog-third-order Adams-Moulton (LF-AM3) time-splitting. The non-hydrostatic implementation of this time-splitting is conveniently obtained by calling into question the Boussinesq assumption and reintroducing acoustic waves together. The resulting compressible processes are treated in the LF-AM3 fast-mode which cannot remain 2D (barotropic) and becomes 3D. The resulting system of equations is given in dimensional form by

where λ is the bulk viscosity. Density components ρ h and δρ are respectively the hydrostatic, Boussinesq and the non-hydrostatic, non-Boussinesq contributions to density.

In the current implementation, the latter compressible density anomaly is linked to the pressure by the first-order relation δρ = c -2 s δp, with δp = p -p h , where p h is the hydrostatic component of pressure and c s is the speed of sound. As in the original ROMS-AGRIF hydrostatic time-splitting, the free-surface anomaly is prognosticated in CROCO fast-mode. The speed of sound, bulk viscosity, and number of fast non-Boussinesq time steps per external Boussinesq time step, N ∆tFast , are given in table 4.

Concerning first the consequences of the explicit modeling of acoustic waves, sensitivity tests have been carried out in the frame of the present study. The propagation speed of acoustic waves has been modified, and no significant consequences have been found on the evolution of KH instabilities. As far as the consequences of the surface-induced processes are concerned, the studied configurations are based on deep pycnocline and shear-stress layers, hence minimizing the interactions between the free-surface and Kelvin-Helmholtz instabilities. Sensitivity studies with a deeper water column (and a deeper pycnocline and shear layer) have shown that the development of the instability is not influenced by the free surface. Only a slight influence in the difference between upper and lower boundary conditions was observed during the pairing and subsequent mixing in the twodimensional case with two billows, as the vertical extent of the final billow becomes large. This leads to an asymmetry of the mixing at the edge of the billow, but this effect does not modify the conclusions presented here.

Appendix B. Net diffusivity estimate