# Diagnosing tracer transport in convective penetration of a stably stratified
layer

Charles W. Powell1<EMAIL_ADDRESS>Peter H. Haynes1  John R. Taylor1  1DAMTP,
University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road,
Cambridge CB3 0WA, UK

###### Abstract

We use large-eddy simulations to study the penetration of a buoyant plume
carrying a passive tracer into a stably stratified layer with constant
buoyancy frequency. Using a buoyancy-tracer volume distribution, we develop a
method for objectively partitioning buoyancy-tracer space into three regions,
each of which corresponds to a coherent region in physical space.
Specifically, we identify a source region where undiluted plume fluid enters
the stratified layer, a transport region where much of the transition from
undiluted to mixed fluid occurs in the plume cap, and an accumulation region
corresponding to the radially spreading intrusion. This method enables
quantification of different measures of turbulence and mixing within each of
the three regions, including potential energy and turbulent kinetic energy
dissipation rates, an activity parameter, and the instantaneous mixing
efficiency. We find that the most intense buoyancy gradients are found in a
thin layer at the cap of the penetrating plume, which forms the primary mixing
stage and exhibits the most efficient mixing. Newly generated mixtures of
environmental and plume fluid joining the intrusion are subjected to
relatively weak turbulence and buoyancy gradients weaken as subsiding fluid is
homogenised. As the intrusion spreads radially, environmental fluid
surrounding the intrusion is efficiently mixed into the intrusion which
dominates the total entrainment of environmental fluid into the plume as
whole. However, the entrainment rate is largest in the plume cap where the
most buoyant fluid is entrained.

###### keywords:

##  1 Introduction

The interaction between active convection and neighbouring stably stratified
regions is relevant to many geophysical flows. An important example is the
tropical upper troposphere and lower stratosphere, where convective plumes
generated by strong thunderstorm complexes can penetrate through the tropical
tropopause layer into the lower stratosphere. Tropospheric air predominantly
enters the stratosphere in the tropics, thereby largely setting stratospheric
composition (Fueglistaler et al., 2009). There is an ongoing debate on the
importance of convective penetration in this transport, particularly with
respect to water vapour (Jensen et al., 2007). Numerical simulations of
convective penetration events have been performed using realistic and complex
meteorological models (Dauhut et al., 2015, 2018) containing many physical
processes but these are computationally expensive and challenging to
interpret. In this study, the working hypothesis is that the process of
convective penetration in this geophysical setting can be greatly simplified
by considering a simple fluid dynamical problem representative of the
environment, in which a region of strong stable stratification is penetrated
by a turbulent buoyant plume generated in a region with weak or zero
stratification. We neglect many complicating factors such as shear in the
strongly stratified layer, the thermodynamics of distinct phases of water, and
associated latent heating in favour of a passive tracer and a quiescent
environment in the lower layer.

This fluid dynamical problem is relevant to many other geophysical processes,
including deep convection in the open ocean. Typically, mixing between the
deep ocean and near-surface water is hindered by the strong vertical density
gradients of the thermocline. In some regions, including several locations in
high latitude oceans and the Mediterranean Sea , intense buoyancy loss from
the ocean surface to the atmosphere results in strong, deep-reaching
convection (Marshall & Schott, 1999; Herrmann et al., 2008). The transport of
surface water into the deep ocean sets and maintains the properties of the
abyss (Marshall & Schott, 1999), both in terms of the general circulation and
also biogeochemical cycles (Ulses et al., 2021). Further examples of naturally
occurring flows involving penetrative convection include modification of
downslope oceanic gravity currents by near-surface convection (Doda et al.,
2023), smaller-scale atmospheric convection below an inversion (Kurbatskii,
2001), as well as the internal structure of many stars where a convectively
unstable layer is bounded above and below by stable layers (Singh et al.,
1994; Masada et al., 2013). Moreover, the fluid dynamical problem of
convective penetration itself is of continuing scientific interest, in
particular the generation of gravity waves by the penetrating plume cap
(Ansong & Sutherland, 2010) and the energetics of the system (Chen &
Bhaganagar, 2023), with applications to many problems (see Hunt & Burridge
(2015) for a discussion of open questions on fountains, i.e. plumes surrounded
by a more buoyant environment). It has only recently become possible to
numerically simulate penetrative convection at mid to high resolution, owing
to constraints such as the domain size needed to control edge effects, among
others.

Here, we are interested in diagnosing the tracer transport that results from
turbulent mixing between plume fluid carrying a passive tracer and the
surrounding environmental fluid in the stratified layer. For a plume entering
a stratified environment, the buoyancy of a fluid parcel will determine its
position once any transient dynamics have relaxed and the fluid settles and
restratifies. In the setup considered here, all tracer-carrying fluid parcels
which penetrate the stratified layer will remain in it, as their buoyancy is
non-zero. In a more complex environmental profile with a lower weakly
stratified layer and an upper strongly stratified layer, there is some
buoyancy threshold above which fluid parcels – carrying tracer – will remain
in the strongly stratified layer. The distribution of tracer in buoyancy space
is therefore the focus of this problem.

Plumb (2007) introduced a tracer-tracer volume distribution – normalised to
form a probability density function (PDF) – to study rapid isentropic mixing
in the stratosphere. The PDF is used to quantitatively represent tracer-tracer
data as a means of analysing model results. They note that this approach is
often unhelpful for analysing observations (or laboratory experiments) since
in-situ observations are often too sparse or do not cover a large enough
spatial domain. However, the method is well suited to numerical simulations
where there is full access to the three-dimensional buoyancy and tracer
fields. Penney et al. (2020) utilise the method to study diapycnal mixing of
passive tracers in a stratified shear flow. Owing to constraints on the
redistribution of fluid in the volume distribution, in the flow they consider
the distribution tends to compact over time. The changing shape of the
distribution may then be used to interpret the mixing process. Here we use a
similar formulation without normalisation, since the volume of regions in
buoyancy-tracer space are of interest in examining the evolution of the plume
and the eventual settling height of tracers carried by the plume.

In this paper we use the buoyancy-tracer volume distribution to diagnose the
transport of a passive tracer in large eddy simulations of a buoyant plume
penetrating a linearly stably stratified layer. The numerical method is
detailed in section 2 and simulation results are presented in section 3. In
section 4 we apply the buoyancy-tracer volume distribution to the stratified
layer in this simple convective penetration setup, in which a source (the
penetrating plume) inputs volume into a confined region of buoyancy-tracer
space. The pre-existing environmental fluid in the stratified layer lies in a
separate confined region of buoyancy-tracer space which can be considered a
‘reservoir’ of environmental fluid that gradually loses volume via mixing with
plume fluid that enters the stratified layer. We show that the flow can be
partitioned into three regions of buoyancy-tracer space: the source region
where volume enters, a transport region through which volume flows during
initial mixing between newly-arrived fluid and the pre-existing reservoir, and
an accumulation region where mixed fluid settles and later homogenises. Each
of these regions of buoyancy-tracer space correspond to coherent regions of
physical space, namely undiluted plume fluid in the rising interior of the
plume, the plume cap, and the radially spreading intrusion, respectively. In
section 5 we analyse diagnostics of the mixing process in each of these
regions.

##  2 Governing equations and numerical model

We consider the penetration of a buoyant plume with source radius $r_{0}$ and
source integral buoyancy flux $F_{0}$ generated in a uniform layer of depth
$H$ into a stably stratified layer with buoyancy frequency $N$. The problem
setup is shown in figure 1. To aid in the examination of the flow evolution
and mixing, we include a passive tracer $\phi$ that satifisfies the same
dynamical equation as buoyancy $b=-g\rho^{\prime}/\rho_{0}$ where
$\rho^{\prime}$ is the density deviation from a reference value $\rho_{0}$.
The tracer is passive in the sense that it has no coupling with the momentum
equation so does not influence the dynamics. The scalar field
$\phi(\boldsymbol{x},t)$ represents the tracer concentration in arbitrary
units, which diffuses according to an associated diffusivity $\kappa_{\phi}$
which for convenience is the same as the diffusivity of buoyancy $\kappa$
here, but may differ in general.

Figure 1: Simulation setup and key quantities & regions for convective
penetration of a buoyant plume into a stratified layer.

We non-dimensionalise using the source integral buoyancy flux $F_{0}$ (with
dimensions $L^{4}T^{-3}$) and buoyancy frequency $N$ in the stable layer. The
length scale is $L=F_{0}^{1/4}N^{-3/4}$ and the time scale is $T=N^{-1}$. We
assume the velocity scale is $L/T$. Owing to the large range of length scales
involved in convective penetration, resolving turbulent scales with direct
numerical simulation would incur a great computational cost. We therefore use
large eddy simulation (LES) with the anisotropic minimum dissipation (AMD)
eddy-viscosity model to represent unresolved scales (Vreugdenhil & Taylor,
2018). The non-dimensional governing equations for velocity $\boldsymbol{u}$
and scalars $b,\phi$ including sub-grid-scale (SGS) contributions are

$\displaystyle\frac{\partial\widehat{u}_{i}}{\partial x_{i}}$
$\displaystyle=0,$ (1)
$\displaystyle\frac{\mathrm{D}\widehat{u}_{i}}{\mathrm{D}t}+\frac{\partial\widehat{p}}{\partial
x_{i}}$
$\displaystyle=\frac{1}{\Rey}\frac{\partial^{2}\widehat{u}_{i}}{\partial
x_{j}\partial x_{j}}+\widehat{b}\delta_{i3}-\frac{\partial\tau_{ij}}{\partial
x_{j}},$ (2) $\displaystyle\frac{\mathrm{D}\widehat{b}}{\mathrm{D}t}$
$\displaystyle=\frac{1}{\Rey\Pr}\frac{\partial^{2}\widehat{b}}{\partial
x_{j}\partial x_{j}}-\frac{\partial\lambda^{(b)}_{j}}{\partial x_{j}}+f_{b},$
(3) $\displaystyle\frac{\mathrm{D}\widehat{\phi}}{\mathrm{D}t}$
$\displaystyle=\frac{1}{\Rey\Pr}\frac{\partial^{2}\widehat{\phi}}{\partial
x_{j}\partial x_{j}}-\frac{\partial\lambda^{(\phi)}_{j}}{\partial
x_{j}}+f_{\phi},$ (4)

where $\widehat{\cdot}$ indicates filtering at the resolved grid scale,
$\delta_{ij}$ is the delta function, $\nu$ is the kinematic viscosity, and
$\kappa$ is the scalar diffusivity of buoyancy $b$ and tracer concentration
$\phi$. The terms $f^{(b)}$ and $f^{(\phi)}$ represent the forcing applied to
the buoyancy and vertical velocity to generate the buoyant plume and are
discussed later. The SGS stress is
$\tau_{ij}=\widehat{u_{i}u_{j}}-\widehat{u}_{i}\widehat{u}_{j}$, the SGS
buoyancy flux is
$\lambda^{(b)}_{j}=\widehat{u_{j}b}-\widehat{u}_{j}\widehat{b}$ and similarly
for the SGS tracer flux $\lambda^{(\phi)}_{j}$. The two dimensionless
parameters are the Reynolds number and Prandtl number,

$\Rey=\frac{F_{0}^{1/2}}{\nu N^{1/2}},\quad\Pr=\frac{\nu}{\kappa},$ (5)

respectively. Note that the molecular diffusivities for $b$ and $\phi$ are
identical here; in general, the Prandtl number for the tracer may differ. The
eddy-viscosity model for the deviatoric component of the SGS stress and the
SGS buoyancy and tracer flux are

$\tau_{ij}^{d}=\tau_{ij}-\frac{1}{3}\delta_{ij}\tau_{kk}=-2\nu_{\text{SGS}}\widehat{S}_{ij},\hskip
20.00003pt\lambda^{(b)}_{j}=-\kappa^{(b)}_{\text{SGS}}\frac{\partial\widehat{b}}{\partial
x_{j}},\hskip
20.00003pt\lambda^{(\phi)}_{j}=-\kappa^{(\phi)}_{\text{SGS}}\frac{\partial\widehat{\phi}}{\partial
x_{j}},$ (6)

where $\nu_{\text{SGS}},\kappa^{(b)}_{\text{SGS}}$ and
$\kappa^{(\phi)}_{\text{SGS}}$ are the non-dimensional SGS viscosity, SGS
buoyancy diffusivity and SGS tracer diffusivity respectively, each determined
by the AMD scheme, and $\widehat{S}_{ij}$ is the non-dimensional shear-rate
tensor. The turbulent diffusivity and viscosity may locally exceed the
prescribed value by several orders of magnitude in regions with intense
turbulence.

We use the codebase DIABLO (Taylor, 2008) to perform three-dimensional LES of
the idealised setup shown in figure 1. DIABLO evolves the Boussinesq Navier-
Stokes equations (1) – (4) discretised using Fourier modes in the two periodic
horizontal directions and second-order finite differences in the vertical
direction with closed top and bottom boundaries. A third-order Runge-Kutta
scheme is used for time-stepping. A $2/3$ dealising rule is applied when
transforming from Fourier to physical space. We use a cubic domain of side
length $L_{\text{domain}}$ with a uniform grid of $512^{2}\times 513$ points.
The side length is chosen large enough that edge effects are not present and
the radially spreading intrusion that forms does not reach the boundary during
the simulation. A sponge layer is added in the top 20% of the domain, where
the velocity is damped towards zero and the buoyancy is damped towards the
initial background stratification $b^{(\text{env})}_{0}(z)=z$ (non-
dimensional), to inhibit the reflection of internal gravity waves from the top
boundary.

