Effect of stratification on the propagation of a cylindrical gravity current

Abstract

Direct numerical simulations (DNSs) of three-dimensional cylindrical release gravity currents in a linearly stratified ambient are presented. The simulations cover a range of stratification strengths $0< S\leq 0.8$ (where $S=(\rho _b^*-\rho _0^*)/(\rho _c^*-\rho _0^*), \rho _b^*, \rho _0^*$ and $\rho _c^*$ are the dimensional density at the bottom of the domain, top of the domain and the dense fluid, respectively) at two different Reynolds numbers. A comparison between the stratified and unstratified cases illustrates the influence of stratification strength on the dynamics of cylindrical gravity currents. Specifically, the front velocity in the slumping phase decreases with increasing stratification strength whereas the duration of the slumping phase increases with increments of $S$. The Froude number calculated in this phase shows a good agreement with models proposed by Ungarish & Huppert (J. Fluid Mech., vol. 458, 2002, pp. 283–301) and Ungarish (J. Fluid Mech., vol. 548, 2006, pp. 49–68), originally developed for planar gravity currents in a stratified ambient. In the inertial phase, the front velocity across cases with different stratification strengths adheres to a power-law scaling with an exponent of $-$1/2. Higher Reynolds numbers led to more frequent lobe splitting and merging, with lobe size diminishing as stratification strength increased. Strong interactions among inner vortex rings occurred during the slumping phase, leading to the early formation of hairpin vortices in weakly stratified cases, while strongly stratified cases exhibited delayed vortex formation and less turbulence.

1. Introduction

Gravity currents or density currents are a horizontal intrusion of different density to an ambient fluid. Gravity currents are observed in many naturally occurring phenomena such as sandstorms (Parsons 2000), powder-snow avalanches (Turnbull & McElwaine 2007) and bushfires (Dold, Zinoviev & Weber 2006). Comprehensive reviews of gravity currents in geophysical flows, laboratory experiments and numerical simulations are given by Simpson (1982) and Meiburg, Radhakrishnan & Nasr-Azadani (2015). We consider a body of heavy fluid initially at rest released into an ambient fluid at $t=0$. The heavy fluid collapses and leads to an intrusion of fluid with a distinct head region. When the gravity current is formed, a well-defined shape is recognised and the flow develops into a highly turbulent head, quasi-steady body and shallower following tail (Cantero et al. 2008; Zordan, Schleiss & Franca 2018).

Experimental (Huppert & Simpson 1980; Alahyari & Longmire 1996; Marino, Thomas & Linden 2005; Patterson, Simpson & Dalziel 2006) and numerical (Blanchette et al. 2005; Cantero, Balachandar & García 2007a; Cantero et al. 2007b, 2008) studies of planar and cylindrical gravity currents have found that propagation of a gravity current undergoes four main stages. As the heavy fluid is released, the gravity current undergoes an acceleration phase (zero to maximum), where the gravitational potential energy of the heavy fluid is converted into kinetic energy. A secondary acceleration is reported by Zhu et al. (2017) and Ooi, Zgheib & Balachandar (2015) who conducted numerical simulations of a circular gravity current on a uniform slope. The secondary acceleration phase is due to the rearrangement and redistribution of the heavy fluid in the current which increases the buoyancy at the downstream end of the current. After the front velocity reaches its maximum, the acceleration phase is followed by the slumping phase, where the front height and speed remain constant. In this phase, the axisymmetric (cylindrical) current head comprises a coherent counter-clockwise (or clockwise, depending on the azimuthal direction being observed) rotating vortex ring. This ring is formed due to the roll-up of the interface shear layer between the heavy and ambient fluids. The effects of stretching and tilting of the vortex ring on the dynamics of axisymmetric (cylindrical) currents have been explored by Patterson et al. (2006) and Cantero et al. (2007a). Next, the gravity current transitions into the self-similar inertial phase in which the buoyancy force is balanced by the inertial force. The gravity current eventually reaches the viscous phase wherein the viscous force dominates the buoyancy force. In the inertial and viscous phases, it has been reported that the front velocity of the gravity current decays as a power law (Fay 1969; Hoult 1972; Cantero et al. 2007b, 2008).

Stratification of the ambient fluid can greatly affect the propagation and dynamics of the gravity current. Additionally, the presence of stratification leads to the generation of internal waves. The specific layer through which the gravity current propagates depends on the relative strength of the current ($\rho _c^*-\rho _0^*$) and the ambient stratification ($\rho _b^*-\rho _0^*$) (Maxworthy et al. 2002; Dai, Huang & Hsieh 2021). In stratified environments, intrusive gravity currents can occur when the density of the heavy fluid matches the density at an intermediate level (see Flynn & Sutherland 2004; Ungarish 2005b; la Forgia et al. 2020; Ottolenghi et al. 2020). Conversely, gravity currents propagate along the bottom boundary when the density of the heavy fluid is equal to or greater than the density at the bottom of the stratified ambient (Maxworthy et al. 2002; Ungarish 2005a; Birman, Meiburg & Ungarish 2007b; White & Helfrich 2008). In this paper, we will only consider the case where the stratification of the ambient fluid is linear for a horizontal domain of uniform depth and the current is denser than the ambient. We will also assume that the release is of full-depth.

Various experimental and numerical studies have been carried out to investigate the propagation of a gravity current into a stratified ambient fluid. Mitsudera & Baines (1992) conducted an experimental study of a downslope gravity current flow in a continuously stratified ambient fluid to model the Bass Strait outflow. An experimental and numerical study was conducted by Maxworthy et al. (2002) to investigate the relationship between the internal Froude number of the gravity current and stratification $S$ in the ambient and

(1.1)\begin{equation} S=\frac{\rho_b^*-\rho_0^*}{\rho_c^*-\rho_0^*}. \end{equation}

The flow regime of the gravity current flow is determined by the Froude number based on the buoyancy frequency which is a dimensionless parameter defined as the ratio of the inertial forces relative to the gravitational forces, i.e. $Fr=u^*_{f,mean}/N^*H^*$, where $u^*_{f,mean}$ is the mean front velocity in the slumping phase, $(N^*)^2 =(g^*/\rho _0^*)(-{\rm d}\rho ^*/{\rm d}z^*)= g^*(\rho _b^* - \rho _0^*)/\rho _0^* H^*$ is the buoyancy frequency, $g^*$ is the gravitational acceleration, $\rho ^*$ is the dimensional fluid density, $z^*$ is the vertical coordinate and $H^*$ is the depth of the domain. Note that in this manuscript, variables with asterisks ($*$) denote dimensional variables. In the subcritical regime ($Fr<1/{\rm \pi}$), the gravity current propagates slower than the maximum speed of the linear internal gravity wave ($N^*H^*/{\rm \pi}$) and minimal (or no) vortices are observed behind the gravity current head, whereas in the supercritical regime ($Fr>1/{\rm \pi}$), the gravity current flow becomes turbulent and strong Kelvin–Helmholtz (K–H) billows can be observed behind the head.

A similar study was conducted by Ungarish & Huppert (2002) using the shallow-water approximation and found excellent agreement with the experimental results presented by Maxworthy et al. (2002). It is worth noting that the propagation speed of the gravity current is reduced by the stratified ambient compared with the homogeneous ambient (where the density of the ambient fluid is constant throughout the domain). A similar phenomenon was previously observed by Lam et al. (2018a,b), Lam, Chan & Ooi (2022b) who conducted direct numerical simulations (DNSs) of a two-dimensional gravity current flow in a stratified ambient to investigate the effects of the initial aspect ratio and the strength of stratification on the dynamics of the gravity current flow. Lam et al. (2018b) reported that in the subcritical regime, the gravity current contains minimal vortices whereas, in the supercritical regime, the gravity current becomes turbulent and strong vortices are observed in the gravity current head.

Planar release gravity currents in stratified (Maxworthy et al. 2002; Ungarish & Huppert 2002; Birman et al. 2007b; Longo et al. 2016; Chiapponi et al. 2018; Lam et al. 2018b; la Forgia et al. 2020; Ottolenghi et al. 2020; Dai et al. 2021; Zahtila et al. 2024), and unstratified ambients (Rottman & Simpson 1983; Shin, Dalziel & Linden 2004; La Rocca et al. 2008; Dai 2015; Pelmard, Norris & Friedrich 2018; Dai & Huang 2020; De Falco, Adduce & Maggi 2021; Maggi, Adduce & Negretti 2022, 2023a; Maggi et al. 2023b), as well as cylindrical release gravity currents in unstratified environment have been widely studied experimentally and numerically (Cantero et al. 2007a,b; Zgheib, Bonometti & Balachandar 2014, 2015a; Dai & Huang 2016). However, experimental and numerical studies on stratified gravity currents with the cylindrical release are limited, and therefore we are interested in cylindrical gravity currents propagating into a linearly stratified ambient.

