Subscribe Now Subscribe Today
Research Article

URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

A. Eghbalzadeh, M.M. Namin, S.A.A. Salehi, B. Firoozabadi and M. Javan
Facebook Twitter Digg Reddit Linkedin StumbleUpon E-mail

This study seeks to explore the ability of unsteady Reynolds-averaged Navier-Stokes (URANS) simulation approach for resolving two-dimensional (2D) gravity currents on fine computational meshes. A 2D URANS equations closed by a buoyancy-modified k-ε turbulence model are integrated in time with a second-order fractional step approach coupled with a direct implicit method and discretized in space on a staggered mesh using a second-order accurate finite volume approach incorporating a high resolution semi-Lagrangian technique for the convective terms. A series of 2D simulations are carried out for gravity currents from both discontinuous and continuous sources. Comparisons with experimental measurements show that 2D URANS simulations in conjunction with sufficiently high grid resolution can capture overall features of the gravity currents including the evolution of the Kelvin-Helmholtz interfacial vortices and the propagation of the energy conserving current head with reasonable accuracy.

Related Articles in ASCI
Similar Articles in this Journal
Search in Google Scholar
View Citation
Report Citation

  How to cite this article:

A. Eghbalzadeh, M.M. Namin, S.A.A. Salehi, B. Firoozabadi and M. Javan, 2008. URANS Simulation of 2D Continuous and Discontinuous Gravity Currents. Journal of Applied Sciences, 8: 2801-2813.

DOI: 10.3923/jas.2008.2801.2813



A gravity current is the flow of one fluid within another driven by the density difference between the fluids. The difference in specific weight that provides the driving force may be due to dissolved or suspended material or to temperature differences. Gravity currents could be continuous or discontinuous, depending on whether or not there is constant supply of heavier material from the source. A comprehensive description of gravity currents and their numerous applications is given in Simpson (1997).

Several laboratory experiments have been performed to study gravity currents (e.g., Huppert and Simpson, 1980; Rottman and Simpson, 1983; Parker et al., 1987; García, 1993; Hallworth et al., 1996; Parsons and García, 1989; Lowe et al., 2002; Ross et al., 2002; Marino et al., 2005). It is difficult to make quantitative measurements during the early stage of the current development. Most experimental measurements are therefore, collected when a current has reached a quasi-equilibrium state.

Numerical modeling has the advantage of capturing the evolution history of a developing current. A number of two-dimensional (2D) and/or three-dimensional (3D) direct numerical simulations (DNS) of gravity currents have been reported in the literature (Härtel et al., 2000; Necker et al., 2002, 2005; Birman et al., 2005; Cantero et al., 2006), which have highlighted the role of large-scale, interfacial instabilities in transporting, entraining and mixing the heavier fluid with the ambient flow. Comparisons of the relative performance of 2D vs. 3D DNS results have led to the following conclusions. 2D DNS of lock-release gravity currents can accurately capture some important features of gravity current flows, including the propagation and the evolution of the current head and the dynamics of the interfacial, Kelvin-Helmholtz vortical structures in the slumping phase. 2D simulations, however, are innately unable to reproduce important three-dimensional instabilities of such flows, including the formation of the lobe-and-cleft structure at the current head and the breakdown of interfacial billows to small structures. DNS is a powerful turbulence simulation tool but due to the strong dependence of the required spatial and temporal resolution to the Reynolds number is inherently restricted to low Reynolds numbers and is not feasible for field-scale simulations. Patterson et al. (2005, 2006) carried out 2D and 3D simulations of 2D and axisymmetric gravity current using an alternative LES approach. Large Eddy Simulation (LES) does relax somewhat the required numerical resolution and permits simulations at somewhat higher Reynolds numbers than possible with DNS, but fully near-wall resolved LES of gravity currents at practical Reynolds numbers is not within the reach of fairly affordable computational resources. For example, in the recent LES study by Ooi et al. (2006) computational grids with approximately 35 million grid nodes were employed to resolve the dynamics of a discontinuous lock-release gravity current in a rectangular domain.

To the best of our knowledge and based on the above review of previous literature, computational approaches that have shown promise as practical engineering simulation tools capable of predicting the rich dynamics of real-life gravity current flows-i.e. in complex geophysical geometries and at full-scale Reynolds numbers-do not exist today. Unsteady Reynolds-Averaged Navier-Stokes (URANS) approach with the statistical turbulence model provides a viable alternative for engineering computations of gravity currents at Practical Reynolds numbers and several such models have been proposed in the literature. The k-ε two-equation model has been widely adopted by researchers to simulate 2D conservative gravity currents (Farrel and Stefan, 1988; Eidsvik and Brors, 1989; Sha et al., 1991; Hurzeler et al., 1995; Bournet et al., 1999; Choi and García, 2002; Imran et al., 2004; Kassem and Imran, 2004; Huang et al., 2005). These studies have yielded results that are in good agreement with the measurements of basic flow properties of the gravity currents (e.g., spreading rate and propagation speed of the front) but tend to produce rather diffused flow patterns at the interface without capturing the interfacial vortex dynamics revealed by experiments and higher level numerical simulations such as LES and DNS.