A buoyant plume is generated in a small forcing region at the bottom of the
domain (see figure 1) with the plume centreline at the middle of the
compuational domain, $x=y=0$. For convenience we define the bottom of the
stratified layer to be $z=0$ and define $t=0$ as the instant at which the
plume penetrates into the stratified layer. We use a volumetric forcing method
in which the vertical velocity, buoyancy and tracer concentration are relaxed
towards a prescribed profile in a thin layer at the bottom of the domain. The
profile is chosen to be the analytic solution to the Morton et al. (1956)
governing equations for a plume with source radius $r_{0}$ and integral source
buoyancy flux $F_{0}=\int_{0}^{\infty}\left.wb\right|_{z=0}\,\mathrm{d}r$,
under the assumption of a Gaussian radial profile. The relaxation is applied
across the entire layer, but the radial profile localises the forcing on $w$
and $b$ to near the centreline. We do not utilise the typical method of
generating a buoyant plume with a simple buoyancy gradient on the bottom
boundary (e.g. van Reeuwijk et al. (2016); Pham et al. (2007)) owing to
pinching of the plume radius close to the bottom boundary where inflow
dominates the diffusive boundary buoyancy flux. Pinching reduces control of
the source radius and results in numerical artefacts due to the horizontal
pseudo-spectral method.

A random 10% perturbation is applied to the radial profile in the forcing
region and to each velocity component in the two grid layers above the forcing
region to initiate turbulence, which develops through the depth of the
unstratified layer. The forcing region is set by two parameters $L_{c}$ and
$L_{p}$ which control a $\tanh$ function modulating the prescribed profiles.
The parameter $L_{c}$ controls the approximate depth of the forcing region and
$L_{p}$ determines how sharply the forcing decays above $z=L_{c}$. To
illustrate, given vertical profiles $r_{m}(z),b_{m}(z)$ consistent with the
desired plume, the forcing $f_{b},f_{\phi}$ (see (3) and (4)) applied to each
scalar $\theta=b,\phi$ is

$\displaystyle f_{\theta}(\boldsymbol{x},t)$
$\displaystyle=\frac{1}{\tau}\left(\theta(\boldsymbol{x},t)-2b_{m}(z)(1+\frac{1}{10}\xi(t))\exp\left[-\frac{1}{2}\frac{x^{2}+y^{2}}{r_{m}(z)^{2}}\right]\right)\left(\frac{1}{2}-\frac{1}{2}\tanh\left(\frac{z-L_{c}}{L_{p}}\right)\right),$
(7)

where $\xi(t)$ is a random number between $-1$ and $1$, used to apply the 10%
perturbation to the analytic profile at every time step. The forcing is the
product of the difference between the dynamical scalar field and the
prescribed vertical profile with a Gaussian radial profile, and a factor
limiting the forcing to $0\leq z\lessapprox L_{c}+L_{p}$. The factor $1/\tau$
controls the coupling strength with the momentum equations. The size of $\tau$
is arbitrary other than being small enough to control against dynamical
variation and large enough to prevent numerical instability.

$\Rey$ | $\Pr$ | $r_{0}$ | $H$ | $L_{\text{domain}}$ | $L_{c}$ | $L_{p}$ | $\tau$
---|---|---|---|---|---|---|---
2.8 $\times 10^{7}$ | 0.7 | 0.12 | 12 | 36 | 1.2 | 0.6 | 1.0
Table 1: Non-dimensional parameters for the large eddy simulation.

In §3–6 we will refer to a simulation with a single set of parameters which
are equivalent to the experimental setup used by Ansong & Sutherland (2010)
except for the source buoyancy flux, which is weaker here. The parameters are
given in table 1 and non-dimensionalised by $N=1\,\mathrm{s}^{-1}$ and
$F_{0}=8\times 10^{-8}\,\mathrm{m}^{4}\mathrm{s}^{-3}$. Henceforth, all values
stated are non-dimensionalised with respect to this choice of $F_{0}$ and $N$.

##  3 Simulation results

###  3.1 Flow structure

Figure 2: Snapshots of the tracer concentration in the $x-z$ plane at the
midpoint of the computational domain where $\phi\geq\phi_{\min}$. Where
$\phi<\phi_{\min}$, the buoyancy field is shown. Snapshots at non-dimensional
times $t=1,5,15$.

Figure 2 shows vertical cross-sections through the centreline of one such
simulation. The plume is identified as regions where the tracer concentration
$\phi\geq\phi_{\min}\equiv 5\times 10^{-4}$ to exclude numerical artefacts,
and where $\phi<\phi_{\min}$ the buoyancy field is shown and referred to as
the environmental buoyancy field $b^{(\text{env})}$. The bottom of the
stratified layer, where $b^{(\text{env})}$ becomes non-zero, is indicated by
the lowest buoyancy contour. Each panel represents a stage in the flow
evolution which proceeds as follows. Figure 2(a) shows initial penetration of
the plume cap. As plume fluid rises through the stratified layer, the upward
acceleration of the plume decreases as the relative buoyancy between the plume
and environment decreases, and once the environmental buoyancy exceeds that of
the plume, the plume decelerates. The kinetic energy of the plume is converted
to potential energy and eventually the rising fluid reverses direction, or
overturns, and begins to subside. At this initial penetration stage, there is
no pre-existing subsiding fluid to inhibit the rising fluid, so the plume
reaches its initial maximum penetration height $z_{\max}$, as shown in figure
2(b). Note that here we define $z_{\max}$ as the distance above the initial
bottom of the stratified layer (see figure 1) rather than the bottom of the
domain. At later times after the initial plume cap has reached $z_{\max}$, the
subsiding fluid inhibits the rise of new plume fluid, resulting in a
(typically) lower quasi-steady state maximum height $z_{ss}$. As fluid
subsides, its buoyancy relative to the surrounding environmental fluid
increases (i.e. becomes less negative) until reaching the level of neutral
buoyancy $z_{n}$. As shown in figure 2(c), fluid that has overturned then
spreads radially as an intrusion. The dynamics as described match
qualitatively with studies of similar set-ups in the literature, for example
Ansong & Sutherland (2010) with an identical setup and similar physical
parameters; Hunt & Burridge (2015) who provide a broad review of the
literature on fountains (i.e. a plume in a uniformly stratified domain); and
Devenish et al. (2010) in which numerical simulations of a fountain are
compared with various theoretical models.

Figure 3: Timeseries of the tracer concentration in the $z-t$ plane at the
centreline of the computational domain. The red dashed contour denotes the
threshold $\phi\geq 10^{-5}$. The maximum penetration height $z_{\max}$ and
the quasi-steady state height $z_{ss}$ are marked.

Figure 3 shows the tracer concentration at the centreline of the computational
domain over time, which visualises the evolution of the maximum height of the
plume cap during the penetration process and the subsequent quasi-steady
state. As has been noted in the literature, the maximum height of the plume
tends to oscillate around $z_{ss}$ (Turner, 1966) but, to our knowledge, the
mechanism setting the frequency of this oscillation is not well understood
(e.g. Ansong & Sutherland (2010)). In the simulation considered here, the
stratification is strong relative to the plume forcing so the quasi-steady
state height $z_{ss}$ is close to the maximum penetration height $z_{\max}$
and the oscillation is weak, as demonstrated in figure 3. For convenience, we
will use $z_{\max}$ to refer to the height of the plume cap.

Internal gravity waves across a range of frequencies are generated during the
penetration process. These waves are visible as small amplitude undulations in
the buoyancy contours in figure 2. The relevance of gravity waves for mixing
and the geophysical problem as a whole depends on context. For example,
gravity waves generated by penetrative convection are considered as important
as turbulence-induced waves for internal angular momentum transport in stars
(Pinçon et al., 2016). Similarly, gravity waves from convection penetrating
from the turbulent surface boundary layer into a weak atmospheric inversion
may carry away a significant fraction of the energy of the convective
overshoot, but a much smaller fraction for a strong inversion (Stull, 1976).
Convectively generated gravity waves from tropospheric thunderstorms modify
the surrounding environment and may promote the organisation of new convection
within the troposphere (Lane & Reeder, 2001), whilst gravity waves in the
stratosphere generated by thunderstorm outflows propagating along the
tropopause are not considered important in the mixing of convective overshoots
owing to the relatively small vertical momentum flux (Flynn & Sutherland,
2004). The influence of gravity waves on mixing in our idealised problem is
beyond the scope of this paper.

It is readily seen in figure 2 that penetration of the plume into the
stratified layer modifies the local environmental buoyancy gradient. In the
initial penetration stage, the bottom of the stratified layer is pushed
upwards, lifting buoyancy contours. Subsequently, subsiding fluid pushes the
buoyancy contours back below the initial position of the stratified layer
interface, where it remains in quasi-steady state. Moreover, the environmental
buoyancy contours surrounding the overturning plume cap oscillate in response
to the rising and falling motion of the plume cap itself. The deformation of
contours in this region locally weakens the buoyancy gradient. The
environmental buoyancy of the initial stratification at the maximum
penetration height $b^{(\text{env})}_{0}(z_{\max})$ represents a plausible
constraint on the maximum environmental buoyancy accessible for mixing with
the plume, and it determines the maximum level of plume fluid detrainment into
the environment (Ansong et al., 2008). However, the details of the
distribution of tracer concentration in buoyancy space depends on the range of
heights at which the most effective mixing between the plume and environment
occurs.

Figure 4: Normalised tracer distribution in buoyancy coordinates in the
stratified layer. Coloured lines are tracer distributions post-penetration at
equal time intervals. The black dashed line is the time-averaged pre-
penetration distribution from the source region indicated in figure 1.
Evolution of total tracer concentration (arbitrary units) in the stratified
layer inset.

###  3.2 Tracer distribution

Figure 4 shows PDFs of tracer concentration $\phi$ in buoyancy coordinates in
the stratified region at fixed time intervals post-penetration. Note that the
total tracer concentration in the stratified region increases with time – see
figure 4 inset – and the normalisation to form a PDF allows comparison of the
distributions at early and late times. A ‘source’ tracer PDF is computed from
a thin layer immediately beneath the stratified layer, referred to as the
source region – see figure 1. The source PDF represents the tracer vs.
buoyancy distribution in the plume _before_ penetrating the stratified region.
Without mixing, the source PDF would persist and tracer would be retained
within a region determined by the buoyancy range of the source PDF. Mixing
during the penetration process manifests as changes in the tracer PDF when
compared with the source PDF. After penetration, but before reaching
$z_{\max}$, the tracer distribution is ‘squeezed’ towards intermediate
buoyancies; the PDF decreases at low buoyancies $b\lesssim 0.3$ and large
buoyancies $b\gtrsim 0.7$ and results in a more concentrated peak in
$0.3\lesssim b\lesssim 0.7$. This occurs due to two simultaneous processes.
Until the plume first overturns at $z_{\max}$, only the edges of the plume –
corresponding to low buoyancy and tracer concentration – are exposed to the
surrounding environment. Fluid at the edge of the plume with low buoyancy and
tracer concentration mixes with more buoyant environmental fluid. As a result,
tracer is moved from low to intermediate buoyancy. Simultaneously, as the
plume rises, turbulent mixing in the undiluted interior of the plume acts to
homogenise the buoyancy and tracer distribution of the plume. This mixes the
large values of buoyancy and tracer concentration on the plume centreline with
smaller values further from the centreline, thus transporting tracer to lower
buoyancies. Once the plume has first overturned at $t\approx 2$, the large
tracer concentrations in the core of the plume are exposed to environmental
fluid of much greater buoyancy. Hence, as this plume fluid mixes with the
environmental fluid, tracer is carried to greater buoyancies, giving the
tracer distribution a growing positive tail as observed in figure 4 from
$t\approx 2$ onwards.

Figure 5: Tracer distribution in buoyancy (left) and height (right)
coordinates in the stratified layer evolving over time post-penetration.

Figure 5(a) shows the tracer concentration (with respect to buoyancy) as a
function of buoyancy and time and figure 5(b) shows the tracer concentration
(with respect to height) as a function of height and time, in which the $z$
and $b$ axes are linked by the initial ambient stratification strength,
$N^{2}$. The height $z$ shown in the right panel then corresponds to the
buoyancy of the initial stratification $b^{(\text{env})}_{0}(z)$ in the left
panel. At late times, most tracer lies in the spreading intrusion at the
neutral buoyancy height $z_{n}$, which coincides with the peak in the tracer
PDF in buoyancy coordinates. Note that the ambient stratification, $N^{2}$,
defines the background buoyancy profile used to link buoyancy and height, but
the ambient buoyancy gradient surrounding the plume tends to be weakened by
the penetrating plume (note the widening buoyancy contours at the bottom of
the stratified layer in figure 2), so tracer can reach much larger heights
than suggested by the corresponding environmental value of buoyancy. This
effect is supplementary to the rise of fluid parcels well above their level of
neutral buoyancy owing to excess kinetic energy. The weakened buoyancy
gradient means that penetration deep into the stratified layer does not
necessarily result in mixing with environmental fluid of significantly larger
buoyancy. Furthermore, undiluted plume fluid parcels which reach
$z_{\max}\approx 3$ can mix with environmental fluid parcels at the
corresponding buoyancy $\left.b\right|_{z_{\max}}\approx 3$ but the resulting
mixture will have a buoyancy that is a volume-weighted averaged of the
buoyancies of the initial fluid parcels. Since the plume buoyancy is
necessarily lower than the environmental buoyancy at $z_{\max}$, significant
tracer concentrations cannot reach $\left.b\right|_{z_{\max}}$.