Birman et al. (2007b) examined the model developed by Ungarish (2006) (where the hydraulic theory from Benjamin (1968) was generalized to a linearly stratified ambient) on the dependency of the front velocity on the stratification strength ($S$). The model was found to work well for weak stratification ($S \leq 0.5$). Ungarish & Huppert (2006) analysed the exchange of energy of gravity currents propagating into an unstratified and linearly stratified ambient with the shallow-water model and two-dimensional (2-D) Navier–Stokes simulations. Good agreement between the shallow-water model and the 2-D simulations was obtained up to the end of the inertial phase as the viscous effects were not included in the shallow-water model. Ungarish & Huppert (2008) then analysed the energy balances of an axisymmetric gravity current in the homogeneous ambient and linearly stratified ambient using the shallow-water model and Navier–Stokes finite difference simulations within a two-dimensional geometry. A fair agreement on the energy changes of the current was also obtained between these two approaches. Similar observations reported by Ungarish & Huppert (2006) were obtained in this study too where the stratification in the ambient enhanced the accumulation of potential energy and reduced the energy dissipation in the system.

Previous literature has explored the dynamics of gravity currents in unstratified ambients with varying aspect ratios and depth ratios of heavy fluid on a horizontal plane (or on a uniform slope). The front location, front velocity and theoretical or empirical models to predict the front velocity during the slumping, inertial and viscous phases have been determined for both planar and cylindrical currents in an unstratified ambient. Some studies have also investigated the effects of stratification on planar currents, but to our knowledge, no studies have investigated fully three-dimensional (3-D) cylindrical gravity currents in a stratified ambient. Such studies are important because cylindrical currents behave differently from planar currents. For example, the spreading rate of a planar current increases linearly, while that of a cylindrical current increases quadratically. Furthermore, the propagation of a current can be affected by a stratified ambient, as reported by Birman et al. (2007b), Lam et al. (2018a,b, 2022a,b) and Dai et al. (2021).

This study aims to investigate the effects of the strength of the stratification, $S$, and Reynolds number, $Re$, on the dynamics of the cylindrical release gravity current flow on the horizontal plane. In this work, we present the results from fully 3-D DNSs of cylindrical gravity currents propagating in a linearly stratified ambient with varying stratification at moderate Reynolds numbers. A detailed study on predicting the front velocity in the different phases as well as the scaling for the front velocity for both unstratified and stratified cases in the slumping phase is explored. Analysis of the density contours is used to compare the flow structures at different Reynolds numbers and stratification strengths. The three-dimensional structure of the advancing front is compared between the unstratified and stratified cases. The paper is structured as follows. Section 2 describes the numerical procedure including the governing equations, initial and boundary conditions, and the formulation of the problem. In § 3, the theoretical and empirical models in predicting the front velocity in the different phases and the impact of stratification on the front velocity in the slumping phase will be analysed. In § 4, we present our simulation results as well as a discussion of these results. Finally, conclusions are drawn in § 5.

2. Computational set-up

The 3-D, cylindrical release gravity currents in a stratified ambient have been simulated using Nek5000, a spectral element, incompressible flow solver (Fischer, Lottes & Kerkemeier 2008). Here, we adopt the Boussinesq approximation where the density difference between two fluids is sufficiently small (less than $5\,\%$) (Turner 1979) to neglect the influence of density differences in the inertial and diffusion terms, and retains only in the buoyancy term (Dai et al. 2021; Cao, Philip & Ooi 2022). The non-dimensional governing equations employed in the study take the form

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

(2.2)\begin{gather}\frac{{\rm D}\boldsymbol{u}}{{\rm D}t}={-}\rho \boldsymbol{\hat{z}} - \boldsymbol{\nabla} p + \frac{1}{Re} \nabla^2 \boldsymbol{u}, \end{gather}

(2.3)\begin{gather}\frac{\partial \rho}{\partial t} + \boldsymbol{\nabla} \boldsymbol{{\cdot}} (\rho \boldsymbol{u}) = \frac{1}{Re Sc}\nabla^2\rho, \end{gather}

where $\rho$ is the density of the fluid, $\boldsymbol {u}=[u,v,w]$ is the velocity for 3-D flow, $p$ is pressure and $\boldsymbol {\hat {z}}$ is the unit vector in the positive $z$-direction. The dimensionless density, $\rho$, is defined as

(2.4)\begin{equation} \rho=\frac{\rho^* - \rho^*_0}{\rho^*_c - \rho^*_0}, \end{equation}

where the symbols $\rho ^*$, $\rho ^*_0$ and $\rho ^*_c$ with asterisks are the dimensional density of the local, top of the domain and heavy fluid, respectively. The value of $\rho$ is bounded between 0 and 1. In the ambient fluid, the dimensionless density in the ambient $\rho _a$ varies linearly with depth $z$ from $\rho _a=\rho _b=S$ at the bottom $(z=0)$ to $\rho _a=\rho _0=0$ at the top $(z=1)$ and

(2.5)\begin{equation} \rho_{a}=\rho_{a}(z,t=0)=S(1-z). \end{equation}

The Schmidt number is $Sc=\nu ^*$/$\kappa ^*$ (where $\nu ^*$ is the kinematic viscosity and $\kappa ^*$ is the molecular diffusivity). Although saline, which is typically used in experiments, has $Sc=700$, it is found that when $Sc$ is of the order of 1 or larger, there is a weak scaling with the dynamics of the gravity current that does not significantly affect the bulk flow results (Härtel, Meiburg & Necker 2000; Necker et al. 2005; Cantero et al. 2007b; Bonometti & Balachandar 2008; Dai 2015). Therefore, $Sc=1$ is used in all simulations in the study.

The height of the domain $H^*$ is taken as the length scale. The velocity scale, $U^*$, time scale, $T^*$, and the Reynolds number, $Re$, are defined as

(2.6a–c)\begin{equation} U^*=\sqrt{g^\prime H^*}, \quad T^*=\frac{H^*}{U^*}, \quad Re=\frac{U^*H^*}{\nu^*}, \end{equation}

where $g^\prime =g^*(\rho _c^*-\rho _0^*)/\rho _0^*$ is the reduced gravity and $g^*$ is the gravitational acceleration acting in the negative $z$ direction.

Figure 1 shows the initial configuration of a full-depth cylindrical-release gravity current in a linearly stratified ambient with $S=0.5$. The horizontal directions are represented by $x$ and $y$, while the wall-normal direction is represented by $z$. The computational domain is a square prism with $L_x = L_y = 30$ ($L_x = L_y = 20$ for the cases with higher Reynolds number to reduce computational cost). A cylindrical lock contained heavy fluid with a density of $\rho _c$ has a height $H=1$ and unit radius ($r_0$) is contained at the centre of the computational domain. The density of the ambient ($\rho _a$) is increased linearly from the top ($\rho _0)$ to the bottom ($\rho _b$). A no-slip boundary condition is employed at the bottom ($z = 0$) of the domain and a slip, impermeable symmetry boundary condition is applied at the top of the domain ($z=1$) and vertical side walls ($x=[-L_x/2,L_x/2]$ and $y=[-L_y/2,L_y/2]$).

Figure 1. Sketch of the computational domain for the 3-D simulation. he horizontal directions are represented by $x$ and $y$, while the wall-normal direction is represented by $z$. The cylindrical region of heavy fluid located in the centre of the domain has a density of $\rho _c$. The heavy and ambient fluid has the same height as the height of the domain $H$. The density of the ambient $\rho _a(z)$ increases linearly from the top $\rho _0$ to the bottom boundary $\rho _b$ as indicated by the lighter grey shading and the $\rho _a(z)$ shown on the top left wall.

In this study, we have systematically investigated the effects of the stratification strength on the dynamics of the cylindrical gravity current. Four stratification strength of $S = 0$, 0.2, 0.5 and 0.8 are considered and simulated at Reynolds numbers of $Re = 3450$ and 10 000. Here, seventh-order polynomial spectral elements are used in this study to reduce numerical error in the simulation. The number of spectral elements employed for the simulations at low and high $Re$ are $N_x \times N_y \times N_z = 190 \times 190 \times 15$ and $270 \times 270 \times 20$, respectively. The grid distribution within the spectral element follows the Gauss–Legendre–Lobatto (GLL) grid spacing. This total number of unique grid points is $1.9\times 10^8$ and $5\times 10^8$ grid points. Grid stretching (geometrical progression with power coefficient) is applied along the wall-normal direction ($z$) where the grid size at the bottom part is denser than the top. The computational grid has a grid spacing of $0.0033 \leq \Delta x = \Delta y \leq 0.0332$ for $Re=3450$ and $0.0021 \leq \Delta x = \Delta y \leq 0.0165$ for $Re=10\,000$. The grid spacing to Kolmogorov scale ratio, $\Delta l/\eta$ (where $\Delta l = (\Delta x \Delta y \Delta z)^{1/3}$ and $\eta$ is the local dissipation) is calculated at different instantaneous time and is always less than 10. This is more conservative than the $\Delta l/\eta \approx 16$ recommended by Lu et al. (2023) and Zahtila et al. (2023), who studied the grid convergence characteristics of spectral element solvers. Therefore, we have ensured that our grid resolution is sufficient to resolve all of the turbulent length scales and also meet the requirement of $\Delta x = \Delta y \approx (ReSc)^{-1/2}$, where $Sc = 1$ (Härtel et al. 2000; Birman, Martin & Meiburg 2005; Dai 2015). A variable time step is used to ensure that the Courant number is always less than 0.5 and the simulation is run for a total time of $T=60$.

3. Theoretical predictions of front velocity

This section describes the theoretical and empirical models for predicting the front velocity of a cylindrical gravity current during the slumping, inertial and viscous phases as it propagates into an unstratified ambient. In the following sections of this paper, we will discuss the effect of stratification on the prediction of the front velocity during the slumping phase.