In this study we seek to demonstrate that URANS model employing a buoyancy-extended statistical turbulence model can indeed capture the rich dynamics of large scale interfacial vortical structures in 2D gravity currents provided that the numerical resolution-both in the sense of the accuracy of the numerical method and the temporal and spatial resolution of the simulation-is sufficiently fine. We develop and apply a robust and second-order-accurate 2D numerical model to reproduce the propagation of gravity currents accompanying the evolution of the interfacial billows. The model employs the Boussinesq form of the URANS equations in conjunction with a buoyancy-extended k-ε closure for the turbulence (Rodi, 1993). We carry out numerical simulations of various continuous and discontinuous gravity currents and compare the results with the experimental visualizations and the measurements published in Huppert and Simpson (1980), Patterson et al. (2005) and Garcia (1993). The grid sensitivity of the numerical solutions is also analyzed for all cases.

Governing Equations: In this study we focus on flows with relatively small density differences for which the usual Boussinesq approximation can be assumed to be valid-i.e., all variations in density can be neglected except for the buoyancy term. Incorporating the Boussinesq hypothesis to relate the Reynolds stresses with the mean rate of strain tensor via an eddy viscosity, 2D URANS equations for incompressible, stratified flow read as:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

where, u and w are the Reynolds-averaged velocity components in x and z directions, respectively; ρO is the density of the ambient fluid; g is the acceleration due to gravity; v and vt are the molecular and eddy viscosity, respectively; ρ is the local mixture density and t is time. P* = P-ρogh is the modified pressure including the gravity terms where P* is total pressure and h is distance from the reference in z direction.

The eddy viscosity is modeled by the buoyancy-modified k-ε turbulence model (Rodi, 1993) as:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

where, k is the turbulent kinetic energy, ε is its dissipation rate and Cμ is a model constant. Quantities of k and ε are obtained from the solution of the transport equations:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

In the Eq. 5-6, D/Dt denotes the material derivative, Pr and B are production terms of turbulent kinetic energy due to the mean velocity gradient and the buoyancy and are estimated, respectively, as:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

In Eq. 7-8σk, σε, and σt = turbulent Prandtl number for k, ε and the source of gravity current (in this study, salinity S), respectively. Cε1, Cε2 and Cε3 are turbulence model constants. The following standard values are specified for the turbulence model constants:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

The turbulent Prandtl number σt, the ratio of the turbulent viscosity and the turbulent diffusivity is set equal to 0.85 which is acceptable in general buoyant shear flow (Rodi, 1987; Farrel and Stefan, 1988).

The empirical coefficient Cε3 is directly related to the buoyancy effect and its optimal value does not have a firm basis, unlike other model coefficients. Rodi (1993) suggested that the coefficient should vary between 0 for buoyant horizontal shear flows and 1 for buoyant vertical shear flow. That is, the optimal value of Cε3 is dependent upon the flow conditions and its choice is controversial. Hossain and Rodi (1982), Rodi (1987), Fukushima and Hayakawa (1990), Fukushima and Watanabe (1990) and Choi and García (2002) suggested Cε3 values ranging from 0 to 0.4 to give a good agreement with experimental results on gravity currents. Since the gravity currents considered in this study are close to the horizontal shear flow, except during the early stage, the choice of a small constant value of Cε3 close to 0 is reasonable based on the recommendation of Rodi (1993). In this study we use a constant value of Cε3 = 0.2.

The URANS and turbulence closure equations are solved simultaneously with the following transport equations for the salinity S, which is used to determine scalar transport and the variation of the fluid density:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Density is assumed to be linearly related to the mean volumetric concentration through the equation of state as ρ = ρo (1.0+mS) where m is a constant.


A fractional step (Yanenko, 1971) method is employed to integrate the governing equations in time coupled with a projection method for satisfying the continuity equation (Namin, 2003). Non-uniform bathymetry of the computational domains is handled by the sigma-type body-fitted grids which fit the vertical direction of the physical domain. The finite-volume integration method is used for different terms in the governing equations. To facilitate the presentation of the algorithm, we present an overview of the entire temporal integration scheme and then the spatial discretization of the various terms and the iterative solution of the resulting discrete equations.

Temporal integration scheme: The fractional step algorithm: The governing equations are discretized in time by fractional-step, or operator splitting, method in which the time advancement is decomposed into a sequence of steps: the advection, the diffusion and the projection. This method allow us deploy the most efficient numerical technique for each of the individual process. Dimensional splitting is also used to reduce the multi-dimensional problem to a series of one-dimensional problems. The various steps are summarized below where the following notation is used for the sake of brevity: Ui = u, w for i = 1 and 2, respectively.

Advection step (i = 1,2)

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Diffusion step (i = 1,2)

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

where, a and b are constant coefficients. If i = 1 then a = 2 , b = 1 and if i = 2 then a =1 , b = 2.

Continuity (or projection) step:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Equation 16 and Eq. 17 are combined by applying the divergence operator to Eq. 17 and enforcing Eq.16 to derive the following Poisson equation for the pressure field:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents


Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

