1. Introduction
Gravity current refers to the phenomenon in which a flow of a given density
$\rho _0$
is released into a lighter ambient fluid of density
$\rho _a\lt \rho _0$
in the presence of a rigid boundary (Huppert Reference Huppert2006; Sher & Woods Reference Sher and Woods2017). The dense current will flow underneath the ambient fluid under the influence of the buoyancy force. The importance of these currents rose due to their occurrence in a vast number of natural and industrial flows. For example, due to intense evaporation, the Mediterranean Sea exhibits higher salinity than the Atlantic Ocean and the Black Sea. This salinity difference drives a dense, gravity-driven outflow from the Mediterranean Sea into these basins (Stashchuk & Hutter Reference Stashchuk and Hutter2003). In the atmosphere, the main examples are the katabatic winds generated when cold air in the mountains flows downhill into warmer ambient air and the sea breezes where marine cold air flows towards the warmer air above the land. These currents are associated with a highly turbulent regime which causes the vertical transport and mixing of mass, heat and pollutants affecting the surrounding climate (Ellison & Turner Reference Ellison and Turner1959; Dallimore, Imberger & Ishikawa Reference Dallimore, Imberger and Ishikawa2001; Miller et al. Reference Miller, Keim, Talbot and Mao2003; Princevac, Fernando & Whiteman Reference Princevac, Fernando and Whiteman2005; Frantz et al. Reference Frantz, Deskos, Laizet and Silvestrini2021). In experimental contexts, gravity currents are often classified by their release type. Instantaneous gravity currents form from the sudden removal of a barrier separating two fluids of different densities. This action creates a current driven by the initial potential energy of the denser fluid, making it a closed system with no external energy exchange. In contrast, continuous gravity currents are generated by a discharge of denser fluid into a resting ambient fluid from outside the domain. These are considered open systems as energy is continuously supplied and can leave the domain.
In order to predict the influence of gravity currents in real life, various models that assess their propagation and dilution properties have been developed. These models are usually accomplished by employing an empirical parameter that relates the known, or easily measured, properties with the required unknown variables. For instance, when dealing with the front location,
$x_f$
, of a continuous, fixed influx gravity current during some intermediate time of propagation,
$t$
, we find the relation (Sher & Woods Reference Sher and Woods2017)
\begin{equation} x_f = \lambda B_0^{1/3}t, \end{equation}
with
$B_0 = g^{\prime}_0 Q_0$
being the source buoyancy flux per unit width, where
$g^{\prime}_0 = g (\rho _0-\rho _a)/\rho _a$
and
$Q_0$
being the buoyancy and the volume flow rate at the source per unit width, respectively,
$g$
being the gravitational acceleration, and
$\lambda (\textit{Fr}_0)$
being the empirical front speed parameter that depends on the source Froude number, defined as
\begin{equation} \textit{Fr}_0 = u_0 / \sqrt {g^{\prime}_0 h_0}, \end{equation}
with
$u_0$
and
$h_0$
being the velocity and the height of the current at the source.
$\lambda(\textit{Fr}_0\lt1) = 0.85+0.17\textit{Fr}_0$
, and
$\lambda(\textit{Fr}_0\gt1) =1.1 \pm 0.1$
(Sher & Woods Reference Sher and Woods2017). Based on
$\textit{Fr}_0$
, two regimes can be determined. In the subcritical regime,
$\textit{Fr}_0\lt 1$
, the supplied dense current accelerates and plunges underneath the ambient fluid and exhibits strong vertical derivatives of density and velocity. For supercritical currents,
$\textit{Fr}_0\gt 1$
, the supplied dense current decelerates and deepens due to ambient fluid entrainment caused by the shearing of a weakly stratified fluid. For further details about the propagation and geometrical characteristics of continuous gravity currents, we refer readers to the experimental study of Sher & Woods (Reference Sher and Woods2017) and the numerical work of Haddad et al. (Reference Haddad, Vaux, Varrall and Vauquelin2022). The propagation characteristics have been studied by several authors, including Britter & Linden (Reference Britter and Linden1980), Shringarpure et al. (Reference Shringarpure, Lee, Ungarish and Balachandar2013) and Hogg et al. (Reference Hogg, Nasr-Azadani, Ungarish and Meiburg2016).
Similarly, the turbulent entrainment and mixing in gravity currents are elementary aspects that require modelling. The inclusion of mixing into a mathematical model usually happens via the equations of conservation of energy. The reason is that the mixing causes a portion of the kinetic energy to be irreversibly transferred into potential energy through the buoyancy flux. In the literature of turbulent stratified fluids, two kinds of energy budget formulae have been considered. First, a formula of the turbulent kinetic energy (TKE) budget type, where the turbulent properties deviating from the mean flow have to be measured. This methodology is quite complex since it requires generating many realisations of the flow in order to define the mean flow. The second, and simpler, energy budget formula is of the available potential energy type, introduced by Lorenz (Reference Lorenz1955) and employed by Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995). The irreversible transfer of energy from kinetic energy to potential energy due to mixing causes the variation of the potential energy of the background state. By examining this variation, the energy flux of mixing can be determined. This methodology explicitly requires knowledge of the density field. A mixing efficiency parameter is used to link the rate of mixing with other measurable or known energy fluxes that lead to mixing. Early analytical and experimental studies on the parametrisation of mixing in various stratified turbulent flows can be found in the literature (Phillips Reference Phillips1972; Posmentier Reference Posmentier1977; Linden Reference Linden1979; Osborn Reference Osborn1980; Ivey & Imberger Reference Ivey and Imberger1991; Park, Whitehead & Gnanadeskian Reference Park, Whitehead and Gnanadeskian1994; Caulfield & Peltier Reference Caulfield and Peltier2000; Strang & Fernando Reference Strang and Fernando2001; Peltier & Caulfield Reference Peltier and Caulfield2003; Caulfield Reference Caulfield2020). In oceanic models, the mixing efficiency, which usually relates the mixing rate to TKE dissipation, is assumed to be a constant value, typically between 0.17 and 0.2 (Osborn Reference Osborn1980). However, a series of studies showed that the mixing efficiency depends on the considered stratified flow configuration, the agent which causes the stratification, the contribution of the vertical derivatives of buoyancy and velocity, and the turbulence intensity. Numerous studies have focused on quantifying mixing in instantaneous gravity currents (Prastowo et al. Reference Prastowo, Griffiths, Hughes and Hogg2008; Fragoso, Patterson & Wettlaufer Reference Fragoso, Patterson and Wettlaufer2013; Ilıcak Reference Ilıcak2014; Hughes & Linden Reference Hughes and Linden2016; Micard et al. Reference Micard, Dossmann, Gostiaux and Venaille2016; Ottolenghi et al. Reference Ottolenghi, Adduce, Inghilesi, Armenio and Roman2016; Mukherjee & Balasubramanian Reference Mukherjee and Balasubramanian2020; Agrawal et al. Reference Agrawal, Ramesh, Zimmerman, Philip and Klewicki2021; Frantz et al. Reference Frantz, Deskos, Laizet and Silvestrini2021; Mukherjee & Balasubramanian Reference Mukherjee and Balasubramanian2021). The mixing efficiency in these studies is usually defined as the ratio of mixing energy to the initial available potential energy and has been computed as a function of the bulk Reynolds number,
$Re_b = \sqrt {g^{\prime}_0 H}H/\nu$
, with
$H$
being the initial height of the column of dense fluid and
$\nu$
being the kinematic viscosity. It was found that the mixing efficiency increases with
$Re_b$
and saturates at an approximate magnitude of 0.11 for
$Re_b\gt 50\,000$
(Prastowo et al. Reference Prastowo, Griffiths, Hughes and Hogg2008; Ilıcak Reference Ilıcak2014; Hughes & Linden Reference Hughes and Linden2016; Agrawal et al. Reference Agrawal, Ramesh, Zimmerman, Philip and Klewicki2021). An experimental examination of turbulent mixing in a statistically steady-state region of a continuous gravity current was conducted by Odier, Chen & Ecke (Reference Odier, Chen and Ecke2014). They found that within the gravity current, the mixing efficiency, defined as the ratio of mixing rate to the rate of turbulence production, is roughly constant and equal to
$0.1 \pm 0.01$
. The mixing and entrainment in continuous gravity currents, under statistically homogeneous or steady-state turbulence conditions, have also been considered in other studies (van Reeuwijk et al. Reference van Reeuwijk, Holzner and Caulfield2019; Salinas et al. Reference Salinas, Zúñiga, Cantero, Shringarpure, Fedele, Hoyal and Balachandar2022).
To our knowledge, no existing study has parametrised the global mixing energy in continuously forced gravity currents. We address this gap by proposing a new framework that relates mixing to the source Froude number
$\textit{Fr}_0$
. Unlike lock-release configurations, the continuous energy input in our open system precludes the use of mixing efficiencies based on fixed initial potential energy. Instead, we relate the mixing energy to the total energy supplied by the source. This approach also differs from oceanic parametrisations, which typically link mixing to measurable quantities such as kinetic energy dissipation. We introduce a dimensionless parameter
$\gamma$
, representing a form of mixing efficiency adapted to open systems, defined as the ratio of mixing energy to total supplied energy. This formulation enables a practical estimation of mixing rates based solely on source properties. The numerical configuration and set-ups are introduced in § 2. The energetics problem of the present configuration is described in § 3.1. The available potential energy framework developed by Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995) is employed for the present configuration in § 3.2. In addition, since continuous gravity currents are open systems, we show how to compute the boundary influx of energy by using source variables in § 3.3. The results of the present study are discussed in § 3.4. A conclusion is presented in § 4.
The 3-D computational domain with the initial and source conditions.
2. Numerical configuration and set-ups
We consider a numerical configuration, shown in figure 1, close to that employed experimentally by Sher & Woods (Reference Sher and Woods2017). The flow fields are obtained by direct numerical simulations (DNS) via the Nek5000 solver (Version 19.0, Argonne National Laboratory, Illinois, https://nek5000.mcs.anl.gov). The three-dimensional (3-D), Newtonian and incompressible Navier–Stokes equations, along with the equation of advection–diffusion of density, are used to resolve the flow field. The Boussinesq approximation is applied. Hence, the governing equations are
\begin{equation} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{\boldsymbol{u}} = 0, \end{equation}
\begin{equation} \frac {\partial \boldsymbol{\boldsymbol{u}}}{\partial t} + (\boldsymbol{\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{\boldsymbol{u}} = - \frac {1}{\rho _a} \boldsymbol{\nabla }p + \nu \boldsymbol{\nabla }^2\boldsymbol{\boldsymbol{u}} + \frac {\rho }{\rho _a} \boldsymbol{{g}}, \end{equation}
\begin{equation} \frac {\partial \rho }{\partial t} + \boldsymbol{\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla }\rho = \kappa \boldsymbol{\nabla }^2\rho, \end{equation}
where
$\boldsymbol{\boldsymbol{u}} = (u_x, u_y,u_z)$
is the velocity vector,
$p$
and
$\rho$
are the pressure and density fields, respectively,
$\nu$
and
$\kappa$
are the kinematic viscosity and mass diffusivity, respectively, and
$\boldsymbol{{g}}=g\boldsymbol{{e_y}}$
being the vector of gravitational acceleration with
$\boldsymbol{{e_y}} = (0,-1,0)$
being a unit vector oriented downward. We assume that the density evolves linearly with the scalar
$S$
, considered here to be the concentration of a dissolved substance, so that
\begin{equation} \rho = \rho _a [1 + \alpha _d(S - S_a)], \end{equation}
with
$\alpha _d = (\rho _0 - \rho _a)/\rho _a$
being the scaled density difference, and
$\rho _a = 1000 \, \mathrm {kg\,m^{-2}}$
being the density of the ambient fluid. The dimensionless concentration of the dissolved substance in the ambient fluid is set to
$S_a = 0$
, making the dimensionless concentration field in a range
$0 \leqslant S \leqslant 1$
where the concentration at the source is represented by
$S_0 = 1$
. All the boundaries, except the inlet and the upper boundaries, satisfy the no-slip and no-penetration conditions for the velocity and an insulated Neumann condition for the scalar. The upper boundary, at
$y = L_y$
, satisfies an open boundary condition with zero stresses. This condition allows the ambient fluid to leave the domain along the upper boundary when new heavier fluid is supplied to the domain. At the inlet, where
$x = 0$
and
$0 \leqslant y \leqslant h_0$
, an influx of velocity
$u_0$
and dimensionless scalar
$S_0 = 1$
are defined. Through the series of simulations, we keep all the source variables constant and we vary only the initial density difference resulting in various
$\textit{Fr}_0$
cases, as shown in table 1. The source Reynolds number,
$Re_0 =Q_0/\nu = u_0h_0/\nu$
, is kept constant at 5000 so that a turbulent regime is held, with
$\nu = 10^{-6}\,\mathrm {m^{-2}\,s^{-1}}$
. The Kolmogorov scale at which the dissipation of the velocity fluctuation occurs is expressed as
$L_K = (\nu ^3/\epsilon)^{1/4}$
, with
$\epsilon$
being the dissipation rate of TKE. Hughes & Linden (Reference Hughes and Linden2016) scaled the dissipation rate in a lock-release gravity current as
$\epsilon = 0.04 (u_f^3/h)$
, with
$u_f$
being the front speed and
$h$
being the characteristic height of the current, approximated for full-depth release as
$H/2$
. Inspired by their scaling, we approximate the dissipation rate in continuous gravity current as
$\epsilon = 0.04 (u_0^3/h_0)$
, using the characteristics source speed
$u_0$
and source height
$h_0$
. We note that our definition serves as an approximation and more refined definitions could be developed for higher accuracy. With this definition, the Kolmogorov scale will be approximately
$0.15\,{\text{mm}}$
. The importance of the advective transport to the diffusive transport of mass is controlled by the product
$Re_0Sc$
with
$Sc = \nu /\kappa$
being the Schmidt number. The diffusion of the dissolved substance takes place at the Batchelor length scale
$L_B = L_K/\sqrt {Sc}$
. The value of Schmidt number for the diffusion of salt in water is
$Sc \sim 600{-}700$
. Therefore, the selection of salt as the dissolved substance will require the governing equation to be solved down to
$L_B \approx L_K/25 = 0.006\,{\text{mm}}$
, which increases the computational cost of the numerical simulation. In order to keep the computational cost manageable, we have employed
$Sc = 1$
so that
$L_B = L_K$
, which is a common practice in the literature (Bonometti & Balachandar Reference Bonometti and Balachandar2008). The 3-D domain is discretised by a total of
$300 \times 58 \times 20$
elements. Each element contains
$8^3$
Gauss–Lobatto–Legendre (GLL) points. As a result, a total of approximately
$1.8 \times 10^8$
grid points are used. The domain was divided vertically into three regions of different vertical spacing. The finer resolution covers the region where the gravity current propagates. The time step is selected based on the Courant–Friedrichs–Lewy (CFL) criterion with a Courant number
$C \lt 0.5$
to ensure the stability of the numerical simulation.
The initial control parameter for the considered simulations. The maximum time
$t_{max}$
represents the physical time of the last extracted file at which the gravity current nearly reaches the right wall of the domain.
3. The analyses of the energetics
In this section we employ the framework of Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995) into the open-system, continuous gravity currents considered in the study. In § 3.1 we explain the initial and source states of the considered energetics problem. We define the potential energy budget in § 3.2. In § 3.3 we use the idea of Prastowo et al. (Reference Prastowo, Griffiths, Hughes and Hogg2008) and Hughes & Linden (Reference Hughes and Linden2016) applied for lock-release gravity current to define a hypothetical state where mixing is neglected. Finally, in § 3.4 the parametrisation of the mixing energy is presented. The simulations have 3-D geometry; however, for simplicity, we choose to work with two-dimensional (2-D) spanwise-averaged fields by assuming a quasi-2-D density field. Therefore,
$\boldsymbol {x}=(x, y, z)$
and
$\boldsymbol{\nabla } = (\partial _x, \partial _y, \partial _z)$
are simply replaced by
$\boldsymbol {x}=(x, y)$
and
$\boldsymbol{\nabla } = (\partial _x, \partial _y)$
, respectively. This selection is commonly employed in the analyses of experimentally and numerically generated gravity currents. The analyses of the full 3-D fields and comparison with the results of the width-averaged fields will be tackled in future work.
3.1. Energetics configuration
The present system, shown in figure 1, consists of a rectangular domain of volume
$\mathcal{V}$
initially filled with ambient fluid of density
$\rho _a$
at rest. The initial total energy inside the domain is negligible, i.e.
$\textit{TE}(t=0 \,\textrm{s}) = 0$
, noting that the potential energy is defined by the density anomaly, making it negligible at
$t=0\,\textrm{s}$
. The heavier fluid, of density
$\rho _0$
, is continuously supplied to the domain and the ambient fluid is allowed to leave the domain at the upper boundary. The total supplementary energy injected to the domain, per unit width, can be expressed as
\begin{equation} \textit{TE}_0(t) = \textit{KE}_0(t) + \textit{PE}_0(t) = \left [\left (\frac {1}{2}\rho _0u_0^2\right)+\left (g\rho ^{\prime}_0\frac {h_0}{2}\right)\right]Q_0t \end{equation}
where
$\textit{KE}_0(t) = ({1}/{2})\rho _0u_0^2Q_0t$
is the source kinetic energy,
$\textit{PE}_0(t) = g\rho ^{\prime}_0({h_0}/{2})Q_0t$
being the source potential energy, and
$\rho ^{\prime}_0 = \rho _0-\rho _a$
being the source density anomaly.
3.2. Total, background and available potential energy reservoirs
At each time, the total,
$\textit{PE}(t)$
, background,
$\textit{BPE}(t)$
, and available,
$\textit{APE}(t)$
, potential energy reservoirs inside the domain are expressed as
\begin{equation} \textit{PE}(t) = g\int _{\mathcal{V}} \rho ^{\prime}(\boldsymbol{x}, t) y {\rm d}\mathcal{V}, \end{equation}
\begin{equation} \textit{BPE}(t) = g\int _{\mathcal{V}} \rho ^{\prime}(\boldsymbol{x},t)y^*(\boldsymbol{x},t) {\rm d}\mathcal{V}, \end{equation}
\begin{equation} \textit{APE}(t) = \textit{PE}(t) - \textit{BPE}(t)= g\int _{\mathcal{V}} \rho ^{\prime}(\boldsymbol{x}, t) [y - y^*(\boldsymbol{x},t)] {\rm d}\mathcal{V}, \end{equation}
with
$\rho ^{\prime}=\rho -\rho _a$
being the density anomaly, and
$y^*(\boldsymbol{x},t)$
being a field that corresponds to the vertical location of a fluid parcel of a certain density
$\rho (\boldsymbol{x},t)$
in the background state. A background density state is a result of an adiabatic rearrangement of fluid parcels in which the gravity drives denser parcels to settle below the lighter ones, establishing a monotonic decrease of density with depth (Peltier & Caulfield Reference Peltier and Caulfield2003). The
$y^*(\boldsymbol{x},t)$
field is usually defined by using a Heaviside function, such as
\begin{equation} y^*(\boldsymbol{x},t) = \frac {1}{L_x}\int _{\mathcal{V}} H([\rho (\boldsymbol{x}^{\prime},t)-\rho (\boldsymbol{x},t)]){\rm d}\mathcal{V}^{\prime}, \end{equation}
where
$L_x$
is the length of the computational domain, and the variables followed by a prime denote dummy variables that sweep through all the points in the domain. An example of the resultant
$y^*(\boldsymbol{x},t)$
field is shown in figure 2(a).
(a) The visualisation of the distribution of
$y^*(\boldsymbol{x}, t_{max})$
with
$\textit{Fr}_0 = 2.7$
. The background density profile (solid blue line) is compared with the hypothetical one (solid red line) for (b)
$\textit{Fr}_0 =0.6$
and (c)
$\textit{Fr}_0 = 2.7$
; here
$\hat {\rho } = (\rho -\rho _a)/(\rho _0-\rho _a)$
. The source energy
$\textit{PE}_0$
,
$\textit{KE}_0$
and
$\textit{BPE}_0$
are compared with the amount of advected energy through the boundaries
$\textit{PE}_{Adv}$
,
$\textit{KE}_{Adv}$
, and
$\textit{BPE}_{Adv}$
, as defined in Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995), for (d)
$\textit{Fr}_0=0.6$
, (e)
$\textit{Fr}_0=1.5$
and (f)
$\textit{Fr}_0=2.7$
.
3.3. Hypothetical background state
Similar to the work of Prastowo et al. (Reference Prastowo, Griffiths, Hughes and Hogg2008) and Hughes & Linden (Reference Hughes and Linden2016) with lock-release gravity currents, we define the source background and available potential energies by considering a hypothetical case where mixing effects are neglected. Let us imagine a density current where mixing processes are ignored. Hence, the background state will take the form of a dense layer, of density
$\rho _0$
, lying underneath the ambient layer of density
$\rho _a$
. In this hypothetical representation, intermediate density values
$\rho _0\lt \rho \lt \rho _a$
do not exist. Therefore, if the volume of the supplied fluid per unit width,
$V_0(t) = Q_0t$
, is extended over the bottom boundary along the length of the domain,
$L_x$
, the depth of the dense layer can be expressed as
$ Y^{*}_0 = Q_0t/L_x$
. The background potential energy of this hypothetical case,
$\textit{BPE}_0(t)$
, is expressed as
\begin{equation} \textit{BPE}_0(t) = \frac {1}{2}g L_x \rho _0^{\prime} (Y^{*}_0)^2 = \frac {1}{2} g \rho _0^{\prime} \frac {Q_0^2}{L_x} t^2. \end{equation}
Therefore, the hypothetical available potential energy is simply
$\textit{APE}_0(t) = \textit{PE}_0(t) - \textit{BPE}_0(t)$
. A comparison between the background density profiles of the mixed and hypothetical non-mixed cases is presented in figure 2(b,c). According to Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995), the advected potential energy,
$\textit{PE}_{\textit{Ad}v}$
, kinetic energy,
$\textit{KE}_{\textit{Ad}v}$
, and background potential energy,
$\textit{BPE}_{\textit{Ad}v}$
, across the surface of an open system are defined by the following equations:
\begin{equation} \textit{PE}_{\textit{Ad}v}(t) = - g \int _{0}^{t} \oint _{\mathcal{S}} \rho ^{\prime} y \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{n} {\rm d}S {\rm d}\tau, \end{equation}
\begin{equation} \textit{KE}_{\textit{Ad}v}(t) = - \frac {1}{2} \rho _a \int _{0}^{t} \oint _{\mathcal{S}} (u_x^2 + u_y^2)\boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{n} {\rm d}S {\rm d}\tau, \end{equation}
\begin{equation} \textit{BPE}_{\textit{Ad}v}(t) = - g \int _{0}^{t} \oint _{\mathcal{S}} \psi (\boldsymbol{x},t) \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{n} {\rm d}S {\rm d}\tau, \end{equation}
where
$\psi (\boldsymbol{x},t)$
is a field obtained from integrating the sorted profiles
$\psi (\rho ^{\prime})=\int ^{\rho ^{\prime}}_0y^*(\rho ^{\prime}_1){\rm d}\rho ^{\prime}_1$
with
$\rho ^{\prime}_1$
being an array of the sorted density anomaly profile,
$\mathcal{S}$
represents the surface of the computational domain,
$\boldsymbol{n}$
being a unit vector normal to the surface
$\mathcal{S}$
and oriented outward from the computational domain, and
$\tau$
being a dummy variable that corresponds to the time. Figure 2(d,e,f) shows a good agreement between our defined source energies (
$\textit{PE}_0$
,
$\textit{KE}_0$
and
$\textit{BPE}_0$
) and the corresponding advective energies (
$\textit{PE}_{\textit{Ad}v}$
,
$\textit{KE}_{\textit{Ad}v}$
and
$\textit{BPE}_{\textit{Ad}v}$
) defined by Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995). This agreement justifies the use of our defined quantities and the assumption of the hypothetical state to represent the amount of energy supplied into the system through the open boundaries.
3.4. Parametrisation of the mixing energy
Once the value of the background potential energy
$\textit{BPE}(t)$
and that of the non-mixed case
$\textit{BPE}_0(t)$
are computed, a global mixing energy can be defined as the difference between the two quantities, i.e.
\begin{equation} E_M(t) = \textit{BPE}(t) - \textit{BPE}_0(t). \end{equation}
This is due to the fact that the background potential energy reservoir does not share its energy with the other reservoirs and only the diapycnal mixing and the source contribution cause its increase. Therefore, we define a global dimensionless mixing rate, normalised by the rate of supply of the total energy at the source
$\textit{TE}_0$
, as
\begin{equation} \mu _M = \frac {\text{d}E_M/\text{d}t}{\text{d}(\textit{TE}_0)/\text{d}t}. \end{equation}
(a) The time variation of the dimensionless mixing rate
$\mu _M$
for different source Froude number cases. The solid lines represent the linear relationship
$\overline {\mu _M} = \gamma \hat {t}$
while the dots represent the measured quantities
$\mu _M$
. (b) The mixing parameter
$\gamma$
as a function of the source Froude number
$\textit{Fr}_0$
. The error bars represent the standard error of the coefficient of the linear curves compared with the associated data. (c) and (d) visualise the distribution of local diapycnal mixing rate
$-\kappa g({\text{d}y^*}/{\text{d}\rho })|\boldsymbol{\nabla }\rho |^2$
, introduced by Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995), for
$\textit{Fr}_0 = 0.6$
and
$\textit{Fr}_0 = 2.7$
, respectively. (e) Shows the ratio of mixing energy to the total energy divided over the length of the current
$x_f$
, as a function of time. (f) The mixing energy in
$\mathrm {J\,m^{-1}}$
as a function of the dimensionless time. This figure serves as a comparison between the measured mixing energy (dots) and that predicted with (3.13) (solid lines).
The computed dimensionless mixing rate of the system
$\mu _M$
, defined by (3.11) is shown as a function of the dimensionless time in figure 3(a), for the considered
$\textit{Fr}_0$
cases with filled circle markers. The result of
$\mu _M$
shows a linear dependency on time. The linear curves
$\overline {\mu _M}(\hat {t})$
are plotted with solid lines, with the overbar representing the linearisation and
$\hat {t}=tB_0^{1/3}/h_0$
being the dimensionless time. We define a dimensionless mixing parameter
$\gamma (\textit{Fr}_0)$
to be the coefficient of the linear curves, i.e.
$\gamma (\textit{Fr}_0) \equiv \Delta \overline {\mu _M}/\Delta \hat {t}$
. This mixing parameter is plotted as a function of source Froude number in figure 3(b).
The results reflect a non-monotonic dependency of the mixing rate on
$\textit{Fr}_0$
in the considered range. The mixing parameter decreases in the subcritical range
$0.6 \leqslant \textit{Fr}_0 \leqslant 1$
as
$\textit{Fr}_0$
increases. However, in the range
$1.5 \leqslant \textit{Fr}_0 \leqslant 2.7$
, which lies in the supercritical regime, the mixing parameter increases as
$\textit{Fr}_0$
increases. The non-monotonic behaviour is a direct consequence of the contribution of the vertical derivatives of streamwise velocity and density. The subcritical currents have higher density difference, e.g. for
$\textit{Fr}_0=0.6$
the density difference is
$\rho ^{\prime}_0 = 110 \, \mathrm {kg\,m^{-3}}$
, which causes the current to propagate faster, increasing the vertical derivative of streamwise velocity and generating a strong shearing rate. The produced strong shear induces the ambient fluid to penetrate through the interface into a mixing layer where intermediate density magnitudes, i.e.
$0\lt \rho ^{\prime}\lt 110 \, \mathrm {kg\,m^{-3}}$
for
$\textit{Fr}_0=0.6$
, appear. The available strong buoyancy does not fully suppress the mixing mechanism, but prevents the mixing layer from deepening and penetrating further into the gravity current; however, mixing the fluid from a magnitude
$\rho ^{\prime} = 110 \, \mathrm {kg\,m^{-3}}$
to an intermediate magnitude
$\rho ^{\prime} = 55 \,\mathrm {kg\,m^{-3}}$
requires the investment of high energy. Figure 3(c) illustrates that for
$\textit{Fr}_0=0.6$
, strong local diapycnal mixing is confined to a shear layer along the upper interface, supporting our argument. Therefore, in the range
$0.6 \leqslant \textit{Fr}_0 \leqslant 1$
, as source Froude number decreases, the density difference increases resulting in stronger shearing rate which drives mixing and invests more of the total energy in generating the intermediate magnitudes of density. On the other hand, the decrease of the density difference
$\rho ^{\prime}_0$
as
$\textit{Fr}_0$
increases in the supercritical range
$1.5 \leqslant \textit{Fr}_0 \leqslant 2.7$
reduces the capability of the buoyancy in suppressing the mixing mechanism. As a result, the mixing layer deepens along the hydraulic jump presented at the near-source region and the mixing penetrates into the current. Even though the density difference is small, e.g.
$\rho ^{\prime}_0 = 5 \, \mathrm {kg\,m^{-3}}$
for
$\textit{Fr}_0=2.7$
, and the shearing process is weaker, the mixing mechanism is more distributed inside the current, which increases the efficiency of mixing of these currents. This is well illustrated in figure 3(d) where the diapycnal mixing for
$\textit{Fr}_0=2.7$
is strong in the hydraulic jump at the near-source region, but downstream its magnitude decreases and it becomes more broadly distributed inside the current. We notice that the critical source Froude number with the minimum mixing parameter exists in the range
$1\lt \textit{Fr}_0\leqslant 1.5$
. The results shown in figure 3(a,b) can be formulated as follows:
\begin{equation} \frac {\text{d}E_M}{\text{d}t} = \gamma \frac {\text{d}(\textit{TE}_0)}{\text{d}t} B_0^{1/3} h_0^{-1} t = \gamma \xi _0 t. \end{equation}
Integrating (3.12) over time gives
\begin{equation} E_M(t) = \frac {1}{2}\gamma \xi _0 t^2, \end{equation}
with
$\xi _0$
being a quantity which is determined purely from the source properties:
\begin{equation} \xi _0 = \frac {\text{d}(\textit{TE}_0)}{\text{d}t}B_0^{1/3}h_0^{-1}=\left [\left (\frac {1}{2}\rho _0u_0^2\right)+\left (g\rho ^{\prime}_0\frac {h_0}{2}\right)\right]Q_0B_0^{1/3}h_0^{-1}. \end{equation}
Here,
$E_M$
is found to be proportional to
$t^2$
. By employing this relationship, the mixing energy can be approximated by knowing only the empirical parameter
$\gamma$
and the source properties included in
$\xi _0$
. The mixing rate,
$\text{d}E_M/\text{d}t$
, can be determined with (3.12) and involved in the models which are based on the energy budget formulae. The mixing rate is a linear function of time with a constant coefficient
$\gamma \xi _0$
.
By dividing (3.13) over the length of the current
$x_f$
, determined from (1.1), we obtain
\begin{equation} \frac {E_M}{\textit{TE}_0}\frac {h_0}{x_f}=\frac {\gamma }{2\lambda }. \end{equation}
This formula tells that the ratio of the mixing energy to the total supplied energy over a horizontal distance
$h_0$
is a constant value
$\gamma /(2\lambda)$
, noting that
$\lambda \approx 1$
as determined by Sher & Woods (Reference Sher and Woods2017). This approximation assumes that the mixing is uniformly distributed along the length of the current at sufficiently late times of propagation. The left-hand side of (3.15), divided by
$h_0$
, is plotted as a function of time in figure 3(e). The curves asymptote at late time towards constant values, justifying (3.15). An amount of 0.7 %–1 % of the total supplied energy is invested in mixing along a distance of 1 m of the length of the current. Equation (3.13) is justified for the considered
$\textit{Fr}_0$
cases with
$h_0 = 4\,{\rm cm}$
and
$Re_0 = 5000$
in figure 3(f).
4. Conclusion
The present work brings new ideas to deal with the energetics of continuous gravity currents. We employed the framework of Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995) to examine the effect of mixing on the density profile of the background state. The results of the present paper can be summarised as follows. First, since continuous gravity currents are open systems, it is necessary to estimate the boundary fluxes associated with the variation of
$\textit{BPE}$
and
$\textit{APE}$
. This was done by the use of the hypothetical, non-mixed density profile to compute
$\textit{BPE}_0$
and
$\textit{APE}_0$
as a source contribution of
$\textit{BPE}$
and
$\textit{APE}$
, respectively. The second, and last, result is the parametrisation of mixing energy
$E_M = \textit{BPE} - \textit{BPE}_0$
. The dependency of the mixing rate on the source Froude number
$\textit{Fr}_0$
is non-monotonic in the range
$0.6\leqslant \textit{Fr}_0 \leqslant 2.7$
, as observed in figure 3. In the subcritical range
$0.6 \leqslant \textit{Fr}_0 \leqslant 1$
, the mixing rate increases as
$\textit{Fr}_0$
decreases. However, in the supercritical range
$1.5 \leqslant \textit{Fr}_0 \leqslant 2.7$
, the mixing rate increases as
$\textit{Fr}_0$
increases. The behaviour within the range
$1\lt \textit{Fr}_0\lt 1.5$
is not yet fully understood and requires additional cases for a comprehensive analysis. A formula that predicts the mixing energy
$E_M$
as a function of the squared time
$t^2$
by including empirical constants,
$\gamma$
and source variables, via
$\xi _0$
variable, was introduced. The theory of Phillips (Reference Phillips1972) and Posmentier (Reference Posmentier1977) suggests that the mixing rate tends to zero as
$\textit{Fr}_0 \rightarrow \infty$
since
$\rho ^{\prime}_0 \rightarrow 0$
. Hence, at some range of
$\textit{Fr}_0\gt 2.7$
, we suggest the existence of the transition to a regime in which
$\gamma$
decreases as
$\textit{Fr}_0 \rightarrow \infty$
. The behaviour of mixing at the other extreme, as
$\textit{Fr}_0 \rightarrow 0$
, is still unclear. The stronger stratification could either completely suppress mixing or create a shear layer with high mixing rates at the upper interface, leading to two distinct outcomes.
Further laboratory experiments and numerical simulation cases, with both 3-D and width-averaged fields, are required to extend the study. It is also of interest to check the effect of different source, boundary and initial conditions on the propagation and diagnostic characteristics of continuous gravity current. For instance, we could explore what happens when we include an initial linear ambient stratification or use a time-varying discharge. These ideas form a starting point for future works in investigating continuous gravity currents.
Acknowledgements
High Performance Computing resources were partially supplied by the EXPLOR centre hosted by the University of Lorraine.
Declaration of interests
The authors report no conflict of interest.
