Diagnosing tracer transport in convective penetration of a stably stratified layer

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 plume fluid in 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 a 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 lie in a thin layer at the cap of the penetrating plume. This provides the primary stage of mixing between plume and environment and exhibits a mixing efficiency around 50 %. Newly generated mixtures of environmental and plume fluid join the intrusion and experience relatively weak turbulence and buoyancy gradients. As the intrusion spreads radially, environmental fluid surrounding the intrusion is mixed into the intrusion with moderate mixing efficiency. This dominates the volume of environmental fluid entrained into the region containing plume fluid. However, the ‘strongest’ entrainment, as measured by the specific entrainment rate, is largest in the plume cap, where the most buoyant environmental fluid is entrained.

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, resulting in vertical transport of trace gases and water vapour (Jensen, Ackerman & Smith 2007; Randel & Jensen 2013). 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.

Another geophysical process where this fluid dynamical problem is relevant is 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), volcanic eruptions that penetrate into the stratosphere (e.g. Textor et al. 2003; Carazzo, Kaminski & Tait 2008) as well as the internal structure of many stars where a convectively unstable layer is bounded above and below by stable layers (Singh, Roxburgh & Chan 1994; Masada, Yamada & Kageyama 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. Much previous work has focused on laboratory studies of buoyant plumes with simple background density profiles, but numerical simulation has recently become feasible (Alfonsi 2011).

Progress towards understanding the contribution of convective penetration to tracer transport in geophysical settings can be made by considering an idealised representation of the problem in which a region of strong stable stratification is penetrated by a turbulent buoyant plume generated in a region with weak or zero stratification. The objective of the study reported in this paper is to diagnose the irreversible diapycnal tracer transport that results from turbulent mixing between plume fluid carrying a passive tracer and the surrounding environmental fluid in the stratified layer. We aim to provide a quantitative description of the mixing involved in this diapycnal transport. Such descriptions are essential in forming parameterisations of convective penetration. Throughout the flow evolution, plume fluid may be distinguished from environmental fluid by the presence of non-zero tracer concentration. Crucially, both the tracer concentration and buoyancy fields are subjected to turbulent mixing, resulting in the entrainment of environmental fluid into the plume and modification of the relationship between buoyancy and tracer within the plume. Plumb (2007) introduced a tracer–tracer probability density function to study rapid isentropic mixing in the stratosphere. Penney et al. (2020) utilised this method to study diapycnal mixing of passive tracers by Kelvin–Helmholtz billows arising in a stratified shear flow. Using buoyancy as one of the tracers, the redistribution of fluid in buoyancy-tracer space was used to interpret the mixing process.

In this paper we use a similar formulation to diagnose the diapycnal transport of a passive tracer in a buoyant plume penetrating a linearly stably stratified layer. The numerical method is detailed in § 2. The evolution of the flow and tracer concentration is presented in § 3. In § 4 we introduce our formulation of the buoyancy-tracer ‘volume distribution’. We show that the flow can be partitioned into three regions of buoyancy-tracer space: the source region where plume fluid enters the stratified layer, a transport region through which volume flows during initial mixing between the plume and environment and an accumulation region where mixed fluid settles and homogenises. Each of these regions of buoyancy-tracer space correspond to coherent regions of physical space that identify the essential structures of the flow, namely the rising plume, plume cap and radially spreading intrusion, respectively. These structures are indicated in figure 1. In § 5 we analyse diagnostics of the mixing process in each of these regions.

Figure 1. Set-up for numerical simulations of a buoyant plume with source integral buoyancy and tracer flux $F_0$ and source radius $r_0$ penetrating from a uniform layer into a linearly stably stratified layer with constant buoyancy frequency $N$. The main plume structures discussed throughout this paper are indicated in red. The pre-penetration region defined in § 3 for diagnostic purposes is shown in grey. The initial buoyancy profile (right) in the stratified environment is $b(\boldsymbol {x},0) = N^2 z$ (dimensional) for $z \geqslant 0$. The maximum penetration height $z_{max}$ and intrusion height $z_n$ above the bottom of the stratified layer are indicated by dotted lines. We also show the forcing region of depth $L_c$ and the forcing modulation profile $f_m(z)$ decaying over a distance $L_p$ in blue, detailed in Appendix A. The (azimuthally averaged) Gaussian profiles of $w, b$ and $\phi$ in the plume rising through the uniform layer are also illustrated in blue.

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 set-up is shown in figure 1. To aid in the examination of the flow evolution and mixing, we include a passive tracer $\phi$ that satisfies the same evolution equation as buoyancy $b = -g\rho '/\rho _0$ where $\rho '$ 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. The scalar field $\phi (\boldsymbol {x}, t)$ represents the (dimensionless) tracer concentration normalised by the tracer concentration on the plume centreline at the source. As illustrated in figure 1, we define the bottom of the initial stratified layer to be $z=0$. We also define $t=0$ as the time at which the plume first penetrates into the stratified layer. The plume source (at the base of the domain) lies at $z=-H \approx -7.97$ for the parameter choices given in table 1. The initial conditions are $\phi (\boldsymbol {x},0) = 0$ throughout the domain whilst $b(\boldsymbol {x},0) = 0$ in the uniform layer $z \leqslant 0$ and $b(\boldsymbol {x},0) = N^2 z$ (dimensional) in the stratified layer $z \geqslant 0$.

Table 1. Non-dimensional parameters for the LES with $N = 1\ \mathrm {s}^{-1}$ and $F_0 = 3.96\times 10^{-7}\ \mathrm {m}^4\ \mathrm {s}^{-3}$ discussed from § 3 onwards.

We generate the buoyant plume by forcing the vertical velocity $w$, buoyancy $b$ and tracer concentration $\phi$ in the shallow forcing region of depth $L_c$ indicated in figure 1. The plume centreline lies at the middle of the computational domain, $x=y=0$. Details on the plume forcing method can be found in Appendix A. At the source, the generated plume has integral buoyancy flux $F_0 = 2\int _0^\infty \bar {w}\bar {b}|_{z=-H}r\,\mathrm {d}r$, where $\bar {\cdot }$ denotes an azimuthal and time average, with no excess momentum flux (i.e. a ‘pure’ plume, see Appendix A for detail). Turbulence is initiated in the plume by applying a 10 % perturbation to the forcing profiles in the forcing region and to all velocity components in the two grid layers above the forcing region. Turbulence develops as the plume rises through the uniform layer and we ensure that, prior to penetrating the stratified layer, the azimuthally averaged vertical velocity, buoyancy and tracer concentration are self-similar with a Gaussian radial profile as expected in a fully developed plume (see Appendices A and B).

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$. The length scale $L$ naturally arises from the Morton, Taylor & Turner (1956) plume equations in a stably stratified environment. Following previous experimental and numerical studies (e.g. Briggs 1965; Devenish, Rooney & Thomson 2010), both the maximum height of the plume $z_{max}$ and the height of the intrusion $z_n$ above the base of the stratification (illustrated in figure 1) scale with $L$ (in the case of a ‘lazy’ plume, when incident momentum is negligible, as considered here).

Owing to the large range of length scales involved in convective penetration, resolving turbulent scales with direct numerical simulation is not feasible due to the computational cost of very high resolution simulations. We therefore use large-eddy simulation (LES) with the anisotropic minimum dissipation (AMD) eddy-viscosity model to represent unresolved scales (Vreugdenhil & Taylor 2018). Large-eddy simulation has been shown to be effective for simulating plumes in previous work in the literature, e.g. Pham, Plourde & Doan (2007). The non-dimensional governing equations for velocity $\boldsymbol {u}$ and scalars $b, \phi$ including sub-grid-scale (SGS) contributions are

(2.1)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\hat{\boldsymbol{u}} = 0, \end{gather}

(2.2)\begin{gather}\frac{\mathrm{D} \hat{\boldsymbol{u}}}{\mathrm{D}t} + \boldsymbol{\nabla} \hat{p} = \frac{1}{\textit{Re}} \nabla^2 \hat{\boldsymbol{u}} + \hat{b} \hat{\boldsymbol{k}} - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau} + f_w \hat{\boldsymbol{k}}, \end{gather}

(2.3)\begin{gather}\frac{\mathrm{D}\hat{b}}{\mathrm{D}t} = \frac{1}{\textit{Re}\textit{Pr}} \nabla^2 \hat{b} - \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\lambda}_b + f_b, \end{gather}

(2.4)\begin{gather}\frac{\mathrm{D}\hat{\phi}}{\mathrm{D}t} = \frac{1}{\textit{Re}\textit{Pr}} \nabla^2 \hat{\phi} - \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\lambda}_\phi + f_{\phi}, \end{gather}

where $\hat {\cdot }$ indicates filtering at the resolved grid scale and $\hat {\boldsymbol {k}}$ is the unit vector in the vertical direction. The terms $f_w$, $f_b$ and $f_\phi$ represent the forcing applied to the vertical velocity, buoyancy and passive tracer to generate the buoyant plume. The details of this forcing are discussed in Appendix A. The SGS stress tensor $\boldsymbol {\tau }$ has components $\tau _{ij} = \widehat {u_i u_j} - \hat {u}_i \hat {u}_j$, the SGS buoyancy flux is $\boldsymbol {\lambda }_b =\widehat {\boldsymbol {u}b} - \hat {\boldsymbol {u}} \hat {b}$ and similarly the SGS tracer flux is $\boldsymbol {\lambda }_{\phi } =\widehat {\boldsymbol {u}\phi } - \hat {\boldsymbol {u}} \hat {\phi }$. The two dimensionless parameters are the Reynolds number and Prandtl number