Following the solution of the pressure equation, the velocity field at the new time step is calculated as follows:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

These fractional step and dimensional splitting approaches save considerable amounts of computational time and preserve overall accuracy of the temporal and spatial discretization schemes (Namin, 2003).

The fractional step approach outlined above for solving the momentum equations is also employed to integrate the transport equations for the turbulence quantities and the salinity. For the later equation only the advection and diffusion steps are required. For the former equations, however, the diffusion step is followed by a step in which the source terms are taken into account explicitly.

Spatial discretization scheme: The discretization of the advection terms in the above iterative scheme is critical for spatial accuracy of the numerical method. We employ the semi-Lagrangian method (Namin, 2003) based on the Fromm scheme (1968) for discretization of the advection terms of the one-dimensional form isolated by the dimensional splitting approach.

The ULTIMATE conservative difference scheme of Leonard (1991) is employed to ensure to maintain monotonicity while preserving spatial accuracy. Using a monotone scheme is essential for numerical stability and the smoothness of the computed solutions in the discontinuities where sharp scalar gradients are present and dominate the dynamics of the flow.

In the second step, the diffusion terms are solved by the semi-implicit second-order Crank-Nicolson scheme to obtain the second intermediate velocity components Ui** as in Eq. 13-15.

Solution of the discrete equations: During the advection step of the fractional step procedure, the discrete equations are integrated in time explicitly. In the diffusion step the Thomas tri-diagonal matrix algorithm is used to solve the discrete equations implicitly. The pressure-Poisson equation is solved using a block tri-diagonal algorithm, which can be formulated by writing the discrete Poisson equation for all the cells.

After calculating the velocity components and pressure field in new time step, the transport equations for the turbulence kinetic energy, energy dissipation rate and the salinity are solved using the same fractional step algorithm used to integrate the momentum equations. Finally the density is calculated using state equation.

Boundary conditions: The boundaries of the computational domain are inlet, outlet, free-surface and solid walls. In the absence of wind shear the net fluxes of horizontal momentum and turbulent kinetic energy are set equal to zero at the free surface. Flat, slip boundary condition of zero velocity gradient normal to the surface is applied and the pressure is set equal to the atmospheric value at the surface. The dissipation rate ε is calculated from the relation given by Rodi (1993):

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

where, D is the local water depth and zs is the coordinate of free surface.

The wall function approach is used to specify boundary conditions at the bottom of the channel in order to avoid the resolution of viscous sublayer. The first grid point off the wall (center of the control volume adjacent to the wall) is placed inside the logarithmic layer based on the buoyancy velocity Ub. In the wall-function approach, the resultant wall shear stress vector τw is related to the velocity vector at the first grid point V2 by the standard, smooth-wall, log law (Wu et al., 2000):

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

where, the subscript 2 refers to the point at the first control volume center at the wall, Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents is the normal distance to the wall; Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currentsis the bed shear velocity; κ is von Karman constant (= 0.41) and E is an integration constant that for smooth beds is equal to 8.43. Here, we employ an algorithm for computing u* in Eq. 26. The near-wall values of turbulent kinetic energy k and dissipation rate ε are given as:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

On the other hand, if the first grid node from the wall locates within the viscous layer Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currentsdue to the evolution of complex vortical structures, these near-wall values are approximated as:

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

In case of the continuous gravity current, known quantities are specified at the inlet for the inflow velocity, the concentration of the gravity current source and the current thickness while turbulent kinetic energy and dissipation rate at the inlet are estimated as follows (Fukushima and Watanabe (1990)):

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

where, uin, kin and hin are the averaged velocity, the turbulent kinetic energy and the thickness of the gravity current at the inlet. At the outlet the normal gradients of all dependent variables are set equal to zero.

In case of lock-release gravity currents, the dense and ambient fluids are completely divided at x = xo ( Fig. 1) at t = 0. With the start of the simulation, a dense fluid propagates rightward along the bottom boundary while the ambient fluid moves leftward along the surface boundary naturally due to the density difference.


Here, we report numerical results for 2D gravity currents from continuous and discontinuous sources along with comparisons with experimental measurements available in the literature. After a brief description of computational details for each case, we first present the comparison of the 2D numerical solutions and available experimental measurements. Subsequently, we discuss the details of the simulated flows and highlight some issues pertaining to the numerical simulation of gravity currents with the 2D URANS model. We have extensively tested the grid sensitivity of numerical solutions for discontinuous source case considered in this study. Numerical solutions computed on a coarse, a fine and three intermediate meshes are presented only for the first test case of discontinuous source case to elucidate the sensitivity of the numerical solution to grid refinement while for other test case a coarse mesh and two essentially grid-converged solutions are included.

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents
Fig. 1: Formation of the lock release gravity current at the initial state and related definitions

For continuous source case, only a coarse mesh and a grid-converged - i.e., the effect of further grid refinement on the solution is not significant - fine mesh solutions are included in this paper due to the space limitation.