The left panel of figure 5 includes a red contour $\phi(b)=10$ which indicates
the extent of the tracer distribution in buoyancy coordinates. Similarly, the
right panel includes a contour $\phi(z)=1$ indicating the extent in height
coordinates. Comparing the evolution of the tracer distribution in buoyancy
and height coordinates, we note that whilst the plume reaches its maximum
penetration height at $t\approx 2$, the greatest tracer distribution extent in
buoyancy space is not reached until $t\approx 6$. If mixing between plume and
environmental fluid during overturning was fast compared with the dynamical
timescale then we would expect tracer to rapidly appear at buoyancies
corresponding with $z_{\max}$. This suggests that the dynamical timescale of
the flow is shorter than the timescale of mixing between plume and
environmental fluid, such that mixing between the largest tracer
concentrations first exposed during overturning and the environment continues
during subsidence, where environmental buoyancy decreases. As a result, we
find that most tracer is transported to intermediate buoyancies between
$\left.b\right|_{z_{n}}$ and $\left.b\right|_{z_{\max}}$, with only very small
amounts of tracer reaching the greatest buoyancies closer to
$\left.b\right|_{z_{\max}}$. Finally, note that the most significant tracer
transport occurs when large tracer concentrations are supplied by the plume
core. This is observed in figure 5, with the most significant changes in the
tracer distribution in buoyancy space – indicated by changes in the position
of the red contour – occurring when when the largest eddies, carrying the most
tracer, reach their maximum height, for example at $t\approx 4$ and $t\approx
6$.

##  4 Diagnosing tracer transport

The turbulent mixing of fluid parcels can be considered a two-step process
(e.g. Davies Wykes & Dalziel (2014)), composed of stirring and molecular
diffusion. Stirring geometrically deforms fluid parcels and stretches surfaces
of constant buoyancy, increasing the area over which molecular diffusion acts.
Moreover, stirring strengthens scalar gradients across buoyancy surfaces,
further promoting diffusion. Stirring is – in principle – a reversible
process, whilst molecular diffusion is irreversible. Therefore, to focus on
_irreversible_ tracer transport, we wish to separate the effects of stirring
and diffusion. Since stirring only deforms fluid parcels, it does not change
the tracer concentration or buoyancy within fluid parcels and therefore only
the diffusive component of mixing results in changes in the distribution of
tracer in buoyancy space.

A limitation of 1-dimensional tracer distributions in buoyancy space is the
lack of information on the volume of fluid parcels at a given buoyancy; a peak
in the tracer PDF may represent relatively few fluid parcels carrying large
tracer concentrations or relatively many fluid parcels carrying small amounts
of tracer. The distinction is important, since the former can result in
stronger gradients upon which diffusion acts and therefore more effective
transport of tracer. Consequently, we move from a 1D tracer PDF in buoyancy
space to considering a 2D volume distribution in buoyancy-tracer space.

The volume distribution formulation presented here builds on previous forms
presented by Plumb (2007) and Penney et al. (2020). In particular, we move
from a probability density formulation to a volume distribution. In this
problem, the total physical volume represented by non-zero parts of the
distribution (i.e. the total plume volume in the stratified layer) changes in
time, so normalising by this volume introduces unnecessary complexity. It is
feasible instead to normalise by the fixed volume of the stratified layer, but
this is arbitrary: we are uninterested in regions with no tracer, and the
normalisation could be changed by increasing the size of the stratified layer
whilst the plume is unchanged. In that case the value of the volume
distribution changes whilst its interpretation does not. Omitting
normalisation also simplfiies the numerical implementation and its
interpretation.

###  4.1 Buoyancy-tracer volume distribution

####  4.1.1 Definition

We consider a volume distribution $W(b^{\prime},\phi^{\prime};t)$ in
$(b,\phi)$-space defined in a fixed volume $V$ as

$W(b^{\prime},\phi^{\prime};t)=\int_{V}\delta(b(\boldsymbol{x},t)-b^{\prime})\delta(\phi(\boldsymbol{x},t)-\phi^{\prime})\,\mathrm{d}V,$
(8)

where $\delta(\cdot)$ is the Dirac delta function with the inverse dimension
of its argument. The definition is such that the volume of fluid in $V$ with
$b^{\prime}<b(\boldsymbol{x},t)<b^{\prime}+\mathrm{d}b^{\prime}$ and
$\phi^{\prime}<\phi(\boldsymbol{x},t)<\phi^{\prime}+\mathrm{d}\phi^{\prime}$
is given by
$W(b^{\prime},\phi^{\prime};t)\mathrm{d}b^{\prime}\mathrm{d}\phi^{\prime}$. An
evolution equation for the volume distribution $W$ can be obtained using the
advection-diffusion equations for $b$ and $\phi$ and is derived in appendix B.
We have

$\frac{\mathrm{d}W}{\mathrm{d}t}=-\bnabla_{(b^{\prime},\phi^{\prime})}\cdot\boldsymbol{F}-S,$
(9)

where $\boldsymbol{F}$ is a volume-weighted average of the diffusive forcing
in (3) and (4), representing the flux of $W$ in $(b,\phi)$-space:

$\boldsymbol{F}(b^{\prime},\phi^{\prime};t)=\int_{V}\kappa(\nabla^{2}b,\nabla^{2}\phi)\,\delta(b-b^{\prime})\delta(\phi-\phi^{\prime})\,\mathrm{d}V.$
(10)

The final term in (9) represents a source or sink of $W$,

$S(b^{\prime},\phi^{\prime};t)=\int_{\partial
V}\boldsymbol{u}\cdot\boldsymbol{n}\,\delta(b-b^{\prime})\delta(\phi-\phi^{\prime})\,\mathrm{d}A,$
(11)

where $\boldsymbol{u}$ is the velocity of the flow in physical space and
$\boldsymbol{n}$ is the outward normal on the boundary $\partial V$ of $V$.
Since we are considering a flow _into_ $V$,
$\boldsymbol{u}\cdot\boldsymbol{n}$ is negative and $S$ acts as a source of
$W$, i.e. $S<0$ is a source. Note that (9) contains no terms in which
advection plays an explicit role except for the source term – which captures
advection through the domain boundary – representing the fact that $W$ remains
unchanged under advection within the domain.

As plume fluid enters the domain $V$ carrying tracer, this is reflected in the
volume distribution $W$ via the source term $S$. As turbulent mixing occurs,
diffusion redistributes $W$ via the flux term
$-\nabla_{(b^{\prime},\phi^{\prime})}\cdot\boldsymbol{F}$. The change in the
volume distribution at $(b^{\prime},\phi^{\prime})$ as a result of turbulent
mixing up to time $t$ is therefore

$M(b^{\prime},\phi^{\prime};t)=-\int_{0}^{t}\nabla_{(b^{\prime},\phi^{\prime})}\cdot\boldsymbol{F}\,\mathrm{d}t^{\prime},$
(12)

such that
$M(b^{\prime},\phi^{\prime};t)\mathrm{d}b^{\prime}\mathrm{d}\phi^{\prime}$ is
the cumulative volume of fluid in V with
$b^{\prime}<b(\boldsymbol{x},t)<b^{\prime}+\mathrm{d}b^{\prime}$ and
$\phi^{\prime}<\phi(\boldsymbol{x},t)<\phi^{\prime}+\mathrm{d}\phi^{\prime}$
that has undergone diffusive mixing up to time $t$. Therefore $M$ represents
the integrated effect of the diffusive flux and we refer to $M$ as the
cumulative mixed volume.

####  4.1.2 Properties

Turbulent stirring acts to bring fluid parcels with different buoyancy and
tracer concentration together in physical space without modifying the volume
distribution since the distribution is invariant under advection within the
domain, as shown above and discussed further in Penney et al. (2020) and Plumb
(2007). This stirring locally intensifies buoyancy and tracer gradients.
Diffusive mixing then acts homogenise the buoyancy and tracer concentration in
a small physical region. Provided the molecular diffusivities of buoyancy and
tracer are equal, a mixture of two fluid parcels lies on a line between the
two parcels in buoyancy-tracer space. The buoyancy-tracer characteristics of
the resulting mixture are a volume-weighted average of the characteristics of
the initial fluid parcels. Therefore, the evolution of the buoyancy-tracer
volume distribution under turbulent mixing is constrained to lie within the
convex envelope of the initial volume distribution in buoyancy-tracer space,
i.e. the smallest envelope that contains all non-zero points of the
distribution. Fluid that enters the domain may join the mixing process, so as
the volume distribution continuously evolves, we conclude that at any given
time, the evolution of the volume distribution lies within the convex envelope
of the current distribution combined with the distribution of any new fluid
that enters the stratified layer.

Figure 6: Schematic diagram of $(b,\phi)$-space illustrating the convex
envelope property of the volume distribution.

A consequence of this constraint on the volume distribution of fluid within
the stratified layer is the tendency to compact over time as turbulent mixing
continuously shifts the buoyancy-tracer characteristics of fluid parcels
towards the mean. Moreover, this implies that compact distributions in
buoyancy-tracer space represent well-mixed fluid. Fluid that enters the domain
may cause the extremes of the distribution to persist, but as mixing with the
environment proceeds, volume will eventually accumulate near the mean of the
entire distribution. We have already seen evidence of this in §3.2 and in
particular figure 4 with large eddies carrying extreme tracer concentrations
eroded via turbulent mixing in the rising plume, resulting in tracer transport
towards the plume core average tracer concentration. Furthermore, note that
the the convex envelope of a line is also a line. Hence when fluid parcels
along a line in buoyancy-tracer space mix, the resulting mixture will remain
on that line. As a result, fluid on a line in buoyancy-tracer space which has
not mixed with fluid parcels elsewhere in buoyancy-tracer space can be
distinguished from fluid which _has_ mixed, since mixed fluid would no longer
be confined to the line.

To illustrate these principles, consider figure 6 which shows a schematic of
the effects of turbulent mixing on a small set of fluid parcels represented in
the volume distribution. Of course, mixing occurs throughout the fluid and it
is not possible to isolate the distinct fluid parcels mixing at any one time.
Nonetheless, applying the principle illustrated here offers a way to link the
dynamical processes in the flow with changes in the buoyancy-tracer volume
distribution by identifying which parts of the flow in physical space
correspond to which regions in buoyancy-tracer space. In figure 6, the blue
dots represent fluid parcels in buoyancy-tracer space that initially lie in
the stratified layer and the black dotted line is the initial convex envelope.
The green dot represents a fluid parcel that enters the stratified layer
during the process of turbulent mixing and must therefore be included in the
convex envelope that the final mixture is constrained to lie within, indicated
by the black dashed line. In the process of turbulent mixing, stirring brings
the blue and green fluid parcels together and diffusive mixing modifies their
properties, leaving the final mixed fluid parcel as indicated by the red dot.
We emphasise that the red dot is contained within the convex envelope of
initial and newly-arriving fluid, but not necessarily within the convex
envelope of the initial fluid only. The grey arrows indicate the flux
$\boldsymbol{F}$ of $W$ associated with the individual mixing event
represented here. Note that we are indicating the direction of
$\boldsymbol{F}$ only; the magnitude would depend on the volume of the parcels
undergoing mixing.

####  4.1.3 Numerical implementation

To examine the numerical simulations considered here, we consider a discrete
formulation of the buoyancy-tracer volume distribution applied to the
stratified layer. Since the bottom of the stratified layer is perturbed below
the initial height $z=0$ post-penetration, we choose the domain $V$ to be
$-1\leq z\leq L_{z}$ to capture the entire stratified region at all times.
Note that in defining the volume distribution, we used a fixed volume $V$ in
physical space. However, by restricting the buoyancy and tracer domains, we
can restrict attention to a subvolume of $V$ which is not necessarily fixed,
but the budget equation (9) still applies. Here, we restrict to $b>0$ so that
regions of $V$ that are not part of the stratified layer are excluded and also
restrict to $\phi>\phi_{\min}$ to isolate the plume within $V$. As a
consequence, the ‘reservoir’ of environmental fluid where $\phi=0$ is
excluded. In interpreting results, we must therefore consider the line
$\phi=\phi_{\min}$ as a source where volume can enter the distribution from
the environment – see discussion on entrainment in §4.5.

We consider a discrete formulation of the buoyancy-tracer volume distribution
in which we subdivide the buoyancy and tracer domains into $N_{b}$ and
$N_{\phi}$ equally sized bins of size

$\delta
b=\frac{b_{\max}-b_{\min}}{N_{b}},\quad\delta\phi=\frac{\phi_{\max}-\phi_{\min}}{N_{\phi}}$
(13)