3.1. Slumping phase

Several theoretical and empirical models have been proposed to predict the front velocity during the slumping phase. First, Benjamin (1968) derived the Froude number of the front using hydraulic theory. The expression of the current speed in terms of Froude number based on the depth $H^*$ of the ambient fluid is expressed as

(3.1)\begin{equation} Fr_B = \frac{u_f^*}{\sqrt{g^\prime H^*}} = \sqrt{\tilde{h}(2-\tilde{h}) \frac{(1-\tilde{h})}{(1+\tilde{h})}}, \end{equation}

where $u_f$ is the constant front velocity and $\tilde {h}=h^*/H^*$ is the dimensionless current height. In the limit of $\tilde {h} \rightarrow 0.5$, $Fr_B = 0.5$. Shin et al. (2004) extended Benjamin's theory by taking into consideration the exchange between the advancing front of the gravity current and the backward disturbance. Both full-depth and partial-depth releases are taken into account in the study, and the speed of the current has the expression of

(3.2)\begin{equation} Fr_S = \tfrac{1}{2}\sqrt{\tilde{D}(2-\tilde{D})}, \end{equation}

where $\tilde {D}=D^*/H^*$ is the dimensionless initial depth of the heavy fluid. In the limit of full-depth release, $\tilde {D}=1$, the above expression yields the same results as (3.1).

Huppert & Simpson (1980), however, reported an empirical expression for the Froude number at the gravity current head based on the experimental measurement as

(3.3)\begin{equation} Fr_{HS} =\begin{cases} \frac{1}{2}\tilde{h}_{HS}^{1/6} & (0.075 \leq \tilde{h}_{HS} \leq 1),\\ 1.19\tilde{h}_{HS}^{1/2} & (0 \leq \tilde{h}_{HS} \leq 0.075), \end{cases} \end{equation}

where $\tilde {h}_{HS}$ is the dimensionless depth of the current just behind the head. In the limit of $\tilde {h}_{HS} \rightarrow 0.5$, the value of $Fr_{HS} \approx 0.45$ is slightly lower than the value of 0.5 predicted by Benjamin (1968) and Shin et al. (2004). Here the subscript $B$, $S$ and $HS$ denotes Benjamin (1968), Shin et al. (2004) and Huppert & Simpson (1980). Previous studies (Bonnecaze et al. 1995; Ungarish & Huppert 1998; Hallworth, Huppert & Ungarish 2001) that study axisymmetric currents have adopted the empirical Froude number reported by Huppert & Simpson (1980) (see (3.3)) in their studies and the results are consistent with the results obtained using different methods such as numerical simulation, experimental measurement and shallow-water approximations. Cantero et al. (2007b), however, reported that the three-dimensionality of the current does not have a strong influence on the speed of the current during the slumping phase but becomes important during the inertial and viscous phases. Numerical results from Cantero et al. (2007b) show a good agreement with the experimental data from Hallworth et al. (2001) when compared with the planar theory. Therefore, the slumping velocity from planar theory will be used for comparison with current simulation results.

The theoretical and empirical models for the mean front velocity in the slumping phase for unstratified gravity currents overestimate the front velocity of gravity currents in a stratified ambient. Ungarish & Huppert (2002) claimed that the mean front velocity in the slumping phase in the Froude condition reported by Huppert & Simpson (1980) (see (3.3)) is a function of $\tilde {h}_{HS}$ with stratification strength ($S$) and height of channel ($H$) as parameters. The proposed mean front velocity in the slumping phase for the gravity current in a stratified ambient is

(3.4)\begin{equation} Fr_{UH}(\tilde{h},S) = Fr_\xi(\tilde{h})(1-S+\tfrac{1}{2}\tilde{h}S)^{1/2}, \end{equation}

where $Fr_\xi (\tilde {h})$ is the Froude number and the subscript $\xi$ denotes $B,HS$ and $C$ (corresponding to Benjamin 1968; Huppert & Simpson 1980 and Cantero et al. 2007b), respectively.

Ungarish (2006) adapts Long's model (Long 1953, 1955; Baines 1995) and the Froude number is defined as

(3.5)\begin{equation} Fr_U(\tilde{h},S)=Fr_\xi(\tilde{h})(1-\tfrac{2}{3}S)^{1/2}. \end{equation}

Here, the subscript $UH$ and $U$ denotes Ungarish & Huppert (2002) and Ungarish (2006). Ungarish (2006) concludes that $Fr$ rescaling is useful for the case with weak stratification but significantly underestimates $Fr$ when $S/\tilde {h}>5$ (larger $S$ and/or small $\tilde {h}$). Birman et al. (2007b) conducted two-dimensional simulations of planar gravity currents propagating in a stratified ambient with initial fractional depth $\tilde {h}_0 = 0.2$, 0.5 and 1 with $S$ ranging from 0.2 to 0.9 at $Re=4000$. The front velocity of the current in the slumping phase as a function of $S$ was compared with the theoretical model by Ungarish (2006). Birman et al. (2007b) reported that the model is found to predict well in determining the Froude number with weak stratification. For strong stratification ($S>0.5$), multiple solutions exist and the solutions deviate from the first theoretical solution but agree closely to the second or third theoretical solutions.

3.2. Inertial phase

After the slumping phase, the gravity current in an unstratified ambient transitions into the first self-similar regime, the inertial phase. Here, the reflected backward propagating wave catches the front with a speed slightly greater than the speed of the front (Rottman & Simpson 1983). During the inertial phase, the buoyancy force is balanced by the inertial force and the front velocity of the current decreases as a power law, $t^{-1/2}$ for the cylindrical current (Cantero et al. 2007b) in the inertial phase. Models for gravity currents in the inertial phase have been proposed by Fay (1969), Fannelop & Waldman (1972), Hoult (1972), Huppert & Simpson (1980) and Rottman & Simpson (1983). The asymptotic behaviour of the cylindrical current in the inertial phase is

(3.6a,b)\begin{equation} r_f = {\rm \pi}^{1/4}\xi_{ip} h_0^{1/4}(r_0t)^{1/2},\quad u_f = \tfrac{1}{2}\xi_{ip} h_0^{1/4}r_0^{1/2}t^{{-}1/2}, \end{equation}

where $r_f$ is the radial front location of the cylindrical current, $\xi _{ip}$ is a constant for the cylindrical current and the subscript $ip$ denotes the inertial phase. The initial height and radius of the lock are $h_0$ and $r_0$. Hoult (1972) and Huppert & Simpson (1980) proposed $\xi _{ip}= 1.3$ and 1.16, respectively. Cantero et al. (2007b) investigates the transition between the different phases on 3-D cylindrical gravity currents at $Re= 895$, 3450 and 8950. A range of experimental data was used to revisit the scaling during the phases of spreading for cylindrical release gravity current. The best-fit value of the constants $\xi _{ip}$ was found to be 1.23 for the cylindrical current.

3.3. Viscous phase

Finally, the gravity current reaches the second self-similar regime, the viscous phase when the viscous force becomes significant compared with the buoyancy force. Fay (1969) and Hoult (1972) predict the self-similar solution for the cylindrical current using a boundary layer approximation

(3.7a,b)\begin{equation} r_f = \xi_{vp, Ht}h_0^{1/3}r_0^{2/3}Re^{1/12}t^{1/4},\quad u_f = \tfrac{1}{4}\xi_{vp, Ht}h_0^{1/3}r_0^{2/3}Re^{1/12}t^{{-}3/4}, \end{equation}

where the subscript $vp$ denotes the viscous phase. The solution from Hoult (1972) has a constant value of $\xi _{vp, Ht}=1.38$. Huppert (1982) revised the analysis by considering the viscous effect over a rigid horizontal surface. The self-similar solutions obtained for the viscous phase differ from those presented in (3.7) and are expressed as

(3.8a,b)\begin{equation} r_f = \xi_{vp, Hp}h_0^{3/8}r_0^{3/4}Re^{1/8}t^{1/8},\quad u_f = \tfrac{1}{8}\xi_{vp, Hp}h_0^{3/8}r_0^{3/4}Re^{1/8}t^{{-}7/8}, \end{equation}

where $\xi _{vp, Hp}=1.197$. Here the subscript $Ht$ and $Hp$ denotes Hoult (1972) and Huppert (1982), respectively.

4. Results and discussion

4.1. Equivalent height of the gravity current

The equivalent height, defined as the azimuthally averaged integrated height of the current in the wall-normal direction ($z$), is defined by

(4.1)\begin{equation} \bar{h}(r,t) =\frac{1}{2{\rm \pi}}\int_{0}^{H}\int_{0}^{2{\rm \pi}}\rho(r,\theta,z,t)\, {\rm d}\theta \,{\rm d}z - \int_{0}^{H}\rho_a(z)\,{\rm d}z. \end{equation}