To facilitate the discussion of the results we adopt a convention similar to that introduced by Lowe et al. (2002) regarding the various regimes of the flow. Namely, we divide the flow in three regimes: 1) the energy conserving head region in which the velocity is roughly uniform with speed same to the head speed; 2) the wake region where the large-scale interfacial vortices dominate the flow and 3) the tail region where interfacial vortices break down into smaller structures and/or blobs of dense fluid are entrained by vortical structures into the lighter fluid. These definitions are slightly different from those in Lowe et al. (2002) but are reasonable to explain the gravity currents from both continuous and discontinuous sources.

Discontinuous gravity currents: A common procedure for generating gravity current in the laboratory is the sudden removal of a vertical gate separating two bodies of dense fluid and less dense ambient fluid on a flat bed in uniform environment, as shown in Fig. 1. The resulting current is referred to as a lock release gravity current whose underlying large-scale dynamics is two-dimensional. We simulate two different experimental configurations investigated by Huppert and Simpson (1980) and Patterson et al. (2005), respectively. For both cases the evolution of the discontinuous gravity current takes place on a flat, horizontal bed.

In the experiments of Huppert and Simpson (1980), a plexiglass channel was filled with fresh water of the depth h0 and divided into two regions by a wooden gate located at a distance x0 from the left end ( Fig. 1). Common salt (NaCl) was added to the water behind the gate until a given uniform density ρs was obtained. The gate was then removed suddenly. Numerical simulations are carried out for two experimental configurations: 1) Case DH1 where h0 = 14.9 cm, x0 = 39 cm and ρs = 1009.3 kg m3 and 2) Case DH2 where h0 = 10 cm, x0 = 30 cm and ρs = 1011.4 kg m3.

Table 1: Configurations for numerical simulations of the discontinuous gravity currents on flat bed experimentally investigated by Huppert and Simpson (1980) and Patterson et al. (2005).
Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents

The gravity current scenario studied experimentally and numerically by Patterson et al. (2005) involves a In the case names, D: Discontinuous gravity current; H: Huppert and Simpson (1980); P: Patterson et al. (2005); 1,2 = Serial number of test cases; cg, mg and fg: Coarse, Medium- and Fine-grid solutions configuration similar to that in Huppert and Simpson (1980) but with the following parameters: h0 = 28.8 cm, x0 = 120 cm and the reduced gravity g` = 4.7 cm s2 (g` = g (ρs0)/ρ0). We refer to this case as Case DP. Computational details for the numerical simulation of these cases are summarized in Table 1 which includes total number of grid nodes (NIxNK), uniform grid spacing Δx and Δz in the horizontal and vertical directions, computational time step Δt and total simulation time T.

Since wall functions are used in all computations reported in this section and in order to facilitate meaningful comparisons between solutions computed on successively refined meshes, the first grid node off the solid wall is placed for all grids at the same distance from the wall. By doing so we ensure that for a given test case wall functions are always applied at the same location within the logarithmic layer and the grid in the wall normal direction is refined only above that location.

To evaluate our numerical solution, the location of the gravity current head computed on the relatively coarse mesh for DH1_cg (Table 1) is compared with the measurement of Huppert and Simpson (1980) in Fig. 2. The result confirms that the 2D simulation can reproduce the head propagation of the gravity current from a fixed volume on a flat bed with good accuracy. These results suggest that the head propagation speed of at least this gravity current essentially does not depend on 3D effects. This conclusion is consistent with and further supported by the previous study of Huppert and Simpson (1980). They showed that the front speed of gravity currents can be predicted with reasonable accuracy using a simple empirical model that expresses the current head propagation speed as a function of the initial amount of dense fluid, the reduced gravity, the channel depth and time (Huppert and Simpson, 1980).

Among a series of experiments reported by Huppert and Simpson (1980), we selected case DH2, for which flow visualizations were reported, to test the ability of our numerical model to capture complex vortical structures and evaluate the sensitivity of the computed solution to mesh refinement.

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents
Fig. 2:
Comparison of the front position of the computed gravity current (solid line) with the experimental measurements (symbols) of Huppert and Simpson (1980)

Numerical solutions obtained on a coarse, a fine and three intermediate meshes, denoted by DH2_cg, DH2_fg22, DH_fg11, DH_fg12 and DH_fg21 (Table 1), are visualized by the density distribution and compared with the experimental visualizations of Huppert and Simpson (1980) in Fig. 3. This Fig. 3 shows snapshots at 4.4, 6.8, 9.7 and 17.5 seconds after the release of dense fluid.

The simulated gravity current is fundamentally in the so-called slumping phase during which the current head propagates at an approximately constant velocity after the current is established in a relatively short period. Figure 3a-b shows initial transient flows after the release of dense fluid when the initial discontinuity disappears quickly and the height of the gravity current becomes uniform. Figure 3c-d, on the other hand, depicts the propagation of the gravity current head of fairly constant height in the slumping regime. Figure 3 leads to a number of important conclusions: (1) with only exception the simulation on the coarsest mesh (DH2_cg), which yields diffused solutions, all other grids capture the onset of Kelvin-Helmholtz billows at the interface; (2) the effect of grid refinement is seen to be far more pronounced in the wake and tail regions than in the head region and during the first two instants; (3) the numerical solutions are more sensitive to the grid refinement in the horizontal direction rather than in the vertical direction.

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents
Fig. 3:
Experimental visualizations (upper) of Huppert and Simpson (1980) and numerical results (lower) visualized by the density showing the time evolution of the gravity current: (a) t = 4.4 sec, (b) t = 6.8 sec, (c) t = 9.7 sec and (d) t = 17.5 sec after removing the gate