where $b_{\min}=0$ and $b_{\max}$ is chosen to be slightly larger than the
initial environmental buoyancy at the maximum penetration height
$b^{(\text{env})}_{0}(z_{\max})$ to maximise the resolution of the
distribution. To accomodate fluctuations in tracer concentration, we choose
$\phi_{\max}$ to be larger than the tracer concentration on the plume
centreline at penetration height $z=0$ as predicted by the MTT plume
equations, which depends on the source integral tracer flux $F_{0}$ (identical
to buoyancy). As discussed above, we set $\phi_{\min}=5\times 10^{-4}$ to
restrict attention to plume fluid and to be consistent with the tracer
concentration threshold used to define the plume in e.g. figure 2. Henceforth
we use $N_{b}=N_{\phi}=128$.

Denoting the centre of a given bin as $(b_{i},\phi_{j})$, the associated value
of the volume distribution is computed as

$W_{ij}(t)=\sum_{V}I_{ij}(\boldsymbol{x},t)\Delta x\Delta y\Delta z$ (14)

where the sum is over all points within the domain, $\Delta x,\Delta y,\Delta
z$ are the grid-cell widths, and the indicator $I_{ij}(\boldsymbol{x},t)$ is
defined as

$I_{ij}(\boldsymbol{x},t)=\begin{cases}1&\left(b(\boldsymbol{x},t)-b_{i},\phi(\boldsymbol{x},t)-\phi_{j}\right)\in\bigl{(}-\frac{1}{2}\delta
b,\frac{1}{2}\delta
b\bigr{]}\times\bigl{(}-\frac{1}{2}\delta\phi,\frac{1}{2}\delta\phi\bigr{]}\\\
0&\text{otherwise}\end{cases}$ (15)

The value of $W_{ij}(t)$ therefore represents the total volume within $V$
where the buoyancy lies within $\delta b/2$ of $b_{i}$ and the tracer
concentration lies within $\delta\phi/2$ of $\phi_{j}$. Note that the in the
continuous formulation, the volume distribution $W(b,\phi;t)$ defined by (8)
must be integrated over $b$ and $\phi$ to yield a volume, whilst $W_{ij}(t)$
itself has dimensions of volume and need only be summed over $i$ and $j$. The
definitions coincide in the limit as $\delta b,\delta\phi\to 0$, such that

$\lim_{\delta b,\delta\phi\to 0}\frac{W_{ij}(t)}{\delta
b\delta\phi}=W(b_{i},\phi_{j};t).$ (16)

The discrete flux term $\boldsymbol{F}_{ij}(t)$ is defined as

$\boldsymbol{F}_{ij}(t)=(F^{b}_{ij}(t),F^{\phi}_{ij}(t))=\sum_{V}I_{ij}(\boldsymbol{x},t)(\dot{b},\dot{\phi})\Delta
x\Delta y\Delta z$ (17)

where $\dot{b}$ and $\dot{\phi}$ are the diffusive forcing terms in the scalar
evolution equations (3), (4) of $b,\phi$ respectively.

The discrete form of the source term $S$ is

$S_{ij}(t)=\sum_{\partial
V}I_{ij}(\left.\boldsymbol{x}\right|_{z=-1},t)w(\left.\boldsymbol{x}\right|_{z=-1},t)\Delta
x\Delta y$ (18)

where $I_{ij}$ is the indicator function as defined in (15) except with
$\boldsymbol{x}=(x,y,-1)$ and $\partial V$ is the surface $z=-1$. The value
$S_{ij}(t)$ represents the total volume of fluid with buoyancy within $\delta
b/2$ of $b_{i}$ and tracer concentration within $\delta\phi/2$ of $\phi_{j}$
which has entered the stratifed layer up to time $t$.

Finally, the discrete form of the cumulative mixed volume $M$ is

$M_{ij}(t)=W_{ij}(t)+\sum_{t^{\prime}}S_{ij}(t^{\prime})\Delta t^{\prime}$
(19)

which represents the change in the total volume of fluid with buoyancy within
$\delta b/2$ of $b_{i}$ and tracer concentration within $\delta\phi/2$ of
$\phi_{j}$ due to mixing. Equivalently, $M_{ij}(t)$ quantifies the integrated
effect of the buoyancy-tracer flux in the given range of buoyancy and tracer
concentration up to time $t$.

###  4.2 Results

The buoyancy-tracer volume distribution $W_{ij}(t)$ is shown in figure 7
(middle) with $x-z$ cross-sections of the tracer field and environmental
buoyancy in the stratified layer (top) to illustrate the dynamics associated
with each snapshot of the volume distribution. The colourscale for $W_{ij}$
indicates the volume of fluid in the stratified layer with the corresponding
buoyancy and tracer concentration. Regions of the volume distribution which
are not coloured indicate that there is no fluid with buoyancy and tracer
concentration in that range. The convex envelope indicated by the red dashed
line is that of $W$ at time $t$ and $S$ at time $t+\Delta t$ where $\Delta t$
is the time interval between data collection in the numerical simulation. At
each time $t$, the volume distribution $W(t+\Delta t)$ is constrained to lie
within the convex envelope shown. The environmental buoyancy at the maximum
height at each time is indicated by the dotted red line. In order to reduce
the influence of numerical artefacts, we exclude points of the volume
distribution where $W_{ij}(t)<10^{-3}$. The same thresholding (based on
$W_{ij}(t)$) is applied to all distributions shown henceforth.

Figure 7: Three instantaneous snapshots of the buoyancy-tracer volume
distribution $W(b,\phi;t)$ (bottom) with $x-z$ cross-section of plume tracer
field $\phi$ in the stratified layer with surrounding environmental buoyancy
contours (top).

The source distribution $S_{ij}(t)$ is also shown on the bottom row of figure
7 and shows that the undiluted plume fluid entering the stratified layer has a
compact distribution in buoyancy-tracer space, i.e. $S$ is non-zero only in
restricted regions. The fluid arriving in the stratified layer lies on a line
since the azimuthally averaged buoyancy and tracer concentrations in a steady-
state plume have, to good approximation, Gaussian radial profiles with the
same width, irrespective of the forcing applied to generate the plume – see
figure 20 in appendix A. We refer to this line as the source or input line.
Departures from the mean cause some spread in the source line but the
resolution of the volume distribution is not always sharp enough to capture
the effect. The core of the plume, where $b$ and $\phi$ are largest,
corresponds to the extreme end of the source distribution, whilst the edges of
the plume where $b$ and $\phi$ are smallest correspond to the region of the
source line closest to the origin in buoyancy-tracer space. Regions where $S$
is non-zero away from the source line (when present) indicate parts of the
plume falling below $z=-1$. Mixing within the undiluted plume as it rises
through the stratified layer redistributes fluid on the source line but does
not move fluid away from the line, by the convex envelope property. In the
buoyancy-tracer volume distribution $W$ shown on the middle row of figure 7,
fluid appearing away from the line found in the source distribution therefore
represents a mixture of plume and environmental fluid. To aid in interpreting
the volume distribution, figure 8 represents this schematically.

Figure 8: Schematic diagram illustrating the correspondence between regions in
physical space and regions in buoyancy-tracer space, and the inaccessibility
of some regions of the convex envelope.

Figure 7(a) shows the plume shortly after penetrating the stratified layer and
before reaching its maximum penetration height. At this stage, only the edges
of the plume are exposed to the environment in the lower stratified layer
which has a relatively small buoyancy. We therefore find volume appearing in a
small region of the convex envelope closest to the origin in buoyancy-tracer
space. In figure 7(b), the plume has reached its maximum penetration height
and undiluted plume fluid in the plume core, which has until this stage been
shielded from the environment by the edges of the plume, overturns and becomes
exposed to environmental fluid near $z_{\max}$ which has relatively large
buoyancy. The volume distribution at this time is now far less compact (see
gray shaded region in figure 8), as environmental fluid with a large range of
buoyancy and a large range of the source line are now able to mix. Note that
we do not find mixing between the extreme end of the source line and the
environment (see hatched region in figure 8) – as undiluted plume fluid in the
core of the plume rises, turbulent mixing continues to erode the buoyancy and
tracer concentration within the source distribution, so that when this fluid
is first exposed to the environment it lies closer to the middle of the source
line.

Figure 7(c) shows the plume in quasi-steady state (QSS) where fluid that
overturns at heights close to $z_{\max}$ has subsided and formed a radially
spreading intrusion. The intrusion dominates the total volume of the plume at
late times and is represented by the region of large $W$ at intermediate
buoyancies $0.5\lesssim b\lesssim 1$ and tracer concentrations
$0.0025\lesssim\phi\lesssim 0.005$. In the QSS plume there is a continuous
input of undiluted plume fluid along a line $\phi\propto b$, mixing between
the plume and environment up to buoyancies corresponding to heights near
$z_{\max}$ and an accumulation of fluid volume at lower buoyancies and tracer
concentrations corresponding to the radially spreading intrusion. Note that
the intrusion contributes to the ‘shielding’ effect by the plume edge, which
prevents undiluted plume fluid in the core from accessing the environment.

Note that a ‘bulge’ in the volume distribution appears to the left of the
source line, which is the result of a numerical artefact: Gibbs ringing causes
small regions of negative $b$ to appear transiently in the core of the plume
near $z_{\max}$, which mix with other fluid in the core of the plume, moving
very small amounts of volume away from the source line. These points are
excluded from the convex envelopes shown and also excluded from calculations
discussed in the following subsection.

###  4.3 Source, transport and accumulation regions

We now restrict attention to the plume in quasi-steady state (discussed in
section 4.4) and consider the cumulative mixed volume distribution
$M(b,\phi;t)$ and the volume flux in $(b,\phi)$-space, $\boldsymbol{F}$. All
buoyancy-tracer distributions defined in section 4.1.1 can be treated as
fields in physical space as well as buoyancy-tracer space, for example
$M(\boldsymbol{x},t)=M(b(\boldsymbol{x},t),\phi(\boldsymbol{x},t);t)$. To
avoid confusion between regions of buoyancy-tracer space and the corresponding
fluid in physical space, we refer to the former as regions and the latter as
classes of fluid: for example, class U refers to undiluted plume fluid
corresponding to the source region $\mathcal{U}=\\{M<0\\}$ of $(b,\phi)$-space
where $M<0$. For convenience, we use the notation $\\{M<0\\}$ as shorthand for
$\\{(b,\phi)|M(b,\phi;t)<0\\}$.

Figure 9: Snapshots at time $t=15$ of the buoyancy-tracer cumulative mixed
volume distribution $M(b,\phi;t)$ in $(b,\phi)$-space (left) and in an $x-z$
cross-section of physical space (right) with surrounding environmental
buoyancy contours.

Figure 9 shows the cumulative mixed volume distribution in both physical and
buoyancy-tracer space at $t=14$ when the plume is in quasi-steady state. The
cumulative mixed volume distribution quantifies the change in volume due to
mixing. Where $M$ is large and positive, there is significant accumulation of
volume in that region of buoyancy-tracer space, whilst small positive values
of $M$ indicate an approximate balance between the volume of fluid parcels
arriving in that region of buoyancy-tracer space and the volume of fluid
parcels leaving. Where $M$ is negative, volume is being lost from the given
region of buoyancy-tracer space and mixed towards other regions.

Figure 10: Cumulative mixed volume distribution $M$ as in figure 9 with flux
vectors $\boldsymbol{F}$ overlaid (left) and flux divergence
$\bnabla\cdot\boldsymbol{F}$ with field lines of $\boldsymbol{F}$ overlaid
(right).

Figure 10 shows a quiver plot of the volume flux $\boldsymbol{F}$ overlaid
onto $M$ (left) and the volume flux divergence $\bnabla\cdot\boldsymbol{F}$
(right). We find that undiluted plume fluid entering the stratified layer
moves primarily along the ’source line’ where $M<0$ before moving into the
$M>0$ region. Environmental fluid is entrained into the plume via mixing and
accumulates where $M$ is maximised. The region of largest $M$ in
$(b,\phi)$-space corresponds in physical space to the intrusion. Since we
expect fluid in physical space to accumulate only in the intrusion, we can
distinguish between the region of $(b,\phi)$-space where $M$ is approximately
constant in time and therefore remains small once reaching quasi-steady state;
and where $M$ is increasing with time.

The regions are distinguished by a time-dependent threshold $m^{*}(t)$ such
that the ‘transport’ region, where much of the transition from undiluted to
mixed plume fluid occurs, is $\mathcal{T}=\\{0<M\leq m^{*}\\}$ and we refer to
fluid in this region as class T. The ‘accumulation’ region is
$\mathcal{A}=\\{M>m^{*}\\}$ with corresponding mixed fluid accumulating in the
intrusion referred to as class A. The numerical method for calculating
$m^{*}(t)$ is based on minimising with respect to $m$ the total flux
divergence,

$\int_{\\{0<M\leq m\\}}\frac{\partial W}{\partial
t}-(-S)\,\mathrm{d}b\mathrm{d}\phi=\int_{\\{0<M\leq
m\\}}\bnabla_{(b,\phi)}\cdot\boldsymbol{F}\,\mathrm{d}b\mathrm{d}\phi,$ (20)

so that the volume that enters the region $\mathcal{T}$ is approximately equal
to the volume leaving $\mathcal{T}$, i.e. the region captures transport of
volume alone. The remaining mixed fluid in class A must then capture the
accumulation of mixed fluid, formed from undiluted plume fluid that has
entered the stratified layer and entrained environmental fluid. The numerical
implementation of this method is discussed further in appendix C.