Note that in this definition of equivalent height, the total mass of the gravity current is $1/R\times \int ^R_0 \bar {h}r\, {\rm d}r$, where $R$ is the radius of the domain. The equivalent height provides a good indication of the distribution of mass in the streamwise direction and the shape of the gravity current (Shin et al. 2004; Marino et al. 2005; Birman et al. 2007a; Cantero et al. 2007b; Dai 2015). However, when stratification is introduced to the ambient fluid, it changes the equivalent height ($\bar {h}$) by an amount proportional to the stratification strength $S$. This is reflected in the second term of (4.1), which is equal to $S/2$ and is zero in an unstratified case where $\rho _a=0$. It should be noted that as the strength of stratification increases, the interface between the intrusion and the ambient becomes smoother and less distinguishable (Sutherland et al. 2022). Before we examine the propagation in detail, we analyse the density contours and equivalent height of the gravity current with different $S$ at $Re=3450$ and 10 000 to characterise the flow in the slumping and inertial phases, the two most persistent and interesting phases.

Figure 2 shows the time evolution of the azimuthally averaged density contour of the gravity current at $Re= 3450$ for $S=0$, 0.2, 0.5 and 0.8, respectively. The equivalent height ($\bar {h}$) of the current is plotted on top of the density contour ($\rho =0.015$). In the slumping phase, strong K–H vortices are formed behind the head of the gravity current. The rolled-up vortices entrain ambient fluid and mix with it, and these structures are reflected in the local maximum of the equivalent height. The height of the current in the head region does not differ significantly during the slumping phase for $S\leq 0.5$. Weaker vortices are formed in the case with strong stratification ($S=0.8$) and the current head results in a smaller head size than the other cases. Also, we can observe the head is slightly lifted from the bottom wall in all cases and in all phases due to the no-slip condition, and a layer of ambient fluid is entering and mixing with the head. A similar observation is also reported by Cantero et al. (2007a). At $t=10.4$, the current with $S=0$ has transitioned into the inertial phase and the K–H billow has mixed with the head. A similar observation is observed for $S=0.2$ and 0.5 in the inertial phase. At a similar time $t=10$ for $S=0.8$, the K–H billow is still presented behind the current and this shows that stratification delayed the transition of the current.

Figure 2. Azimuthal-averaged density contour of the gravity current ($\rho =0.015$) with $S=$ (a) 0, (b) 0.2, (c) 0.5 and (d) 0.8 at $Re= 3450$. The red box shows the head is lifted from the bottom wall. The solid black line shows the equivalent height of the gravity current ($\bar {h}$). SP, slumping phase and IP, inertial phase.

4.2. Front location and front velocity

The front location of the gravity current is typically defined as the location where the equivalent height, $\bar {h}$, is smaller than an arbitrarily small threshold $\delta$ (Cantero et al. 2007b; Dai 2015; Zhu et al. 2017; Chan, Lam & Ooi 2018). However, for a gravity current in a stratified ambient, this method leads to an overestimation in the front location due to a displaced bottom layer. Therefore, the gradient of the equivalent height method is used as it is more robust for flows with strong stratification. This method is explained in more detail and its accuracy is compared with a few different methods in Appendix A.

Figure 3(a) shows the plot of the front location for the cases with stratification strength ranging from 0 to 0.8 at $Re= 3450$ (solid lines) and 10 000 (open circles). The front location of the gravity current is plotted until the front has mixed with the ambient fluid and is not detectable. Increasing the linear stratification strength leads to a reduction in the distance travelled by the front of the gravity current at the end of the slumping phase. At $t=16$, which is approximately when the gravity current for the case $S=0.8$ at $Re=3450$ has dissipated and there is no appreciable or discernible gravity current flow any more, the gravity current propagated for approximately 3.8 lock radii ($r_f \approx 4.8$), which is approximately $22\,\%$ smaller than the distance travelled in the unstratified ambient case ($S=0$), where the current travelled approximately 4.6 lock radii ($r_f \approx 5.6$). Similarly, in the high-$Re$ case, for $S=0.8$, the distance travelled was approximately 4.2 lock radii ($r_f=5.2$), representing a reduction of approximately $17\,\%$ compared to the unstratified case, where the travel distance was 4.9 lock radii ($r_f \approx 5.9$) at $t=17$.

Figure 3. Plot of (a) front location, (b) front location with time offset ($t_0$) and (c) front velocity against time for stratified cylindrical release gravity current with different stratification strength ranging from $S= 0$ to 0.8 at $Re= 3450$ (—) and 10 000 ($\circ$). The stratification is represented with a different colour where red, $S= 0$; green, $S= 0.2$; blue, $S= 0.5$ and cyan, $S= 0.8$. The predicted front location using the theoretical models for $S=0$ at $Re=10\,000$ are included in the plot where ($*$), inertial phase in (3.6); ($\vartriangle$), viscous phase by Hoult (1972) in (3.7); ($\triangledown$), viscous phase by Huppert (1982) in (3.8) with $Re=10\,000$, $h_0=1$ and $r_0=1$. The dashed line shows the maximum speed of the linear internal gravity wave, $N^*H^*/{\rm \pi}$ for different stratification strengths: green, $S= 0.2$; blue, $S= 0.5$ and cyan, $S= 0.8$.

The predicted front location for $S=0$ using the scaling in the inertial and viscous phase is also included in figure 3(b) and is represented with different line styles. It is important to note that the physical simulation time will not correspond to the virtual time origin (Samasiri & Woods 2015; Sher & Woods 2015). Therefore, we have estimated a value of the virtual origin in time, $t_0$, for all the cases in the inertial phase using the method proposed by Sher & Woods (2015) and Samasiri & Woods (2015). Figures 3(b) and 4 show the front velocity profiles collapse in the inertial phase irrespective of the Reynolds number and stratification strength, at least for the range of $Re$ and $S$ considered in this study. The front location of the gravity current for both low and high $Re$ with $0\leq S<0.5$ cases follows the inertial phase scaling of $t^{1/2}$ (represented by $*$) reported by Huppert & Simpson (1980) (defined in (3.6)). For the strongly stratified cases ($S=0.8$), the propagation of the current becomes negligible during the slumping phase due to the insufficient density difference between the current and the ambient at the bottom wall. As a result, there is no longer any gravity current flow.

Figure 4. Plot of the front velocity against time for stratified cylindrical release gravity current with different stratification strengths ranging from $S= 0$ to 0.8 at $Re=$ (a) 3450 and (b) 10 000. The stratification is represented with a different colour where red solid line (circle)$,S=0$; green solid line (circle), $S= 0.2$; blue solid line (circle), $S= 0.5$ and cyan solid line (circle), $S= 0.8$. The theoretical model for $S=0$, (black solid line), slumping phase; ($\cdots \cdots$), inertial phase; (– – –), viscous phase by Hoult (1972) in (3.7); (-$\cdot$-$\cdot$), viscous phase by Huppert (1982) in (3.8).

The proposed scalings in the viscous phase by Hoult (1972) and Huppert (1982) (represented by the $\vartriangle$ and $\triangledown$) underestimate the front locations of the gravity current. However, the trend of the scalings shows a good agreement with the simulation result for both low- and high-$Re$ cases within the range of $0\leq S<0.8$. The front location for $Re$ with $S=0$ in the viscous phase evolves at approximately $t^{0.28}$. This is closer to the viscous scaling $t^{1/4}$ in (3.7), and it is larger than the viscous scaling $t^{1/8}$ in (3.8) by a factor of 2. For the stratified cases, the propagation of the current becomes negligible in the early stage of the viscous phase (or does not transition into the viscous phase). Therefore, no observation is made for those cases. Cantero et al. (2007b) optimised the prefactors of the scalings law and reported the empirical best-fit value of $\xi _{ip}=1.23$, $\xi _{vp, Ht}=2.6$ and $\xi _{vp, Hp}=4.64$ from the front velocity plot. With a slight modification of the prefactors, we found $\xi _{ip}=1.23$, $\xi _{vp, Ht}=2.96$ and $\xi _{vp, Hp}=6.56$ provides a fair agreement with our simulation results. However, these prefactors are not suitable for scaling laws in predicting the front location of the gravity current.

The front velocity of the gravity current is determined by applying a moving-average filter with filter width of $\Delta t=1$ to the front location and then differentiating the filtered front location with respect to the time ($u_f = {\rm d}r_f/{\rm d}t$). The plot of the front velocity against time for all cases is presented in figure 3(c). In the cylindrical release gravity current, similar to a planar release, we observe the propagation of gravity current undergoes four main flow phases as reported by Blanchette et al. (2005), Cantero et al. (2007a,b, 2008) and these phases are discussed in the following.

The first phase is the initial acceleration where the front velocity accelerates from 0 to a maximum value. The maximum front velocity increases with $Re$ but decreases as $S$ increases. As $S$ increases ($S = 0, 0.2, 0.5, 0.8$), the peak value of the front velocity decreases from 0.563 to 0.512, 0.485 and 0.414 at $Re=3450$ and from 0.594 to 0.545, 0.514 and 0.437 at $Re=10\,000$. It is also noted that the duration of the acceleration phase is also extended with increasing $S$ with the gravity current taking a longer time to reach the maximum front velocity. Comparing the duration for the strongly stratified and unstratified cases to reach the maximum front velocity, the case with $S = 0.8$ takes approximately $43\,\%$ longer than $S = 0$ ($t \approx 2.0$ versus $t \approx 1.4$).