As shown in Fig. 3, both DH2_fg11 and DH2_fg12 yield results that are very similar to each other while DH2_fg21 solutions are similar to DH2_fg22 solutions; (4) at the two later instants, the solutions on the finest meshes (fg11, fg12, fg21 and fg22) do not show significant differences from each other, even though the finer meshes do capture somewhat finer scales of motion in the tail regions. These results suggest that the fg22 mesh is adequate for obtaining grid insensitive solutions.

The overall dynamics of the large-scale interfacial vortical structures reproduced by the fine mesh simulations over the gravity current head and in the turbulent wake region and the resulting longitudinal variation of the irregular and unsteady profile of the gravity current heights are in reasonable agreement with the experimental visualizations of Huppert and Simpson (1980). It is observed, however, that the breakdown of the billows reproduced in the rear half of the wake region and in the tail region is too weak to cause intense turbulent mixing of dense and ambient fluids in this region. The billow breakdown is a phenomenon governed by 3D mechanisms and as such it can not be captured by 2D numerical simulations regardless of grid refinement. This is an inherent feature of gravity current flows, which has already been established in the literature in previous 2D DNS studies (Härtel et al., 2000; Necker et al., 2002). The novel finding of present research, however, is that URANS simulations on sufficiently refined grids yield solutions that are comparable with those obtained with 2D DNS both insofar as the onset of the interfacial billows is concerned and the inability of the 2D assumption to capture the rapid breakdown of these intense large-scale structures.

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents
Fig. 4: Variation of the propagation velocity of the computed gravity current nose

The temporal variation of the propagation velocity of the gravity current nose for case DH2_fg22 (Table 1) is shown in Fig. 4 to shows the ability of the URANS model to capture the different regimes of current propagation. As shown in this figure, the head propagation velocity accelerates rapidly to its maximum value during the first 1.5 sec. This is followed by an initially rapid and subsequently more gradual deceleration phase up to approximately t = 4 sec in the initial slumping regime. After the brief acceleration and deceleration periods, the propagation velocity of the computed gravity current head remains nearly constant in the slumping phase. The overall behavior of the current head reproduced by present 2D URANS is consistent with the results of 2D DNS carried out by Birman et al. (2005) and Blanchette et al. (2005) for different lock-release gravity currents in the slumping phase.

Numerical simulations of a full-depth gravity current, which was investigated both experimentally and numerically by Patterson et al. (2005), are carried out and the results obtained on a coarse and two fine meshes (DP_cg, DP_mg and DP_fg, respectively) are compared with the experimental visualizations of Patterson et al. (2005) in Fig. 5. This Fig. 5 shows that, even though smaller scales of the flow in the wake region are captured only on the finer meshes, the overall sensitivity of the solution to mesh refinement is not significant. All numerical solutions reproduce well the variation of the longitudinal profile of the gravity current height and the propagation of the current head without any significant retardation. These results are in better overall agreement with the measurements than the 2D numerical solutions of Patterson et al. (2005) who employed an Implicit Large-Eddy Simulation (ILES), which appeared to suffer from significant retardations of the current head in the slumping phase ( Fig. 1 in Patterson et al., 2005). The reason for this is not entirely clear and can not be easily assessed due to the lack of computational details in the study of Patterson et al. (2005). However, it is worth to point out that Patterson et al. (2005) did not use a turbulence model, which is customary in the ILES approach and carried out calculations on a single very coarse mesh consisting of 256x32 grid nodes, which is significantly coarser than even our coarsest mesh (DP_cg) in the vertical direction. Overall our 2D URANS simulations capture reasonably well the rich dynamics of Kelvin-Helmholtz vortices and the resulting flow patterns are comparable with the DNS results of a similar full-depth gravity current (Härtel et al., 2000). As in the previous case, however, the complete breakdown of the interfacial billows to a 3D state can not be resolved by our simulations due to the inherent 3D nature of the underlying physical processes.

Continuous gravity current: We report simulation results of a gravity current from a continuous source along a flat bed with an abrupt change in slope. The simulation case is one of a series of experiments carried out by Garcia (1993). The geometry of the experimental channel is same as Fig. 6. The current thickness at the inlet was set equal to 3 cm with the help of a sluice gate and the inflow rate per unit width was set equal to 33 cm2 sec1, which gave a layer-averaged inlet velocity of 11 cm sec1. The resulting gravity current reached an equilibrium state transitioning from the supercritical flow along the slope upstream of the breakpoint (x = 5.0 m, Fig. 6 for the location) to the subcritical flow along most of the downstream flat bed after undergoing an internal hydraulic jump. Numerical simulations are carried out for one of the experimental configurations (Garcia, 1993), which involved a conservative saline gravity flow with excess fractional density equal to 0.013. A 24 m long computational domain is employed in order to avoid reflections from the outlet. With reference to Table 2, a uniform mesh is used for the CG_cg and CG_fg simulations. Computational details for numerical simulations are summarized in Table 2.