The cumulative mixed volume distribution with this partitioning is shown in
figure 11, with class U coloured blue, class T coloured green and class A
coloured red. In class U and A, the volume weighted centre-of-mass is shown,
approximately indicating the position in $(b,\phi)$-space towards which mixing
acts to move fluid within each class.

Figure 11: Cumulative mixed volume distribution $M$ at time $t=15$ partitioned
into three classes: U (blue) where $M\leq 0$, T (green) where $0<M\leq
m^{*}(t)$ and A (red) where $M>m^{*}(t)$. Triangles represent the $W$-weighted
centre-of-mass in the respective class.

###  4.4 Quasi-steady state

In this problem, quasi-steady state refers to the long-term behaviour once a
radial intrusion has formed, where undiluted plume fluid continuously arrives
in the stratified layer, mixes with the environment, and accumulates in the
intrusion. This means that in quasi-steady state the volume of undiluted plume
fluid in the stratified layer remains roughly constant whilst the total plume
volume increases monotonically.

The partitioning of $(b,\phi)$-space into regions $\mathcal{U},\mathcal{T}$
and $\mathcal{A}$ only applies in quasi-steady state, and not in the transient
penetration stage. By the convex envelope property, regions with $S=0$ and
$\phi>0$ are necessarily associated with a mixture of plume and environmental
fluid since fluid with these properties (i.e. the associated range of $b$ and
$\phi$) does not arrive in the undiluted plume. Since $W$ is positive by
definition in the plume, these mixtures of plume and environmental fluid have
$M>0$. Regions with $S>0$ are associated with undiluted plume fluid (class U)
on the source line. Here, we must have $M=W+\int S\,\mathrm{d}t<0$. Volume is
continuously lost from the source line by mixing within the undiluted plume
and eventually the environment, but volume is continuously replaced by new
undiluted plume fluid entering the stratified layer. Therefore on the source
line in quasi-steady state, the volume of fluid in the stratified layer with
the given $(b,\phi)$, $W(b,\phi)$, is approximately constant whilst the volume
of fluid that has arrived in the stratified layer with the given ($b,\phi)$,
$\int-S(b,\phi)\,\mathrm{d}t$, increases monotonically with time, so $M<0$
necessarily. Similarly, the possibility of $M=0$ is excluded. In quasi-steady
state we can therefore associate regions of $M<0$ with undiluted plume fluid
_only_ and $M>0$ with mixed fluid _only_. The $M>0$ region may then be further
decomposed to isolate fluid lying in the intrusion as above. Until quasi-
steady state is reached, parts of the source line may have $M>0$ so the
partitioning is not consistent.

Figure 12: Source region, class U and total plume volume used to identify
quasi-steady state in (a) and the decomposition of the plume volume into
undiluted plume fluid, plume cap, and intrusion in (b).

The volume associated with a region $\mathcal{R}$ of $(b,\phi)$-space at time
$t$ is

$V(\mathcal{R};t)=\int_{\mathcal{R}}W(b,\phi;t)\,\mathrm{d}b\mathrm{d}\phi=\sum_{\mathcal{R}}W_{ij}(t)$
(21)

To quantitatively identify QSS, first we define the source line more generally
as the region where the cumulative source volume is positive,
$\mathcal{S}=\\{(b,\phi)|\int_{0}^{t}S(b,\phi;t^{\prime})\,\mathrm{d}t^{\prime}>0\\}$.
We expect that in quasi-steady state
$\mathrm{d}V(\mathcal{S})/\mathrm{d}t\approx 0$ and that the partitioning
discussed above is consistent, in which case the region $\mathcal{S}$ and
$\mathcal{U}$ coincide so $V(\mathcal{S})\approx V(\mathcal{U})$. Note that
there may be some variation in $V(\mathcal{S})$ due to the plume ‘bobbing’, so
to aid in interpretation we also compare with the total plume volume which
grows monotonically. These volumes and the time we identify as the start of
quasi-steady state, $t_{QSS}$, are shown in figure 12(a). The QSS partition of
the total plume volume into the undiluted plume, plume cap, and intrusion
volume is shown in figure 12(b). As expected, we find that the volume of the
intrusion (class A) dominates the plume at late times, as the volume of
undiluted plume fluid and the plume cap remain approximately constant.

###  4.5 Entrainment

We calculate the entrained volume $E$ up to time $t$ as the difference between
the volume of the full plume and the cumulative volume of the source term $S$
up to time $t$, so that $E$ represents the volume of environmental fluid that
has been mixed into the plume. We have

$\displaystyle E$
$\displaystyle\equiv\int_{\\{W>0\\}}W(b,\phi;t)\,\mathrm{d}b\mathrm{d}\phi-\int_{0}^{t}\int_{\\{W>0\\}}-S(b,\phi;t^{\prime})\,\mathrm{d}b\mathrm{d}\phi\mathrm{d}t^{\prime}$
(22)
$\displaystyle=-\int_{0}^{t}\int_{\\{W>0\\}}\bnabla\cdot\boldsymbol{F}(b,\phi;t^{\prime})\,\mathrm{d}b\mathrm{d}\phi\mathrm{d}t^{\prime}$
(23)
$\displaystyle=\int_{0}^{t}\int_{\\{\phi=\phi_{\min}\\}}F_{\phi}(b,0;t^{\prime})\,\mathrm{d}b\mathrm{d}t^{\prime}$
(24)
$\displaystyle\approx\sum_{t^{\prime}}\sum_{i,j=0}F^{\phi}_{ij}(t^{\prime})\Delta
t$ (25)

where the second equality follows from time-integrating (9) and the final
equality follows from Green’s theorem and the fact that the volume flux
$\boldsymbol{F}$ vanishes on the boundary of the $\\{W>0\\}$ region except on
the surface $\phi=\phi_{\min}$ where environmental fluid enters the volume
distribution via entrainment. The final equality shows the numerical
approximation of $E$. Figure 13(a) shows the evolution of the entrainment
profile through the simulation, that is, the entrained volume flux $\dot{E}$
as a function of $b$, calculated by excluding the sum over $i$ in (25) and
then numerically differentiating in time. Figure 13(b) shows the volume
normalised entrainment profile, i.e. the entrained volume flux $\dot{E}$ as a
function of $b$, normalised by the volume of plume fluid with buoyancy within
$\delta b/2$ of $b_{i}$ for each index $i$. We refer to this as the
entrainment rate. We find that the entrained volume flux is largest at
intermediate buoyancies, indicating that most entrained volume is from the
environmental fluid surrounding the intrusion. The largest entrainment rate at
late times is found at large buoyancies. However, during the initial
penetration stage at early times, significantly larger entrainment rates are
found at the largest accessible buoyancy at the given time, indicating strong
entrainment occuring in the plume cap as it pushes further into the stratified
layer.

Figure 13: (a) Entrained volume flux profile and (b) entrainment rate profile.

For a region $\mathcal{R}$ of $(b,\phi)$-space which intersects $\\{\phi=0\\}$
we define the entrainment rate as $\dot{E}(\mathcal{R})/V(\mathcal{R})$, which
provides a clearer comparison of the entrainment between regions with
different physical volumes. The entrained volume can be decomposed into the
volume entrained into the plume cap $E(\mathcal{T})$, the volume entrained
into the intrusion $E(\mathcal{A})$, and the volume entrained just before
penetrating the stratified layer $E(\mathcal{U})$, by splitting the flux
integral (25) between the regions $\mathcal{U},\mathcal{T}$ and $\mathcal{A}$:

$\displaystyle E$
$\displaystyle=\int_{0}^{t}\int_{\\{\phi=0\\}}F_{\phi}(b,0;t^{\prime})\,\mathrm{d}b\mathrm{d}t^{\prime}$
(26)
$\displaystyle=\int_{0}^{t}\left[\int_{\\{\phi=\phi_{\min}\\}\cap\mathcal{U}}+\int_{\\{\phi=\phi_{\min}\\}\cap\mathcal{T}}+\int_{\\{\phi=\phi_{\min}\\}\cap\mathcal{A}}\right]F_{\phi}(b,0;t^{\prime})\,\mathrm{d}b\mathrm{d}t^{\prime}$
(27) $\displaystyle=E(\mathcal{U})+E(\mathcal{T})+E(\mathcal{A}).$ (28)

Figure 14(a) shows the total and partitioned entrained volume and figure 14
shows the entrainment rate for the regions $\mathcal{U},\mathcal{T}$ and
$\mathcal{A}$.

Figure 14: (a) Entrained volume decomposition into classes U, T, A and (b)
entrainment rate in each class.

The relative contribution of ‘fountain-top’ (i.e. plume cap, class T and large
buoyancy) entrainment and ‘lateral’ (i.e. intrusion, class A, low buoyancy)
entrainment to the overall entrainment is noted as an open question in the
study of fountains by Hunt & Burridge (2015). Here we find that the
contribution to the total entrained volume from the plume cap is weak compared
with the intrusion, suggesting that whilst mixing in the overturning process
near $z_{\max}$ is important for setting the maximum accessible buoyancy of
the volume distribution after mixing, entrainment of environmental fluid in
the lower stratified layer is more important for setting the centre-of-mass of
the final distribution. However, the entrainment rate is comparable in class T
and class A despite a significantly smaller volume, indicating stronger mixing
with the environment in the plume cap in quasi-steady-state, as well as during
initial penetration as identified earlier.

###  4.6 Discussion

The results suggest a multi-stage mixing process in quasi-steady state
convective penetration which is illustrated in figure 15. In this figure, we
show schematics of the volume distribution partitioned into the three regions
$\mathcal{U},\mathcal{T}$, and $\mathcal{A}$ in $(b,\phi)$-space and the
corresponding classes of fluid in an $x-z$ cross-section of the plume for each
stage of mixing. An individual fluid parcel moves through each stage from left
to right, corresponding with an increasing value of the cumulative mixed
volume $M$ associated with the parcel, though occasionally the primary mixing
stage may be skipped – see following discussion. In each stage, the convex
envelope of fluid involved is indicated by a gray dashed line in
$(b,\phi)$-space; arrows in $(b,\phi)$-space indicate the movement of
individual fluid parcels due to the mixing process; and arrows in the $x-z$
cross-section indicate the location of mixing in physical space. In the
secondary mixing stage where multiple mixing processes occur simultaneously,
three arrow types are shown which correspond between physical and
$(b,\phi)$-space. The mixing stages in QSS are as follows:

  1. 1.

Mixing within the undiluted plume as fluid rises through the stratified layer
sets the buoyancy-tracer distribution of plume fluid that is first exposed to
environmental fluid near $z_{\max}$. In $(b,\phi)$-space, the plume mixing
stage acts within $\mathcal{U}$ only, homogenising the distrbution and
consequently moving fluid towards the CoM of class U indicated by the blue
triangle in figure 11.

  2. 2.

Primary mixing between the plume and environment occurs in a thin plume cap
region of approximately constant volume near $z_{\max}$. This region may
extend below the plume cap into the intrusion as environmental fluid is
entrained by large eddies that form where subsiding fluid begins to move
radially. In $(b,\phi)$-space, the primary mixing stage acts on the undiluted
plume fluid concentrated near the class U CoM and the surrounding
environmental fluid with buoyancies lying between
$\left.b^{\text{(env)}}\right|_{z_{n}}$ and
$\left.b^{\text{(env)}}\right|_{z_{\max}}$. The resulting mixture is
transported towards the centre of the $\mathcal{T}$ region of $(b,\phi)$-space
as indicated by the field lines in figure 10.

  3. 3.

After environmental fluid is entrained into the plume in the primary mixing
stage, subsiding fluid joins the intrusion and homogenises with its
surroundings within the plume and therefore accesses larger $M$, moving from
class T to class A. Moreover, a number of secondary mixing processes result in
the entrainment of further environmental fluid lower in the stratified layer,
surrounding the intrusion. In $(b,\phi)$-space, mixing occurs between four
regions of fluid as indicated by the convex envelope in figure 15. First, as
subsiding fluid begins to move radially, environmental fluid at the base of
the plume cap is entrained and mixed. This process is indicated by the dashed
arrows. At the edges of the rising plume, undiluted plume fluid mixes with
fluid already in the intrusion, resulting in some fluid parcels moving
directly from class U to class A without entering the primary mixing stage, as
indicated by the dotted arrows in figure 15 and seen in figure 10. Secondary
mixing also occurs at the edges of the intrusion as environmental fluid in the
lower stratified layer is entrained, as indicated by the solid arrows in
figure 15. Finally, turbulent mixing in the interior of the intrusion
homogenises the buoyancy-tracer distribution that accumulates within region
$\mathcal{A}$. As fluid moves radially we therefore find large volumes of
fluid in class A concentrated near the CoM in $(b,\phi)$-space, where $M$ is
largest.

Figure 15: Schematic of the three stage mixing process identified by the
partitioned volume distribution in quasi-steady state convective penetration.
In the initial mixing stage, turbulence in the rising undiluted plume
homogenises the buoyancy-tracer distribution. Primary mixing occurs as fluid
overturns near $z_{\max}$ and first mixes with the surrounding environment.
Secondary mixing occurs as fluid joins the intrusion and spreads radially.