(2.5a,b)\begin{equation} \textit{Re} = \frac{F_0^{1/2}}{\nu N^{1/2}}, \quad \textit{Pr} = \frac{\nu}{\kappa}, \end{equation}

respectively, where $\nu$ is the molecular viscosity and $\kappa$ is the molecular diffusivity for both $b$ and $\phi$. The eddy-viscosity model for the deviatoric component of the SGS stress $\boldsymbol {\tau }^d$ and the SGS buoyancy and tracer flux are

(2.6a–c)\begin{equation} \tau_{ij}^d = \tau_{ij} - \tfrac{1}{3}\delta_{ij}\tau_{kk} ={-}2 \nu_{SGS} \hat{{\mathsf{S}}}_{ij}, \quad\boldsymbol{\lambda}_b ={-}\kappa^{(b)}_{SGS} \boldsymbol{\nabla} \hat{b}, \quad \boldsymbol{\lambda}_\phi ={-}\kappa^{(\phi)}_{SGS} \boldsymbol{\nabla}\hat{\phi}, \end{equation}

where $\nu _{SGS}, \kappa ^{(b)}_{SGS}$ and $\kappa ^{(\phi )}_{SGS}$ are the non-dimensional SGS viscosity, SGS buoyancy diffusivity and SGS tracer diffusivity, respectively, each determined by the AMD scheme. The term $\hat {{\mathsf{S}}}_{ij}$ is the components of the non-dimensional shear-rate tensor. The SGS diffusivities and viscosity may locally exceed the molecular values by several orders of magnitude in regions with intense turbulence.

We use DIABLO (Taylor 2008) to perform three-dimensional LES of the idealised set-up shown in figure 1. DIABLO evolves the Boussinesq Navier–Stokes equations (2.1)–(2.4) discretised using Fourier modes in the two periodic horizontal directions and second-order finite differences in the vertical direction. The boundary conditions on the top and bottom boundary are $\partial _z u = \partial _z v = \partial _z b = \partial _z \phi = 0$ and $w = 0$. A third-order Runge–Kutta scheme is used for time stepping. A $2/3$ dealiasing rule is applied when transforming from Fourier to physical space. We use a cubic domain of side length $L_{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 (which the simulated plume does not reach), where the velocity is damped towards zero and the buoyancy is damped towards the initial background stratification $b(\boldsymbol {x},0) = z$ (non-dimensional), to inhibit the reflection of internal gravity waves from the top boundary. Validation of the numerical method discussed here is detailed in Appendix B.

Henceforth we refer to a single simulation with parameters which are equivalent to the experimental set-up 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 = 3.96\times 10^{-7}\ \mathrm {m}^4\ \mathrm {s}^{-3}$. Henceforth, all values stated are non-dimensionalised with respect to this choice of $F_0$ and $N$. We also drop the hat notation and refer to the resolved variables unless otherwise noted.

3. Flow and tracer structure

The flow evolution is presented in three vertical cross-sections through the plume centreline in figure 2. We identify the plume as regions with tracer concentration $\phi \geqslant \phi _{min} \equiv 10^{-2}$, i.e. we threshold the tracer field at $1\,\%$ of its value on the plume centreline at the source. In the tracer-less environment surrounding the plume we show contours of the buoyancy field. The bottom of the stratified layer, above which the buoyancy of the environment becomes non-zero, is indicated by the lowest buoyancy contour.

Figure 2. Three stages of the flow evolution, shown as $x-z$ cross-sections of the tracer concentration $\phi$ shown where $\phi$ exceeds $1\,\%$ of its value on the plume centreline at the source $z=-H$. Buoyancy contours are shown otherwise. Cross-sections are taken at the plume centreline at non-dimensional times $t = 1, 6.75, 14$. From left to right, the panels show the plume during initial penetration, reaching maximum penetration height and spreading of the intrusion.

Figure 2(a) shows initial penetration of the stratified layer by the plume cap. As the plume rises through the stratified layer, its upward acceleration decreases as the relative buoyancy between the plume and the surrounding environment decreases. Once the environmental buoyancy exceeds that of the plume, the plume decelerates. Eventually, the rising fluid reverses direction, or ‘overturns’, and begins to subside from the maximum penetration height $z_{max}$ (figure 2b). As plume fluid subsides, its buoyancy relative to the surrounding environment increases until reaching the level of neutral buoyancy $z_n$ where the plume fluid forms a radially spreading intrusion – see figure 2(c). The dynamics observed in the simulation agrees qualitatively with studies of similar set-ups in the literature, for example the experiments detailed in Ansong & Sutherland (2010) with an identical set-up and similar physical parameters.

The evolution of the maximum height of the plume cap during penetration and the subsequent quasi-steady state is visualised as a time series of tracer concentration on the plume centreline in figure 3. As has been noted in the literature, the maximum height of the plume tends to oscillate around a quasi-steady state height $z_{ss}$ (Turner 1966) but, to our knowledge, the mechanism setting the frequency of this oscillation (often referred to as ‘plume bobbing’) is not well understood (e.g. Ansong & Sutherland 2010). In the simulation considered here, the quasi-steady state height $z_{ss}$ is close to the maximum penetration height $z_{max}$ and the oscillation is weak. For convenience, we will use $z_{max}$ to refer to the maximum height of the plume. The maximum penetration height $z_{max}$ determines the maximum height at which plume fluid can mix with the environment (Ansong, Kyba & Sutherland 2008), meaning the initial buoyancy at the maximum penetration height, $b = z_{max}$, represents a plausible constraint on the maximum buoyancy accessible for mixing with the plume. However, this constraint can occasionally be exceeded when plume fluid subsiding from the plume cap pulls very buoyant environmental fluid downwards (see figure 2(b) to the left of the plume cap). Here we find $z_{max} = 3.94$ which agrees with experimental estimates of the maximum penetration height in the literature, e.g. $z_{max} \approx 3.8$ for a plume with a round source (List 1982).

Figure 3. Time series of the tracer concentration $\phi (0,0,z,t)$ in the $z\unicode{x2013}t$ plane at the centreline of the computational domain $x=y=0$. The green dashed contour denotes the plume threshold $\phi = \phi _{min}$, i.e. where $\phi$ is $1\,\%$ of its value on the plume centreline. The maximum penetration height $z_{max}$ and the quasi-steady state height $z_{ss}$ are marked.

Internal gravity waves across a range of frequencies are generated during the penetration process. These waves are visible as small amplitude, long wavelength undulations in the buoyancy contours above $z \approx 4$ in figure 2(c). Any influence of internal gravity waves on mixing in this flow will be present in the analyses, but it is beyond the scope of this paper to determine the particular contribution of these waves to mixing.

In the uniform layer, the buoyancy and tracer evolve identically up to a linear factor, i.e. the undiluted plume fluid entering the stratified layer has a linear relationship between $b$ and $\phi$ at each point. This follows from the self-similarity of the buoyancy and tracer concentration profiles in the steady state plume that penetrates the stratified layer (see Appendix B, figure 18). The radial profiles for $b$ and $\phi$ are both Gaussian with the same width but different amplitudes, hence $b \propto \phi$. After penetrating the stratified layer, plume fluid with non-zero buoyancy and tracer concentration mixes with tracer-less environmental fluid and hence the buoyancy and tracer evolve differently. This effect can be quantified using a tracer probability density function (p.d.f.) in buoyancy coordinates $\tilde {\phi }(b; t)$. The p.d.f. is calculated within the stratified layer only. The value of the p.d.f. $\tilde {\phi }(b; t)$ is calculated as the total tracer with buoyancy within a range $b$ to $b + \mathrm {d}b$ in the stratified layer, normalised by the total tracer in the stratified layer $\phi _T(t) = \sum _V \phi (\boldsymbol {x},t)\Delta V$, where $V$ is the stratified layer and $\Delta V$ is the grid-cell volume. The definition of $\tilde {\phi }$ is such that $\sum _B \tilde {\phi }(B; t) = 1$.

Figure 4 shows $\tilde {\phi }(b;t)$ in the stratified layer at fixed time intervals post-penetration. The total tracer in the stratified layer $\phi _T(t)$ is shown inset. The approximately linear increase in $\phi _T$ with time suggests a relatively uniform input of tracer to the stratified layer, carried by the penetrating plume. Owing to the self-similar nature of the penetrating plume, we expect the tracer that enters the stratified layer to have a fixed p.d.f. (with some small variation). This pre-penetration p.d.f. $\tilde {\phi }_0$ can be estimated using a domain $V$ chosen as the pre-penetration region shown in figure 1. This region is a thin layer with (non-dimensional) depth $1/2$ below the bottom of the stratified layer. The pre-penetration p.d.f., shown as a black dashed line in figure 4, represents the tracer p.d.f. in the plume just before it penetrates the stratified layer. Without mixing, $\tilde {\phi }$ in the stratified layer would match the pre-penetration p.d.f. Mixing during the penetration process manifests as changes in the tracer p.d.f. when compared with the pre-penetration p.d.f.

Figure 4. Probability density function $\tilde {\phi }(b, t)$ of tracer as a function of buoyancy $b$ in the stratified layer at fixed time intervals post-penetration shown as coloured lines. The black dashed line shows the time-averaged pre-penetration p.d.f. $\tilde {\phi }_0(b)$, calculated with $V$ chosen as the pre-penetration region indicated in figure 1 and time averaged. The pre-penetration p.d.f. $\tilde {\phi }_0(b)$ shows the tracer-buoyancy relationship within the plume prior to penetrating the stratified layer. Differences between $\tilde {\phi }_0$ and $\tilde {\phi }(b,t)$ represent the effect of mixing. Total tracer in the stratified layer $\phi _T(t)$ shown inset.

Evolution of the post-penetration p.d.f. and changes compared with the pre-penetration tracer p.d.f. highlight two mixing processes during penetration: mixing within the plume during penetration, and mixing between the plume and environment. Where plume fluid carries a large tracer concentration and mixes with the more buoyant surroundings, the positive tail of the tracer p.d.f. increases. This is particularly evident after $t=6.75$ when the plume has reached $z_{max}=3.94$, at which point very large values of buoyancy in the environment become accessible and large tracer concentrations at the centre of the plume are exposed to the environment as plume fluid overturns. At the edges of the plume where tracer concentration is smallest, mixing with the environment again moves tracer from lower to higher values of buoyancy and therefore the p.d.f. decreases where $b$ is small. This effect is supplemented by mixing within the plume, which acts to homogenise the large tracer concentration and buoyancy at the centre of the plume with the lower tracer concentration and buoyancy at the edge of the plume. This acts to narrow the p.d.f. and hence decrease the p.d.f. at large and small values of buoyancy but the effect is only evident before the plume reaches $z_{max}$ at $t=6.75$ and accesses much larger values of buoyancy. 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 p.d.f.

The buoyancy range of the tracer p.d.f. is determined by the maximum penetration height of the plume as well as the rapidity of the mixing between the plume and environment occurring in the plume cap. If fluid quickly reaches $z_{max}$ and subsides before substantial mixing with the environment occurs, only small amounts of the more buoyant environment are entrained and therefore the increase in the p.d.f. at large values of buoyancy is modest compared with a scenario where plume fluid stalls during overturning and significant mixing with the environment occurs. In figure 4, the tracer p.d.f. extent is $b \approx 3$ whilst the environmental buoyancy at $z_{max}$ is approximately $b |_{z_{max}} \approx 3.94$. This suggests the mixing time scale is slow compared with the dynamical time scale, i.e. mixing between the largest tracer concentrations first exposed during overturning and the environment is slow and continues during subsidence, where the buoyancy of the environment decreases.

The tracer p.d.f. hints at competing effects of mixing within the plume and between the plume and the environment. Crucially, the buoyancy and tracer fields are mixed in different ways owing to the linearly increasing buoyancy and vanishing tracer concentration in the linearly stratified environment. Whilst changes in the tracer p.d.f. considered here demonstrate the overall effect on the relation between tracer and buoyancy, it is difficult to extract information on the intensity of mixing between plume and environmental fluid and the specific buoyancy and tracer characteristics of the fluid parcels that mix. Furthermore, the tracer p.d.f. $\tilde {\phi }(b;t)$ does not give information on the volume of fluid parcels at a given buoyancy; a peak in the tracer p.d.f. may represent relatively few fluid parcels carrying large tracer concentrations or 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 diapycnal transport of tracer.

4. Buoyancy-tracer volume distribution

The probability distributions of tracer concentration discussed in § 3 isolate the irreversible transport that results from turbulent mixing. 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. Whilst stirring strengthens tracer gradients across buoyancy surfaces, it is – in principle – a reversible process. However, molecular diffusion results in irreversible changes to the buoyancy and tracer characteristics of fluid parcels and hence changes the tracer distribution.

Here, we use the distribution of volume in buoyancy-tracer space to diagnose mixing in the stratified layer. That is, we map from three-dimensional physical space to a two-dimensional phase space by using the buoyancy and tracer concentration fields to quantify the volume of plume fluid in the stratified layer with each value of $b$ and $\phi$. The total physical volume of plume fluid represented in the distribution changes in time and we do not normalise the distribution to form a p.d.f. Omitting this normalisation simplifies the interpretation of the distribution and its governing equation. The buoyancy-tracer volume distribution formalism presented here builds on previous density–tracer joint p.d.f. formulations presented by Plumb (2007) and Penney et al. (2020).

4.1. Definition and properties

We define the volume distribution $W(B, \varPhi ; t)$ in buoyancy-tracer space such that the volume of fluid in a fixed volume $V$ with $B < b(\boldsymbol {x}, t) < B + \mathrm {d}B$ and $\varPhi < \phi (\boldsymbol {x}, t) < \varPhi + \mathrm {d}\varPhi$ is given by $W(B, \varPhi ; t)\,\mathrm {d}B\,\mathrm {d}\varPhi$. This may be defined as

(4.1)\begin{equation} W(B, \varPhi; t) = \int_V \delta(b(\boldsymbol{x}, t) - B) \delta(\phi(\boldsymbol{x}, t) - \varPhi) \, \mathrm{d}V, \end{equation}

where $\delta ({\cdot })$ is the Dirac delta function with the inverse dimension of its argument. Henceforth, we choose the volume $V$ to be the stratified layer. An evolution equation for $W$ can be obtained using the governing equations for $b$ and $\phi$. See Appendix C for a full derivation. We have

(4.2)\begin{equation} \frac{\partial W}{\partial t} ={-} \boldsymbol{\nabla}_{(B, \varPhi)} \boldsymbol{\cdot} \boldsymbol{F} + S, \end{equation}

where $\boldsymbol {F}(B, \varPhi ; t)$ is the mixing flux distribution and $S(B, \varPhi ; t)$ is the source distribution. The mixing flux distribution $\boldsymbol {F}$ is a vector in buoyancy-tracer space with components formed from the volume-weighted average of the non-advective terms $\dot {b}$ and $\dot {\phi }$ in (2.3) and (2.4) respectively, representing the flux of $W$ in buoyancy-tracer space due to mixing and is defined as

(4.3)\begin{equation} \boldsymbol{F}(B, \varPhi; t) = (F_b, F_\phi) = \int_V (\dot{b}, \dot{\phi})\delta(b(\boldsymbol{x}, t) - B) \delta(\phi(\boldsymbol{x}, t) - \varPhi) \,\mathrm{d}V, \end{equation}

where $\dot {b} = (\textit {Re} \textit {Pr})^{-1} \nabla ^2 b - \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\lambda }_b$ and $\dot {\phi } = (\textit {Re}\textit {Pr})^{-1}\nabla ^2 \phi - \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\lambda }_\phi$. Note that the plume forcing terms $f_b$ and $f_\phi$ are excluded from $\dot {b}$ and $\dot {\phi }$ since the forcing vanishes in the stratified layer. The source distribution $S$ represents a source or sink of $W$ due to boundary fluxes across $\partial V$