The experimental measurements of García (1993) were collected when the current had reached a quasi-equilibrium state. To compare numerical solutions with the measurements, we have continued our unsteady simulations until the solutions in the region of interest (0≤x≤10.0 m) reach a quasi-equilibrium state. The time evolution of the streamwise velocity residual in the region of interest shows that both the CG_cg and CG_fg solutions approach steady state after nearly t = 300 sec.

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents
Fig. 5:
Experimental visualizations (upper picture) and numerical results (lower pictures) visualized by the density showing the time evolution of the gravity current: (a) t = 3.13 sec, (b) t = 5.67 sec, (c) t = 8.14 sec and (d) t = 15.45 sec after removing the gate

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents
Fig. 6: Schematic side view of the 0.3 m wide channel (Garcia, 1993). All dimensions are in meters

Table 2: Configurations for numerical simulations of the continuous gravity current across a break in slope experimentally investigated by Garcia (1993)
Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents
In the case names, C: Continuous gravity current; G: Garcia (1993); cg: Coarse-grid solution and fg: Fine-grid solution

In Fig. 7 and 8 we compare the streamwise velocity profiles at selected locations and iso-contours of density distribution, respectively, computed at t = 300 sec on both meshes with the measurements of García (1993). As shown in Fig. 7, the fine mesh solution (CG_fg) yields more reasonable streamwise velocity profiles except in the vicinity of the jump (4.0 m<6.0 m ). The CG_cg solution underestimates the streamwise velocity especially near the bottom, in the supercritical flow region along the slope (roughly x < 5.0 m). Both numerical solutions (CG_cg and CG_fg) yield streamwise velocity profiles which are significantly more diffused in the vertical direction than the measurements in the vicinity of the breakpoint at x = 5.0 m.

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents
Fig. 7:
Comparison of the vertical profiles of the streamwise velocity measured (symbols) by García (1993) and computed (dashed lines for CG_cg and solid lines for CG_fg ) at selected locations

The CG_cg solution captures reasonably well the velocity profile near the wall (for z < 0.09 m) but partly underestimates the streamwise velocity away from the bottom (for z > 0.09 m), in the subcritical flow region (roughly x > 6.0 m). The variation of the computed velocity profiles along the longitudinal distance is consistent with those observed in the computed density distributions. As shown in Fig. 8, both numerical solutions (CG_cg and CG_fg) reproduce rather diffused distribution of dense fluid without its abrupt longitudinal variation at the internal hydraulic jump region (iso-contour of the density 1002 kg m3 around x = 6.0 m in Fig. 8). This mild variation of dense fluid in the near wall region is responsible for the excessively smooth longitudinal variation of the computed streamwise velocity profiles on the downstream flat bed.

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents
Fig. 8: Comparison of measured (García, 1993) and computed density contours (1,001 and 1,002 kg m1) along the longitudinal direction

The CG_fg solution yields improved prediction of density distribution in the downstream region.

Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents
Fig. 9: Sequential snapshots of density distributions computed for the CG_cg (left) and CG_fg (right) cases at selected time instants

However, the significant discrepancy between the CG_fg predictions and the measurements of high density distribution still remains in the vicinity of the breakpoint (x = 5.0 m). One of possible reasons for this discrepancy could be the 3D nature of the internal hydraulic jump induced in the narrow experimental channel of 30 cm width. For instance, the aspect ratio of the channel width to the gravity current height in the experiments of García (1993) ranges roughly from 1.5 to 3 in the breakpoint region (4.0 m < x = 6.0 m). In this case, the streamwise mean velocity profile is no longer uniform in the lateral direction (Kirkgöz and Ardiçhoğ—lu, 1997). Furthermore, the hydraulic jump is relatively weak (Froude number Image for - URANS Simulation of 2D Continuous and Discontinuous Gravity Currents 1.75), as pointed out by García (1993), which is sufficient to provoke the formation of breaking undular jump (Ohtsu et al., 2003) in the channel flow.

Let us now we turn our attention to the simulated temporal evolution of the continuous gravity current. Figure 9 shows the flow evolution of the gravity current computed on both computational meshes at several instants in time. Density distributions computed at t = 300 sec are also included to depict the predicted gravity current at the quasi-equilibrium state. The CG_fg predicts more energetic vortical structures than CG_cg as the gravity current head travels down the slope and passes the internal hydraulic jump region. Due to the roll-up of large scale vortices over the head, ambient fluid is continuously entrained into the rear part of the head region causing its continuous thinning ( Fig. 9 until t = 60 sec). On the other hand, the current head is seen to grow as it travels downstream as dense fluid is pushed by the vortices from the rear part into the head region. Figure 9 along with our video animations generated using instantaneous solutions further show that the gravity current slightly bounces off the flat bed upon impingement generating a vertically upward flow and a wavy-like structure in the streamwise direction as the gravity current head passes the transient region.