##  5 Mixing diagnostics

The three-stage mixing process identified via the partitioned volume
distribution may be further examined by considering the mixing properties of
fluid lying in each class. In this way, we gain insight into how these mixing
properties may change when the controlling parameters $F_{0}$ and $N$ are
varied.

###  5.1 Characterising mixing

To characterise the mixing properties of fluid in each stage, it is
illustrative to consider the ‘efficiency’ of mixing. The mixing efficency
seeks to relate energy expended in turbulent mixing with the actual mixing
achieved (e.g. Davies Wykes et al. (2015)). The most useful definition depends
on context (e.g. Gregg et al. (2018) for oceanographic contexts and Chemel &
Staquet (2007) for an atmospheric setting) but generally involves two rates of
energy dissipation; the irreversible increase in background potential energy
due to mixing and the kinetic energy dissipated via turbulent motion.

Following Howland et al. (2020) and Holliday & McIntyre (1981), we may treat
the positive semi-definite form of the potential energy,

$E_{p}=\frac{1}{2}\langle b^{\prime 2}\rangle,$ (29)

as a proxy for available potential energy, where $b^{\prime}$ is the departure
from the linear background stratification and $\langle\cdot\rangle$ denotes a
volume average. The irreversible loss of $E_{p}$ that results from the
reduction of buoyancy gradients by mixing is captured by the buoyancy variance
dissipation rate,

$\overline{\chi}=\left\langle\left(\kappa+\kappa_{\text{SGS}}\right)\,|\bnabla
b^{\prime}|^{2}\right\rangle,$ (30)

as the primary sink of potential energy – see appendix C for a derivation of
the full potential energy budget. The energy dissipated via turbulent mixing
is the sum of $\overline{\chi}$ and the dissipation rate of turbulent kinetic
energy,

$\overline{\varepsilon}=\left\langle\left(\nu+\nu_{\text{SGS}}\right)\frac{\partial
u_{i}}{\partial x_{j}}\frac{\partial u_{i}}{\partial x_{j}}\right\rangle,$
(31)

which acts as the primary sink of turbulent kinetic energy. The instantaneous
mixing efficiency $\eta$ is then

$\eta=\frac{\overline{\chi}}{\overline{\chi}+\overline{\varepsilon}}.$ (32)

Here, we use an overbar to denote a volume averaged quantity e.g.
$\overline{\chi}$ which is a scalar, as opposed to the scalar field which is
being averaged e.g. $\chi=\left(\kappa+\kappa_{\text{SGS}}\right)|\bnabla
b^{\prime}|^{2}$. It is instructive to examine the spatial structure of the
scalar fields themselves as the regions which contribute most to the volume
average become clear. Moreover, regions of large $\chi$ indicate intense
buoyancy gradients and regions of large $\varepsilon$ indicate intense
turbulence.

We may also examine the state or ‘activity’ of turbulence in stratified flows
using a pointwise activity parameter,

$I\equiv\frac{\frac{\partial u_{i}}{\partial x_{j}}\frac{\partial
u_{i}}{\partial x_{j}}}{\left|\frac{\partial b}{\partial
z}\right|}=\frac{\varepsilon}{\left(\nu+\nu_{\text{SGS}}\right)\left|\frac{\partial
b}{\partial z}\right|},$ (33)

which can be interpreted as the ratio of the destabilising effects of
turbulent stirring to the stabilising effects of buoyancy and viscosity.
Similarly, $I$ may be treated as the ratio of the (local) buoyancy timescale
$(\partial b/\partial z)^{-1/2}$ to the timescale of development of turbulent
effects $\sqrt{\varepsilon/\left(\nu+\nu_{\text{SGS}}\right)}$ (Ivey et al.,
2008). Regions of large $I$ indicate active turbulence (García-Villalba & del
Álamo, 2011) that is weakly affected by turbulence as it develops quickly. The
parameter $I$ is analogous to the buoyancy Reynolds number $\Rey_{b}$ except
with SGS viscosity included. Moreover, note that the global buoyancy timescale
$N^{-1}$ is often used to define $\Rey_{b}$. However, buoyancy gradients
within the plume differ significantly from the background stratification, so
we instead use a local measure of the buoyancy timescale given by the local
buoyancy gradient. Also note that typically $I$ is defined as a global measure
of the state of turbulence. Here, we consider $I$ as a pointwise field and use
the relative difference in $I$ between distinct regions to indicate the
relative intensity of turbulence in each of those regions. For summary
statistics we define the mean activity parameter as

$\overline{I}=\frac{\langle\varepsilon/\left(\nu+\nu_{\text{SGS}}\right)\rangle}{\langle\left|\frac{\partial
b}{\partial z}\right|\rangle}.$ (34)

Figure 16: $x-z$ cross sections of the $\varepsilon,\chi,I$ and
$\partial_{z}b$ fields at $t=15$.

###  5.2 Results

Figure 16 shows $x-z$ cross sections of the $\varepsilon,\chi,I$ and
$\partial_{z}b$ fields in the plume where $\phi>\phi_{\min}$, with the
surrounding region showing the environmental buoyancy. The results show that
the partition of plume fluid into three regions identifies mixing stages with
different mixing diagnostics: in the rising undiluted plume fluid and the
surrounding subsiding fluid, $\varepsilon$ is largest and roughly an order of
magnitude larger than in the spreading intrusion. $I$ is largest in the rising
component of the plume and comparable between the plume cap and intrusion. The
strongest buoyancy gradients as measured by $\chi$ are found in the primary
mixing region of the overturning and subsiding fluid surrounding the plume
cap. Small regions of intense buoyancy gradient are also found in the rising
column of undiluted plume fluid, but $\chi$ is more sparse here than in the
plume cap. The strongest _vertical_ buoyancy gradients are found at the top
edge of the plume cap, where overturning fluid impinges on the surrounding
environmental fluid. This region coincides with the layer of largest $\chi$.
There is evidence of layering in the intrusion, indicated by the alternating
positive and negative layers of $\partial_{z}b$.

