1. Introduction
Gravity currents, also known as buoyancy or density currents, are flows along a horizontal boundary driven by horizontal density differences, and occur in natural and anthropic environments. To name a few, examples of gravity currents include sea and land breezes, turbidity currents, powder snow avalanches and pyroclastic flows. The reader is referred to Simpson (Reference Simpson1997) and Ungarish (Reference Ungarish2009) for a comprehensive introduction to this topic and review of examples.
Gravity currents have been studied extensively using the lock-exchange set-up in laboratory experiments and numerical simulations (see e.g. Shin, Dalziel & Linden Reference Shin, Dalziel and Linden2004; Cantero, Balachandar & Garcia Reference Cantero, Balachandar and Garcia2007a; Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007b; Cenedese & Adduce Reference Cenedese and Adduce2008; La Rocca et al. Reference La Rocca, Adduce, Sciortino and Pinzon2008, Reference La Rocca, Adduce, Lombardi, Sciortino and Hinkermann2012a,Reference La Rocca, Adduce, Sciortino, Bateman and Bonifortib; Adduce, Sciortino & Proietti Reference Adduce, Sciortino and Proietti2012; Dai Reference Dai2013, Reference Dai2014; Lombardi et al. Reference Lombardi, Adduce, Sciortino and La Rocca2015; Ottolenghi et al. Reference Ottolenghi, Adduce, Inghilesi, Armenio and Roman2016a,Reference Ottolenghi, Adduce, Inghilesi, Roman and Armeniob, Reference Ottolenghi, Prestininzi, Montessori, Adduce and La Rocca2018; Inghilesi et al. Reference Inghilesi, Adduce, Lombardi, Roman and Armenio2018; Lombardi, Adduce & La Rocca Reference Lombardi, Adduce and La Rocca2018; Pelmard, Norris & Friedrich Reference Pelmard, Norris and Friedrich2020). In the classic lock-exchange experiments, the heavy fluid is typically separated from the ambient fluid in the tank by a removable barrier. When the barrier is withdrawn, the fluids of different densities are set into motion. In the classic lock-exchange experiments, the gravity currents propagate on a horizontal boundary without interaction with counterflowing gravity currents.
Collisions between gravity currents are common in nature. Examples of such convergences include sea breeze (Burpee Reference Burpee1979; Noonan & Smith Reference Noonan and Smith1986; Lapworth Reference Lapworth2005), land breeze (Wapler & Lane Reference Wapler and Lane2012; Gille & Llewellyn Smith Reference Gille and Llewellyn Smith2014), thunderstorm outflows (Droegemeier & Wilhelmson Reference Droegemeier and Wilhelmson1985; Intrieri, Bedard & Hardesty Reference Intrieri, Bedard and Hardesty1990), microburst (Orf, Anderson & Straka Reference Orf, Anderson and Straka1996) and atmosphere bores (Clarke, Smith & Reid Reference Clarke, Smith and Reid1981). There is also evidence, from sedimentary deposits, of the collision of turbidity currents on the ocean bed (Simpson Reference Simpson1997). However, few previous studies have considered the collision between two gravity currents. Simpson (Reference Simpson1997) showed the main effect to be the emergence of two undular bores travelling in opposite directions after a collision. Shin (Reference Shin2001) considered the collision between two gravity currents of equal densities but different heights, and developed a theory that predicts the propagation speed of the bores after the collision.
More recently, van der Wiel et al. (Reference van der Wiel, Gille, Llewellyn Smith, Linden and Cenedese2017) have performed laboratory experiments on the collision of two counterflowing gravity currents using the lock-exchange set-up. The effects of differences in gravity current density and height are studied. For symmetric collisions, i.e. same current densities and heights, the interface separating the two currents is vertical and stationary. For asymmetric collisions, i.e. different current densities and heights, the interface is tilted and changes shape in time. Using two-dimensional simulations based on the incompressible Euler equations with the Boussinesq approximation, Cafaro & Rooney (Reference Cafaro and Rooney2018) replicated numerically situations similar to the laboratory experiments of van der Wiel et al. (Reference van der Wiel, Gille, Llewellyn Smith, Linden and Cenedese2017). The interface slope for asymmetric collisions has been shown to be dependent on the buoyancy ratio of the two gravity currents, whereas the maximum height reached by the fluid after the collision has no strong dependence on the buoyancy ratio of the two gravity currents.
To understand turbulent mixing in the collision between two identical counterflowing gravity currents, Zhong, Hussain & Fernando (Reference Zhong, Hussain and Fernando2018) performed laboratory experiments on colliding gravity currents in a Plexiglas tank using a time-resolved particle image velocimetry and planar laser-induced fluorescence system to record instantaneous velocity and density fields synchronously. According to Zhong et al. (Reference Zhong, Hussain and Fernando2018), four flow stages can be identified: independent propagation of gravity currents, their approach while influencing each other, the collision stage, and post-collision slumping of collided fluid. The lifetime of the collision stage between two identical counterflowing gravity currents can be further divided into three phases based on the spatially averaged vertical front velocity in the collision region. During Phase I, the collision produced a rapidly rising density front with intense turbulent kinetic energy production. During Phase II, the density front was flat with negligible vertical velocity. During Phase III, the collided fluid slumped in the collision region. Furthermore, the eddy diffusivity, space–time averaged in the collision region and over the lifetime of collision stage, was recommended as a conditional parametrization for representing collision events in mesoscale meteorological models.
Our computational study on the collision of gravity currents is motivated by and is a continuation of the experimental work by Zhong et al. (Reference Zhong, Hussain and Fernando2018). Our focus is on the collision stage, and this study complements the work by Zhong et al. (Reference Zhong, Hussain and Fernando2018) in several perspectives. In this study, we conduct three-dimensional high-resolution simulations for the collision of two counterflowing gravity currents of equal densities and heights, and we further implement passive tracers in the simulations to highlight the heavy fluids originally contained in the left and right locks. With the help of high-resolution simulations and implementation of passive tracers, the three-dimensional flow structures and the mixing between the heavy fluids of the two colliding gravity currents in the collision region are presented clearly. Energy budgets for the collision of two gravity currents, along with the spatial distribution and temporal evolution of the mean flow and turbulent flow characteristics in the collision region, are now made possible thanks to the detailed flow information provided by the three-dimensional high-resolution simulations. It is our understanding that such complete information on the energy budgets in the collision of gravity currents can be difficult to attain in laboratory experiments. The paper is organized as follows. In § 2, we describe the formulation of the problem and the implementation of numerical methods. The qualitative and quantitative results are presented in § 3. Finally, conclusions are drawn in § 4.
2. Problem formulation
Figure 1 gives a sketch of the initial configuration for the collision of two gravity currents of equal strengths. The two gravity currents are produced from two identical full-depth locks on the left and right ends of the channel. The height of the locks is 
$\tilde {H}$, which is the same as the environment, and the length of each lock is 
$\tilde {L}_0$. The heavy fluid density in the two locks is 
$\tilde {\rho }_1$, and in between the locks is the homogeneous environment of density 
$\tilde {\rho }_0$. The gravitational acceleration is in the 
$-x_3$ direction. Here, we adopt the Boussinesq approximation, in that the density difference is sufficiently small, i.e. 
$(\tilde {\rho }_1-\tilde {\rho }_0) \ll \tilde {\rho }_0$, so that the influence of density variations is retained in the buoyancy term but neglected in the inertia and diffusion terms.
Figure 1. Sketch of the initial condition for the collision of two gravity currents produced from two identical full-depth locks. The heavy fluid inside the left and right locks has density 
$\tilde {\rho }_1$, height 
$\tilde {H}$ and length 
$\tilde {L}_0$. The ambient fluid of the same height has density 
$\tilde {\rho }_0$ and length 
$\tilde {L}_{x_1} - 2 \tilde {L}_0$. The coordinate system follows the right-hand rule, where 
$x_1$, 
$x_2$ and 
$x_3$ represent the streamwise, spanwise and wall-normal directions, respectively. The origin of the streamwise coordinate is placed at the centre of the flow domain. Here, gravity 
$\tilde {g}$ acts in the 
$-x_3$ direction, and the positive spanwise direction points into the page.
The governing equations take the form, using tensor notation,
$$\begin{gather} \frac{\partial u_{j}}{\partial x_j} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial {u}_{i}}{\partial t} + u_j\,\frac{\partial {u}_{i} }{\partial x_j} = \rho e^{g}_{i} - \frac{\partial p}{\partial x_i} + \frac{1}{Re}\,\frac{\partial^2 u_i }{\partial x_j\,\partial x_j}, \end{gather}$$
$$\begin{gather}\frac{\partial {\rho}}{\partial t} + u_j\,\frac{\partial {\rho} }{\partial x_j} = \frac{1}{Re\,Sc}\,\frac{\partial^2 {\rho}}{ \partial x_j\,\partial x_j}. \end{gather}$$
Here, 
${{u}_i}$ denotes the velocity, 
$\rho$ the density, 
${{e}}^{g}_{i}$ the unit vector in the direction of gravity, and 
$p$ the total pressure (including the part created by the density stratification in the flow domain), respectively. The variables without a tilde are dimensionless quantities. The set (2.1)–(2.3) is made dimensionless by the lock height 
$\tilde {H}$ as the length scale and the buoyancy velocity
\begin{equation} \tilde{u}_b=\sqrt{\tilde{g}'_0 \tilde{H}}, \quad \mbox{with}\ \tilde{g}'_0 = \tilde{g}\,\frac{\tilde{\rho}_{1} - \tilde{\rho}_{0}}{\tilde{\rho}_{0}}, \end{equation}as the velocity scale. The dimensionless density is given by
\begin{equation} \rho = \frac{\tilde{\rho} - \tilde{\rho}_0}{\tilde{\rho}_1 - \tilde{\rho}_0}. \end{equation} Other relevant dimensionless parameters are the Reynolds number 
$Re$ and the Schmidt number 
$Sc$, defined by
\begin{equation} Re = \frac{\tilde{u}_b \tilde{H}}{\tilde{\nu}} \quad \mbox{and} \quad Sc = \frac{\tilde{\nu}}{\tilde{\kappa}}. \end{equation}
It is assumed that the heavy fluid and ambient fluid have identical kinematic viscosity 
$\tilde {\nu }$ and diffusion coefficient of density field 
$\tilde {\kappa }$. It has been reported by researchers (e.g. Härtel, Meiburg & Necker Reference Härtel, Meiburg and Necker2000; Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2005; Bonometti & Balachandar Reference Bonometti and Balachandar2008) that the influence of the Schmidt number is weak as long as 
$Sc \approx O(1)$ or larger, and setting the Schmidt number to unity is common practice in the numerical simulations for gravity currents (Cantero et al. Reference Cantero, Balachandar and Garcia2007a,Reference Cantero, Lee, Balachandar and Garciab; Zgheib, Ooi & Balachandar Reference Zgheib, Ooi and Balachandar2016; Dai & Huang Reference Dai and Huang2022). Therefore we use 
$Sc=1$ in all simulations in the study.
The three-dimensional flow domain is chosen as 
$L_{x_1} \times L_{x_2} \times L_{x_3} = 17.5 \times 1.5 \times 1$, where the lock length at each end of the channel is 
$L_0 = 3$. The length in the problem is non-dimensionalized by the lock height 
$\tilde {H}$. This flow domain is chosen to follow the experimental configuration designed by Zhong et al. (Reference Zhong, Hussain and Fernando2018) for the collision of gravity currents (cases C2800, C3500 and C4300) in their study.
The governing equations are solved using the time-splitting method (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988) with the low-storage third-order Runge–Kutta scheme (Williamson Reference Williamson1980) for time advancement. In the time-splitting method, a provisional velocity, which does not satisfy the continuity condition, is calculated in the first step. The pressure field, which satisfies Poisson's equation, is solved in the second step. The pressure field is used to correct the provisional velocity such that the final velocity and the pressure satisfy the complete governing equations. In the computation of the provisional velocity, the convection and buoyancy terms are treated explicitly, and the diffusion terms are treated implicitly with a Crank–Nicolson scheme. For the convection term, divergence and convective forms are used alternately to reduce the aliasing error (Durran Reference Durran1999). In the three-dimensional flow domain, Fourier expansion with periodic boundary conditions is employed in the streamwise and spanwise directions, i.e. 
$x_1$ and 
$x_2$, while Chebyshev expansion with Gauss–Lobatto quadrature points is employed in the wall-normal direction, i.e. 
$x_3$. For the velocity field, we employ the free-slip condition at the top boundary and no-slip condition at the bottom boundary. For the density field, we employ the no-flux condition at both the top and bottom boundaries.
The initial velocity field was set with a quiescent condition in all simulations. The initial density field was prescribed as unity in the heavy fluid regions at both ends of the channel, and zero in the ambient fluid region between the two locks, with an error-function-type transition, of which the thickness depends on the Reynolds number and Schmidt number. The initial density field was seeded with minute, uniform, random disturbances following Härtel, Michaud & Stein (Reference Härtel, Michaud and Stein1997) and Cantero et al. (Reference Cantero, Balachandar, Garcia and Ferry2006). The de-aliased pseudo-spectral code has been employed in previous high-resolution simulations for lock-exchange flows (Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007b; Dai Reference Dai2015; Dai & Wu Reference Dai and Wu2016; Dai & Huang Reference Dai and Huang2022).
In order to distinguish between the heavy fluids originally contained within the left and right locks in the channel, we adopted the technique in Dai, Huang & Hsieh (Reference Dai, Huang and Hsieh2021) and implemented two passive tracers in the simulations, namely 
$C_L$ and 
$C_R$, for the heavy fluids originally contained within the left and right lock regions, respectively. Introducing the passive tracers in the simulations is identical to adding dyes in the locks for visualization in the experiments. The passive tracers were assumed to have identical diffusion coefficient 
$\tilde {\kappa }$ as the density field, and likewise follow the mass transport equation (2.3) as the density field does. The initial concentration for the passive tracer 
$C_L$ (
$C_R$) was set as unity in the left (right) lock region and zero elsewhere, with an error-function-type transition.
In this study, we are interested in the collision of gravity currents, and we consider four Reynolds numbers, i.e. 
$Re=3450$, 
$6450$, 
$8950$, 14 450, which correspond to the four front Reynolds numbers, i.e. 
$Re_f=751$, 
$1398$, 
$1939$, 
$2954$. The front Reynolds number is defined as 
$Re_f = \tilde {u}_f \tilde {d} / \tilde {\nu }$ and is related to the Reynolds number via 
$Re_f = u_f d\,Re$, where 
$u_f$ is the dimensionless front velocity, and 
$d$ is the dimensionless height of the head. To provide sufficient resolution for the simulations of collision of gravity currents, the grid spacing must be of the order of 
$O{(Re\,Sc)}^{-1/2}$ (Härtel et al. Reference Härtel, Michaud and Stein1997, Reference Härtel, Meiburg and Necker2000; Birman, Martin & Meiburg Reference Birman, Martin and Meiburg2005), and the grids employed in the three-dimensional simulations for the preceding four Reynolds numbers follow Dai & Huang (Reference Dai and Huang2022) and are listed in table 1. The time step was chosen such that the Courant number remained less than 
$0.5$.
Table 1. List of simulations for the colliding gravity currents. Four Reynolds numbers, 
$Re=3450$, 
$6450$, 
$8950$ and 14 450, are considered. Quantitative measures based on the simulation results include the dimensionless front velocity 
$u_f$, dimensionless height of the head 
$d$, and the front Reynolds number 
$Re_f$. The domain size is kept the same for all simulations. The grid resolution increases with increasing 
$Re$ to provide sufficient resolution for the colliding gravity currents.
3. Results
3.1. Observations of the collision of gravity currents
We begin by presenting the collision of gravity currents of equal strengths in close proximity to the centre of the flow domain. Here, the origin of the streamwise coordinate 
$x_1=0$ is placed at the centre of the flow domain. Figure 1 shows the side view of the initial configuration for the collision of gravity currents. As the heavy fluids in the two identical locks on both ends of the channel are released, two gravity currents of equal strengths are formed and propagate towards the centre of the flow domain.
In order to more easily visualize the three-dimensional simulations and to take the ensemble (Reynolds) average for the problem in § 3.4, here we use the average in the spanwise direction for the variables in the three-dimensional simulations. The average of a variable 
$f(x_1,x_2,x_3)$ in the spanwise direction is defined as
\begin{equation} \langle \,f \rangle = \frac{1}{L_{x_2}} \int_0^{L_{x_2}} f\left(x_1, x_2, x_3\right) {\rm d}\kern0.06em x_2, \end{equation}
where 
$L_{x_2}=1.5$, and 
$f$ can be a scalar or vector variable, such as the dimensionless density 
$\rho (x_1, x_2, x_3)$, velocity 
$u_i(x_1, x_2, x_3)$, where 
$i=1,2,3$, or pressure 
$p(x_1, x_2, x_3)$ in the three-dimensional simulations.
Figure 2 shows the dimensionless density averaged in the spanwise direction for the collision of gravity currents at 
$Re=6450$ at 
$t=11.60$, 
$13.29$, 
$14.85$, 
$15.27$, 
$16.40$ and 
$19.23$, or equivalently, 
$T=-1.2$, 
$-0.5$, 
$0.16$, 
$0.34$, 
$0.81$ and 
$2.0$. Here, 
$T$ denotes the shifted time, whose definition and time origin will be introduced in § 3.2. Prior to the collision, the two gravity currents of equal strengths propagate towards the centre of the flow domain independently, with the height of the head 
$d \approx 0.516$ and approximately constant front speed 
$u_f \approx 0.420$ for 
$Re=6450$. Based on the height of the gravity current head and the front speed prior to the collision, the front Reynolds number is 
$Re_f = 1398$. The height of the gravity current head, front speed and front Reynolds number prior to the collision in all the cases considered in this study are listed in table 1 for reference.
Figure 2. Spanwise-averaged density fields for the collision of gravity currents of equal strengths at 
$Re=6450$. The two gravity currents are produced from two identical full-depth locks on the left and right ends of the channel. Time instances are chosen at (a) 
$t=11.60$, (b) 
$13.29$, (c) 
$14.85$, (d) 
$15.27$, (e) 
$16.40$ and (f) 
$19.23$, or equivalently, (a) 
$T=-1.2$, (b) 
$-0.5$, (c) 
$0.16$, (d) 
$0.34$, (e) 
$0.81$ and (f) 
$2.0$. The shifted time 
$T$ is introduced in § 3.2.
The two gravity currents propagate towards the centre of the flow domain independently without influencing each other, as shown in figure 2(a) at 
$t=11.60$ (
$T=-1.2$), until the two gravity current fronts are sufficiently close to each other within approximately one lock height, as shown in figure 2(b) at 
$t=13.29$ (
$T=-0.5$). As this stage of flow, the two approaching gravity currents begin to slow down and are about to collide. Afterwards, the two colliding gravity currents make contact, cause updraughts and slumping of collided heavy fluid away from the collision region, as shown in figure 2 at 
$t=14.85$, 
$15.27$, 
$16.40$ and 
$19.23$, or equivalently, 
$T=0.16$, 
$0.34$, 
$0.81$ and 
$2.0$.
In order to visualize the interaction and mixing between the heavy fluids of the two colliding gravity currents, figure 3 shows the spanwise-averaged density fields along with the dynamic pressure 
$p_d$ and concentrations of the passive tracers 
$C_L$ and 
$C_R$ in the collision region. The dynamics pressure 
$\,p_{d}$ is defined such that its gradient is in balance with the material acceleration and the viscous term (Liu & Katz Reference Liu and Katz2006), i.e.
\begin{equation} \frac{\partial p_d}{\partial x_i} =-\frac{{\rm D} {u}_{i}}{{\rm D} t} + \frac{1}{Re}\,\nabla^2 u_i, \end{equation}
Figure 3. Spanwise-averaged (a i–a v) density fields, (b i–b v) dynamic pressure and (c i–c v) concentrations of the passive tracers for the collision of gravity currents of equal strengths at 
$Re=6450$ at five time instances 
$13.29$, 
$14.85$, 
$15.27$, 
$16.40$ and 
$19.23$, or equivalently, 
$T=-0.5$, 
$0.16$, 
$0.34$, 
$0.81$ and 
$2.0$. The shifted time 
$T$ is introduced in § 3.2. Spanwise-averaged velocity (
$U_1$, 
$U_3$) is represented by vectors for reference.
where 
${\rm D}/{\rm D}t$ denotes the material derivative. The dynamic pressure gradient is related to the total pressure gradient via 
$\partial p_d/\partial x_i = \partial p/\partial x_i - \rho e^g_i$, and the spatial average of the dynamic pressure across the field is set to zero. The concentrations of the passive tracers 
$C_L$ and 
$C_R$ are visualized by yellow and cyan colours, respectively, and the mixing of the two passive tracers results in a green colour. As demonstrated by figure 3 at 
$t=13.29$ (
$T=-0.5$), as the two colliding gravity currents are sufficiently close to each other within approximately one lock height, the dynamic pressure between the two approaching fronts rises, and the ambient fluid between the two approaching fronts begins to move vertically upwards. After the two colliding gravity currents make initial contact, the resulting updraughts approach the top boundary, create rising dynamic pressure close to the top boundary, and deflect horizontally, as shown in figure 3 at 
$t=14.85$, 
$15.27$ and 
$16.40$ (
$T=0.16$, 
$0.34$ and 
$0.81$). It is worth noting that as the collision of gravity currents continues, mixing occurs not only between the heavy fluid and ambient fluid – such as the breakdown of Kelvin–Helmholtz billows in the proximity of the density interface in the laboratory experiments by Zhong et al. (Reference Zhong, Hussain and Fernando2018) – but also between the heavy fluids of the two colliding gravity currents, as shown by the green colour in figure 3 at 
$t=16.40$ and 
$19.23$ (
$T=0.81$ and 
$2.0$). The mechanism responsible for the mixing between the heavy fluids of the two colliding gravity currents will be discussed in § 3.3.
3.2. Collided heavy fluid in the 
$H$–$H$ domain
As discussed in Zhong et al. (Reference Zhong, Hussain and Fernando2018), the height of the heavy fluid accumulated in the 
$H$–
$H$ domain is an important parameter in defining the three phases of collision. The 
$H$–
$H$ domain in this study is the region (
$-H/2 \leqslant x_1 \leqslant H/2$; 
$0 \leqslant x_2 \leqslant L_{x_2}$; 
$0 \leqslant x_3 \leqslant H$) placed at the centre of the flow domain. Following Shin et al. (Reference Shin, Dalziel and Linden2004) and Cantero et al. (Reference Cantero, Lee, Balachandar and Garcia2007b), we define an unambiguous metric for the dimensionless height of the heavy fluid in the 
$H$–
$H$ domain as
\begin{equation} h (t) = \frac{1}{L_{x_2} H^2} \int_{-H/2}^{H/2} \int_0^{L_{x_2}} \int_0^{H} \rho ( x_1, x_2, x_3, t )\, {\rm d}\kern0.06em x_3\,{\rm d}\kern0.06em x_2 \,{\rm d}\kern0.06em x_1. \end{equation} Figure 4 shows the time evolution of the dimensionless height of the heavy fluid in the 
$H$–
$H$ domain, where the normalized abscissa is 
$T=(t-t_c)u_f/H$. The time instance 
$t_c$ is chosen at the time instance when 
$h=0.5$, following Zhong et al. (Reference Zhong, Hussain and Fernando2018). Since the gravity current thickness from a full-depth lock-exchange set-up is approximately 
$0.5H$, the time origin 
$T=0$ is the time instant when the two colliding gravity currents have made full contact, i.e. 
$h=0.5$. With such a choice of time origin, the two colliding gravity currents make initial contact at 
$T \approx -0.2$. Prior to the two colliding gravity currents entering the 
$H$–
$H$ domain, the dimensionless height of the heavy fluid in the 
$H$–
$H$ domain remains identically zero, i.e. 
$h=0$. Upon collision, the dimensionless height of the heavy fluid in the 
$H$–
$H$ domain increases rapidly, while the two colliding gravity currents continue to feed into the 
$H$–
$H$ domain. The dimensionless height of heavy fluid in the 
$H$–
$H$ domain reaches a maximum value 
$h \approx 0.9$ at 
$T \approx 1.23$. Afterwards, the dimensionless height of the heavy fluid in the 
$H$–
$H$ domain decreases over time. The maximum height of the heavy fluid in the 
$H$–
$H$ domain and the time evolution of the dimensionless height of the heavy fluid in the 
$H$–
$H$ domain in our study are consistent with the results of Zhong et al. (Reference Zhong, Hussain and Fernando2018), where the experimental set-ups were similar to ours in the simulations.
Figure 4. Time evolution of the dimensionless height of the heavy fluid in the 
$H$–
$H$ domain for the collision of gravity currents at four Reynolds numbers. Symbols: 
$\square$, 
$Re=3450$; 
$\circ$, 
$Re=6450$; 
$\triangle$, 
$Re=8950$; 
$\diamond$, 
$Re=14\,450$. From Zhong et al. (Reference Zhong, Hussain and Fernando2018): solid black line, C2800; dashed red line, C3500; solid blue line C4300. The dimensionless time 
$t$ is shown on the top horizontal axis, and the shifted time 
$T$ is shown on the bottom horizontal axis. The shifted time is defined as 
$T=(t-t_c)u_f/H$, where 
$t_c$ is the instance at which 
$h=0.5$. The vertical dashed lines correspond to 
$T=-0.2$, 
$0.5$, 
$1.2$ and 
$2.8$, where Phase I is 
$-0.2 \leqslant T \leqslant 0.5$, Phase II is 
$0.5 \leqslant T \leqslant 1.2$, and Phase III is 
$1.2 \leqslant T \leqslant 2.8$.
Following Zhong et al. (Reference Zhong, Hussain and Fernando2018), the ‘averaged’ vertical velocity of the density front in the 
$H$–
$H$ domain can be defined as
\begin{equation} w_{f,a} = \frac{{\rm d} h(t)}{{\rm d} t}. \end{equation}
Figure 5 shows the time evolution of the averaged vertical velocity of the density front in the 
$H$–
$H$ domain, and we can identify clearly the three phases in the collision stage as reported in the literature. Here, Phase I is in the time period 
$-0.2 \leqslant T \leqslant 0.5$, during which the two colliding gravity currents make initial contact at 
$T \approx -0.2$, and the averaged vertical velocity of the density front in the 
$H$–
$H$ domain reaches its maximum at 
$T \approx -0.1$. The averaged vertical velocity of the density front in the 
$H$–
$H$ domain then begins to decrease, while the rate of decay of 
$w_{f,a}$ reaches its maximum at 
$T \approx 0.5$. Phase II is in the time period 
$0.5 \leqslant T \leqslant 1.2$, during which the averaged vertical velocity of the density front in the 
$H$–
$H$ domain continues to decrease at a lower rate of decay of 
$w_{f,a}$, but remains positive until 
$T \approx 1.2$, at which 
$w_{f,a}$ clearly changes sign from positive to negative. Phase III is in the time period 
$1.2 \leqslant T \leqslant 2.8$, during which the averaged vertical velocity of the density front in the 
$H$–
$H$ domain remains negative, which indicates that the collided fluid begins to slump away from the collision region.
Figure 5. Averaged vertical velocity of the density front in the 
$H$–
$H$ domain for the collision of gravity currents at four Reynolds numbers. Symbols: 
$\square$, 
$Re=3450$; 
$\circ$, 
$Re=6450$; 
$\triangle$, 
$Re=8950$; 
$\diamond$, 
$Re=14\,450$. The dimensionless time 
$t$ is shown on the top horizontal axis, and the shifted time 
$T$ is shown on the bottom horizontal axis. The shifted time is defined as 
$T=(t-t_c)u_f/H$, where 
$t_c$ is the instance at which 
$h=0.5$. The vertical dashed lines correspond to 
$T=-0.2$, 
$0.5$, 
$1.2$ and 
$2.8$, where Phase I is 
$-0.2 \leqslant T \leqslant 0.5$, Phase II is 
$0.5 \leqslant T \leqslant 1.2$, and Phase III is 
$1.2 \leqslant T \leqslant 2.8$.
3.3. Flow structures in the $H$–$H$ domain
It is now understood that the height of the heavy fluid in the 
$H$–
$H$ domain increases rapidly as the two counterflowing gravity currents collide, and the collision of gravity currents may create geometric distortions of the fronts. Since we are interested in the energy balances for the collision of gravity currents, it is instructive to examine the flow structures in the 
$H$–
$H$ domain as the gravity currents collide before we move on to the mean flow and turbulent flow characteristics in the 
$H$–
$H$ domain.
Figure 6 shows the swirling strength in the 
$H$–
$H$ domain for the collision of gravity currents at 
$Re=6450$ and, for illustrative purposes, the time instance is chosen at 
$T \approx 0.10$ (
$t=14.71$) in Phase I of collision. As shown in figure 6, there exists an array of vertical vortices along the interface separating the two colliding gravity currents. As time proceeds, these vertical vortices are stretched longitudinally and subsequently break up into smaller structures.
Figure 6. Three-dimensional view of the flow structures in the 
$H$–
$H$ domain for the collision of gravity currents at 
$Re=6450$. For illustrative purposes, the isosurface of the swirling strength is 
$\lambda _{ci}=3.68$, and the time instance is chosen at 
$T \approx 0.10$ (
$t=14.71$) in Phase I of collision. The red and blue colours of the isosurfaces of the swirling strength represent the orientations of the vortex in the positive and negative wall-normal directions, respectively.
In order to identify the mechanism responsible for the occurrence of an array of vertical vortices along the interface between the two colliding gravity currents, we study the wall-normal component of the vorticity equation in the 
$H$–
$H$ domain. We take the curl of the momentum equation (2.2), and the wall-normal component of the vorticity equation is
\begin{equation} \frac{{\rm D} \omega_3 }{{\rm D} t} = \underbrace{\frac{\partial u_3}{\partial x_2}\,\frac{\partial u_1}{\partial x_3}}_{S_1} \underbrace{{}-\frac{\partial u_3}{\partial x_1}\,\frac{\partial u_2}{\partial x_3}}_{S_2} \underbrace{{}+\omega_3\,\frac{\partial u_3}{\partial x_3}}_{S_3} \underbrace{{}+\frac{1}{Re_{}}\,\nabla^2 \omega_3}_{S_4}, \end{equation}
where 
$S_1$, 
$S_2$, 
$S_3$ and 
$S_4$ represent the tilting of 
$x_1$ vorticity, twisting of 
$x_2$ vorticity, stretching of 
$x_3$ vorticity, and diffusion of 
$x_3$ vorticity, respectively. Since the gravity is in the negative wall-normal direction, there is no baroclinic production of vorticity in the wall-normal direction. We also remark that the component 
$-(\partial u_2/\partial x_3)(\partial u_1/\partial x_3)$ in the tilting of 
$x_1$ vorticity, and the component 
$(\partial u_1/\partial x_3)(\partial u_2/\partial x_3)$ in the twisting of 
$x_2$ vorticity, have no net effects in 
${{\rm D} \omega _3 }/{{\rm D} t}$ because they cancel exactly.
The contribution of the stretching of 
$x_3$ vorticity, i.e. the 
$S_3$ term, to 
${\rm D} \omega _3 /{\rm D}t$ at 
$T \approx 0.10$ (
$t=14.71$) is shown in figure 7(a). It is worth noting that when the array of vertical vortices develops as the two gravity currents collide, the major contribution to 
${\rm D} \omega _3 /{\rm D}t$ comes from the stretching of 
$x_3$ vorticity, i.e. the 
$S_3$ term. Other contributions to 
${\rm D} \omega _3 /{\rm D}t$ from the tilting of 
$x_1$ vorticity, twisting of 
$x_2$ vorticity and diffusion of 
$x_3$ vorticity, i.e. the 
$S_1$, 
$S_2$ and 
$S_4$ terms, are less than or even of opposite sign to the time rate of change of 
$x_3$ vorticity, and are not shown for brevity. Our results indicate that the occurrence of an array of vertical vortices developing along the interface between the two colliding gravity currents is due to the stretching of pre-existing 
$x_3$ vorticity inside the lobes. Our simulations show that, consistent with previously published reports (Cantero et al. Reference Cantero, Lee, Balachandar and Garcia2007b; Espath et al. Reference Espath, Pinto, Laizet and Silvestrini2015; Dai & Huang Reference Dai and Huang2022), the 
$x_3$ vorticity pre-exists inside the lobes of the two colliding gravity currents even before the two gravity currents enter the 
$H$–
$H$ domain. Inside a lobe moving in the streamwise direction, positive 
$x_3$ vorticity exists in the left portion of the lobe, while negative 
$x_3$ vorticity exists in the right portion of the lobe. As the two gravity currents collide, geometric distortions of the fronts create 
$\partial u_3/\partial x_3$, and consequently, the stretching of pre-existing 
$x_3$ vorticity. The flow structures along with the concentrations of the passive tracers can also be visualized in a horizontal slice at 
$x_3=0.45$ and 
$T \approx 0.10$ (
$t=14.71$) in the 
$H$–
$H$ domain, as shown in figure 7(b). Based on the vorticity and concentrations of the passive tracers, the mixing between the heavy fluids of the two gravity currents is observed in a way that the heavy fluid of one gravity current between a pair of counter-rotating vertical vortices is induced across the interface at 
$x_1 \approx 0$ into the heavy fluid of the other gravity current.
Figure 7. Flow structures visualized with vertical and horizontal slices in the 
$H$–
$H$ domain for the collision of gravity currents at 
$Re=6450$. The time instance is chosen at 
$T \approx 0.10$ (
$t=14.71$) in Phase I of collision. (a) The time rate of change of 
$x_3$ vorticity, 
$D \omega _3/D t$, and the stretching of 
$x_3$ vorticity, i.e. the 
$S_3$ term in (3.5), in the vertical slice at 
$x_1=0.07$. The time rate of change of 
$x_3$ vorticity is visualized by the colour contours, and the positive and negative contributions of the 
$S_3$ term are represented by the thin solid and dashed lines, respectively. (b) Concentrations of the passive tracers 
$C_L$ and 
$C_R$, velocity 
$(u_1,u_2)$ and 
$x_3$ vorticity in the horizontal slice at 
$x_3=0.45$. The concentrations 
$C_L$ and 
$C_R$ are visualized by cyan and yellow, velocity 
$(u_1,u_2)$ is visualized by vectors, and the positive and negative 
$x_3$ vorticities are represented by the thin solid and dashed lines, respectively.
3.4. Energy budgets in the $H$–$H$ domain
From the point of view of energy budgets, the propagation of gravity currents into a homogeneous environment is a process in which the potential energy is converted into kinetic energy and subsequently into dissipation (Dai Reference Dai2015; Dai & Huang Reference Dai and Huang2016; Dai & Wu Reference Dai and Wu2016). The collision of gravity currents, in contrast to the propagation of gravity currents, is a process in which the kinetic energy is converted into potential energy, while dissipation may also occur in the collision process.
It is our understanding that the energy budgets can be difficult to attain in laboratory experiments, and there is no complete information on the energy budgets in the collision of gravity currents. In the following, we will provide a computational analysis of the energy budget based on our three-dimensional high-resolution simulations for the colliding gravity currents. In addition, the flow field is statistically homogeneous in the spanwise direction, and we may derive the ensemble (Reynolds) average in the spanwise direction as the mean flow, and the deviations from the mean flow as the turbulence fluctuations. Therefore, we will also provide an analysis of the mean flow energy and turbulent kinetic energy equations.
For a concise presentation of our analysis, the derivation of the total energy, mean flow energy and turbulent kinetic energy equations is placed in Appendix A, to which readers are referred for more details.
3.4.1. Total energy equation
The evolution equation of the total energy can be written as
\begin{equation} \frac{\partial E_k}{\partial t} + \frac{\partial E_p}{\partial t} = {\mathcal{T}_{total}} +{\mathcal{M}} +{\mathcal{D}_{total}}, \end{equation}
where 
$E_k=\frac {1}{2} u_i u_i$ is the total kinetic energy, 
$E_p= \rho x_3$ is the total potential energy, 
$\mathcal {T}_{total}$ is the rate of transport of total energy, 
$\mathcal {M}$ is the rate of conversion from internal energy due to irreversible diffusion of density, and 
$\mathcal {D}_{total}$ is the rate of dissipation of total energy.
Integration of (3.6) over the entire flow domain 
$\varOmega$ leads to the conservation of energy equation, i.e.
\begin{equation} E^\varOmega_k + E^\varOmega_p + E^\varOmega_d - E^\varOmega_i = E^\varOmega_{p,t=0}, \end{equation}
where superscript 
$\varOmega$ is used to denote integration over the entire flow domain, and 
$E^\varOmega _{p,t=0}$ is the initial potential energy in the system. The total kinetic energy 
$E^\varOmega _k$ and total potential energy 
$E^\varOmega _p$ are
\begin{equation} E^\varOmega_k(t) = \int_{\varOmega} E_k \,{\rm d} V \quad \mbox{and} \quad E^\varOmega_p(t) = \int_{\varOmega} E_p \,{\rm d} V; \end{equation}
the total dissipated energy 
$E^\varOmega _d$ is
\begin{equation} E^\varOmega_d (t) =-\int_{0}^{t} \mathcal{D}^{\varOmega}_{total} (\tau) \,{\rm d} \tau, \quad \mathcal{D}^{\varOmega}_{total} = \int_{\varOmega} \mathcal{D}_{total} \,{\rm d} V, \end{equation}
and the total conversion from internal energy 
$E^\varOmega _i$ is
\begin{equation} E^\varOmega_i (t) = \int_{0}^{t} \mathcal{M}^{\varOmega} (\tau) \,{\rm d} \tau, \quad \mathcal{M}^{\varOmega} = \int_{\varOmega} \mathcal{M} \,{\rm d} V. \end{equation}
Here, the integration of 
$\mathcal {T}_{total}$ over the entire flow domain 
$\varOmega$ vanishes. It is worth noting that 
$E^\varOmega _i$ represents the conversion from internal energy to potential energy due to irreversible diffusion of density (Winters et al. Reference Winters, Lombard, Riley and D'Asaro1995; Dai & Wu Reference Dai and Wu2016), and, as we will show next, the conversion from internal energy can be considered as unimportant for the collision of gravity currents.
For illustrative purposes, figure 8 shows the energy budget in the entire flow domain for the colliding gravity currents at 
$Re=6450$. The energy budgets for the colliding gravity currents at other Reynolds numbers are qualitatively similar and are not shown here for brevity. It is observed that the total energy in the entire flow domain is conserved to a high degree of accuracy. When the two gravity currents propagate towards the centre of the flow domain independently, the potential energy is converted into kinetic energy and into dissipation. As the two gravity currents collide, the potential energy in the entire flow domain begins to increase at the expense of decreasing kinetic energy at 
$T \approx -0.17$ (
$t \approx 14.07$). The potential energy in the entire flow domain continues to increase until 
$T \approx 3.43$ (
$t \approx 22.63$), when the kinetic energy in the entire flow domain begins to increase at the expense of decreasing potential energy. The dissipated energy during the course of propagation, including the collision of gravity currents, comprises up to 
$33\,\%$ of the total energy, while the conversion from internal energy is approximately 
$2.2\,\%$ by the end of the simulation. Therefore, the contribution from the internal energy, when compared with the kinetic energy, potential energy and dissipated energy in the entire flow domain, can be considered as unimportant for the collision of gravity currents.
Figure 8. Energy budget in the entire flow domain for the collision of gravity currents at 
$Re=6450$. All energies are normalized by the initial potential energy in the entire flow domain, 
$E^\varOmega _{p,t=0}$. Symbols: 
$\triangle$, 
$E^\varOmega _k$; 
$\triangledown$, 
$E^\varOmega _p$; 
$\circ$, 
$E^\varOmega _d$; 
$\diamond$, 
$E^\varOmega _i$; 
$\square$, 
$E^\varOmega _k+E^\varOmega _p+E^\varOmega _d-E^\varOmega _i$. The dimensionless time 
$t$ is shown on the top horizontal axis, and the shifted time 
$T$ is shown on the bottom horizontal axis. The vertical dashed lines correspond to 
$T=-0.2$, 
$0.5$, 
$1.2$ and 
$2.8$, where Phase I is 
$-0.2 \leqslant T \leqslant 0.5$, Phase II is 
$0.5 \leqslant T \leqslant 1.2$, and Phase III is 
$1.2 \leqslant T \leqslant 2.8$.
Since we are interested in the collision of gravity currents in the collision region, integration of (3.6) over the 
$H$–
$H$ domain leads to the evolution equation of the total energy in the 
$H$–
$H$ domain:
\begin{equation} \frac{\partial E^{H}_k}{\partial t} + \frac{\partial E^{H}_p}{\partial t} = {\mathcal{T}^{H}_{total}} +{\mathcal{M}^{H}} +{\mathcal{D}^{H}_{total}}, \end{equation}
where superscript 
$H$ is used to denote integration over the 
$H$–
$H$ domain.
Figure 9 shows the temporal evolution of the contributions to the evolution equation of the total energy in the 
$H$–
$H$ domain (3.11) for the colliding gravity currents at 
$Re=6450$. As shown in figure 9(a), that 
$\partial E^H_k/\partial t$ is positive until 
$T \approx 0.01$ and becomes negative afterwards indicates that the total kinetic energy in the 
$H$–
$H$ domain increases in early stages of Phase I of collision, then decreases in the rest of Phase I and Phase II. In Phase III, 
$\partial E^H_k/\partial t \approx 0$ indicates that the total kinetic energy in the 
$H$–
$H$ domain remains nearly unchanged. The potential energy in the 
$H$–
$H$ domain increases in Phase I and Phase II, then decreases in Phase III because 
$\partial E^H_p/\partial t$ is positive until 
$T \approx 1.2$ and becomes negative afterwards. In Phase I, Phase II and early stages of Phase III, energy is transported into the 
$H$–
$H$ domain since 
$\mathcal {T}^H_{total}$ is positive until 
$T \approx 1.26$. In the rest of Phase III, energy is transported out of the 
$H$–
$H$ domain since 
$\mathcal {T}^H_{total}$ becomes negative when 
$T \gtrsim 1.26$. As shown in figure 9(b), the rate of dissipation of total energy in the 
$H$–
$H$ domain, 
$\mathcal {D}^H_{total}$, reaches its maximum magnitude at 
$T \approx 0.22$ in Phase I of collision, and diminishes in Phase II and Phase III. The conversion from internal energy due to irreversible diffusion of density in the 
$H$–
$H$ domain, 
$\mathcal {M}^H$, as discussed above, can be considered as negligible.
Figure 9. Contributions to the evolution equation of the total energy in the 
$H$–
$H$ domain (3.11) for the collision of gravity currents at 
$Re=6450$. All contributions are normalized by 
$u^3_f H^2$. Symbols in (a): 
$\triangle$, 
$\partial E^H_k/\partial t$; 
$\triangledown$, 
$\partial E^H_p/\partial t$; 
$\square$, 
$\mathcal {T}^H_{total}$. Symbols in (b): 
$\diamond$, 
$\mathcal {M}^H$; 
$\circ$, 
$\mathcal {D}^H_{total}$. The dimensionless time 
$t$ is shown on the top horizontal axis, and the shifted time 
$T$ is shown on the bottom horizontal axis. The vertical dashed lines correspond to 
$T=-0.2$, 
$0.5$, 
$1.2$ and 
$2.8$, where Phase I is 
$-0.2 \leqslant T \leqslant 0.5$, Phase II is 
$0.5 \leqslant T \leqslant 1.2$, and Phase III is 
$1.2 \leqslant T \leqslant 2.8$.
The total energy in the 
$H$–
$H$ domain can be further decomposed into the mean flow energy and turbulent kinetic energy, as will be shown in the next subsubsection.
3.4.2. Mean flow and turbulent kinetic energy equations
We are interested in the evolution of mean flow and turbulent flow characteristics, and for this purpose, the ensemble (Reynolds) average for the problem is defined following (3.1), where the brackets 
$\langle \,\rangle$ represent the average performed in the spanwise direction. With the ensemble average, 
$U_i= \langle {u}_i\rangle$, 
$P= \langle p\rangle$ and 
$\bar {\rho }=\langle \rho \rangle$ are the mean flow velocity, pressure and density, respectively; and 
$u_i'=u_i-U_i$, 
$p'=p-P$ and 
$\rho '=\rho -\bar {\rho }$ are the velocity, pressure and density fluctuations, respectively.
The evolution equation of the mean flow energy can be written as
\begin{equation} \frac{\partial E}{\partial t} + \frac{\partial E_p}{\partial t} = {\mathcal{T}_{mean}} + {\mathcal{M}} - {\mathcal{B}_{turb}} - {\mathcal{P}_{turb}} + {\mathcal{D}_{mean}}, \end{equation}
where 
$E=\frac {1}{2} U_i U_i$ is the mean flow kinetic energy, 
$E_p=\bar {\rho } x_3$ is the mean flow potential energy, 
$\mathcal {T}_{mean}$ is the rate of transport of mean flow energy, 
$\mathcal {M}$ is the rate of conversion from internal energy due to irreversible diffusion of density, 
$-\mathcal {B}_{turb}$ is the rate of transfer of energy away from the mean flow due to turbulent buoyancy flux, 
$-\mathcal {P}_{turb}$ is the rate of transfer of energy away from the mean flow due to production of turbulent kinetic energy, and 
$\mathcal {D}_{mean}$ is the rate of dissipation of mean flow energy.
Figure 10 shows the spatial distribution of the mean flow characteristics in the 
$H$–
$H$ domain at five instances (
$T=-0.5$, 
$0.16$, 
$0.34$, 
$0.81$ and 
$2.0$) during the three phases of collision at 
$Re=6450$. The mean flow kinetic energy demonstrates that strong shears are generated in the low-velocity strips during Phase I and Phase II, as shown in figures 10(a i–a iii), and the rate of dissipation of mean flow energy is most evident in the strong shears of the mean flow, as shown in figures 10(g i–g iii). The rate of transfer of energy away from the mean flow due to turbulent buoyancy flux is apparent, as shown by the red and blue strips close to the top boundary in figures 10(e ii–e iv). The turbulent buoyancy flux close to the top boundary is created by a number of horizontal streamwise vortices aligned close to the top boundary. Detailed inspections show that the presence of the top boundary deflects the updraughts of collision into the horizontal direction. The vorticity in these horizontal streamwise vortices aligned close to the top boundary intensifies further due to stretching in the horizontal direction in Phase I of collision, and these horizontal streamwise vortices break up into smaller structures in Phase II of collision. The rate of transfer of energy away from the mean flow due to production of turbulent kinetic energy is apparent in the rising density front and throughout the updraught fluid column where the vertical vortices are located, as shown in figures 10(f ii–f iii). The rate of transfer of energy away from the mean flow due to production of turbulent kinetic energy then weakens in Phase II and Phase III of collision.
Figure 10. Mean flow and turbulent flow characteristics for the collision of gravity currents at 
$Re=6450$. Panels show: (a i–a v) 
$E$, (b i–b v) 
$k$, (c i–c v) 
$\mathcal {T}_{mean}$, (d i–d v) 
$\mathcal {T}_{turb}$, (e i–e v) 
$-\mathcal {B}_{turb}$, (f i–f v) 
$-\mathcal {P}_{turb}$, (g i–g v) 
$\mathcal {D}_{mean}$, (h i–h v) 
$\mathcal {D}_{turb}$. Time instances are chosen at 
$T=-0.5$, 
$0.16$, 
$0.34$, 
$0.81$ and 
$2.0$. Mean flow velocity (
$U_1$, 
$U_3$) is represented by vectors for reference.
Integration of (3.12) over the 
$H$–
$H$ domain leads to the evolution equation of the mean flow energy in the 
$H$–
$H$ domain:
\begin{equation} \frac{\partial E^H}{\partial t} + \frac{\partial E^H_p}{\partial t} = {\mathcal{T}^H_{mean}} + {\mathcal{M}^H} - {\mathcal{B}^H_{turb}} - {\mathcal{P}^H_{turb}} + {\mathcal{D}^H_{mean}}, \end{equation}
where superscript 
$H$ is used to denote integration over the 
$H$–
$H$ domain.
Figures 11(a) and 11(b) show the temporal evolution of the contributions to the evolution equation of the mean flow energy in the 
$H$–
$H$ domain (3.13) for the collision at 
$Re=6450$. The potential energy and conversion from internal energy of the mean flow, when integrated over the 
$H$–
$H$ domain, are identical to their counterparts in the total energy equation in the 
$H$–
$H$ domain (cf. (3.11)), therefore the temporal evolutions of 
$\partial E^H_p/\partial t$ and 
$\mathcal {M}^H$ are included in figure 11 for reference but not discussed for brevity. As shown in figure 11(a), the mean flow kinetic energy in the 
$H$–
$H$ domain increases with positive 
$\partial E^H/\partial t$ in Phase I of collision until 
$T \approx -0.02$, then decreases with negative 
$\partial E^H/\partial t$ in the rest of Phase I and Phase II. In Phase III, the mean flow kinetic energy in the 
$H$–
$H$ domain remains nearly unchanged. The mean flow energy is transported into the 
$H$–
$H$ domain, with positive 
$\mathcal {T}^H_{mean}$ in Phase I, Phase II and early stages of Phase III until 
$T \approx 1.41$., and turns to be transported out of the 
$H$–
$H$ domain with negative 
$\mathcal {T}^H_{mean}$ in the rest of Phase III. As shown in figure 11(b), the rate of dissipation of mean flow energy in the 
$H$–
$H$ domain, 
$\mathcal {D}^H_{mean}$, reaches its maximum magnitude at 
$T \approx -0.17$, which occurs earlier than the time 
$T \approx 0.22$ when 
$\mathcal {D}^H_{total}$ reaches its maximum magnitude. In addition, at 
$T \approx 0.22$ when 
$\mathcal {D}^H_{total}$ reaches its maximum magnitude, the contribution to the total rate of dissipation of energy from 
$\mathcal {D}^H_{mean}$ is approximately 
$22\,\%$ of the maximum magnitude of 
$\mathcal {D}^H_{total}$, which suggests that the dissipation of energy in the 
$H$–
$H$ domain is not predominated by the strong shears of the mean flow as shown in figures 10(g i–g iii). The rate of transfer of energy away from the mean flow due to production of turbulent kinetic energy in the 
$H$–
$H$ domain, 
$-\mathcal {P}^H_{turb}$, reaches its maximum magnitude at 
$T \approx 0.16$, at which the vertical vortices extending throughout the updraught fluid column begin to become distorted and then break up into smaller structures. The rate of transfer of energy away from the mean flow due to turbulent buoyancy flux in the 
$H$–
$H$ domain, 
$-\mathcal {B}^H_{turb}$, reaches its maximum magnitude at 
$T \approx 0.78$, at which the horizontal streamwise vortices close to the top boundary have broken up into smaller structures. It is interesting to note that the rate of transfer of energy away from the mean flow is due primarily to the production of turbulent kinetic energy in Phase I of collision and is due primarily to the turbulent buoyancy flux in Phase II of collision.
Figure 11. Contributions to the evolution equations of the mean flow energy and the turbulent kinetic energy in the 
$H$–
$H$ domain, i.e. (3.13) and (3.15), for the collision of gravity currents at 
$Re=6450$. All contributions are normalized by 
$u^3_f H^2$. Symbols in (a): 
$\triangle$, 
$\partial E^H/\partial t$; 
$\triangledown$, 
$\partial E^H_p/\partial t$; 
$\square$, 
$\mathcal {T}^H_{mean}$. Symbols in (b): 
$\diamond$, 
$\mathcal {M}^H$; 
$\blacksquare$, 
$-\mathcal {B}^H_{turb}$; 
$\bullet$, 
$-\mathcal {P}^H_{turb}$; 
$\circ$, 
$\mathcal {D}^H_{mean}$. Symbols in (c): 
$\triangle$, 
$\partial k^H/\partial t$; 
$\square$, 
$\mathcal {T}^H_{turb}$. Symbols in (d): 
$\blacksquare$, 
$\mathcal {B}^H_{turb}$; 
$\bullet$, 
$\mathcal {P}^H_{turb}$; 
$\circ$, 
$\mathcal {D}^H_{turb}$. The dimensionless time 
$t$ is shown on the top horizontal axis, and the shifted time 
$T$ is shown on the bottom horizontal axis. The vertical dashed lines correspond to 
$T=-0.2$, 
$0.5$, 
$1.2$ and 
$2.8$, where Phase I is 
$-0.2 \leqslant T \leqslant 0.5$, Phase II is 
$0.5 \leqslant T \leqslant 1.2$, and Phase III is
$1.2 \leqslant T \leqslant 2.8$.
The evolution equation of the turbulent kinetic energy can be written as
\begin{equation} \frac{\partial k}{\partial t} = {\mathcal{T}_{turb}} +{\mathcal{B}_{turb}} +{\mathcal{P}_{turb}} +{\mathcal{D}_{turb}}, \end{equation}
where 
$k=\langle \frac {1}{2} u_i' u_i'\rangle$ is the turbulent kinetic energy, 
$\mathcal {T}_{turb}$ is the rate of transport of turbulent kinetic energy, 
$\mathcal {B}_{turb}$ is the rate of transfer of energy to the turbulent flow due to turbulent buoyancy flux, 
$\mathcal {P}_{turb}$ is the rate of transfer of energy to the turbulent flow due to production of turbulent kinetic energy, and 
$\mathcal {D}_{turb}$ is the rate of dissipation of turbulent kinetic energy. While 
$-\mathcal {B}_{turb}$ and 
$-\mathcal {P}_{turb}$ terms transfer energy away from the mean flow in (3.12), 
$\mathcal {B}_{turb}$ and 
$\mathcal {P}_{turb}$ terms supply energy to the turbulent flow in (3.14). It becomes clear that the turbulent buoyancy flux and the production of turbulent kinetic energy simply transfer the energy away from the mean flow into the turbulent flow without changing the total energy.
Figure 10 also shows the spatial distribution of the turbulent flow characteristics in the 
$H$–
$H$ domain at five instances (
$T=-0.5$, 
$0.16$, 
$0.34$, 
$0.81$ and 
$2.0$) during the three phases of collision at 
$Re=6450$. The turbulent kinetic energy demonstrates that the vertical vortices along the interface separating the two colliding gravity currents create strong turbulence extending throughout the updraught fluid column, as shown in figures 10(b i–b iii), and dissipation in Phase I of collision, as shown in figures 10(h i–h iii).
Integration of (3.14) over the 
$H$–
$H$ domain leads to the evolution equation of the turbulent kinetic energy in the 
$H$–
$H$ domain:
\begin{equation} \frac{\partial k^H}{\partial t} = {\mathcal{T}^H_{turb}} +{\mathcal{B}^H_{turb}} +{\mathcal{P}^H_{turb}} +{\mathcal{D}^H_{turb}}, \end{equation}
where superscript 
$H$ is used to denote integration over the 
$H$–
$H$ domain.
Figures 11(c) and 11(d) show the temporal evolution of the contributions to the evolution equation of the turbulent kinetic energy in the 
$H$–
$H$ domain (3.15) for the collision at 
$Re=6450$. As shown in figure 11(a), the turbulent kinetic energy in the 
$H$–
$H$ domain increases at its maximum rate in Phase I, continues to increase at a lower rate in Phase II until 
$T \approx 1.17$, and decreases in the rest of Phase II and Phase III. The turbulent kinetic energy is transported into the 
$H$–
$H$ domain with positive 
$\mathcal {T}^H_{turb}$ in Phase I and early stages of Phase II until 
$T \approx 0.57$, and turns to be transported out of the 
$H$–
$H$ domain with negative 
$\mathcal {T}^H_{turb}$ in the rest of Phase II and Phase III. As shown in figure 11(d), the rate of dissipation of turbulent kinetic energy in the 
$H$–
$H$ domain, 
$\mathcal {D}^H_{turb}$, reaches its maximum magnitude at 
$T \approx 0.43$, which occurs after both the time 
$T \approx -0.17$ when 
$\mathcal {D}^H_{mean}$ reaches its maximum magnitude, and the time 
$T \approx 0.22$ when 
$\mathcal {D}^H_{total}$ reaches its maximum magnitude. Furthermore, at 
$T \approx 0.22$, when 
$\mathcal {D}^H_{total}$ reaches its maximum magnitude, the contribution to the total rate of dissipation of energy from 
$\mathcal {D}^H_{turb}$ is approximately 
$78\,\%$ of the maximum magnitude of 
$\mathcal {D}^H_{total}$, which suggests that the dissipation of energy in the 
$H$–
$H$ domain is predominated by 
$\mathcal {D}^H_{turb}$, as shown in figures 10(h i–h iii). As discussed above, the energy transferring away from the mean flow is supplied to the turbulent flow. As such, the contribution to the turbulent kinetic energy budget is due primarily to production of turbulent kinetic energy in Phase I of collision and is due primarily to turbulent buoyancy flux in Phase II of collision.
4. Conclusions
In this study, we investigated the collision of two counterflowing gravity currents of equal densities and heights using three-dimensional high-resolution simulations of the incompressible Navier–Stokes equations with the Boussinesq approximation. The goal is to deepen our understanding of the flow structures and energetics, including the spatial distribution and temporal evolution of the mean flow and turbulence characteristics, in the collision region, in more detail.
According to Zhong et al. (Reference Zhong, Hussain and Fernando2018), there are four flow stages in the experiments of colliding gravity currents: independent propagation of gravity currents, their approach while influencing each other, collision stage, and post-collision slumping of collided fluid. The lifetime of the collision stage is approximately 
$3 \tilde {H}/\tilde {u}_f$, where 
$\tilde {H}$ is the depth of heavy and ambient fluids, and 
$\tilde {u}_f$ is the front velocity of the approaching gravity currents. Our focus in this study is on the collision stage and, specifically, the three phases involved in the collision stage. The three phases involved in the collision stage are: Phase I, 
$-0.2 \leqslant (\tilde {t}-\tilde {t}_c) \tilde {u}_f/\tilde {H} \leqslant 0.5$; Phase II, 
$0.5 \leqslant (\tilde {t}-\tilde {t}_c) \tilde {u}_f/\tilde {H} \leqslant 1.2$; Phase III, 
$1.2 \leqslant (\tilde {t}-\tilde {t}_c) \tilde {u}_f/\tilde {H} \leqslant 2.8$, where 
$\tilde {t}$ is the time, and 
$\tilde {t}_c$ is the time instance at which the two colliding gravity currents have fully osculated.
During Phase I of collision, geometric distortions of the gravity current fronts result in stretching of pre-existing wall-normal components of vorticity inside the gravity current fronts, and consequently an array of vertical vortices extending throughout the updraught fluid column develop along the interface separating the two colliding gravity currents. These vertical vortices along the interface, in both positive and negative wall-normal directions, are responsible for the mixing between the heavy fluids of the two colliding gravity currents and for the strong turbulence extending throughout the updraught fluid column. Based on our analysis of the energy budget, the turbulent kinetic energy is supplied primarily by the production of turbulent kinetic energy term in Phase I of collision, and the dissipation of energy in the collision region is predominated by the dissipation of turbulent kinetic energy instead of by the strong shears of the mean flow.
Towards the end of Phase I, the presence of the top boundary deflects the updraughts of collision into the horizontal direction, and a number of horizontal streamwise vortices aligned close to the top boundary are created. These horizontal streamwise vortices close to the top boundary induce turbulent buoyancy flux and, during Phase II of collision, break up into smaller structures. While the production of turbulent kinetic energy weakens in Phase II and in Phase III, the transfer of energy from the mean flow to the turbulent flow due to turbulent buoyancy flux reaches its maximum in Phase II of collision. In contrast to Phase I of collision, in which the production of turbulent kinetic energy term is the primary supply in the turbulent kinetic energy, the transfer of energy from the mean flow to the turbulent flow due to turbulent buoyancy flux becomes the primary supply in the turbulent kinetic energy in Phase II of collision. From the point of view of energetics, the production of turbulent kinetic energy and turbulent buoyancy flux simply transfers energy away from the mean flow into the turbulent flow without changing the total energy. During Phase III of collision, the collided fluid slumps away from the collision region, while the production of turbulent kinetic energy, turbulent buoyancy flux and dissipation of energy attenuate.
The three-dimensional high-resolution simulations reported in this study complement the experimental investigation by Zhong et al. (Reference Zhong, Hussain and Fernando2018). While it was evident in the experiments that mixing occurs at the interface between the heavy fluid and the ambient fluid, such as the breakdown of Kelvin–Helmholtz billows in the proximity of the density interface, the mechanism responsible for the mixing between the heavy fluids of the two gravity currents, and a discernible amount of turbulence extending throughout fluid column in the collision region, remained unresolved in the literature. With the help of high-resolution simulations, the flow structures responsible for the mixing between the heavy fluids of the two gravity currents are now revealed, and the spatial distribution and temporal evolution of the mean flow and turbulent flow characteristics in the collision region are presented. It is our understanding that such complete information on the flow structures and energy budgets in the collision of gravity currents can be difficult to attain in laboratory experiments.
The present study is for the collision of gravity currents of equal strengths from a full-depth lock-exchange configuration. Future extensions may include, for example, the collision of gravity currents of different densities and heights, collision of two cylindrical gravity currents relevant for colliding microburst pairs, collision of inwardly propagating gravity currents in a radially symmetric set-up relevant for the collision of sea breezes over near-circular islands, and the collision of gravity currents under the influence of the Coriolis forces. These other aspects would be worth further study.
Acknowledgements
A.D. is grateful for encouragement from Professors P. Linden and S. Dalziel at the University of Cambridge, S. Balachandar at the University of Florida, and M. Garcia and G. Parker at the University of Illinois at Urbana-Champaign. Computational resources are provided by the Computer and Information Networking Center at National Taiwan University.
Funding
The research was supported financially by Taiwan Ministry of Science and Technology through grants MOST-111-2221-E-002-113-MY3 and MOST-111-2811-E002-029-MY3.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of the energy equations
This appendix contains the derivation of the total energy, mean flow energy and turbulent kinetic energy equations.
A.1. Derivation of the total energy equation
The evolution equation of the total energy is obtained by multiplying the momentum equation (2.2) by 
$u_i$, i.e.
\begin{equation} \frac{{\rm D}}{{\rm D} t} \left( \frac{1}{2}\,u_i u_i \right) =- \rho u_3 -\frac{\partial}{\partial x_j} \left( p u_j \right) + \frac{2}{Re}\, \frac{\partial}{\partial x_j} \left( u_i s_{ij} \right) - \frac{2}{Re}\, s_{ij} s_{ij}, \end{equation}
where 
${\rm D}/{\rm D}t$ denotes the material derivative, 
$s_{ij}$ denotes the strain rate tensor, 
$s_{ij}= \frac {1}{2} (u_{i,j}+u_{j,i})$, and 
$u_3$ denotes the velocity component in the 
$x_3$ direction. Using the relationship
\begin{equation} \rho u_3 = \frac{{\rm D}}{{\rm D} t} ( \rho x_3 ) - \frac{{\rm D} \rho}{{\rm D} t}\,x_3 \end{equation}
and 
${\rm D}/{\rm D}t = \partial / \partial t + u_j\,\partial / \partial x_j$, we may rewrite (A1) as
\begin{align} \frac{\partial}{\partial t} \underbrace{\left( \frac{1}{2}\,u_i u_i \right)}_{E_k}{}+ \frac{\partial}{\partial t} \underbrace{( \rho x_3 )}_{E_p} &= \underbrace{- \frac{\partial}{\partial x_j} \left( p u_j + \rho x_3 u_j + \frac{1}{2}\,u_i u_i u_j - \frac{2}{Re}\,u_i s_{ij} \right) }_{\mathcal{T}_{total}} \nonumber\\ &\quad \underbrace{{}+ \left( \frac{1}{Re \,Sc}\,\frac{\partial^2 \rho}{\partial x_j\,\partial x_j} \right ) x_3 }_{\mathcal{M}} \underbrace{{}- \frac{2}{Re}\,s_{ij} s_{ij} }_{\mathcal{D}_{total}}, \end{align}
where 
$E_k=\frac {1}{2} u_i u_i$ is the total kinetic energy, 
$E_p=\rho x_3$ is the total potential energy, 
$\mathcal {T}_{total}$ is the rate of transport of total energy, 
$\mathcal {M}$ is the rate of conversion from internal energy due to irreversible diffusion of density, and 
$\mathcal {D}_{total}$ is the rate of dissipation of total energy.
A.2. Derivation of the mean flow energy equation
To decompose the flow field into the mean flow and turbulent fluctuations, the mean flow is derived by taking the ensemble (Reynolds) average of the flow field. The ensemble average, denoted by the brackets 
$\langle \,\rangle$, is performed in the spanwise direction following (3.1), and the deviations from the mean flow are taken as the turbulent fluctuations. With the ensemble average, 
$U_i= \langle u_i\rangle$, 
$P= \langle p\rangle$ and 
$\bar {\rho }=\langle \rho \rangle$ are the mean flow velocity, pressure and density, respectively; and 
$u_i'=u_i-U_i$, 
$p'=p-P$ and 
$\rho '=\rho -\bar {\rho }$ are the velocity, pressure and density fluctuations, respectively.
The evolution equation of the mean flow energy can be derived by taking the ensemble average of (2.2) and then multiplying the averaged equation by 
$U_i$, i.e.
\begin{align} \frac{\partial \frac{1}{2} U_i U_i}{\partial t} + U_j\,\frac{\partial \frac{1}{2} U_i U_i}{\partial x_j} &=- \frac{\partial}{\partial x_j} \left( P U_j + U_i\,\langle u_i' u_j'\rangle - \frac{2}{Re}\,U_i {S}_{ij} \right) \nonumber\\ &\quad - \bar{\rho} U_3 + \langle u_i' u_j'\rangle\,\frac{\partial U_i}{\partial x_j} - \frac{2}{Re}\,S_{ij} S_{ij}, \end{align}
where 
$S_{ij}=\frac {1}{2} (U_{i,j}+U_{j,i})$ denotes the mean flow strain rate tensor. Using the relationship
\begin{equation} \bar{\rho} U_3 = \frac{\partial \bar{\rho} x_3 }{\partial t} + \frac{\partial}{\partial x_j} \left( \bar{\rho} U_j x_3 \right) - \left[ \frac{\partial \bar{\rho}}{\partial t} + \frac{\partial}{\partial x_j} \left( \bar{\rho} U_j \right) \right] x_3, \end{equation}and noting that taking the ensemble average of (2.3) gives
\begin{equation} \frac{\partial \bar{\rho}}{\partial t} + \frac{\partial}{\partial x_j} \left( \bar{\rho} U_j \right) = \frac{1}{Re \,Sc}\,\frac{\partial^2 \bar{\rho}}{\partial x_j\,\partial x_j} - \frac{\partial \langle\rho' u_j'\rangle}{\partial x_j}, \end{equation}we may rewrite (A4) as
\begin{align} \frac{\partial }{\partial t} \underbrace{\left( \frac{1}{2}\,U_i U_i \right)}_{E}{}+ \frac{\partial }{\partial t} \underbrace{( \bar{\rho} x_3 )}_{E_p} &= \underbrace{-\frac{\partial}{\partial x_j} \left( P U_j + \bar{\rho} x_3 U_j + \frac{1}{2}\,U_i U_i U_j + U_i\,\langle u_i' u_j'\rangle - \frac{2}{Re}\,U_i {S}_{ij} \right)}_{\mathcal{T}_{mean}} \nonumber\\ &\quad \underbrace{{}+ \left( \frac{1}{Re \,Sc}\,\frac{\partial^2 \bar{\rho}}{\partial x_j\,\partial x_j} \right) x_3}_{\mathcal{M}} \underbrace{{}-\left( \frac{\partial \langle \rho' u_j'\rangle }{\partial x_j} \right) x_3}_{-\mathcal{B}_{turb}} \nonumber\\ &\quad \underbrace{{}+ \langle u_i' u_j'\rangle\,\frac{\partial U_i}{\partial x_j}}_{-\mathcal{P}_{turb}} \underbrace{{}-\frac{2}{Re_{}}\,S_{ij} S_{ij}}_{\mathcal{D}_{mean}}, \end{align}
where 
$E=\frac {1}{2} U_i U_i$ is the mean flow kinetic energy, 
$E_p =\bar {\rho } x_3$ is the mean flow potential energy, 
$\mathcal {T}_{mean}$ is the rate of transport of mean flow energy, 
$\mathcal {M}$ is the rate of conversion from internal energy due to irreversible diffusion of density, 
$-\mathcal {B}_{turb}$ is the rate of transfer of energy away from the mean flow due to turbulent buoyancy flux, 
$-\mathcal {P}_{turb}$ is the rate of transfer of energy away from the mean flow due to production of turbulent kinetic energy, and 
$\mathcal {D}_{mean}$ is the rate of dissipation of mean flow energy.
A.3. Derivation of the turbulent kinetic energy equation
The evolution equation of the turbulent kinetic energy can be derived by subtracting the ensemble average of (2.2) from (2.2) itself, multiplying the resulting equation by 
$u_i'$, and then taking the ensemble average again, i.e.
\begin{align} \frac{\partial \langle\frac{1}{2} u_i' u_i'\rangle}{\partial t} + U_j\,\frac{\partial \langle\frac{1}{2}u_i' u_i'\rangle}{\partial x_j} &=- \frac{\partial}{\partial x_j} \left( \left\langle\frac{1}{2}\,u_j' u_i' u_i'\right\rangle + \langle p' u_j'\rangle - \frac{2}{Re}\, \langle u_i' {s}'_{ij}\rangle \right) \nonumber\\ &\quad - \langle \rho' u_3' \rangle - \langle u_i' u_j'\rangle\,\frac{\partial U_i}{\partial x_j} - \frac{2}{Re}\,\langle s'_{ij} s'_{ij}\rangle. \end{align}Using the relationship
\begin{equation} \langle\rho' u_3'\rangle = \frac{\partial}{\partial x_j} \left( \langle\rho' u_j'\rangle x_3 \right) - \left( \frac{\partial \langle\rho' u_j'\rangle}{\partial x_j} \right) x_3, \end{equation}we may rewrite (A8) as
\begin{align} \frac{\partial }{\partial t} \underbrace{\left( \left\langle\frac{1}{2}\, u_i' u_i'\right\rangle \right)}_{k} &= \underbrace{- \frac{\partial}{\partial x_j} \left( \langle p' u_j'\rangle + \langle\rho' u_j'\rangle x_3 + \left\langle\frac{1}{2}\,u_i' u_i'\right\rangle\,U_j + \left\langle\frac{1}{2}\,u_j' u_i' u_i'\right\rangle - \frac{2}{Re}\,\langle u_i' {s}'_{ij}\rangle \right)}_{\mathcal{T}_{turb}} \nonumber\\ &\quad \underbrace{{}+\left( \frac{\partial \langle\rho' u_j'\rangle}{\partial x_j} \right) x_3}_{\mathcal{B}_{turb}} \underbrace{{}- \langle u_i' u_j'\rangle\, \frac{\partial U_i}{\partial x_j}}_{\mathcal{P}_{turb}} \underbrace{{}-\frac{2}{Re}\,\langle s'_{ij} s'_{ij}\rangle}_{\mathcal{D}_{turb}}, \end{align}
where 
$k=\langle \frac {1}{2} u_i' u_i'\rangle$ is the turbulent kinetic energy, 
$\mathcal {T}_{turb}$ is the rate of transport of turbulent kinetic energy, 
$\mathcal {B}_{turb}$ is the rate of transfer of energy to the turbulent flow due to turbulent buoyancy flux, 
$\mathcal {P}_{turb}$ is the rate of transfer of energy to the turbulent flow due to production of turbulent kinetic energy, and 
$\mathcal {D}_{turb}$ is the rate of dissipation of turbulent kinetic energy.
Appendix B. Influence of the Reynolds number
The collisions between gravity currents on the geophysical scale often have Reynolds numbers of orders of magnitude larger than the Reynolds numbers achievable in laboratory experiments and numerical simulations. When the Reynolds number becomes high enough, the turbulent flow reaches a self-similar stage and the scale effects are insignificant (Barenblatt Reference Barenblatt1996). As investigated by Zhong et al. (Reference Zhong, Hussain and Fernando2018), the collisions between gravity currents appear to enter the self-similar stage when 
$Re_f \gtrsim 3000$, which was proposed by Breidenthal (Reference Breidenthal1981) and discussed by Princevac, Fernando & Whiteman (Reference Princevac, Fernando and Whiteman2005).
Our numerical simulations, like the laboratory experiments, provide a controlled environment for the study of the collisions between gravity currents. The Reynolds numbers considered in this study, 
$751 \leqslant Re_f \leqslant 2954$, are smaller than the Reynolds number for the flow to reach the self-similar stage, but it can be expected that the 
$Re_f=2954$ case (
$Re=14\,450$) approaches the self-similar stage.
For illustrative purposes, in the paper we show the energy budget for the colliding gravity currents at 
$Re=6450$ (
$Re_f=1398$). This appendix contains the energy budget for the colliding gravity currents at 
$Re=14\,450$ (
$Re_f=2954$). Figure 12 shows that the energy budget in the entire flow domain for the colliding gravity currents at 
$Re=14\,450$ (
$Re_f=2954$) exhibits initial propagation of gravity currents, their approach while influencing each other, collision stage and post-collision, and appears to be similar to figure 8. Figure 13 shows the temporal evolution of the contributions to the evolution equation of the mean flow energy and the turbulent kinetic energy in the 
$H$–
$H$ domain for the collision at 
$Re=14\,450$ (
$Re_f=2954$). It is observed that the dissipation of energy in the 
$H$–
$H$ domain is predominated not by the mean flow but by the dissipation of turbulent kinetic energy. In addition, the contribution to the turbulent kinetic energy is due primarily to the production of turbulent kinetic energy in Phase I of collision, and is due primarily to the turbulent buoyancy flux in Phase II of collision. Our findings presented in the paper are evidently supported by the collision of gravity currents at 
$Re=14\,450$ (
$Re_f=2954$). It becomes clear now that the 
$Re_f=2954$ case (
$Re=14\,450$) is qualitatively similar to the 
$Re_f=1398$ case (
$Re=6450$), and the collisions between gravity currents at even higher Reynolds numbers are expected to reach the self-similar stage.
Figure 12. Energy budget in the entire flow domain for the collision of gravity currents at 
$Re=14\,450$ (
$Re_f=2954$). All energies are normalized by the initial potential energy in the entire flow domain 
$E^\varOmega _{p,t=0}$. Symbols: 
$\triangle$, 
$E^\varOmega _k$; 
$\triangledown$, 
$E^\varOmega _p$; 
$\circ$, 
$E^\varOmega _d$; 
$\diamond$, 
$E^\varOmega _i$; 
$\square$, 
$E^\varOmega _k+E^\varOmega _p+E^\varOmega _d-E^\varOmega _i$. The dimensionless time 
$t$ is shown on the top horizontal axis, and the shifted time 
$T$ is shown on the bottom horizontal axis. The vertical dashed lines correspond to 
$T=-0.2$, 
$0.5$, 
$1.2$ and 
$2.8$, where Phase I is 
$-0.2 \leqslant T \leqslant 0.5$, Phase II is 
$0.5 \leqslant T \leqslant 1.2$, and Phase III is 
$1.2 \leqslant T \leqslant 2.8$.
Figure 13. Contributions to the evolution equations of the mean flow energy and the turbulent kinetic energy in the 
$H$–
$H$ domain, i.e. (3.13) and (3.15), for the collision of gravity currents at 
$Re=14450$ (
$Re_f=2954$). All contributions are normalized by 
$u^3_f H^2$. Symbols in (a): 
$\triangle$, 
$\partial E^H/\partial t$; 
$\triangledown$, 
$\partial E^H_p/\partial t$; 
$\square$, 
$\mathcal {T}^H_{mean}$. Symbols in (b): 
$\diamond$, 
$\mathcal {M}^H$; 
$\blacksquare$, 
$-\mathcal {B}^H_{turb}$; 
$\bullet$, 
$-\mathcal {P}^H_{turb}$; 
$\circ$, 
$\mathcal {D}^H_{mean}$. Symbols in (c): 
$\triangle$, 
$\partial k^H/\partial t$; 
$\square$, 
$\mathcal {T}^H_{turb}$. Symbols in (d): 
$\blacksquare$, 
$\mathcal {B}^H_{turb}$; 
$\bullet$, 
$\mathcal {P}^H_{turb}$; 
$\circ$, 
$\mathcal {D}^H_{turb}$. The dimensionless time 
$t$ is shown on the top horizontal axis, and the shifted time 
$T$ is shown on the bottom horizontal axis. The vertical dashed lines correspond to 
$T=-0.2$, 
$0.5$, 
$1.2$ and 
$2.8$, where Phase I is 
$-0.2 \leqslant T \leqslant 0.5$, Phase II is 
$0.5 \leqslant T \leqslant 1.2$, and Phase III is 
$1.2 \leqslant T \leqslant 2.8$.
Appendix C. Eddy diffusivity
In modelling environmental turbulent flows, it is conventional to use the eddy diffusivity 
$k_t$ to describe turbulent buoyancy flux:
\begin{equation} -\frac{\partial \langle\rho' u'_j\rangle}{\partial x_j} = k_t\,\frac{\partial^2 \bar{\rho} }{\partial x_j\,\partial x_j}, \end{equation}
where 
$k_t$ is made dimensionless by the lock height 
$\tilde {H}$ as the length scale, and the buoyancy velocity 
$\tilde {u}_b$ as the velocity scale. Note that 
$k_t$ is a function of space and time. For practical purposes, it would be useful to obtain an averaged eddy diffusivity 
$\bar {k}_t$, spatially averaged over the 
$H$–
$H$ domain, rather than pointwise description of eddy diffusivity and an overall eddy diffusivity 
$K_T$, space–time-averaged over the 
$H$–
$H$ domain and Phases I–III, for a parametrization of collision events in mesoscale models.
As validation for our results and a comparison with the laboratory experiments by Zhong et al. (Reference Zhong, Hussain and Fernando2018), figure 14(a) shows the time evolution of 
$\bar {k}_t$, spatially averaged over the 
$H$–
$H$ domain, and figure 14(b) shows the overall eddy diffusivity 
$K_T$, space–time-averaged over the 
$H$–
$H$ domain and Phases I–III, against 
$Re_f$. It is gratifying to note that our results are quantitatively consistent with the experimental data and supportive of the value 
$K_T/u_f \approx 3.6 \times 10^{-3}$ for 
$Re_f \gtrsim 3000$ proposed by Zhong et al. (Reference Zhong, Hussain and Fernando2018).
Figure 14. (a) Time evolution of the eddy diffusivity spatially averaged over the 
$H$–
$H$ domain. Line colours: black, 
$Re=3450$ (
$Re_f=751$); red, 
$Re=6450$ (
$Re_f=1398$); blue, 
$Re=8950$ (
$Re_f=1939$); green, 
$Re=14\,450$ (
$Re_f=2954$). The dimensionless time 
$t$ is shown on the top horizontal axis, and the shifted time 
$T$ is shown on the bottom horizontal axis. The vertical dashed lines correspond to 
$T=-0.2$, 
$0.5$, 
$1.2$ and 
$2.8$, where Phase I is 
$-0.2 \leqslant T \leqslant 0.5$, Phase II is 
$0.5 \leqslant T \leqslant 1.2$, and Phase III is 
$1.2 \leqslant T \leqslant 2.8$. (b) The overall eddy diffusivity spatially averaged over the 
$H$–
$H$ domain and temporally averaged over Phases I–III against 
$Re_f$. Symbols: 
$\circ$, experimental data by Zhong et al. (Reference Zhong, Hussain and Fernando2018); 
$\bullet$, present study. The error bars are the standard deviations of fluctuations from the temporal averaging.
Open access