A series of 2D URANS simulations are carried out to resolve discontinuous and continuous 2D gravity currents. Comparison of the numerical solutions with available experimental measurements leads to the conclusion that high-resolution, 2D URANS simulations using a buoyancy-extended k-ε turbulence model can resolve the rich dynamics of the interfacial vortical structures as well as essential 2D flow features of the gravity currents in the slumping phase. The propagation and the development of the energy-conserving gravity current head are reproduced with good accuracy even on relatively coarse meshes. This result is reasonable because the behavior of the gravity current head in uniform environment is independent of three-dimensional effects in the slumping phase. The fine mesh simulations appear to reproduce reasonably well unsteady flow features in the slumping phase including the formation of the Kelvin-Helmholtz vortices at the interface of dense and ambient fluids; the evolution of the boundary layer as the energy conserving current head travels forward along the bed and the instantaneous velocity field with the magnitude larger than that of the current head in the wake region of the discontinuous gravity current.

Notable discrepancies between 2D simulation results and the experimental measurements of gravity currents are observed in the tail region of discontinuous gravity currents where the 2D large-scale interfacial vortices do not break down fully into smaller structures and the intensity of the ensuing mixing is underestimated. For continuous gravity currents, discrepancies were also observed in the region of the internal hydraulic jump where the 2D simulations underestimate the abrupt vertical mixing. These discrepancies are to be expected, however, as they are associated with phenomena that are inherently three-dimensional and can not be captured by a 2D model. The discrepancies with the measurements notwithstanding, however, our 2D simulation results are significant because they underscore the potential of URANS methods in simulations of gravity current flows.

1:  Birman, V.K., J.E. Martin and E. Meiburg, 2005. The non-Boussinesq lock-exchange problem. Part 2. High-resolution simulations. J. Fluid Mech., 537: 125-144.
CrossRef  |  

2:  Blanchette, F., M. Strauss, E. Meiburg, B. Kneller and M.E. Glinsky, 2005. High-resolution numerical simulations of resuspending gravity currents: Conditions for self-sustainment. J. Geophys. Res., 110: C12022-C12036.
CrossRef  |  

3:  Bournet, P.E., D. Dartus, B. Tassin and B. Vincon-Leite, 1999. Numerical investigation of plunging density current. J. Hydraul. Eng., 125: 584-594.
CrossRef  |  

4:  Cantero, M.I., S. Balachandar, M.H. Garcia and J.P. Ferry, 2006. Direct numerical simulations of planar and cylindrical density currents. J. Applied Mech., 73: 923-930.
Direct Link  |  

5:  Choi, S. and M.H. Garcia, 2002. k-ε turbulence modeling of density currents developing two dimensionally on a slope. J. Hydraul. Eng., 128: 55-63.
CrossRef  |  

6:  Eidsvik, K.J. and B. Brors, 1989. Self-accelerated turbidity current prediction based upon k-ε model turbulence. Cont. Shelf Res., 9: 617-627.
CrossRef  |  

7:  Farrel, G.J. and H.G. Stefan, 1988. Mathematical modeling of plunging reservoir flows. J. Hydraul. Res., 26: 525-537.
CrossRef  |  

8:  Fromm, J.E., 1968. A method for reducing dispersion in convective difference schemes. J. Comput. Phys., 3: 176-189.
CrossRef  |  

9:  Fukushima, Y. and N. Hayakawa, 1990. Analysis of inclined wall plume by the k-ε turbulence model. J. Applied Mech., 57: 455-465.
CrossRef  |  

10:  Fukushima, Y. and M. Watanabe, 1990. Numerical simulation of density underflow by the k-ε turbulence model. J. Hydrosci. Hydr. Eng., 8: 31-40.

11:  Garcia, M.H., 1993. Hydraulic jumps in sediment-driven bottom currents. J. Hydraul. Eng., 119: 1094-1117.
CrossRef  |  

12:  Hallworth, M.A., H.E. Huppert, H.E. Phillips and R.S.J. Sparks, 1996. Entrainment into tow-dimensional and axisymmetric turbulent gravity currents. J. Fluid Mech., 308: 289-311.
CrossRef  |  

13:  Härtel, C., E. Meiburg and F. Necker, 2000. Analysis and direct numerical simulation of the flow at a gravity-current head. Part 1. Flow topology and front speed for slip and no-slip boundaries. J. Fluid Mech., 418: 189-212.
Direct Link  |  

14:  Hossain, M.S. and W. Rodi, 1982. Turbulent Buoyant Jets and Plumes. 1st Edn. Pergamon, Oxford, England, pp: 184.

15:  Huang, H., J. Imran and C. Pirmez, 2005. Numerical model of turbidity currents with a deforming bottom boundary. J. Hydraul. Eng., 131: 283-293.
CrossRef  |  