(4.4)\begin{equation} S(B, \varPhi; t) = \int_{\partial V} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n}\, \delta(b(\boldsymbol{x}, t) - B) \delta(\phi(\boldsymbol{x}, t) - \varPhi) \, \mathrm{d}A, \end{equation}

where $\boldsymbol {u}$ is the velocity in physical space and $\boldsymbol {n}$ is the inward normal on the boundary $\partial V$ of $V$. Since we are considering a flow upwards into $V$, $\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {n}$ is positive and $S$ acts as a source of $W$. Note that whilst $S$ represents the effect of fluxes across the boundary $\partial V$ in physical space, it is distributed in buoyancy-tracer space. Note that (4.2) 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.

Turbulent mixing redistributes volume in buoyancy-tracer space, which results in changes to $W$ via the mixing flux term $-\boldsymbol {\nabla }_{(B,\varPhi )} \boldsymbol {\cdot }\boldsymbol {F}$. The change in $W$ at a point $(B, \varPhi )$ in buoyancy-tracer space as a result of turbulent mixing up to time $t$ is therefore

(4.5)\begin{equation} M(B, \varPhi; t) ={-}\int_0^t \boldsymbol{\nabla}_{(B, \varPhi)} \boldsymbol{\cdot} \boldsymbol{F}(B, \varPhi; t')\, \mathrm{d}t' = W(B, \varPhi; t) - \int_0^t S(B, \varPhi; t') \, \mathrm{d}t', \end{equation}

such that $M(B, \varPhi ; t)\,\mathrm {d}B \,\mathrm {d}\varPhi$ is the change in volume of fluid with $B < b(\boldsymbol {x}, t) < B + \mathrm {d}B$ and $\varPhi < \phi (\boldsymbol {x}, t) < \varPhi + \mathrm {d}\varPhi$ up to time $t$ due to mixing. Therefore $M$ represents the integrated effect of the mixing flux $\boldsymbol {F}$ and we refer to $M$ as the net mixing effect distribution. The second equality in (4.5) follows from time integrating (4.2) and noting that $W(B, \varPhi ; t=0) = 0$ since there is no tracer in the initial stratified layer. Hence, $M$ can be interpreted as a cumulative measure of the changes to $W$ relative to the time-integrated source distribution, i.e. the changes in the volume distribution that arise solely from mixing. The final term in (4.5), which we refer to as the cumulative source distribution, represents the volume of fluid with buoyancy $B < b < B + \mathrm {d}B$ and tracer concentration $\varPhi < \phi < \varPhi + \mathrm {d}\varPhi$ that has entered the stratified layer up to time $t$. The volume distribution $W \geqslant 0$ and the cumulative source distribution is also positive assuming there is a flow into $V$ only. However, $M$ can be positive or negative depending on the relative sizes of the volume distribution and the cumulative source distribution.

The net mixing effect distribution $M(B, \varPhi ; t)$ is positive in buoyancy-tracer space where more volume is present at time $t$ than has entered the stratified layer up to time $t$, i.e. there is a net gain in the volume of fluid with buoyancy $B$ and tracer concentration $\varPhi$ due to mixing. Correspondingly, $M(B, \varPhi ; t)$ is negative where more volume has entered the stratified layer up to time $t$ with buoyancy $B$ and tracer concentration $\varPhi$ than currently exists at time $t$, i.e. there is a net loss in the volume of fluid with buoyancy $B$ and tracer concentration $\varPhi$ due to mixing. The value of $M$ therefore indicates the transfer of volume within $W$ due to mixing; fluid leaves regions of buoyancy-tracer space with $M < 0$ and enters regions with $M > 0$.

To summarise, the distributions $W, S, \boldsymbol {F}$ and $M$ together describe the flow in terms of its effect on buoyancy-tracer space. The volume distribution $W$ is an instantaneous representation of the amount of fluid within the stratified layer with given ranges of values of buoyancy and tracer concentration. Large values of $W$ indicate large volumes of fluid with a narrow range of $b$ and $\phi$, though the fluid parcels corresponding to this range are not necessarily co-located in physical space. The source distribution $S$ represents the volume distribution of fluid that enters the stratified layer from the uniform layer. In the absence of mixing, $W$ would be equivalent to the time integral of $S$. The mixing flux distribution $\boldsymbol {F}$ represents the redistribution of fluid in buoyancy-tracer space due to mixing. The net mixing flux distribution $M$ captures the change in $W$ relative to time-integrated $S$ via $\boldsymbol {F}$ and indicates where there is accumulation or loss of volume due to mixing.

4.2. Idealised example

The effect of an idealised turbulent mixing event on buoyancy-tracer space is illustrated in figure 5. The top row shows the volume distribution containing three initial fluid parcels (blue points) which have entered the stratified layer. As turbulent stirring brings these fluid parcels together, another fluid parcel (green point) enters the domain and all four fluid parcels mix. The resulting mixed fluid parcel (red point) is a volume-weighted average of the fluid parcels involved in the mixing event. The bottom row shows the distributions after the mixing event. The volume distribution $W$ is non-zero (and positive) only where the final mixture lies in buoyancy-tracer space and the direction of the mixing flux vectors $\boldsymbol {F}$ indicates the redistribution of volume. The cumulative source distribution $\int _0^t S \, \mathrm {d}t'$ is positive at the values of $b$ and $\phi$ where the fluid parcels entered the domain and vanishes elsewhere. The net mixing effect distribution $M$ is negative at these points, as volume has been lost, and positive where the mixed fluid parcel lies as volume has been gained. These principles can be used to understand the mixing processes in the physical flow that result in changes in the distribution in buoyancy-tracer space. Whilst it is not possible to isolate distinct fluid parcels that are mixing at any one time, we can identify physical regions of the flow that are subject to turbulent mixing and isolate the corresponding regions of buoyancy-tracer space.

Figure 5. Schematic diagram of the effect in buoyancy-tracer space of an idealised turbulent mixing event between a set of discrete fluid parcels that have entered the stratified layer and a fluid parcel that later enters the stratified layer. The top row illustrates the convex envelope property of the volume distribution $W$, which implies that a mixture of a set of fluid parcels lies within the smallest envelope that contains the distribution of the fluid parcels that are mixed together. The distributions $W, \int S\, \mathrm {d}t$ and $M$ following the idealised mixing event are shown on the bottom row. Positive (arbitrary) values of each distribution are indicated by a circled $+$, negative (arbitrary) values are indicated by a circled $-$. The mixing flux distribution vectors $\boldsymbol {F}$ are indicated by grey arrows.

Turbulent mixing acts to homogenise the buoyancy and tracer concentration of fluid parcels. 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 (Penney et al. 2020). Therefore, as the buoyancy-tracer volume distribution $W$ evolves, it is constrained to lie within its own past convex envelope, i.e. the smallest convex set that contains all non-zero points of the distribution. As illustrated in figure 5, this convex envelope must include fluid that enters the domain during the mixing process. The convex envelope of the initial volume distribution is indicated by the dotted envelope and the dashed envelope indicates the convex envelope including newly arriving fluid parcels. We emphasise that the final mixed fluid parcel 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 principle of homogenising fluid parcels illustrated in figure 5 can be generalised to continuous mixing of fluid in a flow, in which case the convex envelope constraint applies to the volume distribution as a whole. This implies that in the absence of sources the convex envelope reduces over time and converges towards some compact distribution (Penney et al. 2020). In the set-up we consider, fluid entering the stratified layer causes the extremes of the distribution to persist, whilst turbulent mixing acts to continuously shift the buoyancy-tracer characteristics of fluid towards an accumulation region in buoyancy-tracer space.

4.3. Numerical implementation

To examine the numerical simulation detailed in §§ 2 and 3, we use a discrete formulation of the buoyancy-tracer volume distribution introduced in § 4.1. We choose the domain $V$ to be tracer-containing plume fluid within the stratified layer. The stratified layer initially corresponds to the volume $z \geqslant 0$. However, the plume can perturb the bottom of the stratified layer slightly below $z=0$. We therefore define the domain as the region where $-1 \leqslant z \leqslant L_z$, $\phi > \phi _{min}$, and $b > 0$. As a consequence, the ‘reservoir’ of environmental fluid where $\phi = 0$ is excluded. In interpreting results, we therefore consider the boundary $\phi = \phi _{min}$ as a source where volume can enter the distribution from the environment. The entrainment of environmental fluid across this boundary in buoyancy-tracer space into the volume distribution is discussed in § 4.7 and illustrated in figure 7.

The buoyancy and tracer domains are subdivided into $N_b$ and $N_{\phi }$ equally sized bins of size $\delta B = (b_{max} - b_{min})/N_b$ and $\delta \varPhi = (\phi _{max} - \phi _{min})/N_\phi$ respectively. We choose $b_{min} = 0$ and $b_{max} = 4$ since the largest accessible buoyancy is related to the maximum penetration height (Ansong et al. 2008) which is experimentally predicted to be $z_{max} \leqslant 4$ (List 1982). To accommodate fluctuations in tracer concentration, we choose $\phi _{max}$ to be larger than the tracer concentration on the plume centreline $2\phi _m(0)$ at penetration height $z=0$ using the profile predicted by the Morton et al. (1956) plume equations, $\phi _m(z)$, defined in (A5) (see Appendix A for details). We use $\phi _{min} = 10^{-2}$, consistent with the plume threshold introduced in § 3. Henceforth we use $N_b = N_\phi = 256$.

Denoting the centre of a given bin as $(B_i, \varPhi _j)$, the associated value of the volume distribution is computed as

(4.6)\begin{equation} W_{ij}(t) = \sum_V I_{ij}(\boldsymbol{x},t) \Delta x \Delta y \Delta z , \end{equation}

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

(4.7)\begin{equation} I_{ij}(\boldsymbol{x}, t) = \begin{cases} 1 & \left(b(\boldsymbol{x},t) - B_i, \phi(\boldsymbol{x},t) - \varPhi_j\right) \in \left( -\tfrac{1}{2}\delta B, \tfrac{1}{2}\delta B \right] \times \left( -\tfrac{1}{2}\delta \varPhi, \tfrac{1}{2}\delta \varPhi \right], \\ 0 & \text{otherwise}. \end{cases} \end{equation}

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 \varPhi /2$ of $\varPhi _j$. Note that in the continuous formulation, the volume distribution $W(B, \varPhi ; t)$ defined by (4.1) must be integrated over $B$ and $\varPhi$ to yield a volume, whilst $W_{ij}(t)$ itself has dimensions of volume and need only be summed over $i$ and $j$. The continuous and discrete formulations coincide in the limit $\delta B, \delta \varPhi \to 0$, such that

(4.8)\begin{equation} \lim_{\delta B, \delta \varPhi \to 0} \frac{W_{ij}(t)}{\delta B \delta \varPhi} = W(B_i, \varPhi_j; t). \end{equation}

The equivalence (4.8) also applies to the discrete mixing flux distribution $\boldsymbol {F}_{ij}(t)$, the discrete source distribution $S_{ij}(t)$ and the discrete net mixing effect distribution $M_{ij}(t)$ defined by

(4.9)\begin{gather} \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, \end{gather}

(4.10)\begin{gather}S_{ij}(t) = \sum_{\partial V} I_{ij}(\boldsymbol{x}|_{z={-}1}, t) w(\boldsymbol{x}|_{z={-}1}, t) \Delta x \Delta y, \end{gather}

(4.11)\begin{gather}M_{ij}(t) = W_{ij}(t) - \sum_{t'} S_{ij}(t') \Delta t', \end{gather}

where $\Delta t'$ is the simulation time step. In (4.9), $\dot {b}$ and $\dot {\phi }$ are the non-advective terms in the scalar evolution equations (2.3), (2.4) of $b, \phi$ respectively, as defined in § 4.1. In (4.10) we have used the fact that $\boldsymbol {n} = \hat {\boldsymbol {k}}$ on the bottom boundary of the domain $V$.

4.4. Results

The discrete formulation of the distributions given in § 4.3 provides an approximation to the continuous formulation and is presented in all figures shown below. However, the interpretation is the same as the continuous formulation and we will refer to the continuous formulation in all discussions. Quantities derived from the distributions are given in both continuous and discrete forms for completeness. In defining the discrete and continuous formulations we use the arguments $B$ and $\varPhi$, which represent values of buoyancy and tracer concentration respectively. We treat $W$, $\boldsymbol {F}$, $S$ and $M$ as functions of $b$ and $\phi$ to aid clarity, e.g. $W(b, \phi ; t)$, with the interpretation that $b$ and $\phi$ represent values of buoyancy and tracer concentration found in the flow in the same way as $B$ and $\varPhi$ in §§ 4.1 and 4.3.

Figure 6 shows the buoyancy-tracer volume distribution $W(b, \phi ; t)$ (middle row), the source distribution $S(b, \phi ; t)$ (bottom row), and $x-z$ cross-sections of the tracer concentration field and buoyancy contours (top row). These results are shown at three snapshots corresponding with stages of the flow evolution as in figure 2. The distributions are shown only where non-zero, i.e. regions of buoyancy-tracer space which are not coloured indicate that there is no fluid with buoyancy and tracer concentration in that range. In each snapshot of $W$, the red dashed lines show the convex envelope that constrains the evolution of the volume distribution. As seen in the figure, the source distribution lies within the convex envelope of $W(b, \phi ; t)$. Furthermore, as the plume rises and accesses more buoyant fluid in the surrounding environment, the convex envelope is extended along the $\phi =0$ axis as new environmental fluid becomes accessible via mixing.

Figure 6. Three instantaneous snapshots showing the evolution of the buoyancy-tracer volume distribution $W(b, \phi ; t)$ (middle) and source distribution $S(b, \phi ; t)$ (bottom) at non-dimensional times $t=1, 6.75, 14$ corresponding with figure 2. The convex envelope of the volume distribution $W$ at time $t$ is shown as a red dashed line in the middle panel. To aid interpretation, we also show $x-z$ cross-sections of the tracer concentration and buoyancy contours, as in figure 3 (top).

The results shown in figure 6 illustrate how the volume distribution captures the dynamics and mixing processes at each stage of the flow evolution. We first note that the source distribution $S(b, \phi ; t)$ takes positive values only, since there is only a flow into the stratified layer. Furthermore, $S$ is non-zero only on a line through the origin as expected from the linear relationship between $b$ and $\phi$ in the rising plume. We refer to this as the source line. Given that the convex envelope of a set of points on a line segment is the same line segment, mixing of undiluted plume fluid within the plume only redistributes fluid on the source line. When undiluted plume fluid mixes with the surrounding environment, it is moved away from the source line. This offers a clear distinction between undiluted and mixed plume fluid, as illustrated schematically in figure 7. In the buoyancy-tracer volume distribution $W$ shown on the middle row of figure 6, fluid appearing away from the source line therefore represents a mixture of plume and environmental fluid. Further information on the regions of the undiluted plume that mix with the environment is gained by noting that, owing to the Gaussian profiles of the plume pre-penetration, $b$ and $\phi$ are larger near the centreline of the plume and smaller towards the edge of the plume (see figure 6c). Hence fluid near the ‘plume edge’ lies nearest the origin on the source line whilst fluid in the ‘plume core’ lies at the extreme end of the source line.

Figure 7. Schematic diagram illustrating the correspondence between regions in physical space and regions in buoyancy-tracer space. Environmental fluid is represented by the $b$ axis where $\phi =0$ (red), between the bottom of the stratified layer where $b=0$ and the maximum penetration height where $b=z_{max}$. Undiluted plume fluid lies on the source line where $\phi \propto b$ (blue), with large $b$ and $\phi$ in the core of the plume and small $b$ and $\phi$ towards the edges. Mixed fluid lies between these two lines, within their convex envelope (dotted grey). Some regions of the convex envelope (hatched area) are inaccessible due to the shielding effect of the plume edge (where $\phi$ is non-zero but small) and the intrusion surrounding the rising plume. Volume entering the distribution due to entrainment of environmental fluid is indicated by the red arrows.

Figure 6(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 part of the 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 6(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}$ with relatively large buoyancy. The volume distribution at this time covers a wider range of $b$ and $\phi$ as environmental fluid with a large range of buoyancy is able to mix with much of the undiluted plume fluid lying on the source line. However, note that we do not see mixing between the extreme end of the source line and the environment (see hatched region in figure 7). As undiluted plume fluid in the core of the plume rises, turbulent mixing acts to homogenise the buoyancy and tracer concentration within the source distribution, reducing the largest values of $b$ and $\phi$ near the centreline of the plume. Hence, when this fluid is first exposed to the environment in the plume cap, it lies closer to the middle of the source line. We refer to this as the shielding effect. Note that the large values of $b$ and $\phi$ on the source line appear to persist, and may increase, due to new undiluted plume fluid entering the stratified layer.

Figure 6(c) shows the quasi-steady state plume where there is a continuous input of undiluted plume fluid along the source line $\phi \propto b$, mixing between the plume and environment up to values of buoyancy corresponding to heights near $z_{max}$ and an accumulation of fluid at lower values of buoyancy and tracer concentration corresponding to the 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 values of buoyancy and tracer concentration. Note that the intrusion enhances the shielding effect by preventing undiluted plume fluid from accessing the environment before reaching the plume cap.

4.5. Quasi-steady state

In this flow, quasi-steady state refers to the long-term behaviour established once an intrusion has formed. In this state, 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 full plume volume (i.e. the sum of the rising undiluted plume, plume cap and intrusion), and in particular the volume of the intrusion, increases monotonically.

To quantitatively identify the time at which quasi-steady state (QSS) starts, first we define the source line more generally as the region where the cumulative source volume is positive, $\mathcal {S} = \{ (b, \phi ) \vert \int _0^t S(b, \phi ; t') \,\mathrm {d}t' > 0\}$. Next, we can define the volume associated with any region $\mathcal {R}$ of buoyancy-tracer space at time $t$ as

(4.12)\begin{equation} V(\mathcal{R}; t) = \int_{\mathcal{R}} W(b, \phi; t) \, \mathrm{d}b \,\mathrm{d}\phi = \sum_\mathcal{R} W_{ij} (t). \end{equation}

We expect that in QSS $\mathrm {d}V(\mathcal {S})/\mathrm {d}t \approx 0$. However, ‘plume bobbing’ (i.e. the up-and-down motion of the plume cap noted in § 3) results in some variation of the volume of the source line. The quasi-steady volume of undiluted plume fluid $V(\mathcal {S})$ can also gradually increase over time owing to the shielding effect. We therefore introduce an alternative constraint for identifying QSS which utilises the net mixing effect distribution $M$. As illustrated in figure 7, the source line represents undiluted plume fluid arriving in the stratified layer, which introduces volume into the distribution that is eventually mixed away from the source line. In the transient penetration stage, turbulent mixing redistributes fluid on the source line before mixing with the environment. Hence, there is some accumulation on parts of the source line and $M>0$. However, once the plume reaches QSS and mixing with the environment continuously removes volume from the source line, $M$ must become negative. Away from the source line, $S$ vanishes so $M$ is necessarily positive according to (4.5). We define the region $\mathcal {U} = \{ (b, \phi ) \vert M(b, \phi ; t) < 0\}$ and identify QSS as the time when the volume associated with $\mathcal {U}$, $V(\mathcal {U})$, is within $10\,\%$ of the volume of the source line $V(\mathcal {S})$. These volumes and the time we identify as the start of QSS, $t_{QSS} \approx 3.5$, are shown in figure 8(a).

Figure 8. (a) Volume of the source line $\mathcal {S} = \{ \int _0^t S(b, \phi ; t') \, \mathrm {d}t' > 0\}$ (dashed line) and the region $\mathcal {U} = \{ M(b, \phi ; t) < 0\}$ (solid line). The start of QSS (vertical dotted line) is identified as the time when these volumes agree to within $10\,\%$. (b) Decomposition of the full plume volume into undiluted plume fluid, plume cap and intrusion (solid coloured lines) using the partitioning introduced in § 4.6. The total plume input volume up to time $t$ (black dashed line) is shown for reference.

4.6. Source, transport and accumulation regions

We now restrict attention to QSS $t > t_{QSS}$, i.e. ignoring any transient dynamics during initial penetration. Here, we discuss the results for $M$ and show that this distribution can be used to partition plume fluid into three classes which identify coherent regions of the plume.

Figure 9 shows the net mixing effect distribution in both physical and buoyancy-tracer space at $t = 14$, with the mixing flux distribution vectors $\boldsymbol {F}$ overlaid in buoyancy-tracer space. The distribution is represented in physical space by using the buoyancy and tracer concentration to map between buoyancy-tracer space and physical space, i.e. we plot $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. Recall that the net mixing effect distribution, $M$, quantifies the integrated effect of the mixing flux distribution, $\boldsymbol {F}$, or, equivalently, the volume difference between the volume distribution and the cumulative source distribution, representing the volume change due to mixing. As expected, we find $M < 0$ on the source line where undiluted plume fluid is continuously supplied before being mixed away into the $M > 0$ region. Environmental fluid is entrained into the plume via mixing and accumulates in the intrusion where $M$ is maximised. We define class U as undiluted plume fluid corresponding to the source region $\mathcal {U} = \{ M < 0 \}$, introduced in the definition of QSS in § 4.5. For convenience, we use the notation $\{ M\,{<}\,0 \}$ as shorthand for $\{ (b, \phi ) \vert M(b, \phi ; t)\,{<}\,0\}$ henceforth. The mixing flux vectors point along the source line, indicating that mixing within $\mathcal {U}$ is mostly within the plume rather than between the plume and environment, owing to the shielding effect. Once undiluted plume fluid is exposed to the environment, there is a strong mixing flux between $\mathcal {U}$ and the $\phi =\phi _{min}$ axis where environmental fluid joins the distribution.

Figure 9. Snapshots at non-dimensional time $t=14$ of the buoyancy-tracer net mixing effect distribution, $M(b, \phi ; t)$, in buoyancy-tracer space (a) and (b) in an $x-z$ cross-section of physical space. Buoyancy contours are shown in the surrounding environment. The mixing flux distribution, $\boldsymbol {F}$, is overlaid in the top panel and the vectors are uniformly scaled to be visible.

As discussed in § 4.1, we expect volume to accumulate in some region of buoyancy-tracer space. This is clearly demonstrated in physical space, where fluid collects in the intrusion after mixing with the environment. We can distinguish the accumulation region from the ‘transport’ region through which volume moves to reach the accumulation region by identifying a region in which $M$ is small and approximately constant. In this region, plume fluid is actively mixing with the environment and transporting volume away from the source line, but fluid does not accumulate in this region. Then, fluid that has undergone significant mixing and accumulates in the intrusion corresponds with a region 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 \leqslant m^*\}$ and we refer to fluid corresponding to 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 threshold $m^*(t)$ is found by identifying the value $m$ which minimises the total mixing flux convergence $f(m; t)$ associated with a region $\{ 0 < M \leqslant m(t)\}$ at each time $t$, where

(4.13)\begin{equation} f(m; t) = \int_{\{ 0 < M \leqslant m \}} \frac{\mathrm{d} W}{\mathrm{d} t} - S \, \mathrm{d}b\,\mathrm{d}\phi ={-}\int_{\{0 < M \leqslant m\}} \boldsymbol{\nabla}_{(b,\phi)} \boldsymbol{\cdot} \boldsymbol{F} \, \mathrm{d}b\,\mathrm{d}\phi, \end{equation}

so that the volume that enters the region $\mathcal {T}$ is approximately equal to the volume leaving $\mathcal {T}$. Then, remaining mixed fluid lies in the region $\mathcal {A}$ which must 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 using the discrete form of the volume distributions is discussed further in Appendix D. The time variation of the threshold $m^*(t)$ is shown in figure 19 in Appendix D.

The net mixing effect distribution with this partitioning is shown in figure 10, with class $U$ coloured blue, class $T$ coloured green and class $A$ coloured red. Within classes $U$ and $A$, the volume-weighted centre of mass is shown by a coloured triangle, approximately indicating the position in buoyancy-tracer space towards which mixing acts to move fluid within each class. Internal plume mixing of undiluted plume fluid redistributes volume on the source line towards the blue triangle and homogenisation of mixed fluid in the intrusion accumulates volume near the red triangle. Figure 10 demonstrates the correspondence between the regions $\mathcal {U}$, $\mathcal {T}$ and $\mathcal {A}$ and coherent structures of the plume. The source region $\mathcal {U}$ identifies the rising undiluted plume. The accumulation region $\mathcal {A}$ identifies the radially spreading intrusion. The transport region $\mathcal {T}$ corresponds with newly generated mixtures in the plume cap and subsiding fluid joining the intrusion. The partition of the full plume volume into the undiluted plume, plume cap, and intrusion volume is shown in figure 8(b) for $t > t_{QSS}$. 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 each remain approximately constant.

Figure 10. As in figure 9, with the distribution partitioned into three regions: $\mathcal {U}$ (blue) where $M \leqslant 0$, $\mathcal {T}$ (green) where $0 < M \leqslant m^*(t)$ and $\mathcal {A}$ (red) where $M > m^*(t)$. The threshold $m*(t)$ minimises the total mixing flux convergence defined by (4.13). Corresponding fluid classes $U$, $T$ and $A$ (respectively) shown in physical space. Triangles represent the volume-weighted centre of mass in class $U$ (blue) and class $A$ (red). Undiluted plume fluid in $\mathcal {U}$ is mixed towards the blue triangle whilst fluid accumulating in $\mathcal {A}$ is mixed towards the red triangle.

4.7. Entrainment

We calculate the entrained volume $E(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 up to time $t$. We have

(4.14)\begin{align} E(t) &\equiv \int_{\{W > 0\}} W(b, \phi; t) \, \mathrm{d}b \,\mathrm{d}\phi - \int_0^t\int_{\{W > 0\}} S(b, \phi; t') \, \mathrm{d}b \,\mathrm{d}\phi \,\mathrm{d}t' \end{align}

(4.15)\begin{align} &={-}\int_0^t \int_{\{W>0\}} \boldsymbol{\nabla}_{(b, \phi)} \boldsymbol{\cdot} \boldsymbol{F}(b, \phi; t') \, \mathrm{d}b \,\mathrm{d}\phi\, \mathrm{d}t' \end{align}

(4.16)\begin{align} &= \int_0^t \int_{\{\phi=\phi_{min}\}} F_{\phi}(b, 0; t') \,\mathrm{d}b \,\mathrm{d}t' \approx \sum_{t'}\sum_{i} \frac{F^\phi_{i,j=0}(t')}{\Delta \phi} \Delta t', \end{align}

where the second equality follows from time integrating (4.2) and the final equality follows from Green's theorem and the fact that the mixing 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 numerical calculation of $E$ using the discrete form of the mixing flux distribution $F_{ij}^{\phi }$ is given by (4.16).

A common definition of an ‘entrainment profile’ with respect to height is the fractional volume (or mass) increase with height (e.g. de Rooy et al. 2013). This is not a useful definition in the case where the plume overturns, since the rate of change with height captures multiple components of the plume which act to entrain fluid in (potentially) different ways. Exploiting the linear increase of buoyancy with height in the initially linearly stratified environment, we treat buoyancy as a rough proxy for height and define an entrainment profile with respect to buoyancy

(4.17)\begin{equation} e(b_i, t) = \int_0^t F_\phi(b_i, 0; t') \, \mathrm{d}t' \approx \sum_{t'} \frac{F_{i,j=0}^\phi(t')}{\Delta \phi \Delta b} \Delta t', \end{equation}

which represents the volume of environmental fluid entrained up to time $t$ per unit buoyancy. We can then define the volume entrained into a region $\mathcal {R}$ of buoyancy-tracer space which intersects the $\phi = \phi _{min}$ boundary as

(4.18)\begin{equation} E(\mathcal{R}, t) = \int_0^t\int_{\{\phi=\phi_{min}\}\cap\mathcal{R}} F_\phi(b, 0; t') \, \mathrm{d}b \,\mathrm{d}t' = \int_{\{\phi=\phi_{min}\}\cap\mathcal{R}} e(b, t) \, \mathrm{d}b, \end{equation}

which is numerically calculated by summing over the indices $i$ in (4.17) which belong to the intersection of $\mathcal {R}$ with the boundary $\phi =\phi _{min}$. We define the entrainment rate as the time rate of change of the entrained volume into a region $\mathcal {R}$, i.e. $\dot {E}(\mathcal {R}) = \partial _t E(\mathcal {R})$. This quantifies the rate at which volume is entrained into a physical volume represented by a region $\mathcal {R}$ in buoyancy-tracer space. Whilst we expect vigorous mixing at the boundary between a sub-volume of the plume and the environment to result in entrainment, the entrainment rate $\dot {E}$ does not necessarily quantify this since larger volumes would be expected to entrain more volume over time even if the ‘strength’ of the mixing is weaker. To quantify the strength of the entrainment into each sub-volume of physical space corresponding to a partitioning of buoyancy-tracer space, we define the specific entrainment rate as the ratio of the entrainment rate with the volume of each sub-region itself, i.e. $\dot {E}(\mathcal {R}) / V(\mathcal {R})$ for each region $\mathcal {R} = \mathcal {U}, \mathcal {T}, \mathcal {A}$.

Figure 11 shows the evolution of the entrainment profile through the simulation. Significant entrainment occurs over a large range of intermediate buoyancy values, indicating that most entrained volume is from the environmental fluid surrounding the intrusion. Using the partitioning introduced in § 4.6, 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})$. Figure 12(a) shows the entrained volume as well its decomposition. Here, a correction has been made to $E(t)$ to account for numerical artefacts – see Appendix B for details. Figure 12(b) shows the entrainment rate in each of the regions $\mathcal {U}, \mathcal {A}$ and $\mathcal {T}$. The relative contribution of plume cap entrainment and ‘lateral’ entrainment in the intrusion to the overall entrainment is noted as an open question in the study of fountains by Hunt & Burridge (2015), where the plume cap is analogous to the ‘fountain top’. Here, we find that the contribution to the entrained volume from the plume cap (class $T$) is weak compared with the intrusion (class $A$) when in QSS. Under the definition of QSS for this flow given in § 4.5, volume continuously accumulates in the intrusion and hence the volume of the intrusion dominates the volume of the full plume at late times. Since the intrusion spreads radially, there is a greater contact area between the intrusion and the surrounding environmental fluid compared with the plume cap and the rising plume. This allows a greater volume of environmental fluid to be entrained into the intrusion. This suggests that entrainment of environmental fluid from the lower part of the stratified layer into the intrusion is important for setting the centre of mass of the QSS buoyancy-tracer distribution. However, mixing during the overturning process in the plume cap near $z_{max}$ is important for setting the maximum accessible buoyancy of the volume distribution after mixing, and therefore the extent of the accumulation region in buoyancy-tracer space. The specific entrainment rate is larger in class $T$ than in class $A$, indicating stronger mixing with the environment in the plume cap in QSS compared with the intrusion.

Figure 11. Entrainment profile $e(b, t)$ defined by (4.17) at fixed time intervals post-penetration. The value $e(b,t)$ is the volume of environmental fluid entrained up to time $t$, per unit buoyancy, as a function of buoyancy $b$.

Figure 12. (a) Volume of environmental fluid entrained up to time $t$ (black line) and its decomposition into entrainment into the regions $\mathcal {U}, \mathcal {T}$ and $\mathcal {A}$ (coloured lines). These regions are defined in § 4.6 and the entrained volumes $E(\mathcal {U}), E(\mathcal {T})$ and $E(\mathcal {A})$ are computed using (4.18). (b) Specific entrainment rate in each class, defined as the ratio of the rate of change of the entrained volume and the volume, e.g. $\dot {E}(\mathcal {U})/V(\mathcal {U})$. In both plots, the vertical dotted line indicates the start of the QSS.

4.8. Three-stage mixing process

Together, the results presented in this section suggest a multi-stage mixing process in QSS convective penetration. This is summarised in figure 13, in which we show schematics of the volume distribution partitioned into three regions $\mathcal {U}, \mathcal {T}$ and $\mathcal {A}$ of buoyancy-tracer space, and the corresponding classes of fluid in cross-sections of the plume. Fluid moves through each stage of mixing from left to right, corresponding with an increasing value of the net mixing effect $M$, though occasionally the primary mixing stage may be skipped – see following discussion. In each stage, the convex envelope of fluid involved in mixing is indicated by a grey dashed line in $(b,\phi )$-space; arrows in $(b,\phi )$-space indicate the movement of individual fluid parcels due to the mixing process; and circular arrows in the $x-z$ cross-section indicate the location of mixing in physical space. These circular arrows are illustrative and are not intended to indicate the physical nature of the mixing mechanism in each stage. In the secondary mixing stage where multiple mixing processes occur simultaneously, three arrow types are shown which correspond between physical and buoyancy-tracer space. The mixing stages in QSS convective penetration are as follows:

1. (i) Mixing within the rising plume as undiluted plume fluid penetrates into the stratified layer. This fluid is shielded from the surrounding environment. Thus, in buoyancy-tracer space, internal mixing of undiluted plume fluid acts within $\mathcal {U}$ only, homogenising the distribution and consequently moving fluid towards the centre of mass of class $U$ indicated by the blue triangle in figure 10. This sets the buoyancy-tracer characteristics of plume fluid that is first exposed to environmental fluid near $z_{max}$.

2. (ii) Primary mixing between the plume and environment occurs in the plume cap. This mixing may extend below the plume cap into the intrusion as the mixing time scale is slow compared with the dynamical time scale. In buoyancy-tracer space, the primary mixing stage acts on the undiluted plume fluid concentrated near the class $U$ centre of mass and the surrounding environmental fluid with values of buoyancy close to $b = z_{max}$. The resulting mixtures are transported towards intermediate values of buoyancy as indicated by the mixing flux distribution in figure 9(a).

3. (iii) After primary mixing with the environmental fluid, subsiding mixed fluid joins the intrusion and homogenises with its surroundings, moving from class $T$ to class $A$ (dashed arrows). A number of mixing processes act on the fluid accumulating in class $A$. In buoyancy-tracer space, mixing occurs between four regions of fluid as indicated by the convex envelope in figure 13. Secondary mixing with the environment occurs as environmental fluid lower in the stratified layer is entrained into the intrusion as it spreads (solid arrows), dominating the volume entrained by the full plume at late times. 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 (dotted arrows) as seen in figure 9(a). Finally, mixing in the interior of the intrusion homogenises the buoyancy-tracer distribution of fluid that accumulates within region $\mathcal {A}$. As fluid moves radially, large volumes of fluid in class $A$ are concentrated near the centre of mass in $(b,\phi )$-space, where $M$ (and $W$) are largest.

Figure 13. Schematic of the three-stage mixing process in QSS convective penetration of a buoyant plume into a stably stratified layer, identified by the partitioned buoyancy-tracer volume distribution. Buoyancy-tracer space is shown on the top row, with arrows indicating the movement of volume within each mixing stage. The region of buoyancy-tracer space affected by mixing in each stage is indicated by a grey dashed envelope. Physical space is shown on the bottom row with circular arrows indicating where mixing is located. In the secondary mixing stage the three distinct mixing processes are shown as dotted, dashed and solid arrows, both in physical and buoyancy-tracer space.

5. Mixing diagnostics

In this section we use the partitioning introduced in § 4 to examine the statistics of mixing in each stage of the plume evolution.

5.1. Characterising mixing

To characterise the physical nature of the mixing in each stage, we consider the mixing efficiency which relates the total energy expended in turbulent mixing with the actual mixing achieved (e.g. Davies Wykes, Hughes & Dalziel 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. In buoyancy-driven stratified flows, the mixing efficiency is quantified by utilising the partitioning of potential energy into available potential energy (APE) and background potential energy (BPE). For an incompressible, Boussinesq flow BPE is the potential energy that is not available to do work whilst APE represents the energy stored in the buoyancy field if the flow is not in gravitational equilibrium (Lorenz 1955; Davies Wykes et al. 2015). Turbulent mixing irreversibly converts APE into BPE and results in dissipation of turbulent kinetic energy (TKE). The conversion of APE into BPE quantifies the energy expenditure that results in mixing of the buoyancy field and its sum with the TKE dissipation represents the total energy expended by turbulent mixing. The ratio of these two quantities forms the mixing efficiency. Following Howland, Taylor & Caulfield (2020) and Holliday & McIntyre (1981), in the case where $\partial _z b$ is constant in the initial stratified environment we may treat the quadratic form of the potential energy

(5.1)\begin{equation} E_p = \tfrac{1}{2} \langle b'^2 \rangle, \end{equation}

as a proxy for APE. Here, $b'(\boldsymbol {x},t) = b(\boldsymbol {x},t) - z$ is the departure from the linear initial stratification and $\langle {\cdot } \rangle$ denotes a volume average. We refer to $E_p$ as the perturbation potential energy (PE). A full derivation of the perturbation PE budget is given in Appendix E, which follows the derivation as described in Howland et al. (2020) except with SGS terms included. The irreversible conversion of perturbation PE to BPE that results from the reduction of buoyancy gradients by mixing is captured by the buoyancy variance dissipation rate

(5.2)\begin{equation} \bar{\chi} = \left\langle \kappa_{tot} \,|\boldsymbol{\nabla} b' |^2 \right\rangle, \end{equation}

where $\kappa _{tot} = (\textit {Re} \textit {Pr})^{-1} + \kappa _{SGS}^{(b)}$ is the total diffusivity of buoyancy. The buoyancy variance dissipation rate represents the primary sink of $E_p$. The total energy dissipated via turbulent mixing is the sum of $\bar {\chi }$ and the dissipation rate of TKE

(5.3)\begin{equation} \bar{\varepsilon} = \left\langle \nu_{tot} \frac{\partial u_i}{\partial x_j} \frac{\partial u_i}{\partial x_j} \right\rangle, \end{equation}

where $\nu _{tot} = \textit {Re}^{-1} + \nu _{SGS}$ is the total viscosity. The TKE dissipation rate acts as the primary sink of TKE. The instantaneous mixing efficiency $\eta$ is then

(5.4)\begin{equation} \eta = \frac{\bar{\chi}}{\bar{\chi}+\bar{\varepsilon}}, \end{equation}

as in e.g. Howland et al. (2020) and Peltier & Caulfield (2003). Here, we use an overbar to denote a volume averaged quantity. It is instructive to examine the spatial structure of the dissipation rates. In particular, regions of large $\chi$ indicate intense buoyancy gradients and regions of large $\varepsilon$ indicate intense turbulent motion.

Further information on the state of turbulence in stratified flows is often drawn from the buoyancy Reynolds number $\textit {Re}_b \equiv \langle \varepsilon \rangle / \nu N^2$. We define a pointwise activity parameter $I$ and its mean $\bar {I}$

(5.5a,b)\begin{equation} I = \frac{\dfrac{\partial u_i}{\partial x_j} \dfrac{\partial u_i}{\partial x_j}}{\left|\dfrac{\partial b}{\partial z}\right|} = \frac{\varepsilon}{\nu_{tot} \left|\dfrac{\partial b} {\partial z}\right|}, \quad \bar{I} = \frac{\langle \varepsilon /\nu_{tot} \rangle}{\left\langle \left| \dfrac{\partial b}{\partial z} \right| \right\rangle}. \end{equation}

The bulk property $\bar {I}$ is analogous to $\textit {Re}_b$ except with SGS contributions to viscosity included. Also, we replace the global buoyancy time scale $N^{-1}$ with a local measure of the buoyancy time scale $\lvert \partial _z b\rvert ^{-1}$ given by the local buoyancy gradient. This is a more appropriate measure since buoyancy gradients within the plume differ significantly from the background linear stratification and are more representative of the regime in which mixing occurs inside the plume. As with $\textit {Re}_b$, the mean activity parameter $\bar {I}$ 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 time scale $(\partial b / \partial z)^{-1/2}$ to the time scale of development of turbulent effects $(\varepsilon / \nu _{tot})^{-1/2}$ (Ivey, Winters & Koseff 2008). Regions of large $I$ indicate active turbulence (García-Villalba & del Álamo 2011) that is weakly affected by stratification.

5.2. Results and discussion

Cross-sections in physical space of the dissipation rates $\varepsilon$ and $\chi$, activity parameter $I$ and local buoyancy gradient $\partial _z b$ are shown in figure 14. These diagnostics are only shown in the plume, where $\phi > \phi _{min}$. In the surrounding environment, buoyancy contours are shown. The mixing diagnostics within the plume structures identified by classes $U$, $T$ and $A$ are quantified by histograms separated into each class in figure 15. The colours for each class correspond with those used in § 3 (e.g. figure 13). The black dashed lines show the histograms for the full plume, i.e. all fluid within the stratified layer where $\phi > \phi _{min}$. This full histogram is normalised to form a p.d.f. The partitioned histograms are scaled so that the sum of the class $U$, $T$ and $A$ histograms equals the full plume p.d.f. These histograms are summarised by volume averages of the mixing diagnostics within each class, as well as the full plume volume average, in table 2.

Figure 14. The $x-z$ cross-sections of the mixing diagnostics $\varepsilon, \chi, I$ and $\partial _z b$ within the plume, where $\phi > \phi _{min}$ at non-dimensional time $t=14$. Mixing diagnostics are defined in § 5.1. Cross-sections are taken at the plume centreline. Buoyancy contours are shown outside the plume.

Figure 15. (a,c,e,f) Histograms of the mixing diagnostics $\varepsilon, \chi, I$ and $\partial _z b$ defined in § 5.1 (black dashed line) decomposed into histograms within class $U$, $T$, $A$ (coloured lines). The histograms within each class are scaled so that their sum is the histogram of the full plume. In (a,c), the histograms are further separated into those with $\nu _{SGS}=0$ and $\kappa _{SGS}^{(b)}=0$ (light colours) and with $\nu _{SGS}$ and $\kappa _{SGS}$ non-zero (usual colours). (b,d) Histograms of the total viscosity $\nu _{tot}$ and total buoyancy diffusivity $\kappa _{tot}$ respectively. All histograms are computed at $t=15$. All vertical axes show the scaled histogram frequency. The scaling is such that the total area is unity for the full plume histograms and the sum of the partitioned histograms equals the full plume histogram. Similarly in (b,d) the frequency is scaled so that the total area is unity.

Table 2. Volume-averaged mixing quantities $\varepsilon, \chi, I, \partial _z b$ and mixing efficiency $\eta$ defined in § 5.1 at $t=15$. Averages are computed over the full plume and within class $U$, $T$ and $A$. The percentage of the full plume volume associated with each class is given to indicate the relative contribution of each class to the full plume average.

The histograms for TKE dissipation in figure 15(a) and buoyancy variance dissipation in figure 15(c) are further separated based on where $\nu _{SGS}$ and $\kappa _{SGS}$, respectively, are non-zero or vanish. This distinction is made since the total viscosity $\nu _{tot}$ and total dissipation $\kappa _{tot}$ are bimodal with a peak where the SGS contribution vanishes (where the simulation effectively switches from LES to direct numerical simulation, such that $\nu _{tot}$ and $\kappa _{tot}$ reduce to the molecular values) and a peak where the SGS contribution is non-zero – see figures 15(b) and 15(d). The separation of the histograms based on non-zero and zero SGS contribution demonstrates that the bimodality of the histograms for $\varepsilon$ and $\chi$ is a consequence of the bimodal total viscosity and diffusivity alone. The vertical buoyancy gradient histogram is shown on a logarithmic scale since the most extreme values are rare but remain important for mixing.

The results show that the mixing statistics are quantitatively different in each class, suggesting that the mixing regimes differ in each of the three stages of the plume evolution highlighted in § 4.8. The cross-sections in figure 14 show that TKE dissipation $\varepsilon$ is particularly large in the rising undiluted plume and in some regions of the plume cap. This is supported by the histograms which show that TKE dissipation is an order of magnitude larger in class $U$ and T compared with class $A$. The activity parameter is also largest in the rising undiluted plume but comparable in the plume cap and intrusion. The largest values (both positive and negative) of the vertical buoyancy gradient $\partial _z b$ are found in the plume cap (class $T$) with a clear positive bias. The relatively tighter spread of the $\partial _z b$ distribution in class $U$ in figure 15(f) suggests that smaller vertical buoyancy gradients are more common in the undiluted plume fluid as compared with the plume cap and intrusion. This explains the increased magnitude of the activity parameter in class $U$ compared with class $T$ despite similar magnitudes of TKE dissipation; $\varepsilon$ is large in both classes but $\partial _z b$ is generally larger in class $T$ and hence $I$, proportional to the ratio $\varepsilon /\partial _z b$, is smaller in class $T$. The strongest buoyancy gradients as measured by $\chi$ are found in the primary mixing region (class $T$), where fluid overturns and subsides fluid in the plume cap. Intense buoyancy gradients are also found in the rising column of undiluted plume fluid, but $\chi$ is more sparse here than in the plume cap. The strong buoyancy gradients found at the top edge of the plume cap are a result of the relatively less buoyant plume impinging on the more buoyant surrounding environmental fluid. Note that the largest values of $\chi$ at the extreme positive tail of the full plume p.d.f. are almost all from class $T$.

We summarise the mixing regimes described by the mixing diagnostics as follows. In class $U$, we find active turbulence with large dissipation of TKE, consistent with the undiluted plume being unaffected by the surrounding stratified environment due to the shielding effect. Since the undiluted plume fluid becomes well mixed during its rise through the uniform layer, there are relatively small buoyancy gradients. As a result, there is relatively little PE dissipation. In class $T$, rising plume fluid impinges upon the more buoyant surrounding environmental fluid, generating particularly strong vertical buoyancy gradients. Horizontal buoyancy gradients are also generated by overturning. Turbulence advected upwards in class $U$ is carried over into class $T$ but gradually suppressed by the intense buoyancy gradients present in the plume cap, thus reducing the activity parameter. As the turbulent motion stirs these buoyancy gradients, significant PE dissipation occurs and consequently a large mixing efficiency is achieved, around $50\,\%$. 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. There is some indication of a weak stratification and layering 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, TKE dissipation is weak as the buoyant forces driving turbulent motion are weakened by earlier mixing. Whilst both TKE and PE dissipation are weak in class $A$, they are of similar magnitude and hence a greater mixing efficiency is achieved compared with class $U$. This could be attributed to the continued entrainment of environmental fluid above and below the intrusion, introducing small-scale buoyancy gradients which are acted upon by the weak turbulent motion.

6. Conclusions

In this paper we have analysed a LES of a buoyant pure plume penetrating into a linearly stably stratified layer. We have outlined the buoyancy-tracer volume distribution formalism to examine tracer transport via turbulent mixing. Using this formalism, we developed a method for objectively partitioning buoyancy-tracer space into three regions based on the net change in volume due to mixing. Each of these regions identifies corresponds with a class of fluid lying in coherent regions of the plume in physical space. The method distinguishes undiluted plume fluid (class $U$) from mixtures of plume and environmental fluid. Mixed fluid is further partitioned into newly generated mixtures in the plume cap that are actively mixing with the environment (class $T$), and fluid in the radially spreading intrusion that has already undergone significant mixing (class $A$). In buoyancy-tracer space, the intrusion corresponds with an accumulation region (corresponding with class $A$) where volume collects and homogenises. Active mixing with the environment in the plume cap moves volume from the source region, through a transport region (corresponding with class $T$), into the accumulation region. The accumulation region represents the majority of the plume volume at late times. To quantify the mixing regime in each class of fluid, we use the buoyancy variance dissipation rate, TKE dissipation rate, vertical buoyancy gradient and an activity parameter as diagnostic variables in each sub-volume of the plume. The instantaneous mixing efficiency is also calculated by treating the buoyancy variance dissipation rate as a proxy for the energy dissipated in turbulent mixing that actually results in mixing of buoyancy.

Our results demonstrate a three-stage mixing process in QSS penetration of a plume into a stably stratified layer. In the first ‘plume mixing’ stage (class $U$), turbulence within the undiluted rising plume homogenises the buoyancy-tracer distribution as fluid rises into the stratified layer. The turbulent motion near the centreline of the plume is relatively unaffected by the surrounding stratification owing to a shielding effect from the plume edge and surrounding intrusion. This homogenisation process sets the range of buoyancy and tracer concentration which is first exposed to the environment when fluid overturns near the maximum penetration height. The ‘primary mixing’ stage occurs as rising fluid impinges on the more buoyant environment, establishing intense buoyancy gradients in the plume cap (class $T$). The mixing of undiluted plume fluid with the surrounding environment near $z_{max}$ has a particularly large mixing efficiency. The maximum penetration height approximately determines the buoyancy of the environmental fluid involved in the primary mixing stage and sets the extent of the buoyancy-tracer convex envelope for the remainder of the mixing process. 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. Several secondary mixing processes occur 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 plume cap. The volume of environmental fluid entrained into the intrusion during QSS dominates the volume into the plume as a whole at late times, but entrainment in the plume cap is ‘strongest’ in the sense of the largest fractional rate of increase in volume.

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 TKE dissipation significantly larger than PE dissipation. In the plume cap, the intense buoyancy gradients result in large PE dissipation and small TKE dissipation and hence the entrainment of the surrounding environment achieves a large mixing efficiency, with around 50 % of the total energy dissipated by turbulence resulting in mixing. As mixed fluid homogenises in the intrusion and further environmental fluid is entrained, weak buoyancy gradients are continually introduced and eroded by weak turbulence with low TKE dissipation. The mixing efficiency in the intrusion is moderately large, though smaller than in the plume cap.

Models of convective penetration 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 proportion of the full plume volume in each sub-region. 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 the mixing efficiency is small. At late times as fluid accumulates in the intrusion, the plume volume is dominated by the intrusion with weaker buoyancy gradients and turbulence and a greater mixing efficiency. Therefore the turbulent statistics associated with the full plume volume must change over time.

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 volume distribution budget equation (4.2) highlight the mixing processes that occur and the resulting tracer transport.