1. Introduction
Gravity currents are flows driven by differences in density, where a denser fluid moves horizontally under the influence of gravity into a lighter ambient fluid, typically across a boundary surface (Simpson Reference Smagorinsky1999). Such currents occur naturally in a variety of settings, including oceanic, atmospheric and geological flows, as well as in engineered systems such as industrial discharge, pollutant dispersion or oil spill dynamics (Simpson Reference Smagorinsky1999). Since the pioneering insights of eminent scholars like von Karman (Reference Klaassen and Peltier1940), Benjamin (Reference Benjamin1968) and Simpson (Reference Simpson and Britter1972), the complex and captivating dynamics of gravity currents has been the subject of extensive investigation, with numerous researchers striving to unravel their intricate mechanisms and underlying physical processes. The foundational theoretical formulations of the gravity current were laid by Benjamin (Reference Benjamin1968), who provided a theoretical description of their propagation and energy balance. This work was followed by laboratory experiments, notably by Simpson (Reference Simpson and Britter1972), Britter & Simpson (Reference Britter and Simpson1978), Huppert (Reference Huppert and Simpson1982), Rottman & Simpson (Reference Shin, Dalziel and Linden1983), Hallworth et al. (Reference Hallworth, Phillips, Huppert and Sparks1993), Klemp, Rotunno & Skamarock (Reference Kokkinos and Prinos1994), Maxworthy et al. (Reference Meiburg, Radhakrishnan and Nasr-Azadani2002), Shin, Dalziel & Linden (Reference Simpson2004), Lowe, Rottman & Linden (Reference Martin, Rakotomalala, Talon and Salin2005), La Rocca et al. (Reference Van Leer2008) and Fragoso, Patterson & Wettlaufer (Reference Fragoso, Patterson and Wettlaufer2013), which collectively established the fundamental dynamics of gravity currents, exploring their propagation, interaction with boundaries, head dynamics and mixing following the sudden release of dense fluids in various environments.
Over time, gravity currents propagating along smooth surfaces have been extensively studied, particularly in idealised configurations like lock-exchange experiments, where a finite volume of higher density fluid, called the lock fluid, and a lower density fluid, called the ambient fluid, are separated by a removable lock that is instantaneously removed to study the dynamics of these buoyancy-driven currents (Shin et al. Reference Simpson2004; Birman, Martin & Meiburg Reference Birman, Martin and Meiburg2005; Etienne, Hopfinger & Saramito Reference Etienne, Hopfinger and Saramito2005; Lowe et al. Reference Martin, Rakotomalala, Talon and Salin2005; Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007; Martin et al. Reference Maxworthy, Leilich, Simpson and Meiburg2011; Balasubramanian & Zhong Reference Balasubramanian and Zhong2018). Huppert & Simpson (Reference von Karman1980) and Rottman & Simpson (Reference Shin, Dalziel and Linden1983) discovered that, after an initial slumping phase, gravity currents from lock-exchange experiments propagate at a constant speed, characterised by a reduced front velocity, which is controlled by the density difference between fluids. Shin et al. (Reference Simpson2004), in their experiments and theory, showed that, at high Reynolds number, dissipation is not important, as previously thought, and developed theories for current speed based on energy conservation of gravity current flows. Lowe et al. (Reference Martin, Rakotomalala, Talon and Salin2005) demonstrated that, in non-Boussinesq lock-exchange flows with two fluids of equal volume, the lighter current behaves as an energy conserving current and the heavier current as the dissipating current having a lower depth with increasing density contrast.
Complementing these experimental insights, significant advancements in numerical simulations have provided a deeper understanding of the detailed flow structures and mixing processes within gravity currents. Härtel et al. (Reference Härtel, Meiburg and Necker2000a , Reference Huppertb ) used direct numerical simulations (DNS) to offer a highly detailed representation of the turbulent interface between the gravity current and the ambient fluid, shedding light on key aspects of the structure of the gravity current, front velocity, flow instability and mixing. Building on this, Cantero et al. (Reference Cantero, Lee, Balachandar and Garcia2007) performed two-dimensional (2-D) and 3-D simulations, and concluded that the 2-D approximations yield satisfactory predictions for the front position and velocity during the slumping phase. However, as the gravity current transitions into the inertial or viscous phase, the flow exhibits inherently 3-D characteristics. In another study using 2-D and 3-D DNS, Cantero et al. (Reference Cantero, Balachandar, García and Bock2008) demonstrated dynamic deformation at the interface between dense and light fluids due to baroclinic torque, leading to the formation of Kelvin–Helmholtz (KH) vortices. These vortices subsequently experience rapid destabilisation that drives their breakdown into smaller-scale turbulent structures, enhancing mixing and energy dissipation. La Rocca et al. (Reference Van Leer2008) obtained satisfactory agreement on their numerical simulation compared with the experimental results and identified two phases of frontal velocities of gravity current. Ooi, Constantinescu & Weber (Reference Pelmard, Norris and Friedrich2007) performed 2-D energy budget analysis of lock-exchange gravity currents and revealed that the increase in front velocities at higher Grashof number, a dimensionless parameter quantifying buoyancy forces relative to viscous forces, was attributable to reduced dissipation and a more efficient conversion of potential energy to kinetic energy. Hallez & Magnaudet (Reference Hallez and Magnaudet2009) employed a force balance analysis to estimate the onset time of the viscous regime for horizontal lock-exchange gravity currents, with strong agreement between theoretical prediction and computational results. More recently, Dai & Huang (Reference Dai and Huang2024) used a high-resolution computational fluid dynamics (CFD) model to simulate lock-exchange gravity currents and provided new insights into the 3-D flow field within the head of the gravity current.
While numerous studies have been dedicated to understanding the dynamics of single-lock-fluid gravity currents in diverse ambient conditions, the investigation of sequentially released multiple-lock-fluid gravity currents remain largely unexplored. The complex interactions arising from the sequential release of multiple lock fluids introduce additional dynamical features, such as interfacial shear, variable entrainment and differential propagation speeds, which have not been thoroughly addressed in the literature. He et al. (Reference Hirt2021) conducted experimental studies on two-layered lock-fluid gravity currents propagating into a lighter ambient and demonstrated that, following the overtaking of the lower layer by the upper layer, the current exhibits a nearly constant-velocity phase before transitioning to a decelerating phase. Dai (Reference Dai2017) conducted experiments on gravity currents where a two-layered dense lock fluid propagated along a horizontal channel filled with ambient fluid of lighter density. The study demonstrated that the inertial phase Froude number in two-layer density-stratified gravity currents varies systematically with the density difference ratio and buoyancy distribution parameter. Despite these insights, the broader implications of multi-phase gravity currents, particularly their evolution in natural and engineering contexts, remain insufficiently explored. In recent years, the study of colliding gravity currents has gained significant attention, with a strong focus on high-fidelity numerical approaches such as large eddy simulation (LES), which provides a robust framework for capturing the intricate turbulence dynamics that complements laboratory-scale observations (van der Wiel et al. Reference Wu2017; Cafaro & Rooney Reference Cafaro and Rooney2018; Zhong, Hussain & Fernando Reference Zhou and Venayagamoorthy2018; Okon, Zhong & He Reference Ooi, Constantinescu and Weber2021; Kokkinos & Prinos Reference Kokkinos and Prinos2022, Reference Kyrousi, Leonardi, Roman, Armenio, Zanello, Zordan, Juez and Falcomer2023; Dai, Huang & Wu Reference Dai, Huang and Wu2023; Wu Reference Zhong, Hussain and Fernando2024). In addition, several studies have demonstrated that 2-D LES can successfully capture the bulk propagation, frontal structure and dominant KH billow dynamics of gravity currents, particularly during the slumping phase and at moderate Reynolds numbers, thereby providing a useful complementary framework to fully 3-D approaches (Patterson et al. Reference Peltler, Hallé and Clark2005; Ooi et al. Reference Pelmard, Norris and Friedrich2007; Nourazar & Safavi Reference Okon, Zhong and He2017; Pelmard, Norris & Friedrich Reference Peng and Lee2018). However, a comprehensive examination of sequentially released multiple-lock-fluid gravity currents is still lacking, necessitating investigations to unravel their governing mechanisms and broader environmental and industrial significance. In most natural and industrial systems, such as estuaries, coastal regions, reservoirs, cold pools, industrial pollutants, pyroclastic flows, oil and chemical spills and wastewater discharge into rivers and oceans, multiple gravity currents often occur in succession. These sequential currents can interact in complex ways, influencing their speed, structure and mixing efficiency. This study seeks to address this critical gap by investigating the interaction between two sequentially released gravity currents, and developing insights which can later be applied to larger-scale, real-world environmental and industrial applications.
Figure 1. Schematic of two gravity currents generated by sequential release of two lock fluids having
$(\rho _2\gt \rho _1\gt \rho _a)$
(not to scale), where
$\rho_1$
,
$\rho_2$
, and
$\rho_a$
are the densities of lighter lock fluid, heavier lock fluid, and ambient, respectively. The acceleration due to gravity is denoted by
$g$
.
The schematic representation of the two-lock-fluid gravity current configuration in the x-z plane is shown in figure 1. In this set-up, two lock fluids with densities
$\rho _1$
and
$\rho _2$
are sequentially released into an ambient fluid of density
$\rho _a$
with
$\rho _2\gt \rho _1\gt \rho _a$
. This results in the formation of two gravity currents, designated as
$\textit{GC}_1$
and
$\textit{GC}_2$
, propagating with velocities
$U_1$
and
$U_2$
, respectively. For small changes between
$\rho _1$
,
$\rho _2$
and
$\rho _a$
$({\lt}10\,\%)$
and a sufficiently shorter delay in the release time between the two lock fluids,
$t_R$
, it is observed that
$\textit{GC}_2$
eventually overtakes
$\textit{GC}_1$
at a specific time, called the overtaking time,
$t_O$
. In this study, overtaking is considered to have occurred when the non-dimensional thickness of the mixed fluid layer between the ambient fluid and the denser lock fluid at the nose of the advancing current becomes less than 0.125 (
$\Delta x/H \lt 0.125$
, where
$H$
is the total height of the ambient fluid). This criterion reflects the point at which the heavier fluid has effectively accumulated at the head of the current and has overtaken the lighter fluid. As gravity currents propagate, the interaction between fluids induces interfacial instabilities as a result of the shear produced by velocity gradients at the interfaces of the three fluids. These instabilities lead to the formation of KH billows, which play a significant role in enhancing mixing between the fluids by facilitating momentum and mass transfer (Thorpe Reference van der Wiel, Gille, Llewellyn Smith, Linden and Cenedese1973; Britter & Simpson Reference Britter and Simpson1978; Sha, Kawamura & Ueda Reference Simpson1991; Klemp et al. Reference Kokkinos and Prinos1994; Peng & Lee Reference Pokrajac, Venuleo and Franca2010; Dai & Huang Reference Dai and Huang2024). The generated billows disrupt the initial density stratification, leading to a complex entrainment and mixing dynamics. The main focus of this study is to compare the dynamic behaviour of two lock fluids with conventional single-lock-fluid gravity currents. This numerical set-up is particularly valuable for investigating the fundamental mechanisms governing the propagation, interaction and mixing dynamics of gravity currents involving multiple lock fluids released sequentially.
The model used for validation study, as shown in figure 2(a), consists of two fluids of densities
$\rho _1$
and
$\rho _2$
(
$\rho _1\gt \rho _2$
), each having a full depth,
$H$
, of 0.2 m and a length of 0.91 m. In the present study, for comparison of three cases of gravity currents, two models are used, as shown in figures 2(b) and 2(c). The total depth of the fluid,
$H$
, is 0.2 m and the length of the domain with ambient fluid,
$L_a$
, is 4 m. In the first case, the lock fluid consists of two fluids of densities
$\rho _1$
and
$\rho _2$
(
$\rho _2\gt \rho _1$
), each having a streamwise length of 0.1 m and height of 0.2 m. The lock fluids are released simultaneously into the ambient fluid of density
$\rho _a$
. The second case consists of the same configuration as that of the first case, but the lock fluids are released sequentially with a delay in the release time of
$t_R$
. The third case consists of a single full depth lock fluid having length of 0.2 m and the average density of the previous two cases,
$(\rho _1+\rho _2)/2$
. The idea here is to compare the two-lock-fluid cases of simultaneous release (case I) and sequential release (case II) with the well-explored single full depth lock release case (case III). The density values used in all simulations are presented in table 1, and an initial kinematic viscosity,
$\nu$
, of
$1 \times 10^{-6}$
$\rm m^2\,s^{-1}$
was applied for all fluids. Dimensionless fluid densities are expressed as
$\gamma _1={\rho _a}/{\rho _1}$
,
$\gamma _2={\rho _a}/{\rho _2}$
and
$\overline {\gamma }=2\rho _a/(\rho _1+\rho _2)$
. Reduced gravity, which is the driving force of gravity currents, is defined as
$g'=g(\overline {\gamma }^{-1}-1)$
. The average density ratio,
$\overline {\gamma }$
, is used to consistently capture the combined influence of
$\rho _1$
and
$\rho _2$
in all three cases, providing a coherent measure of reduced gravity. The characteristic length and velocity scales used in the study are the initial depth of the lock fluid,
$H$
, and the buoyancy velocity of the flow,
$u_b=\sqrt {g'H}$
, respectively. Time, position and velocity are non-dimensionalised as
$t^*=t(u_b/H)$
,
$x^*=x/H$
and
$u^*=u/u_b$
, respectively. The non-dimensionalised release delay time between the two lock fluids and the non-dimensionalised overtaking time are defined as
$t_R^*=t_R(u_b/H)$
and
$t_O^*=t_O(u_b/H)$
, respectively. The bulk Reynolds number, the local Reynolds number and the Froude number of the flow field are defined as
$Re_b=u_bH/\nu$
,
$Re_l=uh/\nu$
and
$Fr=u/u_b$
, respectively, where
$\nu$
is the kinematic viscosity,
$h$
is the height of the current head and
$u$
represents the instantaneous velocity of the current head.
Table 1. List of all the numerical simulation performed in this study. The study comprises 5 simulations for validation of the model, three simulations for comparison of the three cases of gravity current and 100 simulations for the parametric study
$(\textit{PS})$
of overtaking time and delayed release time (50 each for
$\textit{PS}_1$
and
$\textit{PS}_2$
).
Figure 2. Schematic of the numerical models used in the study for (a) validation, (b) case studies I and II and (c) case study III (not to scale).
In what follows, the governing equations of the gravity current flows, the methodologies adopted in the study and the validation of the numerical model used in the study are presented in § 2. An analysis of the three cases outlined, consisting of two lock fluids that are released simultaneously (case I), sequentially (case II) and a single lock fluid having the average density of the two lock fluids of the previous cases (case III) is presented in § 3 to compare the behaviour of gravity currents generated in these three cases. In § 4, an analysis of several cases of sequential release of two lock fluids is presented, exploring the relationship between the non-dimensionalised release delay time,
$t^*_R$
, non-dimensionalised time of overtaking by the heavier lock fluid,
$t^*_O$
, and the ratio of densities of ambient and lock fluids,
$\gamma _1$
and
$\gamma _2$
. Finally, concluding remarks are presented in § 5.
2. Methodology and model set-up
2.1. Governing equations
Natural gravity currents are buoyancy driven flows that largely comply with the Boussinesq approximation of neglecting density variations except in the buoyancy term (Rottman & Simpson Reference Shin, Dalziel and Linden1983; Härtel et al. Reference Huppert2000b
; Ooi et al. Reference Pelmard, Norris and Friedrich2007; Ungarish Reference Wildt, Hauer, Habersack and Tritthart2009). In this study, we employed the LES technique, which has been extensively utilised in gravity current research to achieve a balance between flow resolution accuracy and computational efficiency (Patterson et al. Reference Peltler, Hallé and Clark2005; Ooi, Constantinescu & Weber Reference Patterson, Simpson, Dalziel and Nikiforakis2009; An, Julien & Venayagamoorthy Reference An, Julien and Venayagamoorthy2012; Constantinescu Reference Constantinescu2014; Meiburg, Radhakrishnan & Nasr-Azadani Reference Meneveau and Katz2015; Nourazar & Safavi Reference Okon, Zhong and He2017; Zhou & Venayagamoorthy Reference Zhou, Cenedese, Williams, Ball, Venayagamoorthy and Nokes2017; Kyrousi et al. Reference La Rocca, Adduce, Sciortino and Pinzon2018; Emami et al. Reference Emami, Mousavi, Hosseini, Fouladfar and Mohammadian2020). Large eddy simulation effectively captures the unsteady energy-containing motions of large eddies while using a turbulence closure model to represent the effects of unresolved subgrid-scale (SGS) motions (Meneveau & Katz Reference Moitro, Dammati and Poludnenko2000; Pope Reference Sha, Kawamura and Ueda2001; Moitro, Dammati & Poludnenko Reference Nasr-Azadani and Meiburg2024). This method involves applying a spatial decomposition technique known as ‘filtering’, which separates the velocity field into a filtered component,
$\overline {u}$
, and a residual component
$u'$
. The effect of all scales of motion below the filter width manifests as an additional SGS stress term,
$\tau _{\textit{ij}} ^s$
, in the momentum equation and as an SGS flux term,
$\chi ^s$
, in the density transport equation. The LES-filtered, incompressible Navier–Stokes equations with the Boussinesq approximation can be expressed in tensor notation as follows:
\begin{equation} \dfrac {\partial \overline {u_i}}{\partial t}+\dfrac {\partial \overline {(u_iu_{\!j})}}{\partial x_{\!j}}=-\dfrac {1}{\rho _o} \dfrac {\partial \overline {p}}{\partial x_i}-\dfrac {\overline {\rho }}{\rho _o}g\delta _{i3}+\nu \dfrac {\partial ^2 \overline {u_i}}{\partial x_{\!j}\partial x_{\!j}}-\dfrac {\partial \tau _{\textit{ij}} ^s}{\partial x_{\!j}}, \end{equation}
with the filtered continuity equation
\begin{equation} \dfrac {\partial \overline {u_i}}{\partial x_i}=0, \end{equation}
and the filtered density transport equation
\begin{equation} \dfrac {\partial \overline {\rho }}{\partial t}+\dfrac {\partial \overline {(\rho u_{\!j})}}{\partial x_{\!j}}=\kappa \dfrac {\partial ^2 \overline {\rho }}{\partial x^2_{\!j}}-\dfrac {\partial \chi ^s}{\partial x_{\!j}}, \end{equation}
where
$\overline {u_i}$
and
$\overline {u_{\!j}}$
are the filtered velocity components,
$t$
is time,
$\rho _o$
is the reference density of the fluid,
$\overline {p}$
is the filtered pressure field,
$g$
is acceleration due to gravity,
$\overline {\rho }$
is the filtered density field,
$\delta _{i3}$
is the Kronecker delta with gravity in the vertical direction (
$j=3$
),
$\nu$
is the kinematic viscosity,
$\kappa$
is the molecular diffusivity and
$\tau _{\textit{ij}} ^s$
and
$\chi ^s$
are the SGS stress tensor and SGS scalar flux, respectively.
The accuracy of LES models relies on the effectiveness of SGS closure models in capturing SGS motions. The widely used Smagorinsky LES model (Smagorinsky Reference Thorpe1963) approximates SGS effects using a linear turbulent eddy viscosity,
$\nu _t$
. The SGS tensor,
$\tau _{\textit{ij}} ^s$
, and the SGS scalar flux,
$\chi ^s$
, are given by
\begin{align} \tau _{\textit{ij}}^s &= \overline {u_i u_{\!j}} - \overline {u_i} \, \overline {u_{\!j}}, \end{align}
\begin{align} \chi ^s &= \overline {\rho u_{\!j}} - \overline {\rho } \, \overline {u_{\!j}}, \end{align}
where the overbar denotes the filtering operation. To close these terms, the Smagorinsky model assumes a linear relationship between the SGS stress tensor,
$\tau _{\textit{ij}} ^s$
, and the filtered strain rate tensor,
$S_{\textit{ij}}$
,
\begin{equation} \tau _{\textit{ij}}^s = -2\nu _tS_{\textit{ij}}. \end{equation}
The turbulent eddy viscosity is defined as
\begin{equation} \nu _t=(C_s\varDelta )^2|S|, \end{equation}
where
$C_s$
is the Smagorinsky constant, typically ranging from 0.1 to 0.2,
$\Delta$
is the filter width (often taken as the mean grid size,
$\varDelta =(\delta x \delta y \delta z)^{1/3}$
) and
$|S|$
is the characteristic strain rate given as
\begin{equation} |S|=\sqrt {2 S_{\textit{ij}} S_{\textit{ij}}}. \end{equation}
Similarly, the SGS scalar flux is modelled using the gradient diffusion hypothesis, defined as
\begin{equation} \chi ^s=-\frac {\nu _t}{Sc_t} \frac {\partial \rho }{\partial x_{\!j}}, \end{equation}
where
$Sc_t$
is the turbulent Schmidt number, typically assumed to be unity for turbulent gravity currents (Härtel et al. Reference Huppert2000b
; Maxworthy et al. Reference Meiburg, Radhakrishnan and Nasr-Azadani2002; Necker et al. Reference Nourazar and Safavi2005; Ooi et al. Reference Patterson, Simpson, Dalziel and Nikiforakis2009; Nasr-Azadani & Meiburg Reference Necker, Härtel, Kleiser and Meiburg2011, Reference Necker, Härtel, Kleiser and Meiburg2014; Pelmard et al. Reference Peng and Lee2018). This formulation provides closure for the SGS terms, allowing for the computation of unresolved turbulent fluxes in the LES framework.
2.2. Model set-up
In this study, 2-D simulations were performed using FLOW-3D with a grid resolution of 1–2 mm. Due to the 3-D nature of FLOW-3D simulations, only one computational cell was used in the spanwise direction, along which free slip boundary condition was imposed. All remaining boundaries were treated as no slip to accurately represent the physical configuration. FLOW-3D, a finite-volume CFD software developed by Flow Science, Inc., has been widely utilised in laboratory-scale gravity current investigations due to its demonstrated capability in resolving essential turbulent flow structures and accurately capturing the entrainment dynamics (An et al. Reference An, Julien and Venayagamoorthy2012; Zhou et al. Reference Zhou and Venayagamoorthy2017; Zhou & Venayagamoorthy Reference Zhou, Cenedese, Williams, Ball, Venayagamoorthy and Nokes2017; Stancanelli, Musumeci & Foti Reference Ungarish2018; Zhou & Venayagamoorthy Reference Zhou and Venayagamoorthy2020; Wildt et al. Reference Yao2020). The volume-of-fluid method of FLOW-3D is used to track fluid–fluid interfaces and free surfaces (Hirt & Nichols Reference Härtel, Carlsson and Thunblom1981). Additionally, its FAVOR (fractional area-volume obstacle representation) technique allows for flexible geometric configuration without changing the mesh (Hirt Reference Hirt and Nichols1993). With second-order monotonicity-preserving upwind differencing (Van Leer Reference Lowe, Rottman and Linden1977) for advective terms and efficient implicit pressure solution using the generalised minimum residual algorithm (Yao Reference Zhou, Cenedese, Williams, Ball, Venayagamoorthy and Nokes2004), FLOW-3D is well suited for high-fidelity LES studies of gravity currents, capturing the key propagation dynamics. The absence of a full 3-D domain may influence the energy cascade and mixing dynamics. Furthermore, 2-D LES inherently suppresses spanwise variability and vortex stretching, thereby limiting the growth of secondary instabilities and allowing KH billows to remain more coherent and longer lived than in fully 3-D flows (Thorpe Reference van der Wiel, Gille, Llewellyn Smith, Linden and Cenedese1973; Peltler, Hallé & Clark Reference Piomelli1978; Corcos & Sherman Reference Corcos and Sherman1984; Klaassen & Peltier Reference Klemp, Rotunno and Skamarock1991; Caulfield Reference Caulfield2021). These effects are further reinforced by the Smagorinsky SGS model, which introduces additional eddy viscosity and is known to damp small-scale fluctuations in stratified shear layers (Piomelli Reference Pope1999; Meneveau & Katz Reference Moitro, Dammati and Poludnenko2000). Despite these limitations, previous studies show that 2-D LES can reliably capture the dominant frontal dynamics and large-scale billow structures of gravity currents during the slumping phase, particularly under moderate Reynolds numbers where small-scale turbulence plays a secondary role (Necker et al. Reference Nogueira, Adduce, Alves and Franca2002; Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007). Thus, while not fully resolving the 3-D instability cascade, the present 2-D framework remains well suited for examining the bulk propagation characteristics of the currents. Furthermore, to reinforce the validity of the 2-D simulations within the scope of this study, a small number of 3-D simulations were performed and compared with the corresponding 2-D results, revealing close agreement in terms of bulk propagation and the overall flow dynamics (results shown in Appendix A).
While the present study does not rely on a shallow-water theoretical framework, it is important to recognise that the initial lock lengths used in the simulations
$(L_1=L_2=0.5H)$
are relatively short and do not satisfy the
$L/H\gt \gt 1$
criterion typically required for hydrostatic balance in shallow-water studies. Consequently, vertical accelerations may arise shortly after the release of the lock fluids, leading to transient departures from a hydrostatic pressure distribution. These vertical motions, although short lived, can influence the early-stage evolution of the gravity currents by inducing interfacial deformation and localised mixing prior to the onset of quasi-steady propagation. Table 1 shows the details of the simulations carried out in this study.
Figure 3. Qualitative validation of the numerical model of this study: coloured snapshots (b), (d), ( f), (h), (j) and (l)) with the experimental results of Lowe et al. (Reference Martin, Rakotomalala, Talon and Salin2005); grey coloured snapshots (a), (c), (e), (g), (i) and (k)) at non-dimensionalised time
$t^*=1.2, 2.3, 3.9, 4.7, 5.9, 7.0$
. The gravity current was generated by the instantaneous release of two lock fluids of densities
$\rho _1=1005.10$
and
$\rho _2=998.06$
$\rm kg\,m^{-3}$
.
2.3. Model validation
The experimental work of Lowe et al. (Reference Martin, Rakotomalala, Talon and Salin2005) was used to validate the numerical model in this study. The validation model, as shown in figure 2(a), was a replica of the experimental model of Lowe et al. (Reference Martin, Rakotomalala, Talon and Salin2005). First, the LES results of the evolution of the gravity current with densities
$\rho _1=1005.10$
and
$\rho _2=998.06$
$\rm kg\,m^{-3}$
, were qualitatively compared with the experimental results of Lowe et al. (Reference Martin, Rakotomalala, Talon and Salin2005), as shown in figure 3. Six density plots were compared at non-dimensionalised times
$t^*$
= 1.2, 2.3, 3.9, 4.7, 5.9 and 7.0. The results show excellent agreement between the two studies based on the positions of the heavier gravity current head that propagates at the bottom and the lighter gravity current head that propagates at the top. The results of this study were able to clearly show the formation of KH billows at the interface of the two fluids due to velocity gradients, as shown in figure 3. Five sets of numerical simulations were performed for the same ambient fluid of density
$\rho _2=998.06$
$\rm kg\,m^{-3}$
, as shown in table 2. The front velocities for the heavier current,
$u_H$
, and the lighter current,
$u_L$
, were evaluated for each simulation and demonstrated a high degree of accuracy, exhibiting strong agreement with the corresponding experimental results by Lowe et al. (Reference Martin, Rakotomalala, Talon and Salin2005), as shown in table 2, having a less than 5 % error in the heavier current velocity,
$E_H$
, and the lighter current velocity,
$E_L$
. These results demonstrate the accuracy and reliability of the numerical model of this study in capturing the key dynamics of the gravity current, including front propagation and KH billow formation.
Table 2. Comparison of the non-dimensionalised front velocities of heavier density current
$(u_H^*=u_H/\sqrt {\textit{gH}({\rho} _1-{\rho} _2)/\rho _2})$
and lighter density current
$(u_L^*=u_L/\sqrt {\textit{gH}({\rho} _1-\rho _2)/{\rho} _1})$
for LES results of this study and experimental results by Lowe et al. (Reference Martin, Rakotomalala, Talon and Salin2005). The error percentages in non-dimensionalised front velocities for heavier density current,
$E_H$
, and lighter density current,
$E_L$
, are based on the current LES study results and experimental results by Lowe et al. (Reference Martin, Rakotomalala, Talon and Salin2005).
3. Lock-fluid release case study
3.1. Interaction of lock fluid with ambient
The three cases of lock-exchange gravity currents were simulated using consistent fluid properties having
$\rho _1 = 1008.28\ \mathrm{}$
,
$\rho _2 = 1013.35\ \mathrm{}$
and
$\rho _a = 998.2\ \mathrm{\rm kg\,m^{-3}}$
, with a release delay of
$t_R = 2\ \mathrm{s}$
for the heavier lock fluid in case II. In all three cases, upon the instantaneous release of the lock fluid into the ambient medium, a strong velocity shear at the interface initiates the formation of distinct vortex structures behind the head of the gravity current, as shown in the density plots at different times presented in figure 4. This phenomenon is attributed to the KH instability that induces mixing and forms a layer of mixed density fluid (Simpson Reference Simpson and Britter1972; Simpson & Britter Reference Stancanelli, Musumeci and Foti1979). As these KH-induced vortices develop, they undergo a cascade of energy, transitioning from coherent large-scale vortices to finer turbulent eddies. This promotes significant entrainment of the ambient fluid, thereby increasing the effective height of the gravity current and simultaneously decreasing the reduced gravity and, thus, the front velocity has minimal impact (Britter & Simpson Reference Britter and Simpson1978; Nogueira et al. Reference Odier, Chen and Ecke2014).
The intricate interplay between the lock fluid and the ambient fluid induces substantial energy dissipation, concentrated within the vortical structures formed at the interface behind the gravity current head. The dissipation associated with vortical structures contributes to local momentum loss in the head and body of the gravity current (Balasubramanian & Zhong Reference Balasubramanian and Zhong2018). However, the primary mechanism of head deceleration, consistent with Rottman & Simpson (Reference Shin, Dalziel and Linden1983), is attributed to the overtaking of the front by the rear bore during the later stages of propagation. We will next analyse the evolution of the gravity current front velocity for the three cases in this study.
Figure 4. Evolution of gravity current for case I, case II and case III at different non-dimensionalised times
$t^*=1.575, 4.724, 7.873, 11.023, 14.172$
. The vertical dotted lines in white are the locations of density and velocity profiles presented in figure 7.
3.2. Propagation of gravity current head
The temporal evolution of a gravity current is classically characterised by three distinct phases: the slumping phase, the self-similar buoyancy-inertia phase and the self-similar buoyancy-viscosity phase (Simpson Reference Simpson and Britter1972; Huppert & Simpson Reference von Karman1980). In the present study, the simulation duration of 30 s is limited to the slumping phase of the gravity current, as defined by the slumping phase criteria established by Huppert & Simpson (Reference von Karman1980). This phase is defined by a relatively constant front velocity, driven primarily by the initial momentum, generated by the conversion of total initial potential energy of the released lock fluid into kinetic energy, and a strong velocity shear at the interface that develops vortical structures due to KH instabilities, resulting in minimal deceleration (Huppert & Simpson Reference von Karman1980; Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007; Gadal et al. Reference Gadal, Mercier, Rastello and Lacaze2023). Figure 4 displays the density field at various times (
$t^*=1.575, 4.724, 7.873, 11.023, 14.172$
) allowing for a qualitative comparison of the positions of gravity current heads and the development of billows and vortical structures at the interface of the gravity current and the ambient for all three cases. In figures 4(a), 4(b) and 4(c), cases I and II demonstrate that the lighter lock fluid, initially in direct contact with the ambient fluid, is primarily responsible for initial entrainment and mixing through the formation of KH billows, while the denser fluid, represented in red, remains comparatively intact for a prolonged duration. This stratified arrangement influences the propagation dynamics, as the lighter fluid, initially concentrated at the head region, dictates the initial front velocity. Consequently, before overtaking of the lighter lock fluid by the heavier lock fluid occurs, the gravity current in case I exhibits a lower front velocity compared with case III, as indicated by the position of the current heads in figure 4(a). However, as the heavier lock fluid in case I eventually overtakes the lighter lock fluid, its higher density contrast and accumulation of the heavier lock fluid at the head region enhance the gravitational driving force, leading to an acceleration phase that ultimately contributes to a higher sustained front velocity post-overtake, as observed in the gravity current head position shown in figure 5. This distinction suggests that the velocity of the gravity current is governed by the evolving composition of the current head, as variations in density and momentum distribution influence the gravitational driving force and entrainment dynamics (Odier, Chen & Ecke Reference Ooi, Constantinescu and Weber2012; Pokrajac, Venuleo & Franca Reference Rottman and Simpson2018).
Figure 5. Evolution of the gravity current head for three case studies. The blue and red diamond shapes indicates the non-dimensionalised time of overtaking by heavier lock fluid,
$t_O^*$
, for cases I and II, respectively.
In comparing case I and case II, the similarity of the vortical structures suggests a similar flow dynamics, yet case II exhibits the lowest front velocity, as shown in figure 4. This outcome is attributed to the delayed release of the heavier fluid, which plays a significant role in the overall front velocity. The analysis reveals that the denser fluid, initially positioned behind the lighter fluid with a higher initial propagation speed, eventually overtakes the lighter fluid (as shown in figures 4
b for case I and 4
d for case II) before it begins to participate in active mixing and entrainment with the ambient. The delayed overtaking event in case II, resulting from the delayed release of the denser lock fluid, emphasises the importance of the release delay time (
$t_R$
) and the subsequent overtaking time (
$t_O$
) in determining the front head velocity. Therefore,
$t_R$
and
$t_O$
are critical parameters, particularly for delayed release cases such as case II, as they influence the propagation and ultimately the entrainment and mixing occurring in a multi-lock-fluid gravity current.
The temporal progression of the gravity current head for three cases, along with the overtaking time,
$t^*_O$
, for case I and case II are presented in figure 5. The results illustrate that the head position is generally most advanced in case I after the overtaking by heavier lock fluid takes place, and the head position is lagging in case II, which is also evident in figure 4. This variation in head position can be primarily attributed to the parameters
$t_R$
and
$t_O$
, as discussed previously. A slight increase in the current head velocity can be observed from the change in slope after the overtaking event for case I and case II, as shown in figure 5. Such a shift was not reported in the previous study of multiple gravity currents (Dai Reference Dai2017), which considered two layers of lock fluids. This further reinforces the critical role of the evolving composition of the gravity current head in governing the propagation dynamics of sequentially released multi-lock-fluid gravity currents.
Figure 6. Comparison of the dynamics of the gravity current for case II (a, c, e) and case III (b, d, f) before (a, b,
$t^*=8.03$
), during (c, d,
$t^*=9.92$
) and after (e, f,
$t^*=11.60$
) overtaking by the heavier lock fluid. The zoomed in plots of the gravity currents show the density fields (in coloured plots) with the superimposed velocity vectors (shown by the arrows).
In figure 6, we present the velocity and density fields at three crucial moments: before (
$t^*$
= 8.03), during (
$t^*$
= 9.92) and after (
$t^*$
= 11.60) the overtaking phase. The propagation dynamics of gravity currents in a two-lock-fluid configuration with sequential release (case II) and the classical single-lock-fluid release (case III) reveal fundamental differences, primarily driven by the overtaking of the lighter lock fluid by the heavier lock fluid in case II. This overtaking is marked by an apparent concentration of the heavier fluid at the current’s head, with a thin layer of the lighter lock fluid situated between the heavier lock fluid and ambient fluids, as shown in figure 6(c). This shift in the composition of the current head leads to notable changes in the dynamics of current propagation, as discussed previously in § 3.2, which affects the general evolution of the flow, the efficiency of entrainment and the density-driven turbulence of the system. Before overtaking occurs, the front in case II advances more slowly compared with case III, as shown in figures 4 and 6, as the lighter lock fluid initially occupies the leading position. However, upon overtaking, the heavier lock fluid accumulates at the head, increasing the gravitational driving force and leading to a subsequent acceleration of the current. This change in the density field in the gravity current’s head also modifies the velocity field at the current’s head as the heavier lock fluid reaches the head region, as depicted in figures 6(a), 6(c) and 6(e). This acceleration after the accumulation of heavier lock fluid at the head of the current after overtaking is made evident by the decrease in the difference of the head front position between figures 6(a) and 6(b) (before overtaking) and figures 6(e) and 6( f) (after overtaking). The velocity vector plots in figure 6 show that the streamwise velocity dominates the head region of the current, enhancing propagation of the current, vortical structures with significant components of streamwise and wall-normal velocities dominate the body of the current enhancing KH instabilities and entrainment in this region and diminished velocities are observed at the tail.
Furthermore, from figures 4 and 6, it is observed that case III exhibits more pronounced KH instabilities compared with case I and case II. The density and velocity fields in figure 6 reveal that case III features more well-developed billows and vortical structures at the interface between the gravity current and the ambient fluid. These instabilities are particularly prominent along the shear layer, indicating stronger turbulence and mixing. In contrast, case II displays relatively fewer and less coherent KH billows, suggesting a more stable interface with reduced turbulent entrainment compared with case III. This difference can be attributed to the density distribution and velocity gradients in the two cases. In case III, the single lock release introduces a stronger and more uniform shear layer at the current–ambient interface, which promotes the growth of KH instabilities. Meanwhile, in case II, the interaction between the two lock fluids before and during overtaking modifies the velocity gradient, potentially suppressing the development of large-scale KH billows in comparison with case III. The enhanced instabilities in case III contribute to higher mixing and entrainment, which significantly influence the evolution and propagation of the gravity current.
Figure 7. Density (left), streamwise velocity (centre) and wall-normal velocity (right) profiles at various sections as shown in figure 4 (
$P_a$
to
$P_h$
are shown as plots (a) to (h)). In all the plots, blue line is for case I, red line indicates case II and the green line represents case III. Plots (a), (d) and ( f) are located at the head, plots (b), (e) and (h) are located at the body and plots (c) and (g) are located at the tail of the gravity current.
3.3. Density and velocity profiles
The head, body and tail regions of the gravity current exhibit distinct dynamical behaviours that collectively influence the entrainment and mixing processes (Klemp et al. Reference Kokkinos and Prinos1994; Nogueira et al. Reference Odier, Chen and Ecke2014). The head, characterised by a distinct bulging structure at the front of the current, is the most dynamically active region, accumulating kinetic energy that drives the current forward (Nogueira et al. Reference Odier, Chen and Ecke2014; Castro-Folker, Grace & Stastna Reference Castro-Folker, Grace and Stastna2023). At the interface between the head and the ambient fluid, sharp density and velocity gradients generate intense shear-induced turbulence, enhancing entrainment through KH instabilities. These instabilities promote mixing by breaking the interface into smaller turbulent eddies, facilitating the exchange of mass and momentum between the current and the surrounding fluid (Benjamin Reference Benjamin1968). Figure 7 shows the profiles of density, streamwise and wall-normal velocities for the three cases at different times and locations (shown in figure 4). The vertical density profiles shown in figure 7 provide valuable insights into the structural evolution and mixing dynamics of the gravity currents in the three release scenarios. These profiles effectively delineate the height of the gravity current and thickness of mixing layer by capturing the vertical stratification and deformation of fluid density. A comparative analysis reveals that, with the exception of figures 7(b) and 7(c), case II (red) consistently exhibits the lowest height of the gravity current across all locations. This reduced height results from the delayed propagation of the denser lock fluid in the sequential release configuration, which limits the vertical development of the current and suppresses early interfacial growth. The delayed release reduces the immediate kinetic energy input and delays the onset of dynamic interaction with the ambient fluid, thereby inhibiting vertical expansion. In contrast, cases I (blue) and III (green) generally exhibit greater current heights, attributed to the more rapid initial propagation and more immediate entrainment of ambient fluid. This is a consequence of the simultaneous release of the entire lock-fluid volume, which generates a stronger gravity-driven front and promotes the development of a thicker shear interface from the outset.
In figure 7, the head region of gravity current is visualised in profiles (a), (d) and ( f), the body of the gravity current is presented in profiles (b), (e) and (h) and the tail of the gravity current is represented by profiles (c) and (g). The head region, as shown in figures 7(a), 7(d) and 7( f), displays sharp changes in densities for all three cases, representing the distinct layers of three fluids. This indicates that the layers of fluids are not actively mixing, involving minimum entrainment in the head region of the current. The streamwise velocity plots at the head region, as shown in figures 7(a), 7(d) and 7( f), indicate a sharp transition between the current velocity to the ambient velocity with nearly horizontal profiles for all three cases. Similarly, the wall-normal velocities in the head region, which are approximately one order magnitude less than the streamwise velocity in all cases, show relatively quick transitions at the interface of the lock fluid and the ambient, indicating minimal shear instabilities, as shown in figures 7(a), 7(d) and 7( f). Thus, the velocity and density profiles depicted in figures 7(a), 7(d) and 7( f) substantiate that the head region of the gravity current is characterised by the presence of thin shear layers and exhibits minimal interfacial mixing, as further corroborated by figure 6. This observation is consistent with findings reported in previous studies, which have similarly identified reduced mixing in the head region of gravity currents (Hallworth et al. Reference Hallworth, Phillips, Huppert and Sparks1993, Reference Hallworth, Huppert, Phillips and Sparks1996; Fragoso et al. Reference Fragoso, Patterson and Wettlaufer2013). The wall-normal velocity profiles at the interface of the gravity current and ambient fluid in figures 7(d) and 7( f) exhibit notable differences across the three cases. In particular, case II shows significantly higher wall-normal velocity magnitudes compared with cases I and III. This enhancement of vertical motion increases interfacial shear and suggests a greater potential for KH instability formation. This behaviour arises due to the slower propagation of the current in case II, resulting from its delayed release of heavier lock fluid. Consequently, the velocity profile for case II is captured behind the nose of the current, where vertical fluctuations are more intense. In contrast, the profiles in cases I and III are extracted at the crown of the gravity current head, as shown in figure 4, where vertical fluctuations are weaker. This positional difference explains the observed variation in wall-normal velocity and the corresponding instability characteristics among the cases.
Following the head, the body of the gravity current transitions into a fully developed turbulent regime dominated by KH billows that further enhance mixing and entrainment (Härtel et al. Reference Huppert2000b ). These billows, evident in the fluctuating density profiles in figures 7(b), 7(e) and 7(h), increase the interfacial area between dense and ambient fluids. The gradual transition of the gravity current streamwise velocity to the ambient fluid streamwise velocity, as shown in figures 7(b), 7(e) and 7(h), indicates significantly thicker shear layers than in the head region, which is also evident in the wall-normal velocity profiles in the same figures. This promotes mixing and entrainment in the body of the gravity current through the formation of KH billows, as evident in figures 4 and 7. In addition to the primary KH-driven mixing, figures 7(b), 7(e) and 7(h) also reveal a localised characteristic of secondary convective instability. Notably, strong non-monotonic variations in the vertical density profiles and corresponding spikes in the wall-normal velocity component suggest transient regions where heavier fluid overlies lighter ambient fluid. This density inversion, in combination with persistent vertical shear, establishes favourable conditions for convective overturning. Such instability enhances vertical momentum exchange and facilitates additional entrainment of the ambient fluid. As a result, the mixing layer in the body of the current not only thickens due to KH billows, but also grows through convective overturning, further promoting interfacial mixing and entrainment. The enhanced mixing and entrainment within the body of the gravity current have also been documented by previous researchers (Härtel et al. Reference Huppert2000b ; Shin et al. Reference Simpson2004; Fragoso et al. Reference Fragoso, Patterson and Wettlaufer2013).
The tail, while exhibiting a lower turbulence intensity, remains significant in dissipating the current’s energy, controlling the mixing rate as the gravity current dissipates and integrates into the ambient medium (Hallworth et al. Reference Hallworth, Phillips, Huppert and Sparks1993). The persistence of small-scale eddies in this region influences the final stages of entrainment, but to a lesser extent than the head and body of gravity current (Hallworth et al. Reference Hallworth, Phillips, Huppert and Sparks1993). Although the turbulent effects are diminished in the tail region, as shown by the reduced streamwise and wall-normal velocities in figures 7(c) and 7(g), the viscous dissipation actively takes place. The combined effects of these regions dictate the evolution of the current, making their study crucial for understanding large-scale environmental and geophysical flows. It should be noted, however, that the simulations presented here are two-dimensional in nature and do not fully capture the 3-D dynamics of KH billow roll-up and breakdown. While the simulation domain resolves the wall-normal and streamwise evolution of the flow, the spanwise structure and secondary instabilities associated with billow instability and turbulent mixing are suppressed. Therefore, interpretations regarding the onset and evolution of KH-induced mixing should be viewed as qualitative, and future work involving fully 3-D simulations is necessary to examine these processes in greater depth.
4. Parametric study of sequential lock-fluid release
4.1. Time of overtake
The sequential release of two lock fluids of varying densities resulted in a combined gravity current with complex interactions. The denser lock fluid, introduced later, propagates faster, eventually overtaking the initially released, lighter lock fluid. This overtaking event changes the composition of the current head, initially dominated by the presence of lighter fluid, to become predominantly the denser fluid after overtaking occurs. This shift has significant effects on the propagation dynamics, entrainment processes and mixing of the combined gravity current.
A parametric study was conducted to examine the non-dimensional times required for the denser lock fluid to overtake,
$t_O^*$
, in terms of critical non-dimensional parameters, specifically, the density ratios
$\gamma _1$
and
$\gamma _2$
, and the non-dimensional time of release,
$t^*_R$
. The density ratios,
$\gamma _1$
and
$\gamma _2$
, have a profound impact on the velocity of the gravity current head and the degree of entrainment, both of which are critical in determining
$t_O^*$
. In particular, the influence of the lighter lock fluid in influencing the overtaking dynamics is substantial. Therefore, the characteristic time scales,
$t^*_O$
and
$t^*_R$
, were non-dimensionalised using
$\gamma _1$
in this parametric study to ensure a consistent frame of reference:
$t_O^*=t_O\sqrt {g(\gamma _1^{-1}-1)/H}$
and
$t_R^*=t_R\sqrt {g(\gamma _1^{-1}-1)/H}$
. This approach not only facilitates a systematic and unbiased comparison of cases with varying
$\gamma _2$
while maintaining a constant
$\gamma _1$
, but also eliminates potential ambiguities in interpreting the results. The effect of these density variations on the overtaking time is systematically analysed and presented in figures 8(a) and 8(b).
Figure 8(a) shows the relationship between the
$t^*_O$
and
$\gamma _2$
values across a range of
$t^*_R$
values for a fixed
$\gamma _1$
(
$\gamma _1$
= 0.995), while figure 8(b) provides analogous results for a lower
$\gamma _1$
value (
$\gamma _1$
= 0.985). The ensemble-averaged front velocity for the simulation duration of 30 s, which is rendered non-dimensional with the buoyancy velocity,
$\overline {U}^*=\overline {U}/u_b$
, is plotted as a colour plot in figures 8(a) and 8(b). As
$t^*_R$
increases, the
$t^*_O$
curves consistently shift upward, indicating a corresponding increase in overtaking time, a trend that aligns with expectations. For fixed
$t^*_R$
, the
$t^*_O$
values rise at an increasing rate with higher values of
$\gamma _2$
, as shown in figures 8(a) and 8(b). The increase in
$\gamma _2$
corresponds to a reduction in
$\rho _2$
, thereby diminishing the available potential energy, which is the main driving force of the current. This reduction in driving buoyancy force directly influences the propagation dynamics of the denser lock fluid, leading to a deceleration in its front velocity, as indicated by the colour plots in figures 8(a) and 8(b). As a result, the time required for heavier lock fluid to overtake the preceding lighter lock fluid increases. Moreover, the prolonged overtaking duration results from the reduced buoyancy forcing of the heavier fluid, which propagates in the wake of the lighter current. As
$\gamma _2$
increases, the density contrast decreases, diminishing the heavier current’s ability to accelerate and overtake. This mechanism gives rise to the observed nonlinear dependence between
$t_O^*$
and
$\gamma _2$
, as shown in figures 8(a) and 8(b). The nonlinearity becomes more pronounced at higher
$t_R^*$
due to increased initial separation and the sustained weaker momentum of the heavier fluid. The overtaking time,
$t_O^*$
, for
$\gamma _1=0.985$
is comparatively higher than that for
$\gamma _1=0.995$
, as shown in figures 8(a) and 8(b), since the weaker buoyancy force in the lighter current in the latter case promotes slower propagation and a shorter overtaking time.
Figure 8. Plot of non-dimensionalised overtaking time by heavier lock fluid,
$t_O^*$
, for different density ratios of heavier lock fluid
$(\gamma _2)$
and non-dimensionalised time of delayed release,
$t_R^*$
, for (a)
$\gamma _1=0.995$
and (b)
$\gamma _1=0.985$
. The colour bar shows the non-dimensionalised ensemble-averaged front velocities
$\overline {U}^*=\overline {U}/u_b$
. The dotted line in black
$(Re_l\approx 10\,000)$
separates the linear and nonlinear regimes of
$t_O^*$
.
The density and velocity profiles during overtaking events provide critical insights into the interplay between stratification, shear and entrainment in multi-lock gravity currents. Figure 9 presents the density, streamwise velocity and wall-normal velocity profiles for four distinct overtaking scenarios in the two-lock-fluid gravity current system. In all cases, the profiles are extracted at a fixed location, 0.1 m upstream of the gravity current head, ensuring a consistent comparison of the dynamics and entrainment characteristics across different combinations of parameters of
$\gamma _1$
,
$\gamma _2$
,
$t_O^*$
and
$t_R^*$
, as detailed in table 3. These cases provide critical insights into the fundamental differences in the behaviour of the current head during the overtaking process. In figure 9, cases (a) and (b) correspond to overtaking events where the lighter lock fluid (
$\gamma _1=0.995$
) propagates into the ambient fluid, with
$\gamma _2=0.985$
and
$0.950$
, and release delay times of
$t_R^*=0.993$
and
$4.965$
, respectively, leading to overtaking times of
$t_O^*=4.310$
and
$7.924$
, respectively, as shown in table 3. Similarly, cases (c) and (d) represent overtaking events for a relatively heavier lock fluid (
$\gamma _1=0.985$
), with
$\gamma _2=0.940$
and
$0.975$
, corresponding to
$t_R^*=1.728$
and
$8.643$
, respectively, and overtaking times of
$t_O^*=4.923$
and
$28.523$
, respectively. Cases (b) and (c) in figure 9 exhibit a larger density contrast and velocity gradient between the gravity current and ambient, resulting in a faster front velocity. This is particularly evident in the streamwise velocity profiles in figure 9(b), where a stronger velocity gradient between the gravity current and the ambient enhances the front velocity of the gravity current, leading to a lower
$t_O^*$
and inducing a linear dependence of
$t_O^*$
on
$\gamma _2$
. For cases (a) and (d) with larger
$\gamma _2$
and longer
$t_O^*$
, the reduced density contrast leads to a weaker buoyancy force for the heavier current, which follows in the wake of the lighter current. This delayed progression results in a larger
$t_O^*$
and a nonlinear dependence of
$t_O^*$
on
$\gamma _2$
, as shown in figures 8(a) and 8(b). Stronger buoyancy forcing and shorter
$t_O^*$
lead to linear behaviour, whereas reduced forcing and prolonged overtaking produce nonlinear trends.
Table 3. Comparison of various parameter values of the four cases of overtaking by the heavier lock fluid. Density, streamwise velocity and wall-normal velocities for these four cases are presented in figure 9. Cases b and c fall in the linear regime while cases a and d fall in the nonlinear regime of
$t_O^*$
, as shown in figure 8.
Figure 9. Plots of (a) density, (b) streamwise velocity (c) and wall-normal velocity during four cases of overtaking, as illustrated in table 3. All the plots were taken at 0.1 m upstream of the gravity current front at the time of overtaking,
$t_O^*$
. Cases (b) and (c) lie in the linear regime while cases (a) and (d) fall in the nonlinear regime of
$t_O^*$
, as shown in figure 8.
4.2. Ensemble-averaged front velocity
In the present parametric investigation comprising 100 simulation cases of two-lock-fluid gravity currents, the ensemble-averaged front velocity,
$\overline {U}$
, was systematically evaluated over the total simulation duration of 30 s. To facilitate a dimensionless analysis, the ensemble-averaged front velocity was normalised using the characteristic buoyancy velocity,
$u_b$
. The resultant dimensionless front velocity,
$\overline {U}^*$
, is presented in a colour-coded regime plot in figure 8, illustrating its variation in all simulated cases. The results indicate a monotonic decrease in
$\overline {U}^*$
with increasing
$\gamma _2$
, which is an expected outcome, as higher values of
$\gamma _2$
correspond to a reduction in the density of the heavier lock fluid
$(\rho _2)$
, consequently weakening the driving force of the gravity current. Furthermore, an increase in the dimensionless release time,
$t_R^*$
, leads to a reduction in
$\overline {U}^*$
, attributable to the delayed participation of the denser lock fluid in the current propagation.
Additionally, the relationship between dimensionless overtaking time,
$t_O^*$
, and
$\gamma _2$
exhibits a linear trend for lower values of
$\gamma _2$
and higher values of
$\overline {U}^*$
. However, as
$\gamma _2$
increases,
$\overline {U}^*$
decreases and the overtaking time gradually deviates from linearity, transitioning into a nonlinear regime. This deviation is primarily driven by the weaker buoyancy forcing of the heavier fluid, which travels behind the lighter current and takes longer to overtake it, as discussed in § 4.1. A threshold behaviour is observed for both
$\gamma _1=0.995$
and
$0.985$
, where overtaking events follow a linear regime when the local Reynolds number,
$Re_l$
, remains approximately above 10 000, as demarcated by the dotted black line in figures 8(a) and 8(b). Beyond this threshold, the heavier current, having weaker buoyancy force and travelling in the wake of the lighter current, experiences a substantial increase in overtaking time, giving rise to the observed nonlinear regime.
5. Conclusions
This study investigates the propagation dynamics of sequentially released gravity currents in a two-lock configuration and compares them with the classical single-lock release gravity current. Through high-fidelity LES simulations, we analysed the dynamic evolution of the gravity current head, revealing the fundamental differences arising from multi-layered lock-fluid interactions. The overtaking dynamics in the two-lock-fluid cases was found to significantly influence the front propagation, with the heavier lock fluid accelerating the current upon overtaking the lighter fluid. A parametric study was carried out to investigate the influence on the non-dimensionalised overtaking time,
$t_O^*$
, through variations in the lock-fluid ratios,
$\gamma _1$
and
$\gamma _2$
, and the non-dimensionalised delay of heavier lock-fluid release,
$t_R^*$
. Based on the findings of this study, the following key conclusions can be drawn:
-
(i) In a two-lock-fluid gravity current released with non-dimensionalised delay time,
$t_R^*$
, initially, lighter lock fluid is responsible for entrainment and mixing provided that
$t_R^*$
is small enough to avoid separation of the two lock fluids. The initial concentration of the lighter lock fluid at the head of the gravity current results in the initial slower speed of the current. -
(ii) As heavier lock fluid accumulates at the head region of the current, acceleration of the current is observed.
-
(iii) The head of the gravity current is dominated by streamwise velocity vectors influencing the front speed of the current. The body of the gravity current consists of components of both the streamwise and wall-normal velocities, mixed density layer and KH instabilities at the interface promoting entrainment and mixing. The tail of the current consisted of small-scale eddies that promote long-term viscous dissipation.
-
(iv) The time of overtake,
$t_O^*$
, increases with higher
$t_R^*$
and
$\gamma _2$
, and lower
$\gamma _1$
. As
$\gamma _2$
increases, the weaker buoyancy force of the heavier fluid trailing behind the lighter current becomes increasingly important, resulting in a nonlinear relationship between
$t_O^*$
and
$\gamma _2$
. A transition from linear to nonlinear behaviour was observed near a threshold local Reynolds number of
$Re_l\approx 10\,000$
. -
(v) The non-dimensionalised ensemble-averaged front velocity,
$\overline {U}^*$
, of the two-lock-fluid gravity currents decreases for increasing
$\gamma _2$
and
$t_R^*$
, which is attributed to a decreased gravitational force and delayed participation of heavier lock fluid in gravity current propagation.
The insights gained from this study have important implications for geophysical and environmental fluid mechanics, particularly in improving the understanding and predictive modelling of transport and mixing processes associated with multiple gravity-driven flows in natural and engineered systems. In particular, the delay time (
$t_R^*$
) and relative density ratios (
$\gamma _1$
,
$\gamma _2$
) investigated here are representative of sequential inflow events or buoyancy-driven discharges of varying densities, such as sediment-laden river plumes, wastewater outfalls or turbidity currents, entering stratified basins, reservoirs or coastal environments.
It is important to note that the simulations in this study are two-dimensional in nature, which prevents the representation of key 3-D processes including the onset of secondary instabilities, spanwise undulations or the rapid breakdown of KH billows that arise from 3-D vortex stretching. These processes typically accelerate the transition to turbulence and enhance small-scale mixing, and their suppression in 2-D simulations may lead to more coherent interfacial structures and delayed billow breakdown. While the results in this study remain physically insightful and broadly applicable for understanding the flow propagation dynamics of multi-lock gravity currents, future studies employing fully 3-D configurations will be necessary to more comprehensively generalise the conclusions, particularly with respect to the dynamics of interfacial instability and mixing.
Figure 10. Iso-surfaces of density from the 3-D Smagorinsky LES at the nose of the gravity current front at non-dimensional times: (a)
${t^*=7.717}$
and (b)
${t^*=13.886}$
, corresponding to densities of (a) 1008
${\rm kg\,m^{-3}}$
and (b) 1012
${\rm kg\,m^{-3}}$
.
Acknowledgements
The authors thank the three anonymous referees for their useful comments and recommendations that have helped us to substantially improve the manuscript.
Funding
We gratefully acknowledge funding from the Office of Naval Research (N00014-22-1-2043).
Declaration of interests
The authors report no conflicts of interest.
Data availability statement
The data supporting the findings of this study are available upon reasonable and scientifically justified request.
Appendix A. Comparison of 2-D and 3-D LES results
To assess the applicability of the 2-D LES framework used in this study, a supplementary 3-D Smagorinsky LES was performed for case II with a uniform grid spacing of 2 mm. Figure 10 shows iso-surfaces from the 3-D simulation in the vicinity of the gravity current head at two representative non-dimensional times,
$t^*=7.717$
and
$t^*=13.886$
, corresponding to density values of
$1008$
and
$1012$
$\rm kg\,m^{-3}$
, respectively.
Figure 11. Qualitative comparison of density fields between (a, d) the 2-D simulation, (b, e) the mid-span slice of the 3-D simulation and (c, f) the spanwise-averaged 3-D simulation at non-dimensional times: (a, b, c)
$\boldsymbol{t}^*=7.717$
and (d, e, f)
$\boldsymbol{t}^*=13.886$
. The vertical white dotted lines indicate the locations of the vertical profiles presented in figure 12.
Figure 12. Profiles of density, streamwise velocity and wall-normal velocity at six streamwise locations (denoted as profiles (a)–( f) in figure 11). The red lines represent the 2-D simulations, the blue lines correspond to the mid-span slice of the 3-D simulations and the green lines denote the spanwise-averaged 3-D simulations.
Figure 11 presents a qualitative comparison of the density field between (a, d) the 2-D simulation, (b, e) a mid-span slice of the 3-D simulation and (c, f) the spanwise-averaged 3-D field at the same times. The comparisons indicate similar front positions, interfacial structure and large-scale flow features across the three representations, suggesting that the dominant propagation characteristics and primary interfacial dynamics are comparable.
A more quantitative comparison is provided in figure 12, which shows vertical profiles of density, streamwise velocity and wall-normal velocity extracted at six streamwise locations (labelled (a)–( f) in figure 11). The 2-D results are shown in red, the mid-span 3-D slice in blue and the spanwise-averaged 3-D profiles in green. Good agreement is observed at locations (a)–(d), while modest differences appear at locations (e) and ( f), particularly near the front region. Overall, these results indicate that, for the present parameter range, the 2-D Smagorinsky LES reproduces the principal large-scale features of the corresponding 3-D simulation.
Open access