16:  Huppert, H.E. and J.E. Simpson, 1980. The slumping of gravity currents. J. Fluid Mech., 99: 785-799.
CrossRef  |  

17:  Hurzeler, B.E., G.N. Ivey and J. Imberger, 1995. Spreading model for a turbidity current with reversing buoyancy from a constant-volume release. Mar. Freshwater Res., 46: 393-408.

18:  Imran, J., A. Kassem and S.M. Khan, 2004. Three-dimensional modeling of density current. I. Flow in straight confined and unconfined channels. J. Hydraul. Res., 42: 578-590.
CrossRef  |  

19:  Kassem, A. and J. Imran, 2004. Three-dimensional modeling of density current.II. Flow in sinuous confined and unconfined channels. J. Hydraul. Res., 42: 591-602.

20:  Kirkgoz, M.S. and M. Ardiclioglu, 1997. Velocity profiles of developing and developed open channel flow. J. Hydraul. Eng., 123: 1099-1105.
CrossRef  |  

21:  Lowe, R.J., P.F. Linden and J.W. Rottman, 2002. A laboratory study of the velocity structure in an intrusive gravity current. J. Fluid Mech., 456: 33-48.
CrossRef  |  

22:  Leonard, B.P., 1991. The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection. Comput. Meth. Applied Mech. Eng., 88: 17-74.
CrossRef  |  

23:  Marino, B.M., L.P. Thomas and P.F. Linden, 2005. The front condition for gravity currents. J. Fluid Mech., 536: 49-78.
CrossRef  |  

24:  Namin, M.M., 2003. A fully three-dimensional non-hydrostatic free surface flow model for hydro-environmental predictions. Ph.D. Thesis. Cardiff School of Engineering, Cardiff University, Cardiff, UK.

25:  Necker, F., C. Hartel, L. Kleiser and E. Meiburg, 2002. High-resolution simulations of particle-driven gravity currents. Int. J. Multiphase Flow, 28: 279-300.
CrossRef  |  

26:  Necker, F., C. Hartel, L. Kleiser and E. Meiburg, 2005. Mixing and dissipation in particle-driven gravity currents. J. Fluid Mech., 545: 339-372.
CrossRef  |  

27:  Ohtsu, I., Y. Yasuda and H. Gotoh, 2003. Flow conditions of undular hydraulic jumps in horizontal rectangular channels. J. Hydraul. Eng., 129: 948-955.
CrossRef  |  

28:  Ooi, S.K., S.G. Constantinescu and L.J. Weber, 2006. Numerical simulation of lock-exchange gravity driven flows. IIHR Technical Report No. 450, Univ. of Iowa, USA.

29:  Patterson, M.D., J.E. Simpson, S.B. Dalziel and N. Nikiforakis, 2005. Numerical modeling of two-dimensional and axisymmetric gravity currents. Int. J. Numer. Meth. Fluids., 47: 1221-1227.
CrossRef  |  

30:  Patterson, M.D., J.E. Simpson, S.B. Dalziel and G.J.F. van Heijst, 2006. Vortical motion in the head of an axisymmetric gravity current. Phys. Fluids., 18: 046601-046601.
CrossRef  |  

31:  Parker, G., M. Garcia, Y. Fukushima and W. Yu, 1987. Experiments on turbidity currents over an erodible bed. J. Hydraul. Res., 25: 123-147.
CrossRef  |  Direct Link  |  

32:  Parsons, J.D. and M.H. Garcia, 1989. Similarity of gravity current fronts. Phys. Fluids., 10: 3209-3213.
CrossRef  |  

33:  Rodi, W., 1987. Examples of calculation methods for flow and mixing in stratified fluids. J. Geophys. Res., 92: 5305-5328.
CrossRef  |  

34:  Rodi, W., 1993. Turbulence Models and Their Applications in Hydraulics-a State-of-Arts Review. 3dr Edn. International Association for Hydraulic Research Monograph, Balkema, Rotterdam, The Netherlands.

35:  Ross, A.N., P.F. Linden and S.B. Dalziel, 2002. A study of three-dimensional gravity currents on a uniform slope. J. Fluid Mech., 453: 239-261.
Direct Link  |  

36:  Rottman, J.W. and J.E. Simpson, 1983. Gravity currents produced by instantaneous releases of a heavy fluid in a rectangular channel. J. Fluid Mech., 135: 95-110.
CrossRef  |  

37:  Sha, W., T. Kawamura and H. Ueda, 1991. A numerical study on sea/land breezes as a gravity current: Kelvin-Helmholtz billows and inland penetration of the sea-breeze front. J. Atmospheric Sci., 48: 1649-1665.

38:  Simpson, J.E., 1997. Gravity Currents in the Environment and the Laboratory. 2nd Edn. Cambridge Univ. Press, New York.

39:  Wu, W., W. Rodi and T. Wenka, 2000. 3D numerical modeling of flow and sediment transport in open channels. J. Hydraul. Eng., 126: 4-15.
CrossRef  |  

40:  Yanenko, N.N., 1971. The Method of Fractional Steps. Springer-Verlag, Berlin.

©  2021 Science Alert. All Rights Reserved