Figures 17, 18 and 19 show PDFs of $\varepsilon,\chi$ and $I,\partial_{z}b$
respsectively. In each figure, the total PDF of the field within the plume
(i.e. where $\phi>10^{-5}$ is shown with a black dashed line. The PDFs of the
field partitioned between classes U, T, and A are shown with solid coloured
lines. The partitioned PDFs are scaled such that their sum matches the total
PDF. We may therefore use figures 17, 18 and 19 to identify the relative sizes
of each field between classes. Note that figures 17 and 18 show a bimodal
distribution for $\varepsilon$ and $\chi$, repsectively. This arises from the
$\nu+\nu_{\text{SGS}}$ and $\kappa+\kappa_{\text{SGS}}$ factor respectively,
each of which is bimodal with a peak where the SGS viscosity or diffusivity
vanishes and a peak where the SGS viscosity or diffusivity is non-zero (see
right panel of figures 17 and 18). To demonstrate this, the partitioned PDFs
in figures 17 and 18 are further separated into PDFs where
$\nu_{\text{SGS}}=0$ and where $\nu_{\text{SGS}}>0$. Now, the sum of all 6
partitioned PDFs equals the total PDF. Evidently, regions where the SGS
viscosity vanishes – and therefore where $\nu+\nu_{\text{SGS}}$ reduces to the
molecular viscosity only – correspond to the smaller peak in $\varepsilon$ and
similarly for $\chi$. These points where the SGS viscosity or diffusivity
vanishes are effectively points where the simulation has switched from LES to
DNS.

Figure 17: PDF of the TKE dissipation rate $\varepsilon$ (left) for the full
plume (black dashed line) and partitioned into classes U, T, A (coloured
lines) at $t=15$. PDF of viscosity $\nu+\nu_{\text{SGS}}$ (right). Figure 18:
PDF of the buoyancy variance dissipation rate $\chi$ (left) for the full plume
(black dashed line) and partitioned into classes U, T, A (coloured lines) at
$t=15$. PDF of diffusivity $\kappa+\kappa_{\text{SGS}}$ (right). Figure 19:
PDFs of activity parameter $I$ (left) and vertical buoyancy gradient
$\partial_{z}b$ (right) for the full plume (black dashed line) and partitioned
into classes U, T, A (coloured lines) at $t=15$.

The PDFs show that $\varepsilon$ is largest in classes U and T and
approximately one order of magnitude smaller in class A which makes up the
majority of the volume of the plume. Similarly, $\chi$ is largest in class T;
an order of magnitude larger than in class A and U. The largest values of
$\chi$ are almost all in class T, with a clear positive bias in the PDF. The
activity parameter $I$ is largest in class U and smaller in classes T and A
with similar distributions. The vertical buoyancy gradient $\partial_{z}b$ is
shown on a log-scale since the most extreme values are rare but most important
for mixing. The extreme values are mostly in class T, particularly large
positive gradients, with the remainder contained primarily in class A. The
relatively tighter spread of the $\partial_{z}b$ PDF in class U indicates that
typically smaller vertical buoyancy gradients are present compared with
classes U and A. The PDF in class T has a clear positive bias, whilst the PDFs
in classes S and A have similar shapes. Note that the generally smaller $I$ in
class T is subject to interpretation: $\varepsilon$ is large but since
$\partial_{z}b$ is also generally large, the ratio is smaller than in class U
where $\partial_{z}b$ is smaller. The volume averages of each mixing quantity
are given in table 2 for the full plume and within each class. To aid in
interpreting the results, the percentage of the total plume volume in each
class is also given.

| Full plume | Class U | Class T | Class A
---|---|---|---|---
Volume % | $100$ | $4.27$ | $11.5$ | $84.2$
Activity parameter $\overline{I}$ | $12.7$ | $81$ | $16$ | $8.37$
Vertical buoyancy gradient $\langle\left|\partial_{z}b\right|\rangle$ | $1.23$ | $0.994$ | $2.49$ | $1.07$
TKE dissipation rate $\overline{\varepsilon}$ | $0.0414$ | $0.292$ | $0.116$ | $0.0185$
Buoyancy variance dissipation rate $\overline{\chi}$ | $0.026$ | $0.0343$ | $0.128$ | $0.0115$
Instantaneous mixing efficiency $\eta$ | $0.386$ | $0.105$ | $0.526$ | $0.384$
Table 2: Volume-averaged mixing quantities defined in §5.1 at $t=15$.

In class U, we find active turbulence with large dissipation of KE, consistent
with the undiluted plume being unaffected by the surrounding stratified
environment in quasi-steady state. As the undiluted plume fluid is well-mixed,
there are relatively small buoyancy gradients. As fluid rises, KE is converted
into PE, resulting in large PE accumulation and relatively little PE
dissipation due to the weak buoyancy gradients. Some PE dissipation is found
in large eddies which have not yet broken down. The mixing efficiency in this
region is low as a result. In class T, undiluted plume fluid encounters the
surrounding environmental fluid with a much larger buoyancy, resulting in
particularly strong vertical buoyancy gradients. Horizontal buoyancy gradients
are also generated by overturning. Turbulent motion advected upwards in class
U is carried over into class T but this motion is gradually suppressed by
buoyancy gradients, reducing the activity parameter. However, the intense
buoyancy gradients are nonetheless enhanced by turbulent stirring which
results in strong PE dissipation and consequently the mixing is particularly
efficient. Finally in class A, which eventually dominates the volume of the
full plume as the intrusion grows, the interior of the intrusion becomes well-
mixed resulting in weak buoyancy gradients, though there is some indication of
a weak stratification forming within the intrusion. Secondary mixing processes
between the intrusion and environmental fluid at the bottom of the stratified
layer can introduce larger vertical buoyancy gradients and result in some PE
dissipation. Overall, KE dissipation is weak as the buoyant forces driving
turbulent motion in the remainder of the plume are weakened. However, whilst
both KE and PE dissipation are weak, they are almost equal and hence the
mixing is moderately efficient (compared with class U). This could be
attributed to the fact that despite relatively weaker turbulent motion,
reduced buoyancy gradients are not able to suppress the turbulent stirring
that does arise.

##  6 Conclusions

In this paper we have analysed a large eddy simulation of a buoyant plume
penetrating into a stably stratified layer and used a buoyancy-tracer volume
distribution to examine tracer transport via turbulent mixing. We developed a
method for objectively partitioning buoyancy-tracer space into three classes,
which correspond with coherent regions of physical space. The method
identifies undiluted plume fluid (class U) and mixtures of plume and
environmental fluid. Further, the partitioning separates the mixed fluid into
newly-generated mixtures in the plume cap and older mixed fluid in the
radially spreading intrusion. In buoyancy-tracer space, the newly-generated
mixtures lie in a transport region (class T) where volume does not accumulate,
whilst the intrusion corresponds with an accumulation region (class A) where
fluid collects and homogenises, which represents the majority of the plume
volume at late times. We use the buoyancy variance dissipation rate, TKE
dissipation rate, vertical buoyancy gradient and an activity parameter as
diagnostics of the mixing process in each sub-volume of the plume. We find
that very intense buoyancy gradients in the plume cap (class T) result in the
efficient mixing of undiluted plume fluid with the surrounding environment. As
newly-generated mixed fluid joins the intrusion and homogenises with fluid
already in the intrusion (class A), the intensity of turbulence decreases and
buoyancy gradients weaken. Nonetheless, we find moderately efficient
entrainment of environmental fluid surrounding the intrusion, which dominates
the overall volume of fluid entrained into the plume. However, the largest
entrainment rates occur intermittently in the plume cap.

Our results demonstrate a three-stage mixing process in quasi-steady state
penetration of a plume into a stratified layer. In the first stage, turbulence
within the buoyant plume homogenises the buoyancy-tracer distribution inside
the plume as fluid rises through the stratified layer, unaffected by the
surrounding stratification. This process erodes extreme values of buoyancy and
tracer concentration in the plume core and sets the range of buoyancy and
tracer concentration which is first exposed to the environment in the
following stage. In the second or ‘primary mixing’ stage, plume and
environmental fluids are efficiently mixed as undiluted plume fluid overturns
near the maximum penetration height. This height approximately determines the
buoyancy of the environmental fluid involved in the mixing, and therefore sets
the extent of the buoyancy-tracer convex envelope for the remainder of the
mixing process. The third or ‘secondary mixing’ stage involves several mixing
processes that act to homogenise mixed fluid in the intrusion. This includes
the entrainment of environmental fluid surrounding the intrusion and mixing
with small amounts of undiluted plume fluid at the edge of the rising plume
that immediately join the intrusion without entering the primary mixing stage.
Also, environmental fluid from near the maximum penetration height can be
pulled down by subsiding plume fluid and mixed with the intrusion.

The statistical properties of turbulence are different in each of the three
stages. The undiluted plume core is the most turbulent (as measured by the
activity parameter) with KE dissipation significantly larger than PE
dissipation. In the plume cap, intense buoyancy gradients are established as
plume fluid impinges on the buoyant surroundings, producing large PE
dissipation compared with small KE dissipation and hence efficient entrainment
of the surrounding environment. As mixed fluid homogenises in the intrusion
and further environmental fluid is entrained, buoyancy gradients are eroded
resulting in some weak PE dissipation. Here, the KE dissipation is also weak
so the mixing is moderately efficient, though less efficient than in the plume
cap.

Models which cannot resolve the processes responsible for mixing and
entrainment must parameterise the effects of mixing on the flow. The markedly
different statistics in each class suggest that each stage of mixing should be
parameterised separately. Parameterisations of mixing in convective
penetration could exploit the changing dynamical balances as the flow evolves.
For example, at early times before the formation of an intrusion, the plume is
dominated by undiluted plume fluid in the plume core with intense turbulence
but inefficient mixing. At late times as fluid accumulates in the intrusion,
the plume volume is dominated by the intrusion with weaker buoyancy gradients
and turbulence. Therefore the fraction of plume fluid in each class changes
with time, as do the turbulent statistics associated with the full plume
volume.

The partitioning method presented here, as well as the buoyancy-tracer volume
distribution formalism as a whole, offers a way to analyse mixing in numerical
simulations of stratified flows. Physical arguments can be made that restrict
the regions of buoyancy-tracer space accessible via mixing and consideration
of terms in the distribution budget equation (9) highlight the mixing
processes that occur and the resulting tracer transport.

## References

  * Ansong et al. (2008) Ansong, J. K., Kyba, P. J. & Sutherland, B. R. 2008 Fountains impinging on a density interface. J. Fluid Mech. 595, 115–139.
  * Ansong & Sutherland (2010) Ansong, J. K. & Sutherland, B. R. 2010 Internal gravity waves generated by convective plumes. J. Fluid Mech. 648, 405–434.
  * Chemel & Staquet (2007) Chemel, C. & Staquet, C. 2007 A formulation of convective entrainment in terms of mixing efficiency. J. Fluid Mech. 580, 169–178.
  * Chen & Bhaganagar (2023) Chen, C. H. & Bhaganagar, K. 2023 Energetics of buoyancy-generated turbulent flows with active scalar: pure buoyant plume. J. Fluid Mech. 954, A23.
  * Craske & van Reeuwijk (2015) Craske, J. & van Reeuwijk, M. 2015 Energy dispersion in turbulent jets. Part 1. Direct simulation of steady and unsteady jets. J. Fluid Mech. 763, 500–537.
  * Dauhut et al. (2015) Dauhut, T., Chaboureau, J.-P., Escobar, J. & Mascart, P. 2015 Large-eddy simulations of hector the convector making the stratosphere wetter. Atmos. Sci. Lett. 16 (2), 135–140.
  * Dauhut et al. (2018) Dauhut, T., Chaboureau, J-P., Haynes, P. H. & Lane, T. P. 2018 The mechanisms leading to a stratospheric hydration by overshooting convection. J. Atmos. Sci. 75 (12), 4383 – 4398\.
  * Davies Wykes & Dalziel (2014) Davies Wykes, M. S. & Dalziel, S. B. 2014 Efficient mixing in stratified flows: Experimental study of a rayleigh–taylor unstable interface within an otherwise stable stratification. J. Fluid Mech. 756, 1027–57.
  * Davies Wykes et al. (2015) Davies Wykes, M. S., Hughes, G. O. & Dalziel, S. B. 2015 On the meaning of mixing efficiency for buoyancy-driven mixing in stratified turbulent flows. J. Fluid Mech. 781, 261–275.
  * Devenish et al. (2010) Devenish, B. J., Rooney, G. G. & Thomson, D. J. 2010 Large-eddy simulation of a buoyant plume in uniform and stably stratified environments. J. Fluid Mech. 652, 75–103.
  * Doda et al. (2023) Doda, T., Ulloa, H. N., Ramón, C. L., Wüest, A. & Bouffard, D. 2023 Penetrative Convection Modifies the Dynamics of Downslope Gravity Currents. Geophysical Research Letters 50 (2), e2022GL100633.
  * Flynn & Sutherland (2004) Flynn, M. R. & Sutherland, B. R. 2004 Intrusive gravity currents and internal gravity wave generation in stratified fluid. J. Fluid Mech. 514, 355–383.
  * Fueglistaler et al. (2009) Fueglistaler, S., Dessler, A. E., Dunkerton, T. J., Folkins, I., Fu, Q. & Mote, P. W. 2009 Tropical tropopause layer. Rev. Geophys. 47 (1).
  * García-Villalba & del Álamo (2011) García-Villalba, M. & del Álamo, J. C. 2011 Turbulence modification by stable stratification in channel flow. Physics of Fluids 23 (4), 045104.
  * Gregg et al. (2018) Gregg, M.C., D’Asaro, E.A., Riley, J.J. & Kunze, E. 2018 Mixing Efficiency in the Ocean. Annu. Rev. Mar. Sci. 10, 443–473.
  * Herrmann et al. (2008) Herrmann, Marine, Somot, Samuel, Sevault, Florence, Estournel, Claude & Déqué, Michel 2008 Modeling the deep convection in the northwestern Mediterranean Sea using an eddy-permitting and an eddy-resolving model: Case study of winter 1986–1987. J. Geophys. Res. 113 (C4), C04011.
  * Holliday & McIntyre (1981) Holliday, D. & McIntyre, M. E. 1981 On potential energy density in an incompressible, stratified fluid. J. Fluid Mech. 107, 221–225.
  * Howland et al. (2020) Howland, C. J., Taylor, J. R. & Caulfield, C. P. 2020 Mixing in forced stratified turbulence and its dependence on large-scale forcing. J. Fluid Mech. 898, A7.
  * Hunt & Burridge (2015) Hunt, G.R. & Burridge, H.C. 2015 Fountains in Industry and Nature. Annu. Rev. Fluid Mech. 47 (1), 195–220.
  * Ivey et al. (2008) Ivey, G.N., Winters, K.B. & Koseff, J.R. 2008 Density Stratification, Turbulence, but How Much Mixing? Annu. Rev. Fluid Mech. 40 (1), 169–184.
  * Jensen et al. (2007) Jensen, E. J., Ackerman, A. S. & Smith, J. A. 2007 Can overshooting convection dehydrate the tropical tropopause layer? J. Geophys. Res. Atmos. 112 (D11).
  * Kurbatskii (2001) Kurbatskii, Albert F. 2001 Computational Modeling of the Turbulent Penetrative Convection above the Urban Heat Island in a Stably Stratified Environment. J. Appl. Meteor. 40 (10), 1748–1761.
  * Lane & Reeder (2001) Lane, T. P. & Reeder, M. J. 2001 Convectively Generated Gravity Waves and Their Effect on the Cloud Environment. J. Atmos. Sci. 58 (16), 2427–2440.
  * Lumley & Newman (1977) Lumley, J. L. & Newman, G. R. 1977 The return to isotropy of homogeneous turbulence. J. Fluid Mech. 82 (1), 161–178.
  * Marshall & Schott (1999) Marshall, John & Schott, Friedrich 1999 Open-ocean convection: Observations, theory, and models. Rev. Geophys. 37 (1), 1–64.
  * Masada et al. (2013) Masada, Youhei, Yamada, Kohei & Kageyama, Akira 2013 Effects of Penetrative Convection on Solar Dynamo. ApJ 778 (1), 11.
  * Morton (1959) Morton, B. R. 1959 Forced plumes. J. Fluid Mech. 5 (1), 151–163.
  * Morton et al. (1956) Morton, B. R., Taylor, G. I. & Turner, J. S. 1956 Turbulent gravitational convection from maintained and instantaneous sources. Proc. R. S. London. Series A. 234 (1196), 1–23.
  * Penney et al. (2020) Penney, J., Morel, Y., Haynes, P. H., Auclair, F. & Nguyen, C. 2020 Diapycnal mixing of passive tracers by Kelvin–Helmholtz instabilities. J. Fluid Mech. 900, A26.
  * Pham et al. (2007) Pham, M. V., Plourde, F. & Doan, K. S. 2007 Direct and large-eddy simulations of a pure thermal plume. Phys. Fluids 19 (12), 125103.
  * Pinçon et al. (2016) Pinçon, C., Belkacem, K. & Goupil, M. J. 2016 Generation of internal gravity waves by penetrative convection. A&A 588, A122.
  * Plumb (2007) Plumb, R. A. 2007 Tracer interrelationships in the stratosphere. Rev. Geophys. 45 (4).
  * van Reeuwijk et al. (2016) van Reeuwijk, M., Salizzoni, P., Hunt, G. R. & Craske, J. 2016 Turbulent transport and entrainment in jets and plumes: A DNS study. Phys. Rev. Fluids 1, 074301.
  * Singh et al. (1994) Singh, Harinder P., Roxburgh, Ian W. & Chan, Kwing L. 1994 Three-dimensional simulation of penetrative convection-penetration above a convection zone. A&A 281 (2), L73–L76.
  * Stull (1976) Stull, Roland B. 1976 Internal gravity waves generated by penetrative convection. J. Atmos. Sci. 33 (7), 1279 – 1286\.
  * Taylor (2008) Taylor, J. R. 2008 Numerical simulations of the stratified oceanic bottom layer. PhD thesis, University of California, San Diego.
  * Turner (1966) Turner, J. S. 1966 Jets and plumes with negative or reversing buoyancy. J. Fluid Mech. 26 (4), 779–792.
  * Ulses et al. (2021) Ulses, Caroline, Estournel, Claude, Fourrier, Marine, Coppola, Laurent, Kessouri, Fayçal, Lefèvre, Dominique & Marsaleix, Patrick 2021 Oxygen budget of the north-western Mediterranean deep- convection region. Biogeosciences 18 (3), 937–960.
  * Vreugdenhil & Taylor (2018) Vreugdenhil, C. A. & Taylor, J. R. 2018 Large-eddy simulations of stratified plane couette flow using the anistropic minimum-dissipation model. Phys. Fluids 30 (8), 085104.

##  Appendix A Validating the numerical method

The capability for our numerical method to reliably simulate a buoyant plume
was validated by detailed comparison between our own LES of a plume in a
uniform environment with direct numerical simulations in the literature of the
same setup. Individual simulations are validated using the framework of van
Reeuwijk et al. (2016), in which the integral theory of Morton et al. (1956)
(henceforth MTT) is used to verify the self-similarity of key variables when
normalised with respect to characteristic scales $w_{m},b_{m}$ as defined by
MTT. Turbulent corrections are included in the form of profile coefficients
which, in a radially averaged sense, account for the eddy terms neglected when
deriving the MTT plume equations and act to modify the shape of self-similar
profiles for each variable (Craske & van Reeuwijk, 2015). An example of this
validation method for a single simulation is shown in figure 20 which shows
azimuthally averaged radial profiles in an interval of heights close to the
top of the uniform layer, where the turbulence and plume structure is
considered fully developed. The vertical velocity $\overline{w}$ and buoyancy
$\overline{b}$ exhibit well-converged Gaussian radial profiles. The normalised
turbulent radial momentum and buoyancy fluxes
$\overline{u^{\prime}w^{\prime}}$ and $\overline{u^{\prime}b^{\prime}}$ also
exhibit the expected self-similarity and suggests that the essential turbulent
plume dynamics are well represented by the numerical scheme. Moreover, the
normalised mean radial specific volume flux $r\overline{u}$ demonstrates that
the MTT entrainment hypothesis $\left[ru\right]_{r=\infty}=-\alpha r_{m}w_{m}$
is supported, indicating the turbulent entrainment in the numerical
simulations broadly matches that found in laboratory experiments. The mean
radial velocity $\overline{u}$ exhibits greater spread in the radial profiles
especially at large radii, likely due to the periodic (i.e. effectively
closed) boundaries compared with the open boundaries used by van Reeuwijk et
al. (2016). Nonetheless, the reasonable convergence of the profiles suggests
that edge effects are not influencing the behaviour of the plume.

Figure 20: Self-similar radial profiles of $\overline{w},\overline{b}$ (top
left), radial momentum and buoyancy fluxes $\overline{u^{\prime}w^{\prime}}$
and $\overline{u^{\prime}b^{\prime}}$ (top right), mean radial velocity
$\overline{u}$ (bottom left) and normalised mean radial specific volume flux
$r\overline{u}$ (bottom right) in the interval $35<z/r_{0}<70$. Dotted line
indicates value of the entrainment coefficient $\alpha$ calculated from a
linear fit of the plume radius and $z$.

The validity of the volumetric forcing method is further supported by
calculating the flux-balance parameter $\Gamma$ introduced by Morton (1959) –
analogous to a plume Richardson number – which is found to be close to $1/2$,
indicating a ‘forced’ plume that is inertia-dominated. We also verify the
assumption of a Gaussian radial profile for the tracer concentration by
testing multiple forcing profiles, all of which result in a Gaussian
concentration profile once the plume structure is fully developed. Finally,
the precise characteristics of the turbulence are tested by calculating the
invariants of the anisotropy tensor as detailed in Lumley & Newman (1977). Our
results indicate turbulence with weak anistropy and axisymmetry, agreeing with
the direct numerical simulations of van Reeuwijk et al. (2016).

##  Appendix B Tracer-tracer volume distribution budget

Consider two tracers with concentrations $b,\phi$ satisfying

$\frac{\partial b}{\partial t}+\boldsymbol{u}\cdot\bnabla
b=\dot{b},\frac{\partial\phi}{\partial
t}+\boldsymbol{u}\cdot\bnabla\phi=\dot{\phi},$ (35)

with $\bnabla\cdot\boldsymbol{u}=0$. Consider a fixed volume $V$ containing
fluid and define

$g_{F}=\int_{V}F(b,\phi)\mathrm{d}V.$ (36)

Then

$\displaystyle\frac{\mathrm{d}g_{F}}{\mathrm{d}t}$
$\displaystyle=\int_{V}\frac{\partial F}{\partial b}\frac{\partial b}{\partial
t}+\frac{\partial F}{\partial\phi}\frac{\partial\phi}{\partial
t}\,\mathrm{d}V$ (37) $\displaystyle=\int_{V}\frac{\partial F}{\partial
b}\left(\dot{b}-\boldsymbol{u}\cdot\bnabla b\right)+\frac{\partial
F}{\partial\phi}\left(\dot{\phi}-\boldsymbol{u}\cdot\bnabla\phi\right)\,\mathrm{d}V$
$\displaystyle=\int_{V}\frac{\partial F}{\partial b}\dot{b}+\frac{\partial
F}{\partial\phi}\dot{\phi}\,\mathrm{d}V-\int_{V}\left(\boldsymbol{u}\cdot\bnabla
b\right)\frac{\partial
F}{\partial\phi}+\left(\boldsymbol{u}\cdot\bnabla\phi\right)\frac{\partial
F}{\partial\phi}\,\mathrm{d}V$ $\displaystyle=\int_{V}\frac{\partial
F}{\partial b}\dot{b}+\frac{\partial
F}{\partial\phi}\dot{\phi}\,\mathrm{d}V-\int_{V}\bnabla\cdot\left(\boldsymbol{u}F(b,\phi)\right)\,\mathrm{d}V,$

since $\bnabla\cdot\boldsymbol{u}=0$. By the divergence theorem we have

$\frac{\mathrm{d}g_{F}}{\mathrm{d}t}=\int_{V}\frac{\partial F}{\partial
b}\dot{b}+\frac{\partial
F}{\partial\phi}\dot{\phi}\,\mathrm{d}V-\int_{\partial
V}\boldsymbol{u}\cdot\boldsymbol{n}F(b,\phi)\,\mathrm{d}S,$ (38)

where $\partial V$ is the boundary of $V$ and $\boldsymbol{n}$ is the outward
normal. Now, consider $F(b,\phi)=I(b;b_{0})I(\phi;\phi_{0})$ where $I$ is the
indicator function for the subset of $V$ where $b>b_{0}$, i.e. $I(b;b_{0})=1$
if $b>b_{0}$ and $I(b;b_{0})=0$ otherwise. Then, $g_{F}$ is the volume of
fluid in $V$ with $b>b_{0}$ and $\phi>\phi_{0}$. Furthermore
$\frac{\partial^{2}}{\partial b_{0}\partial\phi_{0}}g_{F}$ is the ‘volume
density’, i.e. $\frac{\partial^{2}}{\partial b_{0}\partial\phi_{0}}g_{F}\delta
b\delta\phi$ is the volume of fluid in $V$ with $b_{0}<b<b_{0}+\delta b$ and
$\phi_{0}<\phi<\phi_{0}+\delta\phi$.

Now, for $F(b,\phi)=I(b;b_{0})I(\phi;\phi_{0})$ we have

$\displaystyle\int_{V}\frac{\partial F}{\partial b}\dot{b}\,\mathrm{d}V$
$\displaystyle=\int_{V}\dot{b}\delta(b-b_{0})I(\phi;\phi_{0})\,\mathrm{d}V$
(39) $\displaystyle=\int_{S(b_{0},\phi_{0})}\dot{b}\frac{\mathrm{d}S}{\partial
b/\partial n}$ $\displaystyle=-\frac{\partial}{\partial
b_{0}}\int_{V}\dot{b}I(b;b_{0})I(\phi;\phi_{0})\mathrm{d}V,$

where $S(b_{0},\phi_{0})$ is the surface in $V$ where $b=b_{0}$ and
$\phi>\phi_{0}$. Similarly,

$\int_{V}\frac{\partial
F}{\partial\phi}\dot{\phi}\,\mathrm{d}V=-\frac{\partial}{\partial\phi_{0}}\int_{V}\dot{\phi}I(b;b_{0})I(\phi;\phi_{0})\,\mathrm{d}V.$
(40)

Returning to (38) we have the integral form of the budget equation which
governs the buoyancy-tracer cumulative volume distribution,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\int_{V}I(b;b_{0})I(\phi;\phi_{0})\,\mathrm{d}V$
$\displaystyle=-\frac{\partial}{\partial
b_{0}}\left[\int_{V}I(b;b_{0})I(\phi;\phi_{0})\dot{b}\,\mathrm{d}V\right]$
(41)
$\displaystyle-\frac{\partial}{\partial\phi_{0}}\left[\int_{V}I(b;b_{0})I(\phi;\phi_{0})\dot{\phi}\,\mathrm{d}V\right]$
$\displaystyle-\int_{\partial
V}\boldsymbol{u}\cdot\boldsymbol{n}\,I(b;b_{0})I(\phi;\phi_{0})\,\mathrm{d}S.$

The budget equation for the buoyancy-tracer volume distribution is obtained by
taking $\partial^{2}/\partial b_{0}\partial\phi_{0}$ to get

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\int_{V}\delta(b-b_{0})\delta(\phi-\phi_{0})\,\mathrm{d}V$
$\displaystyle=-\frac{\partial}{\partial
b_{0}}\left[\int_{V}\delta(b-b_{0})\delta(\phi-\phi_{0})\dot{b}\,\mathrm{d}V\right]$
(42)
$\displaystyle-\frac{\partial}{\partial\phi_{0}}\left[\int_{V}\delta(b-b_{0})\delta(\phi-\phi_{0})\dot{\phi}\,\mathrm{d}V\right]$
$\displaystyle-\int_{\partial
V}\boldsymbol{u}\cdot\boldsymbol{n}\,\delta(b-b_{0})\delta(\phi-\phi_{0})\,\mathrm{d}S.$

##  Appendix C Numerical method for identifying the accumulation region of
$(b,\phi)$-space

Mixed fluid in the stratified layer is identified by reigons of
$(b,\phi)$-space with $M(b,\phi;t)>0$. Mixed fluid is further partitioned by a
threshold value $m^{*}(t)$ such that:

  * •

$0<M\leq m^{*}$ identifies newly-generated mixed fluid, where undiluted plume
fluid is in the process of primary mixing with the environment;

  * •

$M>m^{*}$ identifies plume fluid which has mixed with the environment and is
accumulating in the radially spreading intrusion.

The threshold $m^{*}(t)$ is chosen to minimise the total flux divergence of
the transport region $\mathcal{T}=\\{0<M\leq m^{*}\\}$, i.e. $m^{*}(t)$ is the
value $m$ that minimises $\left|f(m;t)\right|$ where

$f(m;t)=\int_{\\{0<M\leq m\\}}\frac{\partial W}{\partial
t}-(-S)\,\mathrm{d}b\mathrm{d}\phi=\int_{\\{0<M\leq
m\\}}\bnabla_{(b,\phi)}\cdot\boldsymbol{F}\,\mathrm{d}b\mathrm{d}\phi.$ (43)

The numerical implementation of this method involves, at each time step $k$,
the following steps:

  1. 1.

identify the current maximum of $M$, denoted
$\mathcal{M}=\max_{i,j}M(b_{i},\phi_{j};t_{k})$;

  2. 2.

for each of $N_{M}$ test values of $M$, $m=0,\dots,\mathcal{M}$, calculate

$f(m;t_{k})=\sum_{i,j|0<M_{ij}\leq
m}\left[\frac{W_{ij}(t_{k+1})-W_{ij}(t_{k})}{\Delta t}+S_{ij}(t_{k})\right];$
(44)

  3. 3.

identify the test value $\tilde{m}(t_{k})$ which minimises $f(m;t_{k})$.

Then, once a threshold $\tilde{m}(t_{k})$ has been chosen for each timestep
$t_{k}$, the final threshold $m^{*}(t)$ is chosen by applying a rolling
average with an appropriate number of timesteps. We choose $N_{M}=200$ and a
rolling average width of 10 timesteps. The motivation for applying a rolling
average is to smooth the threshold $m^{*}(t)$ so that the regions
$\mathcal{T}$ and $\mathcal{A}$ do not grow and shrink dramatically in
response to short-term changes in the flux divergence.

Figure 21(a) shows $f(m;t)$ for all times post-penetration and figure 21(b)
shows the corresponding preliminary thresholds $\tilde{m}(t)$ and the final
thresholds $m^{*}(t)$. The total flux divergence of class T for the
preliminary and final thresholds over time is shown in figure 21(c),
indicating the (generally small) error introduced by smoothing the thresholds.

Figure 21: (a) Total flux divergence $f(m;t)$ defined in (43) and the
resulting preliminary and smoothed threshold choice in (b). Error in $f$ from
smoothing process shown in (c).

##  Appendix D Potential energy budget

We start with the governing equations (1) – (3) including the SGS terms and
drop the hat notation. Let $b=z+b^{\prime}$ so that $b^{\prime}$ is the
perturbation from the stratified environmental buoyancy profile. The buoyancy
evolution equation becomes

$\frac{\partial b^{\prime}}{\partial
t}+w+(\boldsymbol{u\cdot}\bnabla)b^{\prime}=\frac{1}{\Rey\Pr}\left[\bnabla^{2}b^{\prime}+\frac{1}{\kappa}\left(\bnabla\boldsymbol{\cdot}(\kappa_{\text{SGS}}\bnabla
b^{\prime})+\frac{\partial\kappa_{\text{SGS}}}{\partial z}\right)\right].$
(45)

When $b^{\prime}$ is the departure from a _linear_ stratification as
considered here, the potential energy (per unit mass) $E_{P}$ takes the
positive semi-definite form

$E_{P}=\frac{1}{2}\langle b^{\prime 2}\rangle$ (46)

where $\langle\cdot\rangle$ denotes a volume average and moreover, $E_{p}$
acts as a proxy for the available potential energy (Holliday & McIntyre,
1981). Multiplying (45) by the perturbation buoyancy $b^{\prime}$ and volume
averaging gives

$\frac{\mathrm{d}E_{p}}{\mathrm{d}t}=\int_{\partial
V}\kappa_{\text{eff}}\,b^{\prime}\bnabla
b^{\prime}\cdot\boldsymbol{\mathrm{d}S}-\overline{\chi}-\overline{J_{b}}-\frac{1}{\Rey\Pr}\frac{1}{\kappa}\langle
b^{\prime}\frac{\partial\kappa_{\text{SGS}}}{\partial z}\rangle$ (47)

Here, we choose the entire stratified layer $z\geq 0$ as the control volume
$V$. The first term on the LHS represents the diffusive buoyancy flux across
the boundary $\partial V$, which is non-zero only where the plume penetrates
the stratified layer. The second term is the volume averaged buoyancy variance
dissipation rate

$\overline{\chi}=\left\langle(\kappa+\kappa_{\text{SGS}})\,|\bnabla
b^{\prime}|^{2}\right\rangle,$ (48)

which represents the primary sink of potential energy. The third term is the
volume averaged buoyancy flux

$\overline{J_{b}}=\langle b^{\prime}w\rangle$ (49)

which represents an exchange between kinetic and potential energy. The last
term captures the effect of the spatially-varying turbulent diffusivity acting
on the background stratification.