After the initial acceleration, the gravity current undergoes a slumping phase. In this regime, the front velocity remains relatively steady. Cantero et al. (2007b) investigated the transition between different phases of spreading by collecting experimental data from the literature and revisited the scaling of the slumping phase. They reported the best-fit value of the mean front velocity in the slumping phase is $Fr_{C}=0.42\tilde {h}_0^{-1/2}$ (where the subscript $C$ denotes Cantero et al. (2007b) and $\tilde {h}_0=h_0/H$ is the depth of initial release). A period of near-constant velocity was observed in the planar currents, but was not clearly identified for the cylindrical current. However, there is a brief period of slower variation before the current exhibits a more pronounced decay (Cantero et al. 2007b). Our simulation results are in good agreement with Cantero et al. (2007b) where the front velocity in the slumping phase decreases at a very slow pace instead of approaching a period of near-constant velocity. The slumping phase is found to occur at a later time at both $Re$ when $S$ increases where the cases with $S=0$ begin at $t \approx 2$, and $t \approx 2.2$, 3 and 4 for the cases with $S= 0.2$, 0.5 and 0.8, respectively (see figure 3c). For a better comparison, the near constant velocity region is used to calculate the mean front velocity. Here, $Fr_{sim} = u_{f,mean}\tilde {h}_0^{-1/2}$ is the Froude number of the gravity current based on the mean front velocity in the slumping phase (table 1). The front velocity remains constant for a longer period as the stratification strength increases for both Reynolds numbers.

Table 1. Mean front velocity in the slumping phase of present simulations expressed as $Fr_{sim}$ with different stratification strength at $Re= 3450$ and 10 000. $Fr_{UH}$ and $Fr_U$ denote the Froude number reported by Ungarish (2006) ((3.5)) and Ungarish & Huppert (2002) ((3.4)), respectively. The flow regime of present simulation are determined by the buoyancy Froude number $Fr=u^*_{f,mean}/N^*H^*$.

For the cases with $S=0$, the mean front velocity in the slumping phase has a value of 0.38 $(Re=3450)$ and 0.41 $(Re=10\,000)$, and gradually decreases to 0.26 $(Re=3450)$ and 0.27 $(Re=10\,000)$ when the stratification strength increases to 0.8. The mean front velocity in the slumping phase has dropped approximately $32\,\%$ in strong stratification and thus $Fr_{sim}$ decreases as $S$ increases.

The flow regime of the gravity current propagating in a stratified ambient is determined by the buoyancy Froude number calculated using the slumping velocity. When the current is moving with a value of $Fr > 1/{\rm \pi}$, it is referred as a supercritical flow and in this regime, the current travels faster than the internal gravity waves. When gravity current has a $Fr$ value smaller than $1/{\rm \pi}$, the gravity current flow is in the subcritical regime. However, the release of the cylindrical gravity current is a developing flow and the local Froude number, based on the instantaneous velocity, changes with time. In the initial acceleration phase, the front velocity increases rapidly from zero and for all cases, is greater than the speed of the internal gravity wave (see figure 3c with the dashed line being the velocity of the linear, mode-one, long internal gravity wave). The gravity current then transitions into the slumping phase after the front velocity reaches a local maximum. In the slumping phase, the front velocity of the current with $S=0.8$ is less than the speed of the internal gravity wave ($N^*H^*/{\rm \pi}$) and is in the subcritical regime. As the flow decelerates in the inertial regime, the flow of the other less stratified cases eventually enters the subcritical regime based on the local Froude number. In the subcritical regime, the oscillation of forward- and backward-propagating internal waves behind the gravity current is more pronounced compared with the supercritical flow. Further discussion on this phenomenon is presented in the following sections.

