INTRODUCTION
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 quasiequilibrium
state.
Numerical modeling has the advantage of capturing the evolution history
of a developing current. A number of twodimensional (2D) and/or threedimensional
(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 largescale, 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 lockrelease 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, KelvinHelmholtz vortical structures
in the slumping phase. 2D simulations, however, are innately unable to
reproduce important threedimensional instabilities of such flows, including
the formation of the lobeandcleft 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 fieldscale 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 nearwall 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 lockrelease
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
reallife gravity current flowsi.e. in complex geophysical geometries
and at fullscale Reynolds numbersdo not exist today. Unsteady ReynoldsAveraged
NavierStokes (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ε twoequation 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 buoyancyextended
statistical turbulence model can indeed capture the rich dynamics of large
scale interfacial vortical structures in 2D gravity currents provided
that the numerical resolutionboth in the sense of the accuracy of the
numerical method and the temporal and spatial resolution of the simulationis
sufficiently fine. We develop and apply a robust and secondorderaccurate
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 buoyancyextended 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 validi.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:
where, u and w are the Reynoldsaveraged 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 v_{t} are the
molecular and eddy viscosity, respectively; ρ is the local mixture
density and t is time. P* = Pρ_{o}gh 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 buoyancymodified kε turbulence
model (Rodi, 1993) as:
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:
In the Eq. 56, D/Dt denotes the
material derivative, P_{r} and B are production terms of turbulent
kinetic energy due to the mean velocity gradient and the buoyancy and
are estimated, respectively, as:
In Eq. 78Ïƒ_{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:
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:
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.
NUMERICAL METHOD
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). Nonuniform bathymetry of the computational
domains is handled by the sigmatype bodyfitted grids which fit the vertical direction of the physical domain. The finitevolume 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 fractionalstep, 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 multidimensional
problem to a series of onedimensional problems. The various steps are
summarized below where the following notation is used for the sake of
brevity: U_{i} = u, w for i = 1 and 2, respectively.
Advection step (i = 1,2)
Diffusion step (i = 1,2)
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:
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:
where
Following the solution of the pressure equation, the velocity field at
the new time step is calculated as follows:
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 semiLagrangian method (Namin, 2003)
based on the Fromm scheme (1968) for discretization of the advection terms
of the onedimensional 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 semiimplicit
secondorder CrankNicolson scheme to obtain the second intermediate velocity
components U_{i}** as in Eq. 1315.
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 tridiagonal matrix
algorithm is used to solve the discrete equations implicitly. The pressurePoisson
equation is solved using a block tridiagonal 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, freesurface 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):
where, D is the local water depth and z_{s} 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 U_{b}. In the wallfunction approach, the
resultant wall shear stress vector Ï„_{w} is related to the
velocity vector at the first grid point V_{2} by the standard,
smoothwall, log law (Wu et al., 2000):
where, the subscript 2 refers to the point at the first control volume
center at the wall, is the normal distance to the wall; is
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 nearwall
values of turbulent kinetic energy k and dissipation rate ε are given
as:
On the other hand, if the first grid node from the wall locates within
the viscous layer due
to the evolution of complex vortical structures, these nearwall values
are approximated as:
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)):
where, u_{in}, k_{in} and h_{in} 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 lockrelease gravity currents, the dense and ambient fluids
are completely divided at x = x_{o} ( 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.
RESULTS AND DISCUSSION
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 gridconverged solutions
are included.

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 gridconverged  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 largescale
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
largescale dynamics is twodimensional. 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 h_{0} and divided into
two regions by a wooden gate located at a distance x_{0} 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 h_{0}
= 14.9 cm, x_{0} = 39 cm and ρ_{s} = 1009.3 kg m3
and 2) Case DH2 where h_{0} = 10 cm, x_{0} = 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). 

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 Finegrid solutions
configuration similar to that in Huppert and Simpson (1980) but with
the following parameters: h_{0} = 28.8 cm, x_{0} = 120
cm and the reduced gravity g` = 4.7 cm s2 (g` = g (ρ_{s}ρ_{0})/ρ_{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.

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 socalled slumping
phase during which the current head propagates at an approximately constant
velocity after the current is established in a relatively short period.
Figure 3ab 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 3cd, 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 KelvinHelmholtz 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.

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 largescale 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 largescale structures.

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
lockrelease gravity currents in the slumping phase.
Numerical simulations of a fulldepth 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 LargeEddy 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 KelvinHelmholtz vortices and the resulting flow patterns
are comparable with the DNS results of a similar fulldepth 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 cm^{2} sec1, which gave a layeraveraged 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 quasiequilibrium 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 quasiequilibrium 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.

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 

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) 

In the case names, C: Continuous gravity current; G: Garcia (1993); cg:
Coarsegrid solution and fg: Finegrid solution 
In Fig. 7 and 8 we compare the streamwise
velocity profiles at selected locations and isocontours 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.

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 (isocontour 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.

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.

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 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 quasiequilibrium
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 rollup 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 wavylike structure in the streamwise
direction as the gravity current head passes the transient region.
CONCLUSIONS
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
highresolution, 2D URANS simulations using a buoyancyextended 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 energyconserving
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 threedimensional
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 KelvinHelmholtz 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 largescale 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 threedimensional 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.