For the full-depth release and with the limit of $\tilde {h}$ and $\tilde {h}_{HS} \rightarrow 0.5$, the Froude number of the gravity current ($Fr_{sim}$) for $S= 0.2$, 0.5 and 0.8 is calculated using (3.4) and (3.5). The scaling proposed by Ungarish & Huppert (2002) and Ungarish (2006) takes into consideration the effects of stratification, where the front velocity decreases with increments of $S$. We recall the best-fit value of the Froude number $Fr_C=0.42\tilde {h}_0^{-1/2}$ reported by Cantero et al. (2007b) and adapt it to the scaling. The results are tabulated in table 1. Average differences of $6.4\,\%$ and $7.2\,\%$ are observed for the low-$Re$ cases across the stratified and unstratified cases compared with the models reported by Ungarish & Huppert (2002) and Ungarish (2006), respectively. The percentage difference decreases to $1.9\,\%$ and $3\,\%$ when the Reynolds number increases to $Re=10\,000$. Both models can capture the stratification effect and provide good predictions for $Fr_{sim}$, however, the scaling reported by Ungarish & Huppert (2002) shows a better agreement compared with that by Ungarish (2006). However, when using the classical formula of Benjamin (1968), $Fr_B(\tilde {h})$ (the second term on the left-hand side of (3.4) and (3.5), it overestimates $Fr_{sim}$ by approximately $30\,\%$ and $22\,\%$ at $Re= 3450$ and 10 000, respectively.

Figure 4 shows the plot of the front velocity with different stratification strengths and different phases at $Re=$ (a) 3450 and (b) 10 000. The gravity current transitions into the inertial phase when the inertial force balances the buoyancy force. For the low-$Re$ cases, the current leaves the nearly constant velocity phase and transitions into the inertial phase at $t= 4.6$, 5, 6.8 and 11.4 with $S =0$, 0.2, 0.5 and 0.8, respectively. The current enters the inertial phase at $t= 3.8$, 4, 5 and 11 for the high-$Re$ case with $S= 0$, 0.2, 0.5 and 0.8, respectively. The front velocity follows the model proposed by Huppert & Simpson (1980) (as defined in (3.6)) with good agreement up to $t= 17$ and 27 for the case with $S= 0$ and 0.2. The density of the current becomes close to the density at the bottom boundary for $S= 0.2$, 0.5 and 0.8 before the end of the inertial phase, and the propagation of the current is not observed when the heavy fluid has mixed with the ambient fluid.

A comparison of the front velocity for the unstratified and stratified cases with the scalings in the viscous phase for $Re=3450$ and 10 000 is also shown in figure 4. For the $S=0$ case at $Re= 3450$ and 10 000, the front velocity leaves the inertial phase at $t= 17$ and 20 and the front velocity decays faster in the viscous regime. In the $S=0.2$ case, the gravity current transitions from the supercritical regime to subcritical gravity current. Here, the internal gravity waves travel faster than the current, pushing the head of the current which results in a reduction in the decay rate of the front velocity. Consequently, the current in the $S=0.2$ case travels further than the unstratified case at later times.

Both the scaling laws for the viscous phase proposed by Hoult (1972) and Huppert (1982) in (3.7) and (3.8) underpredict the simulation results. The velocity of the front in the viscous phase from the simulations is higher than what the scaling laws predict. This is also true when considering the virtual time origin calculated in the viscous phase. It is important to note that for $S = 0.5$ and 0.8, the gravity current flow does not enter the viscous phase as the dense fluid has fully mixed with the ambient fluid at the bottom of the domain. Because the viscous scaling proposed by Hoult (1972) and Huppert (1982) were derived for 2-D gravity currents, they underpredicted the cylindrical gravity current simulations where the three-dimensionality of the current has a strong influence over the propagation speed (Huppert 1982; Cantero et al. 2007b). Similar observations were reported by Cantero et al. (2007b) and best-fit values of the prefactors $\xi _{vp, Ht}=2.6$ and $\xi _{vp, Hp}=4.64$ (see (3.7) and (3.8)) is proposed without modifying the theoretical power law exponents.

For the strongly stratified cases, the current has a similar density close to the bottom boundary ($\rho _c \approx \rho _b$) before transitioning into the viscous phase. Neither of the cases with $S=0.5$ and 0.8 undergoes a transition into the viscous phase. In the case with $S= 0.5$, the density of the current becomes close to the density at the bottom boundary ($\rho _c \approx \rho _b$) before transitioning into the viscous phase. Additionally, the propagation of the current for $S=0.8$ becomes negligible in the slumping phase where front velocity decreases rapidly after the end of the slumping phase. A similar finding was observed for $Re=10\,000$ cases. The propagation of the current becomes negligible after the slumping phase as there is insufficient density difference between the current and ambient at the bottom wall to continue propagating.

4.3. Space-time diagram of the equivalent height $\bar {h}$

The evolution of the azimuthally averaged cylindrical gravity current in a stratified ambient at $Re=3450$ and 10 000 as visualised by the contour of the equivalent height ($\bar {h}$) of the gravity current is shown in figures 5 and 6. The heavy fluid is coloured red while the ambient fluid is represented in blue. The front location of the gravity current is plotted with a solid black line. The red dashed line has a slope corresponding to the maximum speed of the linear internal gravity wave and is plotted for stratified cases only. The crossing white dotted lines, forming a ‘hash’ like pattern, observed on the left side of figure 5(b–d) represent the presence of the forward- and backward-propagating internal gravity waves (refer to the inset plot in figure 5(d) for better visualisation).

Figure 5. Contour of the equivalent height of the gravity current with $S=$ (a) 0, (b) 0.2, (c) 0.5 and (d) 0.8 at $Re=3450$. The heavy and ambient fluid is coloured red and blue, respectively. The solid black line, front location of the current ($r_f$); red dashed line, slope corresponding to the maximum speed of the linear internal gravity wave ($N^*H^*/{\rm \pi}$), front location of the internal gravity wave. The white dotted lines show the movement of the internal gravity wave behind the gravity current. The white arrow indicates the counter-clockwise rotating vortices behind the current head, which propagate in the opposite direction of the current. The red rectangle corresponds to the inset figure and shows the forward- and backward-propagating internal gravity waves behind the gravity current. The colour of the inset in panel (d) is saturated (the colour scale is five times smaller than the original) for better visualisation.

Figure 6. Contour of the equivalent height of the gravity current with $S=$ (a) 0, (b) 0.2, (c) 0.5 and (d) 0.8 at $Re= 10\,000$. The heavy and ambient fluid is coloured red and blue, respectively. The solid black line, front location of the current ($r_f$) and red dashed line, slope corresponding to the maximum speed of the linear internal gravity wave ($N^*H^*/{\rm \pi}$).

The ‘spikes’ observed on the left-hand side of the solid black line during the time interval $5 < t < 15$ indicate the occurrence of the breakdown of the K–H billows behind the gravity current head. These vortices behind the current head transport the heavy fluid towards the head during the early times and eventually break down into smaller-scale turbulent structures at later times. The head of the gravity current continues to propagate forward, while small vortices within the body of the current remain in the same location (or propagate at a slow pace), mix with the surrounding fluid and eventually form the tail of the gravity current at later times. The tail of the current breaks down into multiple sections. Additionally, in figure 5(d), the white arrow indicates that the vortices are rotating counter-clockwise and propagating in the opposite direction of the current.

As the stratification strength $S$ increases to 0.8, the ‘hash’-like pattern continues to grow, indicating the presence of forward- and backward-propagating internal gravity waves behind the gravity current (see figure 5d) and is significant compared with other cases (particularly in the region $0 < r < 5$). In the initial acceleration phase, the gravity current with $S=0.8$ at both $Re$ has a propagation speed similar to the internal gravity wave. The red dashed line in figures 5(d) and 6(d) is almost tangent to the solid black line. Later, the internal gravity wave propagates faster than the gravity current and the red dashed line is plotted below the solid black line. When the gravity current propagates slower than the maximum speed of the linear internal gravity wave, the gravity current is in the subcritical regime. This agrees well with the plots in figure 3(b) where the maximum speed of the linear internal gravity wave is greater than the front velocity in the slumping phase with $S=0.8$. The results are qualitatively similar for the higher Reynolds number cases shown in figure 6.

4.4. Turbulent structures of the gravity current

The 3-D turbulent structure of the cylindrical current is visualised by the isosurface of swirling strength $\lambda _{2}$, which is the imaginary component of the complex eigenvalues of the local velocity gradient tensor (Zhou et al. 1999). In this section, the two extreme cases $S=0$ and 0.8 at $Re=10\,000$ are selected to compare the 3-D flow structures between the unstratified and the strongly stratified cases. Figure 7 shows the time evolution of the isosurface of $\lambda _{2}$ for the unstratified case at $Re=10\,000$.

Figure 7. Time evolution of the isosurface of $\lambda _{2}=-20$ for the case with $S=0$ at $Re=10\,000$. The dotted line represents the location of the spanwise vorticity ($\omega _z$) contour on the $x$–$z$ plane. A closer look at the black and red square region is shown on the left.

At time $t=3.2$ (figure 7a), radially aligned hairpin vortices are formed at the head of the gravity current due to the high shear at the bottom of the no-slip wall. These structures are similar to the tooth-like vortices reported by Dai & Huang (2022) in a planar release gravity current. They found that the merging and splitting of the lobes and clefts of the gravity current are a result of the merging and birth of these vortices.

As the gravity current slumps, inner and outer vortex rings are generated by the gravity current. These structures arise due to the K–H billows developing behind the head of the gravity current and counter-rotating vortex rings formed close to the wall due to the high shear of the near-wall region (see 2-D colour map of vorticity in figure 7). These vortex structures mutually interact with each other causing these structures to bend and tilt with some of these structures getting intertwined (see figure 7a).

Strong interaction between the counter-rotating vortex pairs leads to a chaotic and violent breakdown of the turbulent structures to smaller scales. Smaller secondary filaments perpendicular to the vortex rings are formed (figure 7c,d) and the mechanism leading to the breakdown of these structures is similar to the phenomena observed in a vortex ring collision. McKeown et al. (2020) described this mechanism to be an iterative cascade of the elliptical instability where energy from the larger scales is transferred to smaller scales by the formation of perpendicular structures before finally being dissipated by the smallest scales.

At time $t=5.4$ (figure 7d), where the gravity current is at the end of the slumping phase, the flow at the head of the gravity current appears to be fully turbulent.

When strong stratification is present in the ambient fluid ($S = 0.8$), there are subtle differences in the evolution of the vortical structure of the gravity current compared with the unstratified case. Time evolution of the isosurface of $\lambda _2$ for the cylindrical release gravity current propagating in a strongly stratified ambient at $Re = 10\,000$ is shown in figure 8. Due to the lower slumping velocity, the vortical structures generated by the gravity current are located closer to each other and have a lower circulation strength as the current propagates radially outward. Because of this, the breakdown of the inner vortical rings occurs first compared with the unstratified case where the breakdown of the inner and outer rings occur at roughly the same time (see figures 7c and 8c). The turbulent structures in the flow also appear to be larger than the unstratified case and this is due to the lower local Reynolds number of the head of the gravity current.

Figure 8. Time evolution of the isosurface of $\lambda _{2}=-5$ for the case with $S=0.8$ at $Re=10\,000$. The dotted line represents the location of the spanwise vorticity ($\omega _z$) contour on the $x$–$z$ plane. A closer look at the black and red circled region is shown on the left.

Figures 9 and 10 shows a quadrant of the total isosurface of density and swirling strength $\lambda _{2}$ for the high-$Re$ case with $S= 0$, 0.2, 0.5 and 0.8. The isosurface of density is plotted on the left and the isosurface of $\lambda _{2}$ is on the right and is coloured by the velocity magnitude of the current. For consistent comparison between the stratified cases, a constant non-dimensionalised density contour value of $\tilde {\rho }=(\rho -S)/(1-S)=0.015$ is selected. All the figures are plotted in the same phase (slumping and inertial phases) for a better comparison of the evolution of the gravity current with different stratification strengths. It is important to note that none of the stratified cases, whether at low or high Reynolds numbers ($Re$), transition into the viscous phase.

Figure 9. Time evolution of the isosurface of density (on the left) and $\lambda _{2}$ (on the right) for the case with $S=$ (a) 0 and (b) 0.2 at $Re=10\,000$. The lobe-and-cleft instability is formed at the gravity current head in SP, and the lobe and cleft can be seen clearly in IP. The bulging front indicates a lobe and the indentation between two lobes is a cleft circled with a black circle. SP, slumping phase and IP, inertial phase.

Figure 10. Time evolution of the isosurface of density (on the left) and $\lambda _{2}$ (on the right) for the case with $S=$ (a) 0.5 and (b) 0.8 at $Re=10\,000$. The lobe-and-cleft instability is formed at the gravity current head in SP, and the lobe and cleft can be seen clearly in IP. SP, slumping phase and IP, inertial phase.

The propagation of the current for the case with strong stratification ($S=0.8$) becomes negligible at the later time in the slumping phase and does not transition into the inertial phase (or viscous phase). The figures shown in 10(b) are the current in the slumping phase of spreading.

In the slumping phase, undulations at the front of the gravity current due to the lobe-and-cleft instability are clearly visible for cases with $S=0$, 0.2 and 0.5 but are less prominent for $S=0.8$. Radially aligned hairpin vortices formed at the head of the gravity current due to the high shear at the bottom of the no-slip wall can be observed clearly for $0\leq S<0.5$. However, for $S=0.8$, these hairpin structures are formed in the middle of the slumping phase ($t \approx 6$) and therefore we did not observe any in figure 10(b).

A K–H billow is formed behind the gravity current head for all the cases, corresponding to the rolled up vortices. For $S=0.5$ and 0.8, the vortices around the ring are less significant compared with at $S=0$ and 0.2. The interaction between the inner and outer vortex rings can be observed in the isosurface of $\lambda _2$ for all the cases in figures 9 and 10, but the interaction reduces with increasing stratification strength. The wrinkling of the isosurface of density is due to the turbulent structures developing within the gravity current and the morphology of the isosurface of density is related to the local entrainment.

During the inertial phase, as the current moves radially outward with a decelerating velocity, the turbulent structures undergo a breakdown into smaller scales for values of $0 \leq S \leq 0.5$. This process continues until the heavy fluid becomes completely blended with the surrounding fluid. For the case with $S=0.8$, the hairpin structures are forming in the head at $t=9.6$. The tail of the current begins to separate from the body as the heavy fluid has mixed with the ambient fluid and the head continues to propagate. The lobes and clefts can be seen clearly at the front of the head of the gravity current head when plotting the isosurface of density. The evolution of the radially advancing lobes and clefts structure of the current is discussed further in the next section.

4.5. Lobes and clefts structure

In the last five decades, there have been many experimental and numerical studies conducted to elucidate the formation of the lobes and clefts instability in planar (Simpson 1972; Härtel et al. 2000; Dai & Huang 2022) and cylindrical gravity currents (Cantero et al. 2007a) in the absence of stratification. Formation of the lobes and clefts is due to the no-slip impermeable boundary condition of the bottom wall (Simpson 1972). Simpson (1982) postulated that these structures occur due to the convective instability as the head of the gravity current flows over the less heavy fluid. Stability analysis of the flow at the head of a gravity current conducted by Härtel et al. (2000) showed that the development of the 3-D structures are due to a local instability at the leading edge of the front. Numerical simulation of a planar gravity current conducted by Dai & Huang (2022) noted the importance of tooth-like vortices (also observed in figures 7a,b and 8b,c) that have a very similar appearance to hairpin vortices generated in a turbulent boundary layer, on the merging and splitting of the lobes and clefts structures. When merging and splitting processes are presented, Simpson (1972) reported the empirical relationship for the lobe size for a planar current to have an expression of

(4.2)\begin{equation} \frac{\tilde{\lambda}_l}{\bar{h}_H}=7.4Re_F^{{-}0.39\pm0.02}, \end{equation}

where $\tilde {\lambda }_l$ is the mean wavelength of the lobe, $\bar {h}_H$ is the azimuthally averaged height of the gravity current head, $Re_F = Re \bar {h}_H u_f$ is the local instantaneous front Reynolds number and $u_f$ is the front velocity of the gravity current.

The evolution of the lobes and clefts are visualised by plotting the 2-D contour of density $\tilde {\rho }=0.015$ on the $x$–$y$ plane at height $z = 0.01$ (close to the bottom boundary) with a constant time interval of $\Delta t=0.8$. The cases with $S=0$ at $Re = 3450$ and 10 000 are plotted in figure 11 to investigate the effects of Reynolds numbers on these structures. The colour of lines represents the different phases of the gravity current flow. The density contour lines are initially circular and smooth; however, undulations are observed in between at the end of the acceleration phase and the start of the slumping phase (solid green lines). These undulations grow and develop into lobes and clefts. In the inertial (solid blue lines) and viscous phases (solid black lines), the splitting and merging of lobes are observed. At lower Reynolds number (figure 11a), the density contour in the viscous phase looks more symmetrical than in the high-$Re$ case which contains a larger range of turbulent scales.

Figure 11. Evolution of the front visualised by the contour of $\rho =0.015$ close to the bottom boundary at $Re=$ (a) 3450 and (b) 10 000 with $S=0$. The time separation between contours is $\Delta t = 0.8$. The colour represents the transition between phase: red, initial acceleration phase; green, slumping phase; blue, inertial phase and black, viscous phase.

Figure 12 shows the lobes and clefts structure of the advancing front with $S=0.2$ and 0.5 at $Re=10\,000$. Each contour line represents $\tilde {\rho }=0.015$. The dashed lines show the splitting and merging of the lobes and clefts. A cleft merges with the neighbouring cleft producing a new cleft and larger lobes. The splitting and merging process is observed more frequently in the self-similar inertial and viscous phase compared with the slumping phase. As the current propagates, a lobe splits into two smaller lobes (the distance between dashed lines becomes wider) and multiple clefts merge (the distance between dash lines becomes narrower) which is shown in figure 12(a). The splitting and merging of existing lobes-and-clefts with $S=0.2$ are significantly greater than in the unstratified case shown in figure 11(b).

Figure 12. Evolution of the front visualised by the contour of $\rho =0.015$ close to the bottom boundary at $Re=10\,000$ with $S=$ (a) 0.2 and (b) 0.5. The time separation between contours is $\Delta t = 0.8$. The colour represents the transition between phases: red, initial acceleration phase; green, slumping phase and blue, inertial phase. The dashed lines in panel (a) show the splitting and merging of the lobes and clefts.

The lobes and clefts structure for $S=0.8$ is plotted in figure 13. As the stratification strength increases, the 2-D density contour plot becomes smaller due to the density difference between the heavy fluid and the density at the bottom wall being smaller with increasing $S$. Also, the stratification strength slows down the transitions of the current, and gravity current with $S=0.8$ has a longer duration in the slumping phase compared with the other cases. The splitting and merging process is the slowest for the $S=0.8$ case and has the smallest lobe size. After the end of the slumping phase, the current continues to propagate for a short period before the gravity current has dissipated and there is no gravity current flow any more. The summary of the quantitative information of the lobe-and-cleft structure with different $Re$ and $S$ is tabulated in table 2.

Figure 13. Evolution of the front visualised by the contour of $\rho =0.015$ close to the bottom boundary at $Re=10\,000$ with $S=0.8$. The time separation between contours is $\Delta t = 0.8$. The colour represents the transition between phases: red, initial acceleration phase; green, slumping phase and blue, inertial phase.

Table 2. Quantitative information on the lobe-and-cleft structure. Here, $r_f$ is the radial location of the front, $N_l$ is the total number of lobes in 360$^\circ$ and $\tilde {\lambda }_l=2{\rm \pi} r_f/N_l$ is the mean wavelength of the lobe. SP, slumping phase; IP, inertial phase; VP, viscous phase.

The number of lobes, $N_l$, against time for different $S$ at different $Re$ is plotted in figure 14. The number of lobes is calculated by counting the number of local maxima in the density contour. Spectral analysis of the profile was also conducted but did not result in a distinct peak in the frequency domain due to the diverse distribution of lobe sizes, thus leading to a misleading lobe size. The number of lobes for the case with $S=0$ and 0.2 at low $Re$ does not differ significantly until the current reaches $r_f \approx 5$ and starts to increase after that. For $S=0.5$, $N_l$ increases as the current propagates radially out. The number of lobes for $S=0.8$ remains nearly constant where the gravity current has dissipated and there is no appreciable or discernible gravity current flow any more. For higher Reynolds number cases, the number of lobes for $S=0$ and 0.2 first decreases until $r_f \approx 5$ (at the end and in the middle of the inertial phase). The increase in number $N_l$ as the currents propagate is observed for $S=0.5$ and 0.8 until the gravity current has dissipated. The gravity current at a higher Reynolds number becomes strongly turbulent due to the rapid growth of 3-D disturbances. In addition, the mixing within the current head is significantly greater compared with low-$Re$ cases. As a result, the number of lobes observed in the gravity current increases compared with cases at low Reynolds numbers.

Figure 14. Plot of the number of lobes ($N_l$) against radial front for stratified cylindrical release gravity current with different stratification strength at $Re=$ (a) 3450 and (b) 10 000. The stratification is represented with a different symbol and a cubic spline curve is fitted through the data, $\circ$(—), $S= 0$; *(- - -), $S= 0.2$; $\square$(-$\cdot$-$\cdot$), $S= 0.5$; $\triangle$($\cdots \cdots$), $S= 0.8$. The colour of the symbols represents the transition between phases: red, initial acceleration phase; green, slumping phase; blue, inertial phase and black, viscous phase.

In the slumping phase, the head of the gravity current accumulates the majority of the heavy fluid, resulting in intense mixing. This is due to the elevated head being more energetic than the thinner trailing body. As a result, the entrainment of ambient fluid predominantly takes place in the head of the gravity current (Beghin, Hopfinger & Britter 1981; Ross, Linden & Dalziel 2002; Ross, Dalziel & Linden 2006). Dai & Huang (2022) investigates the merging and splitting processes in the lobe-and-cleft structure at a gravity current head in the slumping phase. They found that the merging process involves the interaction of three hairpin vortices, where the middle hairpin vortex breaks up and reconnects with the two neighbouring hairpin vortices. During this process, a cleft can merge with neighbouring clefts but may never fully disappear. However, in the splitting process, a new streamwise vortex is formed by the parent vortex of opposite orientation. This parent vortex can be either the left part or the right part of the existing hairpin vortex within the splitting lobe. However, as the stratification strength increases, the density difference between the heavy fluid and the ambient fluid becomes smaller, resulting in less mixing within the current head and a decrease in the number of lobes during the slumping phase. When the current transitions into the inertial phase, the lobes and clefts evolve very dynamically, and the splitting and merging processes occur more frequently.

The average size of the lobes can be determined by calculating the mean wavelength of the lobe $\tilde {\lambda }_l=2{\rm \pi} r_f/N_l$ (where $r_f$ is the radial location of the front and $N_l$ is the total number around the radial current). Figure 15 presents the mean wavelength of the lobe as a function of the radial front of the current for both Reynolds numbers. From figure 15(b), it can be observed that the mean wavelength of the lobe decreases as the stratification strength increases at later times where the average size of the lobe decreases. For low-$Re$ cases, a similar trend is observed for $S=0$, 0.2 and 0.5. For $S=0.8$, the lobe size is larger than in the other cases. It is because the gravity current still transitions from the slumping phase to the inertial phase. The gravity current has dissipated and there is no gravity current flow any more before it completely transitions into the inertial phase. The splitting and merging process occurs more frequently in the self-similar phase, and therefore the lobe size for $S=0,0.2$ and 0.5 is smaller than for $S=0.8$.

Figure 15. Plot of the mean wavelength of lobes ($\tilde {\lambda }_l$) against radial front for stratified cylindrical release gravity current with different stratification strength at $Re=$ (a) 3450 and (b) 10 000. The stratification is represented with a different symbol and a cubic spline curve is fitted through the data, $\circ$(—), $S= 0$; $*$(- - -), $S= 0.2$; $\square$(-$\cdot$-$\cdot$), $S= 0.5$; $\triangle$($\cdots \cdots$), $S= 0.8$. The colour of the symbols represents the transition between phases: red, initial acceleration phase; green, slumping phase; blue, inertial phase and black, viscous phase.

Figure 16 shows the dimensionless lobe size, $\tilde {\lambda }_l/\bar {h}_H$ against the front Reynolds number. Discrepancies were observed for both $Re$ cases and this may be due to the model developed by Simpson (1972) being based on planar current, but the current simulation results are for a cylindrical release. The result in the slumping phase for $S=0$ in both $Re$ is in reasonable agreement with Simpson (1972). However, our simulation results show the wavelength of lobes is larger in both the inertial and viscous phases. A similar observation is reported by Cantero et al. (2007a). Despite differences observed between the present study and that of Simpson (1972), our unstratified ($S=0$, see figure 16a,b) simulation results still agree well with those of Cantero et al. (2007a) in both $Re$ cases. This discrepancy is due to the mixing in cylindrical current being greater than a planar current, and the cylindrical current decays more rapidly as the current propagates and mixes radially.

Figure 16. Ratio of the dimensionless mean lobe size, $\tilde {\lambda }_l/\bar {h}_H$ against the front Reynolds number, $Re_F=Re \bar {h}_H u_f$. The simulation results at $Re=$ (a,c) 3450 and (b,d) 10 000 with $S=$ (a,b) 0 and (c,d) 0.8 are plotted at different phases. The solid black line represents the empirical prediction by Simpson (1972): $\tilde {\lambda }_l/\bar {h}_H=7.4Re_F^{-0.39\pm 0.02}$. The figure includes experimental data by Simpson (1972) ($\times$), Dai & Huang (2022) ($+$) and Cantero et al. (2007a) (open symbols for $Re=3450$ and solid symbols for $Re=8950$). The coloured symbol ($\triangle$ and $\triangleleft$ for $S=0$ and 0.8 at $Re=3450$; $\triangledown$ and $\triangleright$ for $S=0$ and 0.8 at $Re=10\,000$): green, slumping phase; blue, inertial phase and cyan, viscous phase.

For $S=0.8$ in figure 16(c,d), the wavelength of lobes in the slumping phase is observed to be similar to the wavelength of $S=0$ in the inertial phase. This is because the transition of the current for $S=0.8$ occurs later than in unstratified cases. A similar observation is found for the case with a higher Reynolds number.

5. Conclusion

Three-dimensional simulations of cylindrical release gravity current are performed with stratification strength varying from 0 to 0.8 at two different $Re$ of 3450 and 10 000. The main objectives of the study are to investigate: (1) the effects of the stratification strength on the cylindrical release gravity current in terms of the propagation of the front and (2) the scaling in the slumping phase for the stratified current. To the best of our knowledge, direct numerical simulations of fully 3-D, full-depth lock release of cylindrical gravity currents in a stratified ambient at a moderate Reynolds number have not been conducted before.

The front location and front velocity of the gravity current in both unstratified and stratified ambient were measured. It was observed that in all cases, the front location of the current decreases with increasing stratification strength, which is consistent with the findings reported in previous studies (Lam et al. 2018a,b; Dai et al. 2021). Similarly, the peak front velocity during the acceleration phase and the constant velocity during the slumping phase decrease with increasing stratification strength and increase with $Re$.

The transition of the current slows by the stratification strength and the gravity current with $S=0.8$ has a longer period of the slumping phase compared with that for the other cases. In the slumping phase, the front velocity of the present simulations is lower than the theoretical prediction of $1/2$ (Benjamin 1968; Shin et al. 2004) for both low- and high-$Re$ cases with $S=0$. The mean front velocity of the unstratified case during the slumping phase is approximately 0.38 for the low-$Re$ case and 0.41 for the high-$Re$ case. These values are in good agreement with the best-fit value of approximately 0.42 reported by Cantero et al. (2007b) for the cylindrical current at large $Re$. However, when the stratification is introduced to the ambient fluid, the models proposed by Benjamin (1968), Shin et al. (2004), Huppert & Simpson (1980) are not suitable as they do not consider the effects of stratification. Consequently, the constant front velocity in the slumping phase of the stratified cases is inconsistent with the constant front velocity predicted by the models. Instead, the scaling proposed by Ungarish & Huppert (2002) and Ungarish (2006), modified with the empirical fit to the front velocity in the slumping phase by Cantero et al. (2007b), effectively captures the effect of stratification, where the front velocity decreases as $S$ increases. Remarkably, the present simulations agree well with the modified scaling, with only a $7\,\%$ difference observed for the low-Reynolds-number cases.

Interestingly, for the high-Reynolds-number case, the difference decreases to $3\,\%$. These results suggest that the modified scaling is valid across a range of Reynolds numbers and can be used to predict the front velocity in the slumping phase for gravity currents propagating in a stratified ambient. The current transitions into the inertial phase after the end of the slumping phase which is observed in cases with weak stratification strength ($S\leq 0.5$). The front velocity in these cases shows very good agreement with the empirical model proposed by Huppert & Simpson (1980) in the inertial phase. In the case with $S=0.8$, the propagation of the current in the inertial phase becomes negligible due to an insufficient density difference between the current and the ambient at the bottom wall. For the cases with $S=0$ and 0.2, a transition to the viscous phase is observed. However, both the viscous scaling proposed by Hoult (1972) and Huppert (1982) underpredict the simulation results even when considering the virtual time origin.

Subcritical flow is obtained for the cases with $S=0.8$ where the propagation of the linear internal gravity wave ($N^*H^*/{\rm \pi}$) is faster than that of the gravity current. In the subcritical regime, the gravity current is smooth and contains minimal vortices. This behaviour is clearly in the contour of the equivalent height in figure 5(d). Gravity currents in the supercritical regime were observed for the remaining cases where the internal gravity wave is locked together with the gravity current head. In the supercritical regime, the gravity current becomes turbulent and strong K–H billows are observed behind the gravity current head. As the gravity current propagates, it gradually slows down due to the entrainment of the linearly stratified ambient fluid. Eventually, the Froude number of the current reaches a critical value ($Fr=1/{\rm \pi}$), triggering a transition to a subcritical flow regime. During this transition, internal gravity waves start to separate from the current head and overtake the current, further decelerating it.

The evolution and breakdown of the turbulent structures are analysed by comparing the isosurfaces of $\lambda _{2}$ at high $Re$ are compared when varying $S$. Strong interaction between the inner vortex rings is observed in the slumping phase. Hairpin vortices, aligned radially, can be observed on the gravity current head during the early stage of the slumping phase for the weakly stratified case ($S<0.5$), as shown by the isosurfaces of density. In contrast, for the strongly stratified case ($S=0.8$), these hairpin vortices do not form until the middle of the slumping phase. The gravity current with $S=0.8$ is less turbulent and the interaction between the inner vortex rings is less prominent. During the transition into the inertial phase, the vortex rings are broken down into smaller rotating vortical structures. The lobes and clefts on the gravity current head are clearly visible on the isosurface of density for $S=0.2$ and 0.5.

The lobes and clefts of the advancing front are compared at both low and high $Re$ with $S=0$. At high Reynolds numbers, the splitting and merging processes are more pronounced, resulting in an asymmetrical appearance of the radial front of the gravity current when visualised by the density contour. In contrast, the radial front of the low-Reynolds-number case appears to be more symmetrical. The lobes and clefts of the advancing front at high $Re$ are also compared when varying $S$. The present simulations reveal that the splitting and merging processes of the gravity current are more significant for weak stratification ($S=0.2$), but less significant for moderate to strong stratification ($S=0.5$ and $S=0.8$). Additionally, the radial front of the current appears more asymmetrical at higher Reynolds numbers but becomes more symmetrical and smaller with increasing $S$. These findings provide valuable insights into the dynamics of gravity currents in stratified environments. The number of lobes in the high-$Re$ cases is significantly greater than in the low-$Re$ cases, but does not vary significantly over time. As the stratification strength increases, the mean wavelength of the lobes decreases. This is attributed to the occurrence of more splitting and merging processes during the self-similar inertial and viscous phases. However, it should be noted that the current slows down at an earlier time as the stratification strength increases. The mean lobe size for both low- and high-$Re$ cases with $S=0$ and 0.8 is consistent with the findings reported by Cantero et al. (2007a).

In this study, the front location and front velocity of the gravity current were examined with respect to their dependencies on $S$ and $Re$. It was observed that increasing $Re$ leads to higher front velocities of the gravity current. Conversely, increasing $S$ results in lower front velocities of the gravity current. These findings have practical implications for predicting the propagation of gravity currents in natural phenomena such as bushfires, where the hot smoke from the fire spreads into a stratified atmosphere. The enhanced stratification slows down the transition of the current and alters the evolution of the gravity current compared with a non-stratified ambient. While previous studies by Dai & Huang (2022) and Cantero et al. (2007a) have reported the splitting and merging processes of 3-D planar and cylindrical currents in a non-stratified ambient, there is currently no literature documenting the splitting and merging processes of lobe-and-cleft structures of cylindrical currents in a stratified ambient. Future extensions of this study could explore different lock aspect ratios or investigate particle-laden flows to further understand the dynamics of cylindrical currents in a stratified ambient.