A Finite Volume Implementation of the Shallow Water Equations for Boussinesq Gravity
Currents
by
Thomas Meyer Hatcher Jr.
A thesis submitted to the Graduate Faculty of
Auburn University
in partial fulfillment of the
requirements for the Degree of
Master of Science
Auburn, Alabama
August 4, 2012
Keywords: Gravity Currents, Lock-exchange, Finite Volume Method, Shallow Water Equations
Copyright 2012 by Thomas Meyer Hatcher Jr.
Approved by
Jos?e Goes Vasconcelos Neto, Chair, Assistant Professor of Civil Engineering
Prabakhar Clement, Professor of Civil Engineering
Xing Fang, Associate Professor of Civil Engineering
Abstract
The shallow water equations (SWE) are a powerful tool for modeling the propagation of grav-
ity currents (GC) because of their relative simplicity, computational efficiency and accuracy. Finite
difference solutions, either based on the method of characteristics (MOC) or the implementation
of numerical schemes such as Lax-Wendroff (LxW) have been traditionally used in such flow com-
putations. On the other hand, the finite volume method (FVM) has been gaining popularity in
several other hydraulic applications, being favored in cases when flow discontinuities are antic-
ipated. This work is focused on an implementation of the finite volume method (FVM) to the
solution of Boussinesq GC using the one and two-layer SWE models. The proposed two-layer
mathematical model is a modification of the the work by [Rottman and Simpson, 1983], adapted
to express such equations in a vectorial conservative format, amenable for FVM implementation.
The traditional solution for the GC front boundary condition (BC), using a characteristic equation
and a front condition, is compared to a new formulation that explicitly enforces local mass and
momentum conservation. Linear numerical schemes (LxW and FORCE) and non-linear schemes
based on the approximate solution of the Riemann problem (Roe and HLL) are implemented in
this framework along with various front conditions. The proposed modeling framework is tested
against experimental data collected by this investigation, and is also compared to previous inves-
tigations. Results indicate that this proposed model has a comparably simple and robust imple-
mentation, being flexible enough to be applied in a wide range of GC flow conditions and presents
good accuracy.
ii
Acknowledgments
I would like to thank my advisor, Professor Jose Goes Vasconcelos, for his experience, hard
work and friendship throughout my time here at Auburn. He introduced me to CFD and experi-
mental hydraulics and was instrumental in the development of this thesis. I also received a great
deal of help from fellow students throughout the experimental portion of this work. Parker Ross
assisted in the construction of the scale model of Mobile Bay. Kyle Moynihan and Kyle Strachan
assisted in a large number of experiments along with Paul Simmons and Carmen Chosie. There
contributions and friendship made the time I spent in the lab enjoyable. Bernardo Trindade and
Gabriel Leite helped introduce me to numerical modeling and their friendship played an important
role throughout this process.
I would like to thank my family for all of their support these last couple of years and through-
out my life. My best friend and wife, Jennifer Kelly Hatcher, has been the centerpiece in my life
these last few years and has provided essential love and support. I would like to thank her for
putting up with me throughout my time at Auburn. My parents, Tom and Sharon Hatcher, have
provided incredible support throughout my life. I could not have asked for better role models to
grow up around. I would like to congratulate my brother, Preston Hatcher, for posting a 1.73 ERA
with 10 saves this past spring and to thank him for teaching me self control as a child. Finally,
I would like to thank all of my friends in the civil engineering department who helped make my
time at Auburn enjoyable.
iii
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Description of GC flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Laboratory experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3 GC Flow stages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.4 Modeling approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1 The description of the front . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.1 Theoretical investigations of GC front conditions . . . . . . . . . . . . . . 18
2.1.2 Empirical investigations of GC front conditions . . . . . . . . . . . . . . . 21
2.2 Experimental investigations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Analytical description of GC flows with integral models . . . . . . . . . . . . . . . 24
2.3.1 Deep ambient with constant Froude number . . . . . . . . . . . . . . . . . 24
2.3.2 Generalized solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4 Numerical modeling of gravity currents . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.1 One-Layer SWE model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4.2 Two-Layer SWE model . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.4.3 Navier-Stokes models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.5 Hyperbolic numerical schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.5.1 Linear schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
iv
2.5.2 Nonlinear schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3 Knowledge Gap and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.1 Experimental program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.1.1 Physical apparatus construction . . . . . . . . . . . . . . . . . . . . . . . 51
4.1.2 Initial salinity measurements . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.1.3 Front trajectory measurements . . . . . . . . . . . . . . . . . . . . . . . . 56
4.1.4 MicroADV measurements . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2 One-layer SWE model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.3 Two-layer SWE model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.4 Aspects of the numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.4.1 Numerical schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.4.2 Boundary Conditions: DC and MOC . . . . . . . . . . . . . . . . . . . . 64
4.4.3 Stability and time step calculation . . . . . . . . . . . . . . . . . . . . . . 66
5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.1 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.2 Numerical model results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.3 Comparison between experimental results and numerical predictions . . . . . . . . 76
5.4 Comparison with previous models . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6 Conclusions and future developments . . . . . . . . . . . . . . . . . . . . . . . . . . 97
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
v
List of Figures
1.1 Aerial view of the front that forms between the Connecticut river and the Long Island
Sound estuary. Source: environmentalheadlines.com . . . . . . . . . . . . . . . . . . 2
1.2 A particle-driven GC generated from the terrorist attacks of 9/11. Source: 911re-
search.wtc7.net. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 A haboob or dust cloud smothering a residential area, my3monsters.com. . . . . . . . 5
1.4 High Reynolds number lock-exchange GC for Boussinesq fluids along a horizontal
surface presented in [Slim, 2006]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Schematic diagram of a partial depth lock-exchange experiment in a rectangular chanel.
Adapted from [Rottman and Simpson, 1983]. . . . . . . . . . . . . . . . . . . . . . . 10
2.1 Schematic diagram of Benjamin?s formulation for the case of a heavier fluid propagat-
ing underneath a lighter ambient fluid. Adapted from [Klemp et al., 1994]. . . . . . . . 20
2.2 Schematic diagram of a GC in an initially moving ambient system. Adapted from
[Wright and Paez-Rivadeneira, 1996]. . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Schematic diagram for both fronts that is generated from a GC flowing into a shallow
ambient (f > 0:5). Adapted from [Ungarish and Zemach, 2005]. . . . . . . . . . . . . 36
2.4 Two-wave structure for the HLL scheme presented in [Toro, 2001]. . . . . . . . . . . 45
4.1 The experimental apparatus with stationary ADV probes and a video camera moving
with the GC front. The initial GC depth h0 is equal to the ambient depth H. . . . . . . 52
4.2 The photograph on the left is an aerial view of Mobile Bay, Alabama. Source: Landsat
program. To the right is a snapshot of the depth-averaged velocities (arrows) and
the water levels (contours) predicted by the estuarine and coastal ocean model. The
vertical line in the middle of the bay represents the shipping canal. Presented in [Chen
et al., 2005]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3 Different components of the construction for the acrylic channel. . . . . . . . . . . . . 54
4.4 Photograph of the leading edge of a GC that illustrates the procedure in which the
front trajectories were measured. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
vi
4.5 Schematic diagram of a GC that illustrates the cell structure for the DC BC at a given
time step. The cell LE represents the last computational node that is completely filled
by the front of the GC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.1 Dimensionless trajectory and velocity for the GC front. . . . . . . . . . . . . . . . . . 70
5.2 Velocity hydrograph comparison for 1%, 2%, and 3% salt in the 9.14 m tank. . . . . . 72
5.3 GC propagation for the one and two-layer SWE over a wide range of fractional depth:
f = 0:01, 0:5 and 1:0. The DC BC is used with 400 discretization cells. . . . . . . . . 75
5.4 Front trajectory results for the one-layer SWE model using Roe scheme and DC BC
with various front conditions and from experiments. . . . . . . . . . . . . . . . . . . . 77
5.5 Front trajectory results for the two-layer SWE model using Roe scheme and DC
boundary conditions with various front conditions and from experiments. . . . . . . . 78
5.6 Slumping for the one and two-layer SWE models using the HLL and Roe scheme with
two front conditions: Benjamin (B) and Huppert Simpson (HS). x0 = 1, h0 = 1 and
e = 2% using 800 computational cells. . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.7 Velocity hydrographs with B and HS front conditions and 200 computational cells.
The experimental data (e = 2:0%) is provided from MicroADV devices in the slump-
ing stage: x = 7:68. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.8 Velocity hydrographs with the MRS front condition and 200 computational cells. The
experimental data is provided from MicroADV devices at x = 7:68. . . . . . . . . . . 89
5.9 Comparison between measured ambient velocities and two-layer SWE predictions us-
ing B and HS front conditions with 200 computational cells. The experimental data
(e = 3:0%) is provided from MicroADV devices at x = 7:70. . . . . . . . . . . . . . 90
5.10 Comparison of GC profiles at various time steps between the proposed two-layer
model and [Rottman and Simpson, 1983] model for the critical condition: f = 0:5. . . 92
5.11 Dimensionless slumping comparison between the proposed model and [Ungarish and
Zemach, 2005] model using Roe scheme and modified FORCE scheme. . . . . . . . . 93
5.12 Dimensionless front velocity results for the two-layer model over a large range of
fractional depth. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
5.13 Dimensional slumping comparison between the SWE and EFDC for a) t=1.8sec and
b) t=39sec. The discretization size for the SWE models is 800 cells. . . . . . . . . . . 96
vii
List of Tables
5.1 Experimental initial conditions and results for each experiment. Values for um were
determined from a digital camera tracking the GC front. . . . . . . . . . . . . . . . . 69
5.2 Sensitivity of each modeling approach to the number of discretization nodes. The
error for each model is compared to its computed propagation time at 1,000 nodes.
The Roe scheme is used with B front condition and kDt = 0:1. . . . . . . . . . . . . . 74
5.3 Comparison between lock-exchange experiments and one-layer SWE model predic-
tions using both BC approaches at the GC front, kDt = 0:1 and a 400-cell solution do-
main. The time of propagation for the GC in the experiments was 93.67 s (e = 1:0%)
and 67.63 s (e = 2:0%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.4 Comparison between lock-exchange experiments and two-layer SWE model predic-
tions for different time steps and front conditions using the DC boundary condition.
The GC time of propagation in the experiments was 93.67 s (e = 1:0%) and 67.63 s
(e = 2:0%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.5 Two-layer SWE results using various values for b for the RS front condition with 400-
cell solution domain. When b = 1:215, this front condition is referred to as MRS.
The GC time of propagation in the experiments was 93.67 s (e = 1:0%) and 67.63 s
(e = 2:0%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.6 Computational time comparison between the one-layer (1L) and two-layer (2L) SWE
models and between the DC and MOC front BC strategies for various discretization
sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
viii
List of Notations
The following symbols are used in this thesis:
x0 = Initial length of denser fluid
h0 = Initial depth of denser fluid
h; h1 = Depth of GC
h2 = Depth of ambient fluid
hLE = Depth of leading edge or front of GC
H = Initial ambient depth, total depth of system
f = Initial depth ratio or fractional depth = h0=H
fLE = Depth ratio of the GC leading edge
u;u1 = Velocity of GC
u2 = Velocity of the ambient fluid
uLE = Velocity of the GC leading edge
rc;r1 = Density of bottom GC
ra;r2 = Density of ambient fluid
Dr = Difference between the current and ambient densities =rc ra
e = Relative density difference =(rc ra)=rc
g0 = Reduced gravitational acceleration = ge
um = Mean velocity at the GC front
Re = Reynold?s number
Rem = Mean Reynold?s number
FrH = Total depth densimetric Froude number
FrLE = Froude number at the GC leading edge
b = Scaling factor for the front condition
ix
~U = Vector of conserved variables
~F(U) = Vector of conserved variables fluxes
~S(U) = Vector of source terms
A;B = Source term parameters
Dx = Dimension of space discretization
Dt = Dimension of time discretization
kDt = Scaling factor for time step calculation
Cr = Courant number
MOC = Method of Characteristics
DC = Dual-Cell boundary condition
x
Chapter 1
Introduction
Gravity currents (GC) (i.e. buoyancy or density currents) are driven by a density difference
between two or more fluids and generally travel in a quasi horizontal plane. In a large number of
cases the difference in density is small so that the buoyancy forces are of equal or greater magnitude
than the inertial forces. The primary discrepancies in the density are caused by two factors: the
accumulation of dissolved solids or temperature differences between fluids. Examples of GC in
which an accumulation of solids causes an increase in density include exchange flows in estuaries,
sediment entrainment in rivers and the expansion of dense volcanic dust clouds. When a cold
stream enters a warmer lake or vice versa, the density current is driven by temperature differences.
In this thesis high Reynold?s number GC flows, which are formed from the discharge of saltwater
into freshwater, are analyzed in the context of the shallow water equations (SWE).
In water bodies such as rivers and estuaries, GC flows are a common occurrence and have a
substantial impact on the flow field. When a river with a high discharge rate empties into an estuary
with negligible tidal effects, a lighter GC may form at the free surface with a salt wedge (denser
GC) propagating on the channel bottom into the river (Figure 1.1). In the Connecticut River during
the spring months, there is a large discharge due to snow melt, which forms a salt-wedge regime
with the Long Island Sound estuary (see Figure 1.1). Detailed measurements of the density profile
have been recorded for the lighter GC that forms at the water surface [Simpson, 1997].
Sea-lochs are deep inlets from the ocean in which tidal and river effects are important. In
systems with large amounts of river flow, well-mixed layers of brackish water develop on the
river side of a sill or obstacle. If the tide rises high enough, a salt wedge may rise over this sill
into the brackish region. This process has important effects on the pollution of the inlet and may
determine the amount of dissolved oxygen in an unmixed zone that lies underneath the brackish
1
Figure 1.1: Aerial view of the front that forms between the Connecticut river and the Long Island
Sound estuary. Source: environmentalheadlines.com
region. Overturn due to dense water flowing over the sill may occur at rare intervals. In Loch Eil
in Scotland, the overturn occurs at a frequency of about once a year, and this occurrence is caused
by a GC flow [Simpson, 1997].
One of the first observed GC flows was related to volcanic eruptions in which hot gas and
ash are expelled at high velocities. This ash cloud can travel hundreds of kilometers and is the
most destructive component of the eruption [Baxter, 2000]. Once the buoyancy forces outweigh
the initial momentum, the fluid rises vertically as a buoyant plume, which is outside the scope
of quasi horizontal GC flows [Turner, 1973]. As the ash cloud moves higher and higher into
the atmosphere, entrainment becomes more and more important [Slim, 2006]. Thus, the reduced
gravity is decreased in magnitude until neutral buoyancy is achieved, which is inevitable because of
the ambient stratification. At this depth, the plume begins to move horizontally as an intrusive GC.
2
Figure 1.2: A particle-driven GC generated from the terrorist attacks of 9/11. Source: 911re-
search.wtc7.net.
Similar examples occur at the ocean floor where 75% of the annual volcanic activity occurs. Small-
scale emissions through hydrothermal vents and large-scale eruptions rise as buoyant plumes and
have an important impact on ocean currents [Simpson, 1997].
Instead of forming a buoyant plume, a volcanic eruption can also spread as a pyroclastic GC
when the density of the fluid is sufficiently large. Moreover, when the ash cloud is too dense
to continue rising, the flow collapses and forms a ?ground-hugging? surge [Slim, 2006]. These
pyroclastic flows can travel at speeds of up to 160 ms 1 with temperatures reaching 900 C [Blong,
2000]. Since the late 1700s, the primary cause of death in volcanic eruptions is due to pyroclastic
flows (i.e. approximately 27% or 220,000 deaths) [Baxter, 2000]. A similar type of GC flow
occured due to the terrorist attacks on 9/11/2001 (see Figure 1.2).
A number of other types of GC flows can also occur in the context of volcanoes. When heavy
rains re-mobilize the unconsolidated ash deposited from a volcanic eruption, a lahar or mud flow
3
is formed. This hazardous type of GC can travel at high velocities for several kilometers and
can recur years after the initial eruption [Blong, 2000] and [Slim, 2006]. In 1985, the volcanic
eruption of Nevado del Ruiz in Columbia claimed 25,000 lives. Four major lahars were formed,
which traveled a total distance of 60 km at 10 ms 1. These mud flows over topped the banks a
river channel in Armero and followed the bank into the town, which resulted in the majority of the
casualties for this event [Simpson, 1997].
The GC examples presented above are higher velocity and lower viscosity flows. Thus, the
Reynolds number is large (Re > 1;000). However, the flow of molten rock is an example of a
low velocity and high viscosity flow. In addition, these low Reynolds number flows (Re < 10) are
non-Newtonian meaning that there is a nonlinear relationship between stress and strain. However,
the focus of this thesis is on the simpler Newtonian fluids, which are valid for gases and water
flows. In addition, the buoyancy forces are balanced by inertial forces instead of by viscosity so
that viscous forces are neglected.
Another common atmospheric GC is a haboob or dust cloud, which is a cold atmospheric
front several hundred meters in depth that can travel at over 20 ms 1. These types of GC flows
are generated from thunderstorm outflows over sandy and arid landscapes. The turbulence from
the outflow entrains dust, which has been measured at 40 mg=m3 [Simpson, 1997]. However,
the difference in density due to the entrained dust particles is generally a small fraction of 1%
compared to the effects of temperature. Therefore, it is common to model these haboobs only
based on the temperature differences between the two fluids [Simpson, 1997].
GC flows are often categorized depending on the description of the density difference. When
this difference is caused by temperature or dissolved solids, the GC is referred to as homogeneous.
However, sometimes GC flows are driven by suspended particles that progressively settle out of
the GC, and the density of the current is reduced. These flows are categorized as particle-driven
whereas a combination of these two is known as particle-laden. A haboob is an example of a
particle-laden GC; however, this type of flow is more homogeneous dependent. On the other
4
Figure 1.3: A haboob or dust cloud smothering a residential area, my3monsters.com.
hand, a buoyant plume that occurs due to volcanic eruptions, ocean outfalls, etc. eventually move
horizontally as a GC once neutral buoyancy is reached and is known as an intrusion [Slim, 2006].
In this thesis, the focus is on turbulent, homogeneous GC flows that travel along a horizontal
surface in a Boussinesq system (i.e. ra=rc 1 where ra is the ambient density and ra is the current
density). The thin layer assumption, which is widely used in the context of GC, is utilized in order
to validate the SWE. Two thin-layer numerical models are proposed: one-layer SWE (the ambient
fluid is neglected in the governing equations) and two-layer SWE (the governing equations are
formulated based on the current and the ambient fluids). In addition, lock-exchange experiments
are conducted to compare with the numerical models and to measure the GC velocity structure.
1.1 Description of GC flows
GC flows are driven by a density difference in the vertical direction, which results in a pressure
difference in the horizontal plane. This pressure difference is balanced by the dynamic velocity
field in the horizontal direction. The gravitational effect associated with the density difference is
5
referred to as the reduced acceleration due to gravity, g0: [Ungarish, 2009].
g0=
r
c ra
rc
g =eg (1.1)
in which ra is the ambient fluid density, rc is the density of the current and e is the relative density
difference. The reduced gravity is a fundamental feature in the GC propagation and causes the
decrease in velocity for GC compared to traditional dam break flows. This parameter replaces the
gravity (g) in the SWE and plays a role in the front condition (described in 2.1).
The first quantitative study of GC flows was done by von Karman, who evaluated the spread
of poisonous gas [Huppert, 2006]. Subsequently, [von Karman, 1940] formulated a classical rela-
tionship for the GC depth and velocity using Bernoulli?s theorem in which a heavier fluid advances
into the lighter atmosphere:
u LE
(g0h LE)1=2 = Fr LE (1.2)
in which hLE and uLE are the depth and velocity of the front or leading edge of the GC, respectively.
FrLE is the Froude number of the GC front, which was evaluated by von Karman to equalp2 in the
context of a relatively deep ambient. [Benjamin, 1968] argued that von Karman?s formulation was
invalid because he applied Bernoulli?s theorem across a streamline characterized with head losses.
Subsequently, Benjamin re-derived this expression (1.2) using the momentum integral, and ended
up with the same results. This result surprised Benjamin; however, because both approaches used
a different integral of the Euler equation, they could not have reached different results [Huppert,
2006].
The Froude number for a GC flowing into a deep ambient is actually closer to unity. More-
over, [Huppert and Simpson, 1980] obtained FrLE from their experiments and [Bonnecase et al.,
1992] stated that viscous drag and Reynold?s stresses are the reason for this discrepancy. [Shin
et al., 2004] formulated a new expression for the GC front in which FrLE = 1 for the deep ambient
scenario. Their result differs from Benjamin and von Karman?s value because [Shin et al., 2004]
approach includes the energy transfer between the GC and the ambient fluid, which is important
6
for f!0. In the experiments conducted in this thesis f = 1; however, the numerical model is able
to simulate any initial condition for the fractional depth (0 f 1).
Two problems emerge from equation 1.2. For example, when a GC propagates into a shallow
ambient fluid, the Froude number also becomes dependent on the ambient depth H or in effect,
the depth ratio or fractional depth at the leading edge of the GC (i.e. fLE = hLE=H). In addition,
(1.2) represents one equation for two unknowns (hLE and uLE) [Huppert, 2006] (further discussed in
subsection 4.4.2). [Benjamin, 1968] solved the first potential problem by deriving an expression for
FrLE that is a function of fLE by applying a flow force balance (section 2.1). For the lock-exchange
flow that is simulated in this thesis, two approaches are utilized at the front BC in conjunction with
Benjamin?s condition or an empirical alternative.
The leading edge or front of a GC forms a distinct region that separates the advancing GC
and ambient fluids. Mixing is prevalent at the GC front, but the entrainment between the two fluids
is swept upstream in the form of shear instabilities. Towards the front of the GC a head or nose
generally forms in which the flow depth is greater than that of the fluid upstream. This nose region
is characterized by ??a zone of breaking waves and intense mixing?? [Simpson, 1997]. Initially, the
nose of a GC moving along a horizontal surface is quasi-steady; however, in the case of an incline,
the size of the GC nose increases with slope [Simpson, 1997]. For an inviscid GC flowing on a
rigid and horizontal surface (Figure 1.4), the depth of the GC nose is dependent on the ambient
depth, which varies from the deep ambient f = 0 to the full depth lock-exchange problem f = 1. In
addition, the shape of the nose is dependent on the turbulence in the ambient fluid and the ambient
velocity structure.
This study is not concerned with viscous effects, but if viscosity becomes important, the rate
of advance for the GC is decreased and the turbulence at the nose is lessened. For high velocity
and low viscosity GC (i.e. Re > 1;000), the Reynolds number does not play a significant role in
the propagation. For instance, in a thunderstorm the Reynold?s number is typically in the order of
1;000;000. Moreover, field observations suggest that the flow structure in these highly turbulent
examples is similar for all Re > 1;000 GC flows [Simpson, 1997]. Viscosity can play an important
7
Figure 1.4: High Reynolds number lock-exchange GC for Boussinesq fluids along a horizontal
surface presented in [Slim, 2006].
role for flows in which Re < 1;000; moreover, when Re < 10, viscosity tends to dominant over
inertia in the GC flow [Simpson, 1997]. When a viscous oil slick spreads over top a water body,
surface tension eventually exceeds viscosity in the level of importance [Fannelop and Waldman,
1971].
For high Reynold?s number GC flows, turbulent mixing is prevalent, and the instabilities
travel upstream of the GC nose. This turbulent mixing feature is illustrated in Figure 1.4 in which
most of the information is presented from [Britter and Simpson, 1978] and [Simpson and Britter,
1979]. The mixing in GC flows is a complex process that usually occurs in the following two
forms: billows that roll up the nose of the GC and a pattern of lobes and clefts. The billows
are formed from shear instabilities at the interface of the fluids because the two fluids flow in
opposite directions. The lobes and clefts are formed by the influence of turbulence against the
ground [Simpson, 1997]. Entrainment, the mixing between fluids, plays an important role in the
propagation of the GC, and the importance of entrainment increases with slope. Moreover, the
billows that entrain the ambient fluid play an important role in mass transfer, which decreases the
front velocity.
The velocity of the GC front is one of the key parameters and has received a great deal of
attention from researchers. This value can be approximated by equating the mean potential energy
8
loss to the kinetic energy gain [Simpson, 1997]. For a lock exchange flow where the height of the
lock equals h0, this relationship becomes:
1
2g
0h0 = 1
2u
2 (1.3)
Isolating the velocity (u) in (1.3) yields:
U =
p
g0h0 (1.4)
in which U is an estimated potential velocity of the current and h0 is the initial GC depth. Although
the velocity estimation in (1.4) is an unrealistic simplification, the results are typically of the same
order of magnitude [Ungarish, 2009]. It is customary to normalize the results for the velocity with
this velocity estimation.
As discussed above, one of the key parameters in GC flows is the Reynolds number, which
governs the importance of viscosity. An estimate of the magnitude for the Reynold?s number can
be determined from equation 1.4 [Simpson, 1997]:
Re =Uh0=n (1.5)
in which n is the kinematic viscosity. For this study, we are concerned with high Reynolds number
flows (i.e. Re >> 1;000), and the final viscous stage that may eventually develop is neglected
[Simpson, 1997]. In addition, the present work is focused on thin layer systems (h0=x0 < 1) so that
the SWE are more applicable because there is less interfacial mixing [Ungarish, 2009].
1.2 Laboratory experiments
Experiments are conducted to validate the proposed numerical model using a MicroADV
device as well as digital cameras. One of the advantages of experimental research is that the
GC flows are similar to large-scale environmental examples. Moreover, experiments are typically
9
easier to conduct than field measurements and it is easier to control parameters such as the effects
of the ambient fluid velocity. Experimental investigations from [Britter and Simpson, 1978] and
[Simpson and Britter, 1979] provided important information on the GC structure especially at
the nose that forms near the leading edge. One of the first GC experiments (the lock-exchange
problem) was conducted by removing a vertical barrier that separates freshwater from saltwater
[Keulegan, 1957] and [Barr, 1967]. A similar approach was conducted in this study in which high
definition digital cameras tracked the GC front while MicroADV probes measured the internal GC
velocity distribution.
Constant volume lock-exchange experiments are typically conducted by separating a denser
fluid with a density r+Dr, length x0 and depth h0 from a lighter ambient fluid r of depth H in a
long rectangular tank (see Figure 1.5). For full depth lock-exchange flows, h0 = H where as partial
depth experiments are conducted when h0 < H. For partial depth releases, the ambient fluid of
depth H h0 is added on top of the denser fluid so that the depths on both sides of the gate are
equal (for a more comprehensive discussion, see [Shin, 2001]). Once this vertical gate is removed,
a GC flow begins in which the denser fluid propagates underneath the lighter ambient. When h0
equals H, a GC flow also develops at the free surface moving in the opposite direction of the denser
bottom current. Once the upstream moving depression wave (or hydraulic jump) that develops for
the denser fluid reaches an upstream physical boundary, a reflection occurs, and a nose region is
developed towards the front of the GC [Rottman and Simpson, 1983].
Experiments are often preferred to large-scale field measurements in the context of GC flows
[Shin, 2001]. Laboratory experiments are easy to repeat in order to establish consistency with the
data collection. They are less expensive and can be performed in a more controlled environment
so that the data collection is of higher quality than for field trials. In addition, it is easier to isolate
a certain aspect of the flow for a better understanding [Shin, 2001]. For field measurements, it is
often difficult to eliminate the effects of wind, boating and tidal effects in estuaries.
10
Figure 1.5: Schematic diagram of a partial depth lock-exchange experiment in a rectangular chanel.
Adapted from [Rottman and Simpson, 1983].
The differences between laboratory experiments and environmental flows can occur because
of the difference in length scales. In order to achieve dynamic similarity between these two meth-
ods, the relative importance of the controlling forces (e.g. inertial, viscous, etc.) should be main-
tained [Shin, 2001]. The dimensionless Reynold?s number should be consistent for the laboratory
experiments.
UL
n
Prototype
= ULn
Experiment
(1.6)
where U is the velocity scale and n is the kinematic viscosity. The length scale L within the
Reynold?s number is typically smaller in the laboratory; however, the problem is reduced if water
is used instead of air for atmospheric flows [Shin, 2001]. Moreover, the kinematic viscosity of
water is 0.01 cm2s 1 at 20 C and 0.15 cm2s 1 for air [Simpson, 1997]. Another important factor
is that GC flows are weakly dependent on the Reynold?s number when the Reynold?s number is
greater than 1,000 as previously described.
The Peclet number is another relevant dimensionless parameter that is related to mass trans-
port phenomena.
Pe = ULk (1.7)
11
in which k is the diffusivity of the fluid. The Peclet number is the ratio of the advective and
diffusive terms in the momentum equation [Shin, 2001]. In the majority of GC applications, the
order of magnitude of the Peclet number is very large, sometimes greater than 1 million. Therefore,
diffusion effects are negligible [Shin, 2001]. In laboratory experiments, which are presented in this
thesis, the diffusivity of pure salt is 1.1 x 10 5cm2s 1 so that the Peclet number was also greater
than 1 million [Shin, 2001].
1.3 GC Flow stages
The inviscid GC propagation in a rectangular channel is described by three stages [Ungarish,
2009]: slumping, transition, and self-similar. The transition stage is the least interesting of the three
stages because it more or less serves as a bridge between the other two. However, the beginning of
the transition stage (i.e. the end of the slumping stage) is an important feature that effects the GC
depth and velocity (further discussed in 5.2).
At the initial release that marks the beginning of the slumping stage, the depth decreases
dramatically in comparison with the other stages. Moreover, the shape of the interface and the
velocity of the GC undergo a transformation. Although the depth (h) drops to at least half of the
initial depth (h0) after the initial release due to the energy conserving theory of [Benjamin, 1968],
the velocity and the depth of the GC nose remain almost constant for the remainder of the slumping
stage.
[Ungarish, 2009] describes the slumping stage in two sub-categories starting with a motion
analogous to a dam-break wave, which is characterized by the initial drop in depth and increase
in velocity for the GC. This type of motion occurs at the beginning of the lock exchange flow
when stationary fluid remains in the system (i.e. upstream of the release). Towards the front of the
GC ?the interface has a negative slope, corresponding to a rarefaction (expansion) pressure wave,
under which the fluid is accelerated? [Ungarish, 2009]. Next, the entire depth of the current begins
to descend abruptly sending a disturbance towards the nose. A positive slope for this disturbance
begins to form, which travels towards the front of the GC [Ungarish, 2009].
12
The second stage is transitional in nature in which the GC ?is under the influence of the
initial conditions and the front boundary conditions? [Ungarish, 2009]. Once the downstream
propagating bore or wave (depending on f) that originated at a physical boundary reaches the GC
front, the depth and velocity of the front begin to decrease. Moreover, the quasi steady motion
that forms during the slumping stage is ended and there is an unsteady transition to the self-similar
stage. The importance of the initial conditions begins to decrease as the GC propagates further
downstream and enters this self-similar stage.
In the final stage (i.e. self-similar), the inertial forces begin to balance with the buoyancy
forces to retain an almost steady flow. For the inviscid assumption, the GC flow remains in the self-
similar stage until there is contact with a physical boundary (wall) or the simulation is terminated
after a certain minimum GC thickness (i.e. h!0) otherwise viscosity would have to be included
in the computations. Because the depth and the velocity decrease with time, ??the importance of
the viscous friction relative to the inertial terms increases monotonically during the self-similar
propagation?? [Ungarish, 2009]. The following relationship formulated in [Huppert and Simpson,
1980] describes the time when the inviscid assumption becomes invalid:
x =(x50h50gr=n2)17 (1.8)
According to this relationship, the inviscid assumption is valid for all of the experiments and nu-
merical simulations conducted in this study,
1.4 Modeling approaches
One of the first applications in which GC were studied was oil spreading over top of a water
body. Because of the relatively small thickness of oil compared to water in the context of oil spills
within the ocean, the solution of the governing equations can be simplified to analytical solutions.
First, it is customary to reduce the Navier-Stokes equations into the depth-averaged SWE, which
are referenced throughout this paper.
13
[Fannelop and Waldman, 1971] and [Hoult, 1972] used the SWE in order to derive self-similar
solutions. They focused on the three major phases in the propagation of oil on a water body:
gravity-inertial, gravity-viscous, and surface tension-viscous. [Fannelop and Waldman, 1971] as-
sumed that the SWE are time-independent, which is required in order to derive these analytical
solutions, so that time derivatives vanish in the governing equations.
In order to derive the SWE in the context of oil spreading using a single pair of governing
equations (one-layer SWE), some simplifications are made. For example, the viscosity in the mo-
mentum equation is neglected because the viscosity of crude oil is one to two orders of magnitude
greater than water, so the horizontal spread of oil is considered independent of the vertical direc-
tion. However, if the Reynold?s is greater than 1,000, viscosity is typically neglected for saltwater
intrusions as well. The following derivation follows the formulation in [Fannelop and Waldman,
1971]; however, subtle changes are made in the notation for consistency with this thesis, which
focuses on a heavier GC propagating underneath a lighter ambient fluid (e.g. saltwater into fresh-
water).
[Fannelop and Waldman, 1971] derives the 1-D one-layer SWE from the 2-D Navier-Stokes
equations with the incompressible and Boussinesq assumptions:
?u
?x +
?v
?y + ju=x = 0 (1.9)
?u
?t +u
?u
?x +v
?u
?y =
1
r
?P
?x +
m
r
?2u
?y2 (1.10)
?v
?t +u
?v
?x +v
?v
?y =
1
r
?P
?y g (1.11)
in which u and v represent the longitudinal and vertical velocity components, respectively. P is
the pressure, r is the density and m is the shear stress. This set of equations (1.9-1.11) is valid
for both 2-D and axisymmetric GC flows. For the planar case, j = 0, and for the radial case,
j = 1. In order to derive the 1-D SWE from these 2-D Navier-Stokes equations, the thin layer
assumption is utilized. Thus, equations (1.9-1.11) are integrated across the GC layer depth (h),
14
and v?u?y becomes zero in the x-momentum equation. In addition, the vertical accelerations are
neglected (i.e. v ?v?y = 0). After integration, (1.11) is reduced to the hydrostatic relation between
the pressure in the layer and the depth. Then, an expression is derived by differentiating this
pressure relation with respect to x and assuming the weight of the GC is equal the weight of the
displaced ambient fluid.
?P
?x =rg
0?h
?x (1.12)
in which g0 = (r1 r2)=r1 where r1 is the density of the heavier current and r2 is the density
of the ambient fluid. Equation 1.12 is substituted into equation 1.10 in order to obtain a single
momentum equation. This thesis is concerned with inviscid GC, so the last term in the RHS of the
x-momentum equation (1.10) is omitted and the resulting expression for momentum is:
?u
?t +u
?u
?x +g
0?h
?x = 0 (1.13)
Following the formulation in [Fannelop and Waldman, 1971], equation (1.13) is not sufficient
to solve for the two unknowns (h and u). Therefore, the continuity equation is required alongside
of this momentum equation. In order to formulate the updated continuity equation from (1.9), the
vertical accelerations are neglected. Then, the difference in the velocity (Dv) is equated to the
continuity equation integrated along the depth:
Dv =
Z h
0
?u
?x + ju=x
dh = h?u?x juh=x (1.14)
in which:
Dv = ?h?t +u?h?x
Therefore, the resulting continuity equation is written as:
?h
?t +u
?h
?x +h
?u
?x + ju=x
= 0 (1.15)
15
Together, equations (1.13) and (1.15) form the 1-D SWE, which are written in primitive for-
mat. This set of equations is amenable to the planar and axisymmetric scenarios, which are widely
studied for both analytical models and SWE models. However, this thesis is not concerned with
axisymmetric GC, so the radial term is omitted in the continuity equation. Moreover, the chain
rule is utilized in order to write the resulting SWE in conservative format (i.e. in terms of h and
uh), which is preferred near shocks: [LeVeque, 1992]:
?h
?t +
?uh
?x = 0 (1.16)
?uh
?t +
?
?x
(uh)2
h +
1
2g
0h2
!
= 0 (1.17)
16
Chapter 2
Literature Review
In this chapter, the theoretical and empirical GC front investigations are discussed. Some
of the relavent experimental investigations into GC flows are presented. Subsequently, the use of
analytical box models to simulate these GC flows are discussed before moving on to the more
complex numerical models. The one and two-layer SWE models and the corresponding boundary
conditions are analyzed in detail since this discussion leads to the proposed numerical models
presented in the Methodology, Chapter 4. The advantages and disadvantages between the SWE
models and the Navier-stokes models are presented in section 2.4.3, and this chapter is concluded
with a discussion of hyperbolic numerical schemes that are used in SWE models.
2.1 The description of the front
The description of the GC front (i.e. h and u / uh) has received considerable attention in recent
decades. It is well known that the front velocity is very important to the GC propagation and
interfacial mixing. In the widely used integral models (discussed in 2.3), the description of the
front is one of the primary components along with volume continuity. Moreover, the front of the
GC is typically utilized as a boundary condition in the SWE. The SWE are invalid at the GC front
because of the importance of vertical accelerations, so the addition of a front condition is required
to complete the model (described ahead in 4.4.2).
One of the most traditional problems in unsteady, open channel hydraulics is the 1-D dam
break problem, which is similar to the description of a GC front. In the dam break problem with a
dry, frictionless bed scenario on one side, the front of the fluid (i.e. water) approaches a feather tip
solution [Sturm, 2010]. Through characteristic analysis, Ritter was able to compute the velocity of
the dam break front, which is described by the following relationship: 2u = c0 where c0 =pgh0
17
is the initial wave celerity [Klemp et al., 1994]. In addition, a similar expression was obtained for
the velocity of the upstream moving depression wave: u = co.
[Schoklitsch, 1917] compared these dam break results to the propagation of GC by means of
laboratory experiments. It was apparent that the upstream moving wave of the GC compared well
with the dam break analysis. However, the velocity of the GC front was about 1=2 of the dam break
front velocity. This discrepancy lead to a new relationship (equation 1.2), which was formulated
by [von Karman, 1940]. [Abbott, 1961] altered this relationship by replacing the Froude number
with an empirical coefficient b that was determined from experiments:
uLE =b
p
g0hLE (2.1)
in which the subscript LE denotes the front or leading edge of the gravity current. One notices that
this expression (2.1) equates to the dam break front velocity when b = 2 and the density of the
ambient fluid is much smaller than the GC fluid density [Klemp et al., 1994]. In the description of
GC propagation, the correct value for b or the Froude number has been the focus of a number of
research investigations. Because of the importance and intrinsic complexity of interfacial mixing
(specifically at the GC front), empirical relationships are often applied in SWE models. However,
there have been many attempts to improve the theoretical description of the GC front.
2.1.1 Theoretical investigations of GC front conditions
From experimental observations, it is known that the celerity of a GC front depends on the
density difference (Dr) between the fluids and the depth ratio at the GC front region (i.e. fLE). In
equation 2.1 the density difference is incorporated in the form of the reduced gravity. However, the
depth ratio is not involved in this expression, so it is expected that b is a function of the depth ratio.
One of the main theoretical breakthroughs for the description of the GC front was formulated by
[Benjamin, 1968] in which the Froude-like number is a function of the fractional depth.
18
In order to formulate this expression at the leading edge of the GC, [Benjamin, 1968] equated
the buoyancy forces to the drag forces between the two fluids at the front. For Ritter?s traditional
dam break problem, the ambient fluid density is neglected as well as the shear with the ambient
fluid, which leads to the feather tip front. In GC flows these simplifications are strictly invalid.
Benjamin?s analysis was applied to a system in which an air cavity intruded over top of the water-
filled pipe with a rectangular cross-section undergoing discharging at one end. The following
development follows [Klemp et al., 1994], who applied Benjamin?s analysis in the context of
a heavier fluid traveling underneath a lighter fluid (e.g. the propagation of a saltwater wedge into
freshwater). In order to obtain the relationship for the GC front velocity, the flow is assumed steady
and depth-averaged. The vertical momentum equation (1.11) is simply reduced to the hydrostatic
relationship [Klemp et al., 1994]:
?P
?z = g
0 (2.2)
The x-momentum equation is formulated from equation 1.10, and the effects of shear are
neglected. This expression is integrated from a point far ahead r of the GC front to a point l
far behind the GC across the channel of depth H, which yields the following flow force balance
[Klemp et al., 1994]: Z
H
0
P
r +u2r
dz =Z H
0
P
l +u2l
dz (2.3)
in which ur is the ambient velocity at r and ul is the ambient velocity at l. Pl and Pr represent the
pressures at l and r, respectively (see Figure 2.1).
In order to compute this flow force balance, the continuity of volume is required: ul =
urH=(H h). Moreover, Pl is determined by applying the vertical momentum equation (2.2)
far behind the nose: Pl = 12u2r g0h. When these two expressions are applied to (2.3) along with
the assumptions that Pr = 0 far ahead of the front and Ps = 12u2r at the stagnation point x f , the result
is Benjamin?s formula:
uLEp
g0hLE = Fr(fLE)=
(2 f
LE)(1 fLE)
1+fLE
1=2
(2.4)
19
Figure 2.1: Schematic diagram of Benjamin?s formulation for the case of a heavier fluid propagat-
ing underneath a lighter ambient fluid. Adapted from [Klemp et al., 1994].
This important expression represents one equation for two unknowns: uLE and hLE and is a solution
for equation 2.1. From this relationship (2.4), it is clear that the velocity of the GC depends on
the reduced acceleration due to gravity, the depth and the depth ratio. The implementation of this
expression in numerical models is discussed in subsection 4.4.2.
[Shin et al., 2004] extended Benjamin?s expression to include the effects of mass transfer
between the front of the ambient fluid advance and the front of the heavier fluid. Their results are
exactly the same as Benjamin?s for the lock-exchange problem because the interfacial waves (in
which energy is transferred between the two fluids at the interface) generated from the ambient
front are unable to reach the leading edge of the GC. However, the results from [Shin et al., 2004]
expression improve in accuracy over Benjamin?s solution as the depth ratio decreases to a deep
ambient problem [Shin et al., 2004]. For Boussinesq fluids, the expression formulated by [Shin
et al., 2004] is:
Up
g0h = Fr(fLE)=
1
2
s
h
H
2 hH
(2.5)
This expression was reduced from the original non-Boussinesq relationship; furthermore,
[Shin et al., 2004] expects the non-Boussinesq relationship to break down for large density dif-
ferences. This reasoning is based on their assumption that the current and ambient fronts are
connected by a horizontal interface with uniform flow conditions. This assumption is warranted
20
for Boussinesq systems; however, in non-Boussinesq fluids the interface is more complex and this
assumption is invalid [Keller and Chyou, 1991].
2.1.2 Empirical investigations of GC front conditions
Huppert and Simpson (HS): This front condition was constructed empirically from the ex-
periments of [Huppert and Simpson, 1980] and [Keulegan, 1957]. It performs better for full
depth lock releases but is able to simulate all depth ratios [Ungarish, 2009].
Fr(fLE)=
8>
><
>>:
1
2f
13 (0:075 f
LE < 1)
1:19 (fLE 0:075)
(2.6)
Rottman and Simpson (RS): In this front condition, Benjamin?s theoretical condition is al-
tered with an empirical front condition b. [Rottman and Simpson, 1983] adjusted b based
on their partial depth experiments because their model was unable to simulate large initial
values for f. On the contrary to the HS front condition, the RS front condition performs well
for small fractional depths but not as well for full depth releases [Ungarish, 2009].
Fr(fLE)= bp2
(2 f
LE)(1 fLE)
1+fLE
1=2
(2.7)
Ungarish and Zemach (UZ): [Ungarish and Zemach, 2005] formulated a new front condition
that bridges the HS and RS front conditions, so it is able to perform well for a large range of
fractional depth.
Fr(fLE)=(1+3fLE) 1=2 (2.8)
Kranenburg (K): [Kranenburg, 1978] adjusted Benjamin?s theoretical front condition to al-
low for an initially moving ambient. In addition, they utilized an empirical coefficient k that
they estimated to equal 0.6. While the focus of this thesis is on counter-flows in which the
ambient fluid has an initial condition of zero velocity, in many hydraulic applications the
21
Figure 2.2: Schematic diagram of a GC in an initially moving ambient system. Adapted from
[Wright and Paez-Rivadeneira, 1996].
ambient fluid has non-zero velocity, either coflows or counter-current flows with respect to
GC (e.g. cold discharges in river flows). The expression presented by [Kranenburg, 1978] is:
C1
(grH)1=2 = Fr(h=H)=
h=H(1 h=H)(2 h=H)
1+h=H +k(1 h=H)
1=2
(2.9)
in which C1 = uLE q=H is equal to the velocity uLE of the GC for ambient fluids with initially
zero velocity. Notice that the LHS of Equation 2.9 is a function of the ambient fluid depth
instead of the GC depth. The schematic diagram for this system is illustrated in Figure 2.2.
Wright (W): [Wright et al., 1990] developed a generalized relationship that is able to simu-
late both co-flows and counter-flows. The depth ratio f in this relationship is based on the
depth behind the nose (see Figure 2.2). Several experiments conducted for various ambient
velocities demonstrated the effectiveness in this approach [Wright and Paez-Rivadeneira,
1996].
C1
(grH)1=2 = Fr(h=H)=[1 (1+qr)h=H]
h=H(1 h=H)3
(1 h=H)3 q3r(h=H)3
1=2
(2.10)
where qr is the velocity ratio between the ambient fluid and the GC (qr = q2=q1) in which q = uh.
22
2.2 Experimental investigations
Along with analytical and numerical research, experimental investigations have provided im-
portant insights into GC flows. [Simpson and Britter, 1979] studied a large range of dimensionless
depths for GC flows and determined the resulting dimensionless velocities, rate of mixing, and
the depth of the mixed layer behind the head. Their apparatus consisted of a moving floor and a
downstream weir, which allowed them to halt the GC flow in order to obtain more accurate mea-
surements. A similar procedure was utilized by [Wright and Paez-Rivadeneira, 1996] to analyze
the effects of coflows and counterflows. For lock-exchange flows, [Rottman and Simpson, 1983]
analyzed a large range of initial depth ratios (f) and compared the results their two-layer SWE
model. For large initial depths, [Rottman and Simpson, 1983] concluded experimentally that the
upstream moving jump occurs when f 0:7 and that the reflection off of the upstream boundary
produces a hydraulic drop that eventually overtakes the GC front marking the beginning of the self
similar stage. For smaller initial depths, the upstream moving jump is replaced by a depression
wave [Rottman and Simpson, 1983].
In [Hacker et al., 1996] a digital image technique (DigImage, [Dalziel, 1993]) was used to de-
termine the density structure for lock-release GC flows. They were able to observe the detrainment
of dense fluid at the GC nose by breaking Kelvin-Helmholtz waves, which produced a stratified
region behind the nose. The diluted fluid in this region was replaced by dense fluid behind the
nose, which traveled toward the front of the nose. Thus, this process causes a re-circulation pro-
cess and the entire GC becomes diluted. The work by [Hacker et al., 1996] utilized the following
expressions to compute the total depth densimetric Froude number (FrH) and the mean Reynolds
number (Rem) for the lock-exchange problem:
FrH = umph
0 g0
Rem = 12 um h0n (2.11)
where h0 is the initial depth of the denser fluid, g00 is the initial reduced gravity, um is the mean
velocity of the GC front and nu is the kinematic viscosity. Rem is multiplied by 12h0 because this is
23
the energy conserving depth described by [Benjamin, 1968]. These dimensionless parameters are
provided for the experiments conducted for this thesis in section 4.1.
[Gerber et al., 2011] used a particle image velocimetry (PIV) to measure the Reynolds stress
and shear production of turbulence of a stably stratified GC. Their experimental results compared
well to their Reynold?s averaged Navier-stokes numerical model. [Firoozabadi et al., 2010] used
MicroADV probes to measure the turbulence energy, Reynolds stresses and turbulence intensity
of 3-D GC flows for various discharges, concentrations and slopes. As the discharge and the
concentration increased, the turbulence kinetic energy also increased. Their results also indicate
that the normalized turbulence intensity does not change with bed slope.
A large number of experimental investigations have utilized sodium chloride (NaCl) in their
experiments to generate Boussinesq fluids ([Rottman and Simpson, 1983], [Shin et al., 2004],
etc.). The same approach was implemented in the experiments presented in this thesis. For non-
Boussinesq lock-exchange flows, [Lowe et al., 2005] used sodium iodide (NAI) in addition to
NACL in order to obtain density ratios (r2=r1) between 0.6 and 1.0. From their non-Boussinesq
experiments, they determined that in most cases the lighter current retains its energy-conserving
depth while the denser current is dissipative, so the depth decreases according to a decreasing
density ratio.
2.3 Analytical description of GC flows with integral models
For a quick and simple estimation of the GC propagation, integral models (also referred to as
box models) can be a useful modeling option. In integral models only the resistance of the advance
(the front condition) and the conservation of buoyancy are imposed. Thus, GC flows for an integral
model are described by the spreading of an equal area rectangle, so horizontal variations in the flow
are ignored [Slim, 2006]. While this is an unrealistic assumption; however, these modeling results
can provide important insights into the average flow depth and velocity [Ungarish, 2009].
24
2.3.1 Deep ambient with constant Froude number
The first integral model is derived for a homogeneous GC propagating into a deep ambient
f 0. However, similar models have also been proposed for particle-driven GC for both incom-
pressible [Hogg et al., 1999] and compressible fluids [Timmermans et al., 2001]. In the context
of the deep ambient, the Froude number, which is the ratio between the inertial and gravitational
forces, is assumed constant (valid for the deep ambient problem). For the two-dimensional or pla-
nar case, there are two key variables in the internal cell calculations: GC length xN(t) and height
hN(t). As a result, two equations are required to complete the box model [Ungarish, 2009]. The
first equation is the simple volume continuity equation:
xN(t)hN(t)= V (2.12)
in which V is the volume of the dense GC. For the fixed volume GC (e.g. the lock-exchange
problem) with inertial-buoyancy balance (neglecting viscosity), the derivation begins with the front
condition. This BC, which can be derived through dimensional analysis, with the Froude number
(Fr) notation used in [Benjamin, 1968] is written as:
uN = dxNdt = Fr(g0hN)1=2 (2.13)
in which dxN=dt represents the GC front velocity; the Froude number is defined as Fr = u=pg0hN.
This BC is the same relationship that is implemented for the SWE (described ahead in this paper).
Unlike in Benjamin?s formulation, the Fr in (2.13) is assumed constant for this derivation. The
combination of (2.12) and (2.13) yields the following relationship for xN(t) (or for hN(t)) [Ungar-
ish, 2009]: Z
t2
t1
= 1V 1=2
Z xN
x 0
[Fr(f(x))] 1x 1=2dx (2.14)
25
where f(x)= V=(Hx) through algebraic manipulation.
x N = x 0
1+ 32 Fr(g
0 V)1=2
x 0 3=2 t
!2=3
(2.15)
Once xN is determined from (2.15), h is computed from the continuity equation (2.12). The same
strategy can be implemented to obtain an integral model (2.16) for the axisymmetric scenario:
r = r0
1+2Fr(g
0 V)1=2
ro 2p1=2
!
(2.16)
where r is the radius and r0 is the initial condition.
2.3.2 Generalized solution
[Huppert and Simpson, 1980] developed an innovative integral model for both the two-dimensional
and axisymmetric scenarios. Instead of using a constant Fr, they implemented the empirical front
condition (HS) that was derived from their experiments and from [Keulegan, 1957]. As a result,
they obtain a solution for both the slumping and the self-similar stage, which depends on the depth
ratio (f). For the two-dimensional case, [Huppert and Simpson, 1980] substitutes the continuity
equation (2.12) into the HS front condition. For the slumping stage, this relationship is written as:
Fr = dxNdt =(g0 V=xN)1=2 = 12
V
xNH
1=3
(0:075 hN=H < 1) (2.17)
After integration and isolating xN, the slumping stage relationship is:
xN = x0
1+ 712
g03 V H 2
x0 7
1=6
t
!6=7
(x0 xN xs) (x0 x xs) (2.18)
in which x0 is the initial lock length and xs is the slumping distance, which [Huppert and Simpson,
1980] define as the distance when (hN=H 0:075)
26
For the self-similar stage (i.e. hN=H 0:075), the authors propose that the Froude number is
1:19 according to the HS front condition. Following the same procedure as the slumping stage,
the relationship for the self-similar stage is:
xN =
h
x3=20 +1:78(g0 V)1=2(t ts)
i2=3
(xs < x ) (2.19)
in which x is the distance when viscous effects become larger than inertial effects (1.8) [Huppert
and Simpson, 1980]:
For the axisymmetric problem, which is derived in a similar manner to the two-dimensional
case, the solution in the slumping stage is: [Huppert and Simpson, 1980]:
r =
r4=30 23p 1=6(g02 V H2)1=6t
3=4
(r0 r rs) (2.20)
where rs and r in cylindrical coordinates are equivalent to xs and x for the two-dimensional
problem. The relationship for the self-similar stage is: [Huppert and Simpson, 1980]:
r =[r2s +2:37p 1=2(g0 V)1=2t]1=2 (rs r r ) (2.21)
In summary, integral models are simple to implement and require minimal computational
time as compared to the SWE and Navier-Stokes models since they lend themselves to analytical
solutions. The development of such models does not require numerical methods or schemes, which
can become extremely complex in other models (e.g. higher order Riemann solvers). The cost of
this reduction in complexity is the use of assumptions: that GC spread in a quasi-steady manner
with a constant depth in an exactly rectangular shape. Nevertheless, integral models can provide
important insights into the flow of GC in a channel.
The results for the deep ambient scenario are satisfactory depending on the context. However,
when the initial depth ratio approaches the lock-exchange problem (f > 0:5), the integral model
is unable to accurately capture the depth profile and other flow features although the Huppert and
27
Simpson formulation is a viable alternative. As a result, more complex models (e.g. SWE and
Navier-Stokes) are required for the accurate simulation of these flows.
2.4 Numerical modeling of gravity currents
Experimental studies (described in section 2.2) are an imperative tool for the understanding
of GC flows. However, it is not practical to construct scale models for every type of scenario that
is present in the field. It is very difficult to simulate experimentally several common flow features:
tidal intrusions, complex geometries, non-Boussinesq systems, flows over porous substrates, etc.
Analytical and numerical models are essential for a better understanding of GC propagation in
which experimental data and/or field observations are used to calibrate these numerical models.
Once calibrated, these numerical models can become powerful tools, yielding useful predic-
tions of GC flows that are very important to the scientist or engineer in the field, who may need to
determine the mixing mechanisms and/or flow field for a variety of issues such as the transport of
contaminants. While analytical models are a viable option for these situations, they suffer from se-
vere assumptions as described in the previous text. As a result, many GC flow applications require
the use of numerical modeling. There are two primary methods for the numerical modeling of GC
flows: SWE models and Navier-Stokes models, which are described below.
2.4.1 One-Layer SWE model
The SWE models are generally separated into two groups: one-layer and two-layer. The
two-layer SWE models differ in that they include the effects of the ambient fluid velocity in the
internal cell calculations. Both of these models can provide an accurate estimation of the GC front
and are computationally efficient compared to the more complex Navier-Stokes models. In the
following text the details of the SWE models are analyzed and several of the important research
investigations are included.
The one-layer SWE models (4.2) are widely used for GC flowing into deep ambient fluids.
[Ungarish and Huppert, 2004] used the one-layer SWE to simulate a GC traveling at the base of
28
a stratified ambient and determined that the stratification reduces the GC velocity. The key at-
tractions in the single-layer SWE are the simplicity and the computational efficiency. Because the
velocity of the ambient fluid is neglected, discontinuous solutions are generally not observed for
flows in simplified rectangular geometries without obstacles. Details of this one-layer formula-
tion are found in [Ungarish, 2009] for both Boussinesq and non-Boussinesq fluids. In this study,
Boussinesq fluids are focused; however, it is relatively simple to formulate the one-layer model
without the Boussinesq assumption. Therefore, the more generalized one-layer, non-Boussinesq
SWE model is the focus of this discussion:
?h
?t +
?uh
?x = 0
?uh
?t +
?
?x
(uh)2
h +
1
2g
0h2
!
= 0
For this SWE model, the one-layer assumption (i.e. the ambient fluid velocity is neglected)
poses inaccuracies for larger depth ratios. However, the single-layer SWE model provides fairly
accurate results for the GC trajectory for depth ratios approaching unity f!1 even though the one
layer assumption is least valid [Ungarish, 2007]. For the full depth lock-exchange problem, the ac-
curacy is believed to be caused by a fortunate balance of errors from different simplifications (e.g.
ambient velocity effects, entrainment, etc.) [Ungarish, 2007]. The ambient fluid velocity increases
the GC velocity, but this increase is balanced by the reduction in velocity due to entrainment. The
single-layer SWE model cannot simulate either one of these complex flow features, but the com-
bined effect of neglecting these features equates to an accurate estimation for the GC front velocity.
On the other hand, the flow depth in the initial slumping stage is inaccurately computed by this
SWE model for larger initial depth ratios. Therefore, it is customary to move to the two-layer SWE
model for this scenario.
[Ungarish, 2007] formulated a single-layer SWE model for a large range of depth ratios and
density differences. Furthermore, in this study Boussinesq GC flows for the lock-exchange prob-
lem were analyzed as a special case of their non-Boussinesq model. Although there are several
29
important underlying assumptions in this approach, the generality, absence of empirical coeffi-
cients and computational efficiency are attractive [Ungarish, 2007]. The characteristic equations
for the one-layer model that are implemented in [Ungarish, 2007] are:
g012 dhh1=2 du = 0 on dxdt = u (g0h)1=2 (2.22)
in which the second expression represents the characteristic flow velocities. Although these equa-
tions (2.22) can be used to compute the internal cell calculations, [Ungarish, 2007] implemented
the Lax-Wendroff two step finite difference method. This numerical scheme is well documented
and performs well in the absence of shocks [Chaudhry, 2008], which are generally not present in
the one-layer model.
However, the characteristic approach was still used in the computations by [Ungarish, 2007].
Since Benjamin?s theoretical front condition was used at the GC front BC, a characteristic equation
was also used to solve for both unknowns (i.e. hLE and uLE). Therefore, the traditional solution to this
problem has been to use a characteristic equation along with a front condition at this BC (described
further in 4.4.2). For the other boundary condition that is located upstream of the release, either a
reflective [Toro, 2001] or characteristic [Sturm, 2010] boundary condition may be used.
One of the inconvenient components of SWE models for GC is the need for a boundary condi-
tion at the GC front. Therefore, [Ruo and Chen, 2007] formulated the modified SWE for inviscid
GC in order to eliminate the need for a front condition. Therefore, the only required boundary
conditions are at the physical boundaries or end walls. In order to formulate their mathematical
model, [Ruo and Chen, 2007] included the effects of the ambient resistance, which is normally
accounted for by the front condition.
The continuity equation for the modified SWE is the same as in the traditional SWE . How-
ever, the momentum equation is altered in order to account for the ambient resistance. [Ruo and
Chen, 2007] begins their derivation with equations (1.9) and (1.11). Then, the pressure term is
30
reduced from the following relationship, which is also used in the two-layer formulations:
P(x;y;t)= Pi(x;t)+rg(H h) (2.23)
where Pi(x;t) is the pressure at the fluid interface, and H and h are the ambient and current depth,
respectively. In addition, equation (1.11) is implemented in order to develop a new x-momentum
equation:
?u
?t +u
?u
?x +gr
?h
?x (2.24)
However, a new problem emerges in this mathematical model: How to determine the pressure
at the fluid interface. [Ruo and Chen, 2007] implements the unsteady Bernoulli equation along
the instant streamline in order to overcome this issue. Their solution for Pi(x;t) is the following
relationship:
P i(x;t)= Pt ragh 12ra(eVi)2 ra?F?t
y=h(x;t) (2.25)
in which Vi is the absolute velocity of the ambient fluid at the interface, F is the velocity potential so
that V = F, and Pt is the total pressure of the ambient fluid at a location far from the current [Ruo
and Chen, 2007], which in this case is constant. When equation 2.25 is substituted into equation
2.24, it becomes clear that the new difficulties become computing eVi and f. Furthermore, [Ruo
and Chen, 2007] breaks apart eVi into two components: tangential velocity wt and normal velocity
wn. It is then assumed that ?(w2n)=?x >>?(w2t )=?x, which leads to the following relationships:
Vi = wn = ?h=?tp1+(?h=?x)2 ??x
?F
?t
= ??x
(?h=?t)2
1+(?h=?x)2
(2.26)
so that momentum equation for the modified SWE is [Ruo and Chen, 2007]:
?u
?t +u
?u
?x +gr
?h
?x =
g
2
?
?x
(?h=?t)2
1+(?h=?x)2
(2.27)
31
in which g = ra=r. A look at this highly nonlinear momentum equation suggests that (2.27) is
very similar to (1.17). However, there is an added source term on the RHS of (2.27), so [Ruo and
Chen, 2007] solved this momentum equation using the perturbation method. Although they obtain
accurate results that compare well to the more traditional SWE models for smaller depth ratios,
[Ruo and Chen, 2007] does not believe that their approach can be extended to the lock-exchange
problem.
2.4.2 Two-Layer SWE model
The two-layer SWE have received more attention than the single-layer alternative for GC flow
modeling. For instance, the lock-exchange problem is a relatively easy experiment to perform in
which the two-layer SWE model is preferred. Numerous research investigations have focused on
a wide range of GC flows and the use of the two-layer SWE. For example, [Bonnecaze et al.,
1993] applied the two-layer SWE to particle-driven GC flows. However, this discussion focuses
on homogeneous GC flows on a horizontal surface. As previously described, the two-layer SWE
are preferred when f!1. When f!0, the results converge to the one-layer SWE model results.
[Rottman and Simpson, 1983] developed the first two-layer SWE model for GC flows in the context
of partial-depth releases. For larger fractional depth (i.e. f > 0:5), they encountered difficulties,
which are described below. In the case of a two-layered model, the SWE represent a system with
four PDEs to account for each phase?s mass and linear momentum conservation [Rottman and
Simpson, 1983].
?h1
?t +
?
?x(u1h1)= 0 (2.28)
?h2
?t +
?
?x(u2h2)= 0 (2.29)
?u1
?t +u1
?u1
?x =
1
r1
?P i
?x g
?h1
?x (2.30)
?u2
?t +u2
?u2
?x =
1
r2
?P i
?x g
?h1
?x (2.31)
32
The subscripts (1, 2) represent the lower (heavier) and upper (lighter) fluid layers, respectively.
[Rottman and Simpson, 1983] implemented equation (2.28) as the sole continuity equation for
their model. Then, the ambient variables (i.e. h2 and u2) were computed from the following rela-
tionships:
h2(x;t)= H h1(x;t) (Depth uniformity) (2.32a)
u2 = u1h1h
2
(Local continuity) (2.32b)
where equation 2.32a is deduced from the rigid lid assumption in which the total depth of the
system does not change. The relationship for the ambient velocity (equation 2.32b) is formulated
from the two continuity equations (2.28 and 2.29) along with 2.32a and is valid for constant volume
gravity currents.
In order to utilize one relationship for the equation of motion and to eliminate Pi (the pressure
at the fluid interface), equations (2.30) and (2.31) are substituted together into a single momentum
equation:
(1+ra)?u1?t +
1 raH +h1H h
1
u1?u1?x +
g0 r(1+a)3 u
21
H
?h
1
?x = 0 (2.33)
Furthermore, [Rottman and Simpson, 1983] implemented the Boussinesq simplification, which
eliminates r (i.e. r2=r1), yielding:
?u1
?t +(1 2a)u1
?u1
?x +(1 b)g
0?h1
?x = 0 (2.34)
The mathematical model that is formulated in [Rottman and Simpson, 1983] serves as the
basis of the proposed model, which is presented in section 4.3. The two-layer mathematical model
formulated in [Rottman and Simpson, 1983] is presented below in primitive format. By using
the assumption of local flow continuity and overall depth uniformity, the governing equations are
33
reduced to 2 PDEs that represent the conservation of mass and momentum of the GC:
?h1
?t +
?
?x(u1h1)= 0
?u1
?t +(1 2a)u1
?u1
?x +(1 b)g
0?h1
?x = 0
(2.35)
where
a = h1H h
1
b = h1H + u
21
g0H
1 h1H
2
(2.36)
[Rottman and Simpson, 1983] implemented the first expression in (2.35) as the sole continuity
equation for their model along with the single momentum equation. However, one of the draw
backs is clear from the parameters a and b in which the solution is undefined as the fractional
depth approaches unity h!H. A solution to this problem is proposed in section 4.3.
The context in which these two-layered models are most often implemented is in the release
of dense fluids into relatively small ambient depths (e.g. the lock-exchange problem). The dense
fluid is often originally contained by a barrier that is suddenly removed and a wall is located
upstream of the release, which represents a physical boundary. Once the four PDEs are reduced
to a single pair of PDEs, the numerical solution that is implemented in [Rottman and Simpson,
1983] involves characteristic analysis where these PDEs are transformed into the system of ODEs
presented below:
h1 du1dh
1
(1 2a)u1 +c = 0
c = u1(1 a) u21a2 +g0h1(1 b) 1=2
(2.37)
One of the main benefits in this approach is computational efficiency and simplicity. However,
the method of characteristics is not able to simulate shocks in which the characteristic lines of the
same family intercept one another [Chaudhry, 2008].
In addition to the internal cell calculations, the c+ equation is implemented at the GC front
boundary condition to provide closure to solve for the two local unknowns, hLE and uLE. The results
in [Rottman and Simpson, 1983] compare well with experiments and exhibit minimal oscillations
throughout the slumping stage; however, their front condition is only tested for smaller depth ratios
34
(f 0:5). Furthermore, their front condition (RS) decreases in accuracy for larger depth ratios,
especially for the full depth lock-exchange problem [Ungarish and Zemach, 2005].
When f 0:5, the numerical results for [Rottman and Simpson, 1983] two-layer SWE model
are non-physical (i.e. there is a multi-valued solution for the depth). Moreover, the upstream
propagating shock that is generated for f > 0:5 becomes stronger as it travels at higher speeds, and
the reflection off of a back wall triggers the formation of another sharp interface that eventually
catches up with the GC front. While the numerical modeling suggests that the limit for this sharp
interface to develop is f >0:5, experiments from [Rottman and Simpson, 1983] point that this limit
is in fact closer to f > 0:7, albeit the transition is not exactly clear. The discrepancy, according to
[Rottman and Simpson, 1983], is a result of the smooth profile in the velocity for the experiments
and the mixing between the two fluids, which is neglected in the shallow water model. As a result,
the two layer model developed by [Rottman and Simpson, 1983] is accurate quantitatively for f
0:5; the authors pointed that numerical simulations with f > 0:5 would require new developments
in the theory of flow discontinuities.
The work by [Klemp et al., 1994] extended the two-layer model developed by [Rottman and
Simpson, 1983] in order to solve for the depth ratios greater than 0.5. Instead of using the method
of characteristics, [Klemp et al., 1994] implemented a finite difference model using the second-
order accurate Leapfrog scheme for the solution of the two-layered SWE. They enforce a bound-
ary condition at the upstream propagating discontinuity based on an analytical relationship for
hydraulic jumps. [Klemp et al., 1994] specified the pressure far downstream pr(d) according to
the following expression:
pr(d)= 12 V 2f u022 +g0D (2.38)
in which Vf is the speed of the upstream moving jump, u02 is the velocity of the ambient fluid just
downstream of the jump and D is the dissipation term (see Figure 2.3).
Applying the flow force balance far upstream and downstream and using equation 2.38 and the
continuity relations, [Klemp et al., 1994] obtained the following solution for the upstream moving
35
shock:
V 2f
g0H =
f2 a2 2D0
[2f2 a=a]+(1 f)2(1 2a)=(1 a)2 (2.39)
in which h f is the height of the upstream moving jump, f = h0=H, a =(h0 h f)=H and D0=D=H.
[Klemp et al., 1994] set the dissipation term (D) equal to zero for consistency with Benjamin?s
theory and based on the evaluation of weak hydraulic jumps. Equation 2.39 is not sufficient to
solve for the moving jump BC, so [Klemp et al., 1994] implemented an additional condition for
steady-state jumps based on the formulation in [Stommel and Farmer, 1952]:
u
1 +Vf
2
g0h1 +
u
2 +Vf
2
g0h2 = 1 (2.40)
The schematic diagram for a GC flowing into a shallow ambient fluid is illustrated in Figure 2.3.
This analytical relationship was required for their SWE model for f > 0:5 because the leading
edge of the backward propagation wave travels slower than the trailing region until the front of the
wave steepens into a discontinuity. [Klemp et al., 1994] describes this phenomenon as a hydraulic
jump according to long wave theory; however, this jump differs from the classical hydraulic jump
because it does not form an abrupt discontinuity. Instead, [Klemp et al., 1994] points that the slope
of this jump remains finite due to non-hydrostatic effects, which are outside the scope of long wave
theory. [Klemp et al., 1994] presented their results for a large range of depth ratios; however, they
have focused their model on the initial stages of the current propagation, up to the instant when the
upstream moving front contacts the end wall.
The work by [Ungarish and Zemach, 2005] presents a model that extends the theory of
[Rottman and Simpson, 1983] and [Klemp et al., 1994] in order to include the entire slumping
stage in their formulation, extending the simulation beyond the instant when the backward prop-
agating front contacts the upstream wall boundary. [Ungarish and Zemach, 2005] implemented a
second-order accurate finite difference technique similar to [Klemp et al., 1994], but they used the
Lax-Wendroff (LxW) two-step numerical scheme. The model performs an explicit tracking of the
upstream propagating front before and after the reflection off the back wall, and these fronts are
36
Figure 2.3: Schematic diagram for both fronts that is generated from a GC flowing into a shallow
ambient (f > 0:5). Adapted from [Ungarish and Zemach, 2005].
treated as discontinuities in the depth. In a certain way, the treatment of these fronts resembles a
shock fitting approach used to track open channel bores presented in [Cunge et al., 1980]. When
the shock contacts the end wall, the model by [Ungarish and Zemach, 2005] predicts the propaga-
tion of a discontinuity towards the nose, which is also explicitly tracked. [Ungarish and Zemach,
2005] track this bore until the end of the dam break stage where the depth and the velocity of the
current begin to decrease.
One of the drawbacks in the approach by [Ungarish and Zemach, 2005] is the added model
complexity in introducing two additional boundary conditions to track the flow discontinuities. In
addition, the model requires two solution methods depending on the depth ratio. If f > 0:5, an
explicit feature tracking approach is required; however, if f 0:5, the simpler approach is imple-
mented in which no shock tracking methods are utilized except at the front of the GC [Ungarish,
2009]. In addition, experimental results indicate the discontinuities that form for two-layer GC
flows are not represented by an exactly vertical interface between the two fluids. Therefore, it
is an idealized approximation in which additional experimental confirmation is necessary. This
motivated the development of an unified and more general shock-capturing approach that is com-
putationally consistent for all depth ratios, which is formulated and discussed in section 4.3.
37
Although the BC at the GC front, upstream hydraulic jump and the reflected bore have been
discussed, there is still another feature that needs to be addressed. For gravity currents flowing into
shallow ambients (f 0:8 or larger depending on FrLE), a critical condition emerges at the GC
nose [Ungarish, 2009]. This constraint occurs when b* < 1 (see equation 2.41) so that the GC front
cannot travel faster than the c+ characteristic, which would cause cavitation [Klemp et al., 1994].
Therefore, the depth and velocity constraints are computed by setting b* = 1. This critical change
of behavior is clearly seen from the following dimensionless, two-layer characteristic relationships
[Ungarish, 2009]:
8
>><
>>:
dh*
du* =
1
1 b*
h
a*u* p(a*u*)2 +(1 b*)h*
i
; on
dx*
dt* = c = u
*(1 a*) p(a*u*)2 +(1 b*)h*
(2.41)
in which the asterisk denotes dimensionless variables. The parameters, a* and b*, are equivalent to
the expressions in equation 2.36 in dimensionless form. In equation 2.41 and throughout Chapter
5, the dimensional variables are converted to dimensionless variables, denoted with an asterisk,
according to the following procedure [Ungarish, 2009]:
fx*; h*; t*; u*g=fx=x0; h=h0; t=T; u=Ug (2.42)
where
U = g0h0 1=2 T = x0U
By using the critical nose approach described above, the constraints when implementing Ben-
jamin?s condition are h = 0:347H and u* = 0:527. Thus, the depth and velocity cannot be greater
than these constraints.
More recently, [Adduce et al., 2011] formulated an alternative two-layer SWE model in which
the conserved variables are: r1h1, r2h2, u1 and u2. Their primary objective was to include the
effects of entrainment and shear forces, which decrease the GC front velocity. These features
38
enable [Adduce et al., 2011] to implement the theoretical front condition from [Shin et al., 2004].
However, their approach is more complex and is still largely dependent on empirical parameters.
For the majority of SWE numerical models, the important effects of entrainment and shear
forces require the use of an empirical front condition. However, [Adduce et al., 2011] avoided this
approach by included these effects in the SWE formulation:
?(r1h1)
?t +
?(r1h1u1)
?x =r2ue
?(r2h2)
?t +
?(r2h2u2)
?x = r2ue
?u1
?t +
?
?x
u2
1
2 +
r
2h2 +r1h1
r1
gcosq
= gsinq t 1b +t 21r
1h1
?u2
?t +
?
?x
u2
2
2 +
r
1h1
r1 +
r2h2
r2
gcosq
= gsinq + t 2b +t 12r
2h2
(2.43)
in which r1 and r2 are the densities of the denser GC and the ambient fluid, respectively. ue [L=T]
is an entrainment parameter and q is the angle between the bottom of the tank and the horizontal
[Adduce et al., 2011]. t1b and t2b represent the shear stress at the lower and upper boundaries,
respectively. t12 and t21 are the interfacial shear stresses in each direction. The entrainment
parameter (ue) was computed from an empirical relationship based on the formulation in [Ellison
and Turner, 1959]. [Adduce et al., 2011] altered their expression in order to simulate the full depth
lock-exchange system:
ue
u1 =
k F2r
F2r +5 (2.44)
Then, [Adduce et al., 2011] computed the density of the GC r1(t) according to the following
expression:
r1(t)=
M1 +
Z t
0
dt
Z xn(t)
0
r2ue dx
V1 +
Z t
0
dt
Z xn(t)
0
ue dx
(2.45)
in which M1 is the mass and V1 is the volume of the GC.
The advantage of this approach is that it takes into account the entrainment between the two
fluids, which allows for a more accurate simulation of the GC front. However, the accuracy for this
39
approach is highly dependent on the empirical coefficient k, which was calibrated from experiments
to be equal to 0.48. One questions whether other cases of GC flows would have the same k value. In
addition, the model complexity and the computational time for this model should be substantially
greater compared to the approach from [Ungarish and Zemach, 2005], or the model proposed in
this thesis.
2.4.3 Navier-Stokes models
For GC research, the majority of Navier-Stokes models are solved using Direct Numerical
Simulation (DNS) or Large Eddy Simulation (LES). These models are able to simulate interfacial
mixing and provide high-resolution results at the GC nose. A large number of insights into GC
flows has been gained through high-resolution DNS models. [H ?artel et al., 2000] determined that
the stagnation point is not at the GC front but at a point just upstream. [Necker et al., 2002] com-
pared 2-D and 3-D models for particle-driven GC and concluded that 3-D models are required to
accurately model the spanwise vortices or interfacial instabilities. These features survive over long
times for 2-D simulations but are more rapidly broken down in the 3-D models. [Birman et al.,
2005] investigated non-Boussinesq GC and confirmed that the denser GC front dissipates more for
larger density differences whereas the lighter ambient front remains very close to the energy con-
serving theory seen in Boussinesq systems. Moreover, the denser front is reduced to substantially
less than half the channel height depending on the density difference while the light front remains
near half the channel height. Although DNS models are the best choice theoretically, there are
many drawbacks. The most important limitations are the intensive computational requirements
and the increased complexity. Moreover, for 3-D simulations in present-day computing systems,
the DNS method is limited to GC flows where Re = O(1;000) [Slim, 2006].
Another option for high resolution simulations of GC flows is Large Eddy Simulation (LES).
[Ozgokmen et al., 2009] analyzed several LES models and concluded that a hybrid approach,
which combined the dynamic eddy viscosity and approximate deconvolution approaches, com-
pared best to the DNS results and was 1,000 times more computationally efficient. However, in
40
the first Navier-Stokes models that were used to simulate GC flows the details of the interfacial
instabilities were not the primary objective. [Klemp et al., 1994] solved the two-dimensional non-
hydrostatic incompressible Boussinesq equations in addition to the two-layer SWE to determine
the validity of Benjamin?s front condition. The results from [Klemp et al., 1994] two-dimensional
simulations indicate good agreement in comparison with the two-layer SWE. The main differences
between their two models are the simulation of the upstream moving disturbance (for f > 0:5) and
the ?hydrostatic circulation in the vicinity of the head? [Klemp et al., 1994]. [Ungarish and Hup-
pert, 2004] also found good agreement between the one-layer SWE and Navier-Stokes simulations
in the context of a GC propagating at the base of a stratified ambient fluid.
[Klemp et al., 1994] also determined that Benjamin?s front condition performed well in com-
parison with their model. Moreover, they concluded that the discrepancy between the SWE models
and experiments is not caused by Benjamin?s theory but by other simplifications in the SWE mod-
els the most important being the inability to simulate interfacial mixing. The validity of Benjamin?s
theoretical expression is further discussed in section 5.3 with the use of MicroADV devices.
2.5 Hyperbolic numerical schemes
Since the proposed numerical models for this thesis are implemented for the one and two-layer
SWE, which is an example of a hyperbolic PDE, the focus is on hyperbolic numerical schemes.
The majority of research and implementations for these numerical schemes is based on the more
traditional water-air component of hydraulics. In the context of dam-breaks, an extensive compar-
ison for both linear and nonlinear explicit numerical schemes is presented in [Zoppou and Roberts,
2003] for the SWE. Moreover, [Vasconcelos and Nunes, 2009] presented a detailed analysis for
all of the numerical schemes that are implemented in this thesis: Lax-Wendroff (LxW), Modified
FORCE (MFORCE), HLL, and Roe. They analyzed each of the numerical schemes based on error
analysis (i.e. L1 and L?) and computational efficiency in the context of the dam break problem.
Although the Roe scheme was generally the most accurate of the compared schemes, the HLL
41
scheme was considered the most appropriate choice based on the improved computational perfor-
mance and similar accuracy. In the absence of shocks, the LxW scheme performed very well. The
MFORCE scheme was accurate for each scenario, but the computational time rivaled the nonlin-
ear schemes, which are considered superior [Vasconcelos and Nunes, 2009]. The implementation
of these numerical schemes to inviscid GC using the finite volume method is further discussed in
4.4.1.
Numerical schemes are applied on broad range of applications in hydraulics. More recently,
nonlinear numerical schemes based on the solution of the initial value problem presented by Rie-
mann (i.e. Riemann solvers) have been applied in a large range of hydraulic problems that include:
one and two-dimensional dam-breaks [Zoppou and Roberts, 2003], the structural collapse of wa-
ter supply reservoirs in urban areas [Zoppou and Roberts, 1999], the simulation of shocks in pipe
filling events [Vasconcelos et al., 2009], etc. However, these schemes are rarely utilized in GC
models. The main drawback for numerical schemes based on the Riemann problem is the in-
crease in complexity without a clear increase in accuracy. A method to determine the adequacy
of a numerical scheme is to assess computational efficiency, the possibility of shocks, etc. How-
ever, whenever explicit schemes are used, it is essential that the CFL condition is satisfied. The
maximum time step that is allowed for this condition is computed from the following expression:
Dt =Cr Dx juj+c (2.46)
in which Cr is the Courant number.
In order to sustain model stability when using explicit schemes, Cr should be less than or
equal to one. When the time step is well below this limit and shocks exist, linear schemes suffer
either from numerical diffusion or oscillations depending on the type of scheme [Godunov, 1959].
On the other hand, nonlinear numerical schemes are not affected by this issue. For example, in
the simulation of pipe filling events, it may be possible that the air phase requires small Courant
numbers in certain locations along the pipe. In order to sustain a stable simulation, the Courant
42
condition (2.46) must be satisfied and nonlinear schemes are much more effective [Vasconcelos
et al., 2009]. Implicit schemes can overcome this difficulty; however, if the Courant number is
increased beyond unity, the accuracy of the simulation is compromised when shocks are present.
In this section of the thesis, some of the most commonly used linear and nonlinear explicit
schemes are presented and examples of their implementation are provided. This discussion begins
with linear schemes, which are presented in the Methodology subsection 4.4.1, all in the context
of the Finite Volume Methods (FVM). The linear schemes that are subsequently presented are all
Centred methods. The advantage of Centred schemes with respect to alternatives (e.g. Upwind
schemes) is that ?explicit local information on wave propagation for constructing approximations
to the numerical flux? is not required [Toro, 2001]. In other words for more complex flow events
in which is difficult to determine the location and the time at which the information is generated,
Centered-cell schemes are preferred. Upwind schemes are superior for other problems (e.g. a
contaminant is discharged at an upstream boundary) in which it is clear where the information
generated.
In a large number of recent models, the FVM is used because of its ability to implement
nonlinear schemes, which are outside the scope of the FDM. While the FVM method is imple-
mented in this study (presented in 4.4.1) for both linear and nonlinear schemes, the use of the
FDM in previous GC investigations is more frequent. These numerical schemes are implemented
for mathematical models of the following conservative form:
?~U
?t +
? ~F(U)
?x =
~S(U) (2.47)
where ~U is the vector of the conserved variables, ~F(U) is the flux vector and ~S(U) represents
the source terms. For the FVM, the conserved variables are updated according to the following
expression: The various FVM schemes propose expressions for the interfacial fluxes ~F(U)n+1=2i+1=2 in
order to solve this update equation for the conserved variables.
~Un+1i =~Uni + Dt
Dx
~
F(U)n+1=2i 1=2 ~F(U)n+1=2i+1=2
+Dt ~S(U) (2.48)
43
2.5.1 Linear schemes
For GC flows, the most widely used linear schemes are presented:
Lax-Friedrich (LxF): This 1st order accurate, linear and monotone scheme is an extension
of the unstable FTCS (forward-time centred-space) scheme. The general approach for GC
problems is to move on to second-order accurate schemes due to the problems with numerical
diffusion that are exhibited in first-order linear schemes. In addition, the second-order accu-
rate schemes are more accurate for smaller discretization sizes. The Lax-Friedrich scheme
is provided in [Toro, 2001] in FVM format:
~ULFi+1=2 = 1
2
h ~
F(U)ni + ~F(U)ni+1
i
+ 12 DxDt
h~
Uni ~Uni+1
i
(2.49)
Lax-Wendroff (LxW): The LxW scheme is second-order accurate and is solved with predictor
and corrector steps. This scheme is widely used in hydraulic problems and is one of the best
options if shocks are not present [Vasconcelos and Nunes, 2009]. However, in the presence
of flow discontinuities, the LxW scheme produces spurious non-physical oscillations in the
results. Despite this drawback, this scheme is widely used in the context of GC with the
aid of artificial viscosity [Ungarish and Zemach, 2005], which reduces the magnitude of
the oscillations. One of the drawbacks of this approach is that the amount of diffusion that
needs to be introduced is not easy to estimate a priori in order to avoid excessive smearing
[LeVeque, 1992]. The Lax-Wendroff two-step scheme in FDM format is provided in [Sturm,
2010] for the St. Venant equations. In the following relationship, the Lax-Wendroff scheme
is rewritten in a general format that is also applicable to the 1-D SWE.
8
>><
>>:
~Un+1=2i+1=2 = ~Uni +~Uni+12 + Dt2Dxh ~F(U)ni ~F(U)ni+1i
~Un+1i =~Uni + DtDxh ~F(U)n+1=2i 1=2 ~F(U)n+1=2i+1=2i+Dt ~F(U)ni
(2.50)
44
MacCormick: This second-order accurate linear scheme yields results that are similar to
the Lax-Wendroff scheme. However, [Zoppou and Roberts, 2003] determined that that the
MacCormick scheme is slightly more accurate in the context of the dam break problem. In
that study the error analysis was performed for the following parameters: normalized depth
and normalized velocity. [Zoppou and Roberts, 2003] determined that the error associated
with the combination of the MacCormick scheme and artificial viscosity is almost double
the amount of the MacCormick scheme alone depending on the damping parameter. Unfor-
tunately, the presence of spurious oscillations in the presence of shocks generally requires
artificial viscosity. Despite this drawback, the MacCormick scheme has been utilized in a
large number of GC models. [Adduce et al., 2011] implemented the MacCormack scheme
in FDM format for their two-layer SWE model based on the implementation of [Garcia and
Kahawita, 1986] for the St. Venant equations, which are closely related to the 1-D SWE.
8
>>>
>>><
>>>
>>>
:
~Upi =~Uni DtDxh ~F(U)ni ~F(U)ni 1i+Dt ~S(U)ni predictor
~Uci =~Upi DtDxh ~F(U)pi+1 ~F(U)pi i+Dt ~S(U)pi corrector
~Un+1i = ~Uni +~Uci2 step n+1
(2.51)
Leapfrog: This numerical scheme differs from the MacCormick and LxW schemes in that it
is second-order accurate in time but not in space. In addition, this scheme is not as widely
used as the previous two second-order accurate linear schemes. However, the Leapfrog
scheme was implemented successfully in the work by [Klemp et al., 1994]. However,
because they end their simulation when the upstream moving jump contacts the physical
boundary, the performance of the Leapfrog scheme for the lock-exchange problem is not
completely clear. It is expected that this scheme does not perform well after the reflection
45
Figure 2.4: Two-wave structure for the HLL scheme presented in [Toro, 2001].
from the upstream physical boundary, which is discussed in chapter 5. Details of this numer-
ical scheme is described in [Sturm, 2010] in FDM format for the St. Venant equations:
~Un+1i =~Un 1i Dt
Dx
~
F(U)ni+1 ~F(U)ni 1
(2.52)
2.5.2 Nonlinear schemes
For the nonlinear schemes, the focus is on two of the most widely used implementations of
the Riemann problem: HLL and Roe. The use of nonlinear schemes in GC models is not frequent,
and for the lock-exchange problem, no examples were found. Therefore, the following analysis is
based on more traditional open channels flows and follows the formulation in [Toro, 2001].
HLL: The HLL scheme is widely used in gas dynamics and hydraulics. However, it has
rarely if ever been implemented for the SWE in the context of lock-exchange GC flows. The
nonlinear HLL scheme assumes a two-wave structure that is correct for purely 1-D problems
(see Figure 2.4). The HLL flux (FHLL) occurs in the intermediate or star region between the
two wave speeds, SL and SR. If shear waves or species equations are introduced, the HLL
scheme is inadequate because the middle wave is ignored [Toro, 2001].
46
In order to utilize the FVM update expression (2.48), the interfacial fluxes ( ~F(U)i+1=2) must
be computed. For the HLL scheme, ( ~F(U)i+1=2) is determined from the following relation-
ship:
~F(U)i+1=2 =
8
>>>
>>><
>>>
>>>
:
FL
FHLL = SRFL SLFR +SRSL(uR uL)S
R SL
FR
(2.53)
where SL and SR are the wave speeds, the subscript L denotes the ith cell and the subscript R
denotes the cell i+1. These wave speeds can be computed from the following relationship
that was suggested by [Toro, 2001]:
SL = uL cLqL SR = uR +cRqR (2.54)
and qK is given by:
qK =
8
>>>
<
>>>
:
s
1
2
(h
+hK)h
h2K
if h > hK;
1 if h hK
(2.55)
There are a multiple options available in which to determine the variables in the intermediate
or star region (h and u ). [Toro, 2001] recommended the Two-Rarefaction Riemann Solver:
8
>><
>>:
h = 1g0 12(aL +aR + 14(uL uR) 2;
u = 12(uL +uR)+aL aR
(2.56)
Roe: For the SWE, the first-order accurate Roe scheme is presented in [Toro, 2001], which
forms the basis of the following formulation. The accuracy of the Roe scheme compares
well to the HLL scheme [Zoppou and Roberts, 2003]; however, Roe?s scheme is generally
47
more computationally intensive [Vasconcelos and Nunes, 2009]. The nonlinear Roe scheme
was originally applied to the Euler equations, but was first applied to the SWE by [Glaister,
1988]. For the Roe scheme, the quasi linear form of (2.47) is approximated with a linear
system:
?~U
?t +A
?~U
?x = 0 (2.57)
where A is the Jacobian matrix.
To compute the Roe averages, there are two common methods for SWE models: Roe-Pike
[Roe and Pike, 1984] and Glaister [Glaister, 1988]. For the Roe-Pike approach presented in
[Toro, 2001], the Roe-averages are computed from the following relationships:
(Roe Pike)
8
>>>>
>><
>>>
>>>
:
?u = uL
ph
L +uR
ph
Rp
hL +phL
?h =phLhR
?c =
q
1
2(c
2L +c2R)
(2.58)
where ?u, ?h and ?c are the Roe averages for the velocity, the depth and the celerity, respectively.
For the Glaister approach the only difference is in the calculation of ?u, which is presented
below:
(Glaister)
?u = uL
ph
R +uR
ph
Lp
hL +phR (2.59)
Once the Roe averages are determined, the eigenvalues ( ?lj) and eigenvectors ( ?R(j)) are com-
puted from the following relationships, which are presented in [Toro, 2001]:
?l1 = ?u ?a
?l2 = ?u+ ?a
?R(1) =
h 1
?u ?c
i
?R(2) =
h 1
?u ?c
i (2.60)
48
In addition, the wave strengths are computed from the Roe averages:
wavestrengths
8
>><
>>:
?a1 = 12
h
Dh ?h?cDu
i
?a2 = 12
h
Dh+ ?h?cDu
i (2.61)
The final step in the Roe scheme is to compute the inter-cell fluxes so that the conserved
variables can be updated.
~F(U)i+1=2 = 1
2
~
F(U)ni + ~F(U)ni+1
12
2
j=1
?ajj?ljj?R(j) (2.62)
49
Chapter 3
Knowledge Gap and Objectives
Although there have been significant insights into GC flows in recent years, there are still
several unanswered questions:
The link between the SWE models with theoretical front conditions from [Benjamin, 1968]
and [Shin et al., 2004] and the experimental results is not clear. [Klemp et al., 1994] states
that the error from the model is due to the SWE assumptions and not from the theoretical
front conditions. Since the SWE models do not account for interfacial mixing, which slows
down the GC advance, the velocity is overestimated.
In experimental studies, it is difficult to determine the depth of the nose due to interfacial
instabilities (e.g. Kelvin-Helmholtz). [Shin et al., 2004] presents a solution for this problem
where h is a function of (g, r, H). However, it is difficult to determine the density variation
from experiments.
For state of the art two-layer SWE models, two modeling approaches are required depending
on the depth ratio (f 0:5 and f > 0:5). For larger initial depth ratios, the model is much
more complex and involves cumbersome, explicit tracking of flow features. A simpler and
general model that is valid for all depth ratios is obviously desired.
Nonlinear numerical schemes are becoming more and more popular in gas dynamics and
hydraulics; however, there use in GC flows is still incipient. Therefore, the advantage of
these nonlinear schemes on GC is not well known.
The BC at the front of the GC has been traditionally solved with a front condition and a
characteristic equation. However, for more complex mathematical models (e.g. the two-
layer SWE), the characteristic equations also become more complex [Rottman and Simpson,
50
1983] and [Klemp et al., 1994]. In addition, the continuity errors are oftentimes unaccept-
able (especially when depth constraints are implemented at the GC front). A relevant, robust
solution for GC front calculations that improves continuity errors would be a clear improve-
ment.
This work is linked to the development of a new numerical model based on the FVM to simu-
late the propagation of GC in the environment that can describe all relevant stages of the propaga-
tion up to the self-similar stage. The proposed model should rely on a simple implementation and
on conservative solutions of the SWE to simulate waves and flow discontinuities. The only explicit
tracking that is enforced in the proposed model is at the GC front, and in this work several front
condition alternatives were tested and compared. Details of the numerical model development as
well as the experimental procedure are presented in the Methodology (Chapter 4).
51
Chapter 4
Methodology
The procedure for the experimental investigations and the numerical modeling developments
is discussed in the subsequent text. Both components are constructed in such a way to compare the
front trajectory and to develop velocity hydrographs near the middle of the tank. The development
is based on the full depth lock-exchange problem for Boussinesq fluids in a rectangular tank.
Examples of the notation and a sketch of the experimental apparatus are provided in Figure 4.1.
4.1 Experimental program
Lock-exchange experiments were conducted in this work for three different density differ-
ences e = 1%;2%and3%) with the following initial conditions: x0 = 76:2cm and h0 = 40:6cm.
There was no initial or constant ambient motion in the experiments conducted for this study. In-
stead of analyzing a large range of initial parameters (e.g./ x0, h0, Dr, etc.), the experimental
contributions are based on the analysis of the MicroADV probe results. However, the trajectory
of the nose is also obtained from high-definition digital cameras for a comparison with the SWE
models. Several experiments were conducted in smaller scale tanks in order to develop a consistent
experimental program that involved the mixing procedure and the gate removal. Once consistent
results for the GC trajectory were obtained in the smaller tanks, the experiments were performed
in the 9.14 m tank that was used in this work. Two to three duplicates were performed for each run
to ensure that the results were consistent and properly obtained.
4.1.1 Physical apparatus construction
The 9.144 m long acrylic tank that is utilized in this thesis aims to represent a section of the
navigation channel of Mobile Bay in Alabama, USA (see Figure 4.2. The full scale model was
52
Figure 4.1: The experimental apparatus with stationary ADV probes and a video camera moving
with the GC front. The initial GC depth h0 is equal to the ambient depth H.
used in a separate work to investigate the exchange flows between the Mobile and Tensaw rivers
and the Gulf of Mexico that occur within the bay. The rectangular acrylic tank represents a channel
9.754 m that connects two reservoirs (see Figure 4.3). For this thesis, the channel was separated
from the rest of the physical model with two wooden barriers approximately 2 cm thick with rubber
attachments at each end. After these barriers were locked in place by two metal clamps, the length
of the channel was reduced by 0.305 m at each end. In order to perform the experiments, a third
gate was constructed to separate the two fluids, which initiated the GC flow upon removal.
The acrylic channel was constructed in four segments because the acrylic sheets are only 2.44
m in length. In each component a 15 cm wide sheet that represents the channel bottom was welded
to two side sheets (40.6 cm wide). Therefore, the maximum depth that could be achieved in the
channel is 40.6 cm. After the side pieces were welded on top of the bottom sheet, the effective
channel width was 12.7 cm.
Before the the acrylic sheets were welded together, openings were drilled into the acrylic.
Then the acrylic sheets were tightened together with metal screws in order to increase the welding
53
Figure 4.2: The photograph on the left is an aerial view of Mobile Bay, Alabama. Source: Landsat
program. To the right is a snapshot of the depth-averaged velocities (arrows) and the water levels
(contours) predicted by the estuarine and coastal ocean model. The vertical line in the middle of
the bay represents the shipping canal. Presented in [Chen et al., 2005].
adhesion and to provide structural support. In the front of the physical model on the outside of
the tank, small acrylic pieces were used to connect the adjacent acrylic sheets (Figure 4.3). On
the inside of the tank, adhesive tape was placed between the acrylic pieces with silicon caulking
applied along the edges in order to obtain a strong waterproof seal. As a result, the entire channel
formed a cohesive and waterproof unit.
4.1.2 Initial salinity measurements
The difference in density can be obtained by dissolving salt in one side of the tank CITE or
by using a different temperature CITE. If the primary difference between the fluids is temperature,
then the colder fluid will flow underneath the warmer fluid. In this study salt is dissolved in one
side of the tank, so the saltwater will flow underneath the ambient freshwater.
When dissolving the salt into the denser fluid, it is imperative that the resulting density is
accurately computed. The methods that have been implemented to determine the density are:
salinity tables, polynomial equations, hydrometers and refractive index measurements. Although
54
Figure 4.3: Different components of the construction for the acrylic channel.
55
reference tables are readily available, this method is not the most accurate because the volume of
the fluid is known only to a few percent. In addition, the salt is never pure and may vary from one
brand to another. Moreover, salt may contain certain amounts of water and may not completely
dissolve in the fluid.
An alternative and more accurate approach is to determine the refractive index of the solution,
which depends on the amount of salt. [Shin, 2001] utilized an electronic refractometer, which is
usually designed for medical use, with standard tables in which the density is determined from the
value in Brix. According to [Shin, 2001], the error associated with this approach is 0.1 Brix or
0.005 gcm 3. The smallest value for the density difference in their experiments was 0.1 gcm 3,
so the maximum error was less than 5%. The refractive index also depends on the temperature of
the fluid, which can be calibrated by the instrument. In order to further reduce these errors due
to temperature differences, [Shin, 2001] allowed the fluids to reach room temperature before the
simulation.
However, because of the simplicity and accurate results, the difference in density for this work
was estimated from a polynomial expression in www.csgnetwork.com/h2odenscalc.html, which
was provided by the University of Michigan and NOAA:
r =r0 +Ac+Bc3=2 +0:00048314c2 (4.1)
where
A = 0:824493 0:0040899T +0:000076438T 2 0:00000082467T 3 +0:0000000053675T 4
B = 0:005724+0:00010227T 0:0000016546T 2
r0 = 1000 1 (T +288:9414)508929:2(T +68:12963) (T 3:9863)2
in which A, B and C are polynomial parameters, r is the density of the denser fluid in kgm 3, T is
the temperature in C and c is the concentration in mg=l. The initial salt concentration was deter-
mined from this relationship, and the resulting density difference was re-evaluated using precision
56
Figure 4.4: Photograph of the leading edge of a GC that illustrates the procedure in which the front
trajectories were measured.
hydrometers. These hydrometers determine the density with an error of approximately 5% or less.
In addition to providing reliable data, experiments were performed relatively quickly.
4.1.3 Front trajectory measurements
One of the most common data recording techniques that has been used in past decades is to
record the trajectory of the GC front [Huppert and Simpson, 1980]. Although the position and the
velocity of the front has been accurately obtained in previous studies, newer technology allows for
a more accurate description. In the present study, digital cameras were used to record the position
of the front at 30 frames per second. This, however, is not the first time that this type of quality of
video camera has been used for GC front trajectories [Adduce et al., 2011]. Measurements of the
trajectory of the front are very important in the evaluation of numerical models.
In order to obtain values for the position of the GC nose at various positions along the channel,
a grid (4 in by 1 in) was drawn on the front of the tank to record the longitudinal position of GC
fronts. The video camera was operated in a way that followed the front of the GC. Because the
57
video camera also recorded the time, the velocity distribution of the GC front for each experiment
was obtained.
4.1.4 MicroADV measurements
For the velocity hydrograph measurements, MicroADV probes sampling at 20 Hz were placed
near the middle of the tank in order to determine the longitudinal velocities for the denser GC and
the ambient fluid. The lower ADV probe was placed low enough to measure the GC at 10.2 cm
from the channel bottom. The upper probe was placed at 33.0 cm from the channel bottom in
order to measure the ambient velocities. Although the ADV probes can measure the velocity in
three-dimensions, the main focus of the data analysis was related to the longitudinal velocities.
In this analysis, it was important to position the MicroADV probes as close to vertical as
possible in order to accurately measure the velocity in the longitudinal direction. The probes were
held in place with a wooden support, and adhesive tape was applied to hold the them in the vertical
direction. The wooden support contains a cantilever section that hangs over the channel so that the
ADV probe can be placed in the center of the tank. The data was recorded at 20 Hz for about 250
seconds or when the GC front reflected back and forth twice.
4.2 One-layer SWE model
As mentioned, the objective was to develop a robust and generalized Boussinesq gravity cur-
rent model for the lock-exchange problem. It was decided to minimize the use of explicit track-
ing of flow features (i.e. discontinuities appearing when f > 0:5). Using the one and two-layer
SWE written in physically conservative format, it was possible to implement linear and nonlinear
shock-capturing schemes within a finite volume framework in an attempt to better describe the
propagation of gravity currents, as presented below.
The one-layer SWE models (equation 4.2) are widely used for GC flowing into a deep am-
bient. In an attempt to establish clarity and simplicity, the one-layer SWE are written in vectorial
58
conservative format:
?~U
?t +
? ~F(U)
?x =
~S(U)
~U =
2
64 h
uh
3
75 ~F(U)=
2
66
4
uh
(uh)2
h +
1
2g
0h2
3
77
5 ~S(U)=
2
640
0
3
75
(4.2)
This mathematical model is solved within the finite volume framework by integrating the
previous PDE system, yielding the following expression that performs the updates on the conserved
variables for each computational cell i at the time step n+1:
~Un+1i =~Uni + Dt
Dx
~
F(U)n+1=2i 1=2 ~F(U)n+1=2i+1=2
where ~F(U)n+1=2i 1=2 are inter-cell fluxes (not to confuse with ~F(U)), which are calculated according
to various numerical schemes formulations. Details of such formulations are presented in [Toro,
2001] and are omitted here for brevity.
A major limitation of one-layer model is the inaccuracy when the ambient flow depth is not
large in comparison with the depth of the current [Ungarish, 2009]. Moreover, the one-layer SWE
cannot capture the slumping characteristics of gravity currents with large depth ratios where the
effects of the ambient velocity are more important. In order to appropriately simulate the advance
of gravity current fronts, SWE models (one and two-layered) require the use of a front condition.
4.3 Two-layer SWE model
For the one-layer model, the velocity of the ambient fluid is neglected in the internal cell
calculations. However, if the fractional depth becomes large enough (particularly when f > 0:5),
59
then the ambient velocity becomes an important parameter. The mathematical model for the two-
layer SWE in a Boussinesq system was formulated in [Rottman and Simpson, 1983]:
?h
?t +
?
?x(uh)= 0
?u
?t +(1 2a)u
?u
?x +(1 b)g
0?h
?x = 0
(4.3)
where
a = hH h b = hH + u
2
g0H
1 hH
2
Although equation 4.3 can produce robust results for certain GC applications, the previous
equation has two limitations. First it is not expressed in terms of physically conserved variables (u
instead of uh). Second, while amenable to be implemented in a finite difference framework, it is
not in the conservative format used in finite volume models (i.e. equation 4.2). Thus, in cases of
flow discontinuities, shocks will be simulated with incorrect velocities. In this work we propose to
further manipulate the previous expression using the continuity equation and the chain rule (shown
below) in order to obtain conservative variables.
?q
?t =
?uh
?t =
?u
?t +
?h
?t
The terms that cannot be presented in the fluxes of the conservative variables ( ~F(U)) are isolated
as source terms, which yields the following expression for the linear momentum conservation
equation:
?uh
?t +
?
?x
(uh)2
h +
1
2g
0h2
!
= 2auh?u?x +bg0h?h?x (4.4)
60
Now, combination of the equation (4.4) with the continuity equation for the dense layer yields
a system of PDEs that is expressed in conservative, divergent format:
?~U
?t +
? ~F(U)
?x =
~S(U)
~U =
2
64 h
uh
3
75 ~F(U)=
2
66
4
uh
(uh)2
h +
1
2g
0h2
3
77
5 ~S(U)=
2
66
4
0
A?u?x +g0hB?h?x
3
77
5
(4.5)
where A and B are written as:
A = 2uh
2
H h B =
(1 h=H)2 + u=pg0H 2
(1 h=H)2 (4.6)
In the above equations, A and B are components of the source terms presented by [Rottman and
Simpson, 1983] written in such a way to demonstrate that the numerator and denominator in both
expressions approach 0 as f!1 and h!H, since the current velocity (u) approaches 0 in such
conditions. For h > 0:95H, special handling of the source terms is required.
A straight forward inspection of equation 4.6 indicates that the terms A and B appearing in the
source term of equation 4.5 become undefined as h=H!1 where u!0, which in turn renders the
numerical solution unstable. For the depth gradient term, this causes no problems except at f = 1;
however,~S = 0 at f = 1. In this scenario B = 1 because u = 0 due to the horizontal free surface.
For the velocity gradient component, the fix is not as simple. Moreover, instabilities are generated
for h=H ranging from 0.95 to 1.0. In order to overcome this dilemma, L?Hopital?s rule is applied
as h=H!1 so that a limiting value for A is formulated:
lim
h / H!1 A = 4uh (4.7)
In the proposed two-layer SWE model, this use of L?Hopital?s rule was applied for f > 0:95.
Moreover, the velocity gradient term is negligible for h=H > 0:9, so the possible side effects are
negated. Using this approach, the two-layer SWE model is able to accurately simulate GC flows of
61
all fractional depth. In this implementation, the parameter A for the velocity gradient component
of the source term is provided in the following expression.
A =
8
>><
>>:
2uh2
H h (h=H < 0:95)
4uh (h=H 0:95)
(4.8)
The source terms are evaluated numerically using a finite difference calculation, and in this
implementation a simple, central difference scheme scheme was used. The two-layer SWE may
be solved by implementing a numerical scheme in a FVM update formulation (equation 2.48).
In summary, the advantage of the proposed mathematical model is the simplicity since the same
numerical scheme is applied everywhere in the solution domain at all time steps, except of course
at the boundary conditions. Thus, the only boundary conditions that are used in the computations
are the physical boundaries constraining the flow and the GC front. There is no explicit tracking
of shocks except at the leading edge of the GC.
4.4 Aspects of the numerical solution
A brief discussion on the strategy for the numerical solution implementation for these hyper-
bolic PDE systems is presented, which include the selected numerical schemes; boundary condi-
tions, including a new formulation to compute flow at the gravity current nose; and aspects related
to the numerical stability and time steps that are required in the proposed model.
4.4.1 Numerical schemes
As mentioned, to perform updates in the conserved variables using finite volume method,
it is necessary to calculate the interface fluxes. While the proposed model has been tested with
several linear and nonlinear schemes, four schemes are selected to perform the comparison in this
work: the first-order accurate, linear FORCE scheme; second-order accurate, linear Lax-Wendroff
62
scheme; and two first-order accurate, nonlinear Roe and HLL schemes. All of these schemes are
widely applied in CFD problems, and a detailed analysis is provided in [Toro, 2001].
Lax-Wendroff (LxW): The linear, two-step LxW scheme was presented in subsection 2.5.2
for the FDM. For the FVM, the expression for this second-order accurate scheme is [Toro,
2001]: 8
>>>
>>>
<
>>>
>>>:
~Un+1=2i+1=2 = 12(~Uni +~Uni+1)+ Dt2Dx( ~F(U)ni ~F(U)ni+1)
~F(U)n+1=2
i+1=2 =
2
66
4
uhn+1=2i+1=2
(uh)2
h +
1
2g
0h2
n+1=2
i+1=2
3
77
5
(4.9)
Modified FORCE: The original FORCE scheme discussed in [Toro, 2001] is an average
50%-50% of the interfacial fluxes from the LxW and LxF schemes. This work follows the
idea from [Vasconcelos and Nunes, 2009] and adds more weight to LxW scheme (i.e. 70%)
so that diffusion is further reduced without a noticeable increase in oscillations at shocks.
This results in the modified FORCE scheme adopted here:
~F(U)n+1=2
i+1=2 =
1
2
0:3 ~F(U)n+1=2i+1=2
LxF +0:7 ~F(U)n+1=2i+1=2
LxW
(4.10)
HLL: The nonlinear HLL scheme has performed well in the more common dam break prob-
lem [Vasconcelos and Nunes, 2009] and should provide accurate and computationally effi-
cient solutions for lock-exchange GC flows. The primary difference between this formula-
tion and other SWE models is that the gravity g is replaced with the reduced gravitational
acceleration g0 in the flux calculations. The procedure for the HLL scheme is provided in
subsection 2.5.2. For the HLL scheme, the expression for the interfacial fluxes ( ~F(U)i+1=2)
63
that is used in the update expression (2.48) is:
~F(U)i+1=2 =
8
>>>
>>><
>>>>
>>:
FL
FHLL = SRFL SLFR +SRSL(UR UL)S
R SL
FR
in which the subscript L denotes the ith cell and the subscript R denotes the cell i+1. The
description of wave speeds (S) and the fluxes (F) is provided in subsection 2.5.2. In the
proposed model, the intermediate or star variables are computed with the Two-Rarefaction
Riemann Solver, which was recommended by [Toro, 2001]. However, similar results were
obtained with the Riemann Solver based on exact depth positivity. Therefore, it is assumed
that the choice of approximate-state Riemann solver for the intermediate variables is not
crucial for these GC flows.
Roe: The nonlinear Roe scheme was presented in subsection 2.5.2. The main difference
in this implementation from other SWE models is that the gravity (g) is replaced with the
reduced gravity (g0) due to Boussinesq gravity currents. The proposed numerical model was
tested with the Glaister [Glaister, 1988] and Roe-Pike [Roe and Pike, 1984] approaches to
compute the Roe averages, and the results were almost identical. For the Roe scheme, the
expression for the interfacial fluxes (~F(~U)i+1=2) that is used in the update expression (2.48)
is:
~F(U)i+1=2 = 1
2
~
F(U)ni + ~F(U)ni+1
12
2
j=1
?ajj?ljj?R(j)
where the wave strengths ( ?aj), eigenvalues ( ?lj) and eigenvectors ( ?R(j)) are provided in subsection
2.5.2.
64
4.4.2 Boundary Conditions: DC and MOC
At the boundaries of the solution domain (either between a moving gravity front and a physical
boundary or between two physical boundaries) numerical schemes cannot be used to determine the
two flow variables U=[h;uh]T and boundary condition calculations are required. At physical walls
the alternatives to solve the problem are generally straightforward. One way is to either enforce
zero flow (u = 0) and solve for the depth h using a relevant characteristic equation [Sturm, 2010];
an alternative would be the use of a reflective boundary condition as presented in [Toro, 2001].
The other boundary condition calculation is at the GC front, and two alternatives are im-
plemented for this moving boundary in the case of the one-layer SWE model: the traditional
characteristic approach (MOC) and a proposed approach based on local continuity and momentum
conservation (DC) (described ahead in this subsection). Both alternatives enforce a constraint on
the depth when the ambient depth H decreases and approaches the lock-exchange problem (e.g.
0:347 H for Benjamin?s front condition) [Klemp et al., 1994] and [Ungarish and Zemach, 2005].
The MOC BC that is only implemented here for the one-layer model involves the combination of a
C+ characteristic equation and the enforcement of a front condition (as explained in section 2.1).
This approach provides two equations that allow for the solution of the two unknowns [hLE;qLE]T at
the gravity current cell LE, which denotes the front or leading edge. Assuming that the velocity of
cell 1 (i.e. the wall boundary) will always be zero, the resulting set of equations are:
ucell:1 +2
p
g0hcell:1 = 2
p
g0hcell:1 = uLE +2
p
g0hLE uLE =
p
g0hLEFr LE (4.11)
which are solved iteratively (for this study, the Bisection method was used). One notices that the
characteristic equation transports the Riemann invariant value u+2c from the upstream wall (cell
1) until the GC nose according to the characteristic velocity: u+c. In addition, the solution for
[hLE;uLE]T allows for the determination of uhLE.
The boundary condition strategy for the two-layered system could also have followed the
same strategy, but the more complex structure of the characteristic equation (as shown in [Rottman
65
and Simpson, 1983]), and consequently Riemann invariants, motivated the development of an
alternative approach. This alternative method explicitly enforces continuity, linear momentum,
the kinematic condition and a front condition at the leading edge of the GC. This front region is
assumed to span over cell LE, a finite volume cell fully occupied by the dense fluid, and over cell
LE+1, which undergoes the advance of the current and is not yet completely filled (see Figure 4.5).
Because the front region exists between two computational cells, this method was named Dual-Cell
(DC), and it requires five equations to solve for the unkown parameters: h LE, uh LE, h LE+1, uh LE+1 and
Dx LE+1. At the upstream front BC cell (i.e. LE), continuity is enforced along with the x-momentum
equation. The continuity equation and a front condition is implemented at the downstream cell (i.e.
LE+1). The continuity equations that were used in the SWE models are provided below:
Continuity
8
>><
>>:
dALE
dt =Dx
dhLE
dt = uh LE-1 uh LE+1 cell: LE
dALE+1
dt =
d(Dx LE+1 hLE+1)
dt = uh LE cell: LE+1
(4.12)
in which A is the cell area (A =DxhLE for cell LE and A =DxLE+1hLE+1 for cell LE+1). The continuity
equation is slightly more complex for cell LE+1 because both the depth hLE+1 and the distance DxLE+1
change with time (see Figure 4.5). The x-momentum equation that was used in the DC boundary
condition is provided from the following expression:
Fx =
cs
?mLE+1uLE+1
cs
?mLEuLE (4.13)
in which ?m = rc uh is the mass flow rate. The implementation of the momentum and continuity
equations for the SWE models is provided below along with the other two equations (i.e. front
66
condition and kinematic condition), which make up the DC boundary condition:
DC
8
>>>>
>>>
>>>>
>>>>
<
>>>
>>>>
>>>
>>>>
>:
I r
c uh n+1LE
uh n+1
LE+1
h n+1LE
uh n+1LE-1
h n+1LE
+F 2 F1 = 0
II h n+1
LE = h
n
LE +
Dt
Dx(uh
n+1
LE-1 uh
n+1
LE+1)
III uh n+1
LE+1 = h
n+1
LE+1
pg0h n+1
LE+1 FrLE+1
IV Dx n+1
LE+1 =Dx
n
LE+1 +Dt
uh n+1LE+1
h n+1LE+1
V uh n+1
LE =
1
Dt
h
Dx n+1LE+1(h n+1LE+1 h nLE+1)+h n+1LE Dx n+1LE+1 Dx nLE+1
i
(4.14)
where n is the time step index, F1 is the upstream hydrostatic force at the interface between cells LE-1
and LE, and F 2 is the downstream hydrostatic force at the interface between cells LE and LE+1. These
five equations represent respectively: 4.14 I the linear momentum balance at cell LE; 4.14 II the
local continuity at cell LE; 4.14 III the enforcement of a front condition relating depth and velocity
at cell LE+1; 4.14 IV the kinematic condition, which updates the front position Dx n+1LE+1, at cell LE+1;
and 4.14 V the local continuity at cell LE+1.
These equations are solved iteratively for the pair of cells at the edge of the computational
domain until the advance of the GC front at cell LE+1 (i.e x n+1LE+1, given by u LE+1Dt) exceeds the
cell size Dx. When this happens, the front of the GC enters a new cell and Dx LE+1 is reset to zero.
This alternative to compute the boundary conditions is shown to be accurate for both one-layer and
two-layer SWE models, and have compared well with experimental data collected in this work as
well with previous investigations.
4.4.3 Stability and time step calculation
Due to the choice of applying explicit schemes for the numerical solution, there is an up-
per limit for the computational time step Dt established by the Courant-Friedrichs-Lewy (CFL)
condition:
67
Figure 4.5: Schematic diagram of a GC that illustrates the cell structure for the DC BC at a given
time step. The cell LE represents the last computational node that is completely filled by the front
of the GC.
Dt Dxmax(u+pg0h) (4.15)
The enforcement of this condition is not sufficient to render the model stable. What was
determined is that in some cases, even when equation 4.15 is satisfied, the conditions become
unstable at the front. When the time step was reduced and the number of time steps to traverse cell
LE+1 increased, stability is obtained in the computations. Thus, a second limit for time step is also
enforced in which the velocity scales with the celerity:
Dt kDt Dxpg0h
LE
(4.16)
in which kDt is a constant factor used in the computations ranging between 0.05 and 0.5, depending
on the initial depth ratio and the density difference. However, stability was typically achieved with
kDt = 0:1. Later in the model development, a change was made to the force calculations of the
boundary conditions so that the Courant condition governed the stability. As a general rule, the
68
relationship that was implemented to compute the time steps is:
Dt min
Dx
max(u+pg0h);kDt
Dxp
g0hLE
(4.17)
69
Chapter 5
Results and Discussion
5.1 Experimental results
Front trajectories
It is convenient to display the results for lock-exchange GC flows in terms of dimensionless
units (see section 2.4.2). Furthermore, the distance is normalized by the initial current length
(x0), the depth by the initial GC depth (h0), the velocities by the initial celerity (pg0h0), etc.
(see [Ungarish, 2009]). This procedure allows for a comparison for various initial conditions
and verifies consistent experimental results. In addition, other researchers can easily compare
there data to the dimensionless results. For the experiments conducted in this study, the initial
conditions and some of the dimensionless results (see section 2.2) are presented in Table 5.1. It is
clear from the results that the inviscid assumption is justified (Rem >> 1;000). The consistency for
FrH demonstrate the consistency in which the experiments were performed.
Table 5.1: Experimental initial conditions and results for each experiment. Values for um were
determined from a digital camera tracking the GC front.
Run x0 h0 r1 g0 um Rem FrH
m m kgm -3 ms -2 ms -1
1 0.762 0.406 1010 0.107 0.098 19952 0.468
2 0.762 0.406 1020 0.202 0.135 27657 0.472
3 0.762 0.406 1028 0.277 0.158 32423 0.473
In Figure 5.1 dimensionless GC front trajectory and velocity results are provided for the 9.14
m tank. After an initial reduction in velocity, a near constant velocity region forms denoting the
slumping stage (see section 1.3). The velocity results for the three experiments are almost identical
70
Figure 5.1: Dimensionless trajectory and velocity for the GC front.
(i.e. steady) except at the beginning of the experiment where the removal of the gate is very impor-
tant and could not be absolutely equal between repetitions. Likewise, the change in position with
time reveals that the velocity is almost constant (actually very slowly decreasing) within the slump-
ing stage. These results suggest that the experiments were conducted with sufficient accuracy to
yield consistent results.
MicroADV results
MicroADV results were obtained in the 9.14 m segment of the channel at one location for all
of the experiments: x = x=h0 = 7:68. The depth (h0) for the experiments was 40.6 cm and the
initial lock length (x0) was 76.2 cm. Thus, the thin layer assumption for the shallow water theory
is valid. For h0=x0 < 1, the mixing between the GC and the ambient fluids is not significant, and
the majority of breaking waves are confined to the interface [Hacker et al., 1996]. The results for
71
the MicroADV probe are displayed in terms of dimensionless units in Figure 5.2. Thus, the data
for all three experiments should fall along the same path.
The results for the MicroADV probe provide a clear representation of the internal velocity
structure for lock-exchange experiments. It is apparent from Figure 5.2 that pressure pulses are
consequences of surface waves generated from the gate removal. An oscillatory pattern forms
from this wave celerity; however, there is a clear interface caused by the arrival of the GC front. As
the time increases, the pressure pulses are reduced. Furthermore, a parabolic decrease in velocity
occurs after the GC front moves past the ADV probe until the end wall reflection generates negative
velocities. The one and two-layer SWE models are able to simulate the GC propagation beyond
the end wall reflection so that the velocity magnitudes are compared.
The negative peaks that occur when t 40 represent the front of the reflected GC from the
far wall so that the GC travels back towards its initial location. After the negative velocity front
arrives at the MicroADV probe, one notices that the velocity quickly changes back to positive
values before a negative curve develops. Thus a dense solitary wave was separated from the rest of
the flow at the initial reflection with the far wall. This solitary wave was generated from a hydraulic
jump, which was diminished due to the stably stratified layer that was left behind the GC nose from
turbulent mixing [Simpson, 1997].
72
Figure 5.2: Velocity hydrograph comparison for 1%, 2%, and 3% salt in the 9.14 m tank.
73
5.2 Numerical model results
The results for the one and two-layer SWE models are illustrated in Figure 5.3 for a wide
range of fractional depth. In Figure 5.3 f = 0:01, f = 0:5 and f = 1:0 for the deep ambient,
critical condition and full depth scenario, respectively. The DC BC was implemented at the front
of the GC with Benjamin?s front condition and 400 discretization nodes. Unless otherwise noted,
kDt = 0:1 for the full depth and critical condition. For the deep ambient problem, kDt = 0:2 in order
for the code to remain stable.
It is clear from figure 5.3 that the two-layer SWE model approaches the one-layer alternative
for the deep ambient scenario. The depth profile and the GC trajectory are almost identical at each
point in time. In the deep ambient scenario, the velocity of the ambient fluid is much smaller in
magnitude thus less momentum is transferred to the GC. However, as the depth ratio begins to
increase, the two SWE models begin to diverge from one another. The resistance of the ambient
fluid becomes more important, and the two-layer model is the better alternative [Ungarish, 2009].
The decrease in accuracy of the one-layer SWE model for smaller ambient depths is due to
the increased importance of the ambient velocity. As the fractional depth increases from the deep
ambient to the full depth lock-exchange problem, the velocity increases for both the one and two-
layer models. Moreover, for the deep ambient scenario, a larger pressure from the ambient fluid is
being exerted on the GC, which corresponds to a greater horizontal velocity.
One of the interesting features in Figure 5.3 is that the GC in the two-layer model increases
in velocity relative to the one-layer model as the depth ratio increases to the critical condition
and then to the full depth lock-exchange problem. The reason for this discrepancy is due to the
increased accuracy of the two-layer model in computing the length of the slumping stage, which
is underestimated with the one-layer model. The two-layer model simulates the reflection with
the upstream boundary with a sharp reflected bore, which eventually reaches the front. For the
one-layer model, this reflection generates a wave, which reaches the GC front much faster, so the
depth of the front and the resulting velocity are decreased quicker than in the two-layer model.
74
Although though the GC front in the two-layer model is constrained because of the more
complex characteristics, the depth and velocity of the front are greater than in the one-layer model
for the majority of the slumping stage. This constraint has been validated by high-resolution mod-
els that suggest the GC front rarely moves faster than the c+ characteristic for Boussinesq fluids
[Ungarish, 2009]. Because the two-layer model computes a longer slumping stage for the shallow
ambient condition, the GC propagates at a faster.
In table 5.2 the sensitivity of the DC boundary condition is analyzed in the context of the
number of discretization nodes. For the one-layer model, the DC boundary condition is compared
to the traditional MOC boundary condition. The results indicate that the DC BC is more accurate
for simulations with smaller discretization nodes. However, both models converge in accuracy at
approximately 400-800 nodes.
In previous two-layer models using the shock-tracking approaches, a smaller number of dis-
cretization nodes was often used to resolve the flow discontinuities. The discontinuous interfaces
for these models are assumed vertical, so the numerical scheme does not need compute this feature
even though post shock oscillations are a potential side effect [Ungarish, 2009]. By comparison, the
shock-capturing approach proposed in this thesis requires more sophisticated numerical schemes
and improved grid discretization in order to accurately reproduce the discontinuous interfaces. A
comparison between these two approaches is presented in Figure 5.11.
Table 5.2: Sensitivity of each modeling approach to the number of discretization nodes. The error
for each model is compared to its computed propagation time at 1,000 nodes. The Roe scheme is
used with B front condition and kDt = 0:1.
Model Parameter Number of nodes
50 100 200 400 800 1000
GC propag. time (s) 77.6 71.1 68.1 68.0 67.2 67.31L: MOC
Continuity error (%) +15.3 +5.5 +1.2 +1.0 -0.1 0
GC propag. time (s) 66.9 65.6 64.6 65.5 65.5 65.81L: DC
Continuity error (%) +1.7 -0.3 -1.8 -0.4 -0.4 0
GC propag. time (s) 58.4 56.2 55.1 55.4 55.4 55.52L: DC
Continuity error (%) +5.3 +1.4 -0.6 -0.1 -0.1 0
75
Figure 5.3: GC propagation for the one and two-layer SWE over a wide range of fractional depth:
f = 0:01, 0:5 and 1:0. The DC BC is used with 400 discretization cells.
76
5.3 Comparison between experimental results and numerical predictions
The assessment of the proposed numerical solutions for the one and two-layer SWE with
regards to the experimental results is done in the following comparisons:
1. Experimental measurements of the GC front trajectory are compared to the results obtained
from the one and two-layer models using the Roe scheme;
2. Experimental results for the overall GC propagation time are compared to numerical predic-
tions with different numerical schemes, front conditions, front BC strategies (for the one-
layer model) and time steps (for the two-layer model). Continuity errors in each alternative
are also compared;
3. Velocity measurements obtained with MicroADV probes for both the denser GC and the
upper ambient fluid are compared to the predicted velocity hydrographs obtained with the
one and two-layer SWE models;
Trajectory results
Experimental results for the trajectory of the GC front in the 9.14 m tank were obtained from
high definition digital cameras, recording at 30 frames per second. The advance of the front was
measured against a grid placed on the acrylic walls, and typical results are presented in Figure
5.4. The GC flow was caused by the removal of a gate located at x0 = 0:762m, which separated
saltwater from freshwater. The depth of both fluids prior to the gate removal was 0.406 m, and the
relative density difference was approximately 2%. As described in [Simpson, 1997], the advance
of the GC front is quasi-steady within the slumping stage, as one notices that the slope of the
trajectory varies only slightly with time.
Figure 5.4 presents numerical predictions for the GC front trajectory obtained with the one-
layer SWE model using the Roe scheme to solve flow in internal cells associated with various
front conditions. The DC front boundary condition equations were used in a 400-cell computa-
tional domain. In the chart, the acronyms ?B?, ?HS? and ?UZ? represent for the front conditions
77
Figure 5.4: Front trajectory results for the one-layer SWE model using Roe scheme and DC BC
with various front conditions and from experiments.
proposed by Benjamin [Benjamin, 1968], Huppert and Simpson [Huppert and Simpson, 1980] and
Ungarish and Zemach [Ungarish and Zemach, 2005], respectively. Results presented in this figure
for the coordinate of the front are normalized by the original position of the gate x0, while the time
was normalized by x0=pg0h0. One notices that the best prediction obtained with the one-layer
model was the one obtained with Benjamin?s front condition followed by Huppert and Simpson.
In general, the tested numerical alternatives were able to capture well the general trend of the
experimental trajectory of the GC front.
A similar comparison, this time obtained with the two layer model is presented in Figure 5.5;
the options for the numerical solution are the same as for the one-layer model described above.
Results are now reversed in that the predictions with Benjamin?s front condition were the worst,
while results obtained with HS and UZ front conditions simulated the experimental measurements
very well.
78
Figure 5.5: Front trajectory results for the two-layer SWE model using Roe scheme and DC bound-
ary conditions with various front conditions and from experiments.
Effects of the boundary conditions and numerical schemes
The effect of the numerical scheme selection, boundary condition solution strategy (for one-
layer model), time step size (for two-layer model) and front condition choice to the accuracy of the
model predictions are compared and summarized in Tables 5.3 and 5.4 for the one and two-layer
SWE models, respectively. The assessment of the model accuracy is measured using three criteria:
1) the difference between the experimental propagation time and the respective numerical predic-
tion; and 2) the continuity (mass balance) errors; and 3) the computational efficiency involved in
the calculations of the GC flow. Numerical results, all obtained with a 400-cell computational do-
main, are compared to the experimental condition (x0 = 76:2cm; h0 = 40:6cm;Dr = 0:02gcm 3),
but results for larger density differences follow the same general trend.
One recalls that the application of the one-layer model is not recommended for the conditions
used in experiments due to the large ratio between the GC depth and the ambient depth, particu-
larly at the initial stages of the flow. This condition may be one explanation why the continuity
errors obtained for this model associated with the traditional characteristic equation are too high,
79
regardless of the front condition, as it may be seen in Table 5.3. For this scenario, the constraint
based on the energy conservation theory (f 0:5), causes an increase in the continuity error for
the characteristic BC. The use of the alternative DC boundary condition equations, however, seems
to address most of the issues with continuity errors. As far as the velocity of the GC front, the best
results were obtained with Benjamin?s front condition for the one-layer model.
In general, the worst numerical scheme with regards to the accuracy of the GC propagation
time and continuity error is the second-order LxW scheme, possibly due to the artificial viscos-
ity that must be introduced in the scheme so that oscillations are controlled. The best numerical
scheme for this condition is the modified FORCE, but both non-linear schemes tested also pro-
duced accurate results. The nonlinear Roe and HLL schemes performed well for all conditions that
were tested. In addition, these nonlinear schemes are more stable for larger density differences and
partial depth releases and perform better for smaller time steps than the linear schemes (Table 5.4).
In Figure 5.6 the results with the Roe and HLL schemes are superimposed, and the depth profile
and front trajectory are almost identical.
The conditions tested in the experiments are more appropriately simulated with the use of a
two-layer SWE formulation even though the front trajectory is also simulated accurately with the
one-layer model. In general continuity errors are all under 1.4% for the two-layer model (Table
5.4), and the largest divergence between the measured and predicted propagation time is 18%.
This error was expected because [Shin et al., 2004] stated that Benjamin?s theory over-predicts the
propagation time by about 20%. As discussed in the Methodology chapter, a boundary condition
for the GC front based on the characteristic equation was not implemented for the two-layer SWE
model. Thus it was decided here to compare the effects of changing the time step size by using
two values for the coefficient kDt presented in equation 4.17.
The best agreement for the two-layer model was observed for the case when Huppert and
Simpson?s front condition was used in terms of the GC front propagation time. With regards to the
continuity error, both HS and UZ front conditions presented continuity errors consistently smaller
than 1% while the results using Benjamin?s front condition were consistently the worst. There was
80
Figure 5.6: Slumping for the one and two-layer SWE models using the HLL and Roe scheme with
two front conditions: Benjamin (B) and Huppert Simpson (HS). x0 = 1, h0 = 1 and e = 2% using
800 computational cells.
81
a slight increase in the continuity errors for the smaller kDt tested. All two-layer numerical schemes
tested presented similarly accurate results, with a slight advantage to the non-linear schemes. The
use of the LxW scheme resulted in instabilities that deemed this shock-capturing model unstable,
and thus it was decided not to present the results here.
Based on the error computed for the propagation time for the two-layer model, which was
almost always greater in magnitude than 6%, an additional front condition was implemented.
[Rottman and Simpson, 1983] proposed an approach based on Benjamin?s condition that allowed
for a calibration constant. From their partial-depth experiments, they calibrated this parameter b
to equal 1.0. However, this term changes for the full depth lock-exchange problem. Therefore,
b is re-calibrated in order to determine a value that satisfied the experiments for this work (This
condition is referred to as MRS). The results are summarized in Table 5.5, and b = 1:215 (see
equation 2.7) provided the best results. For the full depth lock-exchange problem, it is impor-
tant to note that the depth and velocity constraints at the GC front change depending on the front
condition (see subsection 2.4.2). The front constraints for the MRS condition are: hLE 0:414H,
uLE=pg0h0 0:448, respectively. Moreover, the error with experiments for the two-layer SWE
model using the MRS front condition is constistently less than 1%.
According to Table 5.3, the one-layer SWE are very accurate in conjunction with Benjamin?s
theoretical front condition. It is speculated that the reason for these good results is that the errors
associated with interfacial mixing/entrainment balance with the errors associated with neglecting
the ambient velocity. In the actual GC propagation, the ambient velocity transfers momentum to
the GC, which increases the GC velocity. On the other hand, turbulent entrainment, which is most
prevalent at the GC front [Simpson, 1997], halts the advance of the GC. [Ungarish, 2007] similarly
obtained accurate results for the full depth release using the one-layer SWE even though these
equations are least valid for this depth ratio. When the two-layer SWE model is implemented,
this balance of errors is changed because the ambient velocity is taken into account within the
internal cell calculations. Therefore, because turbulent entrainment is neglected, the errors are
greater compared to the one-layer model.
82
Table 5.3: Comparison between lock-exchange experiments and one-layer SWE model predictions
using both BC approaches at the GC front, kDt = 0:1 and a 400-cell solution domain. The time of
propagation for the GC in the experiments was 93.67 s (e = 1:0%) and 67.63 s (e = 2:0%).
Experiment: e = 1:0%
DC BC equations Characteristic BC equations
B LxW FORCE HLL Roe LxW FORCE HLL Roe
GC propag. time (sec) 91.54 89.90 90.21 89.87 92.43 91.74 92.02 93.71
Error w/ experim. (%) -2.3 -4.0 -3.7 -4.1 -1.3 -2.1 -1.8 +0.0
Continuity error (%) -5.0 -2.0 -2.1 -2.4 -8.3 -7.4 -7.7 -8.0
HS LxW FORCE HLL Roe LxW FORCE HLL Roe
GC propag. time (sec) 108.39 107.42 107.49 107.49 109.94 109.66 109.76 109.63
Error w/ experim. (%) +15.7 +14.7 +14.8 +14.8 +17.4 +17.1 +17.2 +17.0
Continuity error (%) -7.4 -2.0 -2.3 -2.3 -12.9 -11.4 -11.8 -10.8
UZ LxW FORCE HLL Roe LxW FORCE HLL Roe
GC propag. time (sec) 122.02 120.09 119.95 120.19 125.18 124.68 124.81 124.56
Error w/ experim. (%) +30.3 +28.2 +28.1 +28.3 +33.6 +33.1 +33.2 +33.0
Continuity error (%) -7.3 -2.0 -2.1 -2.3 -12.6 -11.0 -11.6 -10.7
Experiment: e = 2:0%
DC BC equations Characteristic BC equations
B LxW FORCE HLL Roe LxW FORCE HLL Roe
GC propag. time (sec) 66.55 65.45 65.83 65.43 67.23 66.75 66.93 68.01
Error w/ experim. (%) -1.6 -3.2 -2.7 -3.3 -0.1 -1.3 -1.0 +0.6
Continuity error (%) -4.88 -2.03 -3.47 -2.23 -8.28 -7.43 -7.65 -8.18
HS LxW FORCE HLL Roe LxW FORCE HLL Roe
GC propag. time (sec) 78.82 78.12 79.06 78.14 79.96 79.78 79.83 79.76
Error w/ experim. (%) +16.5 +15.5 +16.9 +15.5 +18.2 +18.0 +18.0 +17.9
Continuity error (%) -7.33 -2.04 -2.42 -2.22 -12.89 -11.51 -11.83 -11.02
UZ LxW FORCE HLL Roe LxW FORCE HLL Roe
GC propag. time (sec) 88.73 87.30 87.56 87.32 91.05 90.68 90.77 90.62
Error w/ experim. (%) +31.2 +29.1 +29.5 +29.1 +34.6 +34.1 +34.2 +34.0
Continuity error (%) -7.27 -2.05 -2.10 -2.23 -12.63 -11.12 -11.58 -10.88
83
Table 5.4: Comparison between lock-exchange experiments and two-layer SWE model predictions
for different time steps and front conditions using the DC boundary condition. The GC time of
propagation in the experiments was 93.67 s (e = 1:0%) and 67.63 s (e = 2:0%).
Experiment: e = 1:0%
kDt = 0:1 kDt = 0:05
B FORCE HLL Roe FORCE HLL Roe
GC propag. time (sec) 76.15 76.15 76.14 76.15 76.15 76.15
Error w/ experim. (%) -18.7 -18.7 -18.7 -18.7 -18.7 -18.7
Contintuity error (%) +1.0 +1.3 +1.3 +1.5 +1.3 +1.4
HS FORCE HLL Roe FORCE HLL Roe
GC propag. time (sec) 100.26 99.65 99.63 99.56 99.68 99.65
Error w/ experim. (%) +7.0 +6.4 +6.4 +6.3 +6.4 +6.4
Contintuity error (%) +0.5 +0.6 +0.6 +0.8 +0.7 +0.7
UZ FORCE HLL Roe FORCE HLL Roe
GC propag. time (sec) 107.48 106.36 106.28 106.28 106.39 106.32
Error w/ experim. (%) +14.7 +13.5 +13.5 +13.5 +13.6 +13.5
Contintuity error (%) +0.5 +0.6 +0.6 +0.8 +0.7 +0.7
Experiment: e = 2:0%
kDt = 0:1 kDt = 0:05
B FORCE HLL Roe FORCE HLL Roe
GC propag. time (sec) 55.38 55.38 55.38 55.39 55.39 55.38
Error w/ experim. (%) -18.1 -18.1 -18.1 -18.1 -18.1 -18.1
Contintuity error (%) +1.4 +1.3 +1.3 +1.6 +1.4 +1.4
HS FORCE HLL Roe FORCE HLL Roe
GC propag. time (sec) 72.46 72.46 72.19 71.50 72.48 72.22
Error w/ experim. (%) +7.1 +7.1 +6.7 +0.57 +7.2 +6.8
Contintuity error (%) +0.7 +0.6 +0.7 +1.0 +0.7 +0.8
UZ FORCE HLL Roe FORCE HLL Roe
GC propag. time (sec) 77.36 77.32 76.77 75.58 77.36 76.80
Error w/ experim. (%) +14.4 +14.3 +13.5 +11.8 +14.4 +13.6
Contintuity error (%) +0.7 +0.6 7 +1.0 +0.7 +0.8
84
Table 5.5: Two-layer SWE results using various values for b for the RS front condition with 400-
cell solution domain. When b = 1:215, this front condition is referred to as MRS. The GC time of
propagation in the experiments was 93.67 s (e = 1:0%) and 67.63 s (e = 2:0%).
Experiment: e = 1%
b based on equation 2.7 1.0 (RS) 1.1 1.2 1.215 1.3 p2 (B)
GC propag. time (s) 131.35 112.81 95.73 93.33 83.64 76.16
Error w/ experim. (%) +40.2 +20.4 +2.2 +0.4 -10.7 -18.7
Continuity error (%) +0.51 +0.55 +0.70 +0.73 +0.91 +1.29
Experiment: e = 2%
b based on equation 2.7 1.0 (RS) 1.1 1.2 1.215 1.3 p2 (B)
GC propag. time (s) 95.53 82.08 69.45 67.88 60.83 55.40
Error w/ experim. (%) +41.3 +21.4 +3.0 -0.4 -10.1 -18.1
Continuity error (%) +0.51 +0.55 +0.70 +0.73 +0.91 +1.29
The computational time associated with the for SWE model alternatives is an important fea-
ture that is compared in Table 5.6. This comparison is made between the one and two-layer models
and between the DC and MOC boundary conditions for the one-layer model. The first section of
Table 5.6 displays the computational time for the GC to travel the full length of the tank with the
following initial conditions: x0 = 1m, Ltank = 10m, h0 = 1m and e = 2%. The second section
of Table 5.6 compares the different modeling approaches. The results indicate that the DC BC is
consistently more computationally efficient than the MOC BC. It is apparent from these results
that the computational time is strongly dependent on the boundary condition solution strategy for
the GC front. Also, the source terms that are present in the two-layer model have an important
effect on the computational time (20%-50%, see Table 5.6). The two-layer model combined with
the MOC BC and explicit tracking of the discontinuities is expected to provide even worse com-
putational times as the one-layer model for comparable domain discretizations. Therefore, the
two-layer SWE shock-capturing model proposed in this thesis is potentially attractive in terms of
computational efficiency, fairly low continuity errors and above all conceptual simplicity and ease
of implementation.
It was found that the selection of the numerical scheme has a minor importance on the overall
accuracy. It was determined that the nonlinear schemes are more stable for partial depth releases
85
Table 5.6: Computational time comparison between the one-layer (1L) and two-layer (2L) SWE
models and between the DC and MOC front BC strategies for various discretization sizes.
Num. of cells 100 200 400 800 2000
1L: MOC (s) 1.37 5.20 26.93 88.78 179.00
1L: DC (s) 1.02 2.07 2.62 9.58 36.32
2L: DC (s) 1.40 2.77 3.17 13.28 52.57
1L: MOC/DC 1.3 2.5 10.3 9.3 4.9
DC: 2L/1L 1.4 1.3 1.2 1.4 1.5
and for larger density differences. There was severe diffusion for the first-order Lax-Friedrich
scheme. In addition, the second-order accurate linear schemes tested were unable to remain sta-
ble for the two-layer SWE model in the context of the full depth lock-exchange problem due to
excessive oscillations that form during the reflected bore. In previous models these second-order
accurate schemes performed adequately because a shock-fitting procedure was used. In shock-
capturing models, the choice of numerical scheme is more important since the discontinuity is not
explicitly tracked.
The most important factors in the simulation of GC propagation compared to lock-exchange
experiments were: 1) selection of the appropriate mathematical model for the numerical solu-
tion (best was two layer SWE); 2) selection of the appropriate front condition (modified RS for
two-layer model); and 3) the selection of the nose BC solution strategy (most efficient was DC
equations).
Velocity measurements and predictions
A comparison between the predicted and measured velocities at an intermediate point in the
tank are presented in Figure 5.7. Velocity values are normalized by pg0h0 while time values
are normalized by xo=pg0h0 (the dimensionless variables are denoted with asterisk). Velocity
measurements were performed with the MicroADV probes sampling at a frequency of 20 Hz. The
oscillations observed at the beginning of the velocity measurements were caused by the removal
of the gate separating fresh and salt water. Suddenly, at about t = 15 the jump in the velocity
86
values corresponds to the arrival of the front at the ADV station, at x = x=h0 = 7:68. Velocity
measurements indicate a gradual decline in the velocity after an initial peak that was caused by
the arrival of the GC front. Measurements were continued after the gravity current advanced the
entire length of the tank in order to observe the reflections. As the reflected current reached the
MicroADV probes, the velocity becomes slightly negative, which was an expected outcome.
Figure 5.7a) presents a comparison between the experimental measurements of the velocity
with the predictions from the one-layer SWE model using the Roe scheme, DC boundary condi-
tions and either Benjamin?s or Huppert and Simpson?s front conditions. These front conditions
represent the theoretical implementation and the most accurate empirical condition that was tested
for the full depth release. One notices good agreement between both simulations and the experi-
mental data, with Benjamin?s front more accurately predicting both the arrival of the front and the
reflected wave. The predicted velocity decay during the motion of the gravity current is milder
than the measurements. Moreover, for the one-layer model the use an empirical front condition is
not as important since the theoretical formulation in [Benjamin, 1968] is provides accurate results.
A comparison between the predicted GC velocity using the two-layer SWE model and the
experiments is presented in Figure 5.7b). As in the previous case, Roe scheme and DC boundary
conditions were used along with Benjamin?s and Huppert and Simpson?s front conditions. From
recent studies [Ungarish and Zemach, 2005], it is well known that the HS front condition performs
well for the full depth lock-release, which is reinforced in 5.7b). On the other hand, Benjamin?s
condition overestimates the GC front velocity by about 20% [Shin et al., 2004]. Therefore, one
would expect that the HS front condition would provide a better estimate for the GC velocity
structure. However, the magnitude (specifically at the velocity peaks) is better represented by
Benjamin?s condition even though the front advance is better approximated by the HS condition.
These results suggest that there is intense mixing and entrainment at the GC front (stated in
[Simpson, 1997]), which decrease the velocities from their theoretical value. Although there have
been many explanations regarding the inaccuracy of Benjamin?s theory [Klemp et al., 1994] and
[Shin et al., 2004], Figure 5.7b) provides a clear understanding of the source of these problems in
87
the context of the velocity structure. The results in 5.7 indicate that the best agreement is observed
with the HS front condition and the two-layer model.
The results with the MRS front condition are presented in Figure 5.8 for the two-layer model to
compare with two experiments: e = 1% and e = 2%. Although the HS front condition performed
well in predicting the GC velocities, the results suggest that the MRS front condition performs best
at predicting the arrival of both GC fronts: initial propagation and first reflection. It performs a
little better than the HS front condition at predicting the peak velocities at the front but not as well
as the B front condition. All things considered for the two-layer model, the MRS front condition is
the best method tested for the Boussinesq, full depth lock-exchange problem.
In Figure 5.9 the experimental results (Dr = 3:0%) for the upper MicroADV probe are com-
pared to the two-layer SWE model predictions. The MRS and B front conditions are used with
200 computational cells. As seen in Figure 5.7b), the B front condition overestimates the arrival
of both fronts: initial propagation and first reflection. The peak velocities compare well to both
front conditions tested but are better represented by the MRS front condition. In addition, the MRS
condition performs very well at predicting both GC front locations: t 15 and t 40. After
the arrival of the GC front (t 20 30 in Figure 5.9), the ambient velocities computed with the
two-layer SWE model underpredict the MicroADV results.
88
Figure 5.7: Velocity hydrographs with B and HS front conditions and 200 computational cells.
The experimental data (e = 2:0%) is provided from MicroADV devices in the slumping stage:
x = 7:68.
89
Figure 5.8: Velocity hydrographs with the MRS front condition and 200 computational cells. The
experimental data is provided from MicroADV devices at x = 7:68.
90
Figure 5.9: Comparison between measured ambient velocities and two-layer SWE predictions
using B and HS front conditions with 200 computational cells. The experimental data (e = 3:0%)
is provided from MicroADV devices at x = 7:70.
5.4 Comparison with previous models
This subsection compares the results of the proposed model to the ones obtained by previous
models. Two models were selected for this comparison, both of these using two-layer formulations:
Rottman and Simpson [Rottman and Simpson, 1983] and Ungarish and Zemach [Ungarish and
Zemach, 2005]. The former model successfully simulated gravity current flows for partial depth
releases (for f 0:5), but had limitations in representing the back-propagation of the wave during
the slumping stage of gravity currents. The second model presented by [Ungarish and Zemach,
2005] is able to successfully simulate flows for cases when f > 0:5, but requires explicit tracking
of the upstream moving hydraulic jump and the reflected bore. As it is shown below, the proposed
model is able to simulate flows for any value of f without needing to perform explicit tracking of
any flow feature except the GC front.
91
Figure 5.10 presents a comparison between the proposed model and the Rottman and Simpson
model results for the case when f =0:5. The results presented are the profiles of the GC at different
times along the simulation. All profiles are presented in a non-dimensional fashion, with depth
normalized by h0, length by x0 and time by x0=pg0h0. [Rottman and Simpson, 1983] presented
results using Benjamin?s front condition for a qualitative analysis; they obtained better agreement
with experiments using their front condition. The proposed model results were obtained using
Roe?s scheme and Benjamin?s front condition. One notices an excellent agreement between both
models.
A more challenging simulation is the case when f > 0:5 due to the formation of the backward
moving feature resembling a discontinuity in the depth. A simulation for the case when f = 1 was
performed with the proposed model in the conditions used in the work by [Ungarish and Zemach,
2005], and the results of both model?s GC profiles are presented in Figure 5.11. Two different
modeling alternatives are presented for this comparison: Figure 5.11a) presents results obtained
with the Roe scheme; and Figure 5.11b) presents results obtained with the Modified FORCE
scheme. In both cases the HS front condition was used with the limiting value hLE = 0:427H
at the GC front, following [Klemp et al., 1994]. Both modeling alternatives were successful in
reproducing the model results presented by [Ungarish and Zemach, 2005]. Results obtained with
the MFORCE scheme were slightly sharper (i.e better at simulating the vertical interface of shocks)
than the results obtained with Roe scheme. For the Roe scheme, which performed similarly to the
HLL scheme, Glaister?s approach [Glaister, 1988] was used to compute the Roe averages; however,
the same results were obtained with the Roe-Pike approach (see [Toro, 2001]). One also notices
that the few high frequency spurious oscillations presented in [Ungarish and Zemach, 2005] model
(which uses LxW scheme) are not present in the proposed model results.
In the previous discussion, the two-layer SWE model results were presented for three different
initial depth ratios: f = 0, f = 0:5 and f = 1. In the following comparison, front conditions are
tested for a larger number depth ratios in the context of the dimensionless front velocity. The results
are presented in Figure 5.12 and were obtained using Roe?s scheme and three front conditions: B,
92
Figure 5.10: Comparison of GC profiles at various time steps between the proposed two-layer
model and [Rottman and Simpson, 1983] model for the critical condition: f = 0:5.
93
Figure 5.11: Dimensionless slumping comparison between the proposed model and [Ungarish and
Zemach, 2005] model using Roe scheme and modified FORCE scheme.
94
Figure 5.12: Dimensionless front velocity results for the two-layer model over a large range of
fractional depth.
HS and UZ. One notices that as the GC propagates into shallower ambient depths (thus larger
f values) the velocity normalized by pg0h0 decreases. Results obtained with Benjamin?s front
simulation were the largest, while results obtained with Huppert and Simpson?s front condition
had an intermediate maximum at about f = 0:10. The results obtained with the UZ front yielded
the smallest GC front propagation velocity. These results are compared with the numerical results
presented by [Ungarish and Zemach, 2005] and with the experimental results by [Rottman and
Simpson, 1983], and the best agreement with experimental results from the proposed model was
observed for the UZ front condition.
95
A final comparison is made between the Environmental Fluid Dynamics code (EFDC) and
the two SWE models (with the HS front condition), which is presented in [Hatcher et al., 2012].
If EFDC is used to compare the SWE models (Figure 5.13), the depth profile is much better ap-
proximated by the two-layer model. However, both SWE models are able to accurately simulate
the position of the front if the correct front condition is used. The positive results between the
two-layer SWE models and EFDC indicate the possibilities of the thin layer models to GC flows.
In addition, the position of the nose and the slope of the reflected bore are very similar for both
approaches. It is important to note that the SWE model simulates a smaller magnitude for the
length of the nose region and is unable to model the mixing between the two fluids.
96
Figure 5.13: Dimensional slumping comparison between the SWE and EFDC for a) t=1.8sec and
b) t=39sec. The discretization size for the SWE models is 800 cells.
97
Chapter 6
Conclusions and future developments
In this thesis Boussinesq lock-exchange flows were analyzed experimentally and numerically
with the one and two-layer SWE models. These types of flows are common in canals, rivers, estu-
aries, etc. and play an important role in mixing and contaminant transport. Thus, computationally
efficient and accurate numerical models are important to determine the trajectory and velocities of
the GC. Although the model was constructed for Boussinesq, 1-D lock-exchange flows, the pro-
posed BC at the GC front can be applied to a larger range of flows (e.g./ non-Boussinesq). In order
to test the numerical model, MicroADV probes were used to measure the velocities of the bottom
current and the ambient fluid, and a high definition video camera was used to track the GC front.
This study presents a modeling framework based on the FVM to simulate GC flows using the
one and two-layer SWE. The motivation was the overall success in the application of the FVM
in the context of free surface flow simulation. While the simulation of GC flows with numerical
models has been performed for almost three decades, some conditions associated with GC flows
are more challenging to be represented numerically. This is particularly true for the case when GC
propagate into small ambient depths (f 0:5), as in such cases a backward propagating wave with
a large depth discontinuity poses difficulties to SWE-based models. To date, such flows have been
simulated with models that perform explicit tracking of this hydraulic jump and the reflected bore,
which results in additional model complexity.
The proposed model aims to present a simple, albeit robust, way to perform GC flow com-
putation to any ambient depth. Previous two-layer SWE models required different approaches
depending on the initial fractional depth. In the proposed two-layer SWE model, two main contri-
butions are presented. The first is a numerical model in which a single shock-capturing framework
98
is used to simulate GC flows for any initial depth ratio. The second contribution is a new al-
ternative to compute the GC front boundary condition enforcing explicitly mass and momentum
conservation, referred to as Dual-Cell (DC). This method was generally able to increase model
accuracy, particularly for smaller discretization sizes while being substantially more computation-
ally efficient than the alternative based on the method of characteristics (MOC). For some of the
discretization sizes that were tested, the DC BC was more than 10 times faster than the MOC.
The characteristic backward moving flow discontinuity cannot be properly simulated with
a single-layer model as is discussed in [Ungarish, 2009]. Results obtained with the proposed
two-layer SWE model are much better qualitatively and were more successful in reproducing ex-
perimental results with the MicroADV probes. Moreover, the two-layer SWE compared very well
with previous numerical models both for small and large ambient depths. The development of
this numerical solution involved the application of the SWE equation in conservative variables and
the separation of some terms of the expression derived by Rottman and Simpson [Rottman and
Simpson, 1983] to be handled as source terms. Such source terms required special treatment to
overcome difficulties as f 1, and the use of L?Hopital rule in the two-layer model was successful
for this purpose.
The velocity gradient component of the two-layer source term vector required the use of
L?Hopital?s rule for h=H > 0:95. This velocity gradient term was not important in the source
term calculations for large fractional depth, so the effect on the numerical scheme was negligible.
However, if L?Hopital?s rule is not used the numerical code was unstable. On the other hand, the
depth gradient term played the most important role in establishing the discontinuity instead of the
depression wave for large values of f. When this discontinuity (instead of the depression wave)
reflects off of the upstream boundary, the reflection causes the GC nose to form in case of the
two-layer model. Once the upstream reflection occurs, the source terms are only important at the
interface of the reflected bore (where the depth gradient is large), which eventually reaches the
front of the GC.
99
The resulting one and two-layer SWE models were tested with different front conditions, dif-
ferent strategies to compute the nose boundary condition and for various fractional depth. The
most important factors as far as accuracy is concerned were: 1) proper selection of the mathe-
matical model; 2) selection of the front condition expression; and 3) selection of nose boundary
condition calculation. Unlike the initial expectation, three out of four numerical models performed
comparably well in the simulations. The Lax-Wendroff scheme (even with the use of artificial
viscosity) produced the worst results due to the oscillations observed and the code crashing at the
upstream reflection.
The use of a modified FORCE scheme and the non-linear Roe and HLL schemes was gen-
erally very successful in reproducing the experimental results collected in this investigation, par-
ticularly when associated with the two-layer SWE formulation and the modified version of the RS
front condition (MRS) that was presented in this thesis. The empirical coefficient b in this expres-
sion was calibrated for the full depth lock-exchange experiments, and the best agreements with
experiments was b = 1:215. This value is clearly larger than the value provided in [Rottman and
Simpson, 1983] (b = 1:0) for partial depth releases. It is reiterated that an important feature of
the proposed model is the ability to simulate the upstream moving discontinuity observed during
the slumping stage without the need of explicit tracking. Moreover, the objective was to develop
a shock-capturing two-layer SWE model that can simulate GC flows without the restriction of the
ambient depth.
Although there has been extensive research for GC flows in recent decades, further investi-
gations are required for better understanding. The effects of entrainment and complex geometries
are still poorly understood in the context of the SWE models. For laboratory experiments, it is
difficult to accurately measure the GC depth. For this study, there was an attempt to utilize pres-
sure transducers to overcome this problem; however, surface waves smear the results indicating
that larger-scale GC flows are required if transducers are to be used. In future studies additional
CFD models will be tested (e.g. OpenFOAM or FLOW3D). In addition, the proposed two-layer
modeling approach should be extended to non-Boussinesq problems and axisymmetric GC flows.
100
Bibliography
[Abbott, 1961] Abbott, M. B. (1961). On the spreading of one fluid over another, part II: the wave
front. La Houille Blanche, 6:827?836.
[Adduce et al., 2011] Adduce, C., Sciortino, G., and Proietti, S. (2011). Gravity currents produced
by lock-exchange: experiments and simulations with a two layer shallow-water model with
entrainment. J. Hydraulic Eng.
[Barr, 1967] Barr, D. I. H. (1967). Densimetric exchange flow in rectangular channels. iii large
scale experiments. Houille Blanche, 6:619?631.
[Baxter, 2000] Baxter, P. J. (2000). Impacts of eruptions on human health. In Encyclopedia of
Volcanoes, pages 1035?1043. Academic Press.
[Benjamin, 1968] Benjamin, T. B. (1968). Gravity currents and related phenomena. J. Fluid
Mech., 31:209?248.
[Birman et al., 2005] Birman, V. K., Martin, J. E., and Meiburg, E. (2005). The non-boussinesq
lock-exchange problem. part 2. high-resolution simulations. J. Fluid Mech., 537:125?144.
[Blong, 2000] Blong, R. (2000). Volcanic hazards and risk management. In Encyclopedia of Vol-
canoes (e. H. Sigurdsson), pages 1215?1227. Academic Press.
[Bonnecase et al., 1992] Bonnecase, R. T., Huppert, H. E., and Lister, J. R. (1992). Particle-driven
gravity currents. Journal of Fluid Mechanics, 250:339?369.
[Bonnecaze et al., 1993] Bonnecaze, R. T., Huppert, H. E., and Lister, J. R. (1993). Particle-driven
gravity currents. J. Fluid Mech.
[Britter and Simpson, 1978] Britter, R. E. and Simpson, J. E. (1978). Experiments on the dynam-
ics of a gravity current head. J. Fluid Mech., 88:233?240.
[Chaudhry, 2008] Chaudhry, M. H. (2008). Open-Channel Flow. Springer, second edition.
[Chen et al., 2005] Chen, Q., Zhao, H., Hu, K., and Douglass, S. (2005). Prediction of wind waves
in a shallow estuary. J. Waterw., Port, Coastal, Ocean Eng., pages 137?148.
[Cunge et al., 1980] Cunge, J. A., Holly, F. M. J., and Verwey, A. (1980). Practical Aspects of
Computational River Hydraulics. Pitman Publishing Ltd., London, UK.
[Dalziel, 1993] Dalziel, S. B. (1993). Rayleigh-taylor instability: experiments with image analy-
sis. Dyn. Atmos. Oceans, 20:127?153.
101
[Ellison and Turner, 1959] Ellison, T. H. and Turner, J. S. (1959). Turbulent entrainment in strati-
fied fluids. J. Fluid Mech., 6:259?274.
[Fannelop and Waldman, 1971] Fannelop, T. K. and Waldman, G. D. (1971). Dynamics of oil
slicks. AIAA, 10:506?510.
[Firoozabadi et al., 2010] Firoozabadi, B., Afshin, H., and Bagerpour, A. (2010). Experimental
investigation of turbulence specifications of turbidity currents. J. Applied Fluid Mech., 3(1):63?
73.
[Garcia and Kahawita, 1986] Garcia, R. and Kahawita, R. A. (1986). Numerical solutions of the
de saint venant equations with the mac-cormack finite difference scheme. International Journal
for numerical methods in fluids, 6:259?274.
[Gerber et al., 2011] Gerber, G., Diedericks, G., and Basson, G. R. (2011). Paricle image ve-
locimetry measurements and numerical modeling of a saline density current. J. Hydraulic Eng.,
pages 333?342.
[Glaister, 1988] Glaister, P. (1988). Approximate riemann solutions of the shallow water equa-
tions. J. Hydraulic Res., 26.
[Godunov, 1959] Godunov, S. K. (1959). Finite difference method for numerical computation of
discontinous solution of the equations of fluid dynamics. Matematicheskii Sbornik, 47(3):271?
306.
[Hacker et al., 1996] Hacker, J., Linden, P. F., and Dalziel, S. B. (1996). Mixing in lock-release
gravity currents. Dynamics of Atmospheres and Oceans, 24:183?195.
[Hatcher et al., 2012] Hatcher, T. M., KC, M., Vasconcelos, J. G., and Fang, X. (2012). A compar-
ison between numerical modeling approaches for the simulation of gravity currents. Submitted
to the EWRI conference 1/2012.
[H ?artel et al., 2000] H ?artel, C., Meiburg, E., and Necker, F. (2000). Analysis and direct numerical
simulation of the flow at a gravity-current head. part1. flow topology and front speed for slip
and no-slip boundaries. J. Fluid Mech., 418:189?212.
[Hogg et al., 1999] Hogg, A. J., Ungarish, M., and Huppert, H. E. (1999). A box model for non-
entraining, suspension-driven gravity surges on horizontal surfaces. Eur. J. Mech. B - Fluids,
19:139?165.
[Hoult, 1972] Hoult, D. P. (1972). Oil spreading on the sea. Anna. Rev. Fluid Mech., 4:341?368.
[Huppert, 2006] Huppert, H. E. (2006). Gravity currents: a personal perspective. J. Fluid Mech.,
554:299?322.
[Huppert and Simpson, 1980] Huppert, H. E. and Simpson, J. E. (1980). The slumping of gravity
currents. J. of Fluid Mech., 99:785?799.
[Keller and Chyou, 1991] Keller, J. J. and Chyou, Y. P. (1991). On the hydraulic lock-exchange
problem. Z. Angew. Math. Phys., 42.
102
[Keulegan, 1957] Keulegan, G. H. (1957). An experimental study of the motion of saline water
from locks into fresh water channels. Nat. Bur. Stand. Rept. Technical Report 5168.
[Klemp et al., 1994] Klemp, J. B., Rotunno, R., and Skamarock, W. C. (1994). On the dynamics
of gravity currents in a channel. J. Fluid Mech., 269:169?198.
[Kranenburg, 1978] Kranenburg, C. (1978). Internal fronts in two layer flow. J. Hydraul. Div.
ASCE, 104(HY10):1449?1453.
[LeVeque, 1992] LeVeque, R. J. (1992). Numerical Mehods for Conservation Laws. Birkhauser
Verlag.
[Lowe et al., 2005] Lowe, R. J., Rottman, J. W., and Linden, P. F. (2005). The non-boussinesq
lock-exchange problem. part 1. theory and experiments. J. Fluid Mech., 537:101?124.
[Necker et al., 2002] Necker, F., Hartel, C., Kleiser, L., and Meiburg, E. (2002). High-resolution
simulations of particle-driven gravity currents. Intl J. Multiphase Flow, 28:279?300.
[Ozgokmen et al., 2009] Ozgokmen, T. M., Lliescu, T., and Fischer, P. F. (2009). Large eddy
simulation of stratified mixing in a three-dimensional lock-exchange system. Ocean Modell.,
26:134?155.
[Roe and Pike, 1984] Roe, P. L. and Pike, J. (1984). Efficient construction and utilisation of ap-
proximate riemann solutions. Computing Methods in Applied Scence and Engineering.
[Rottman and Simpson, 1983] Rottman, J. W. and Simpson, J. E. (1983). Gravity currents pro-
duced by instantaneous releases of a heavy fluid in a rectangular channel. J. Fluid Mech.,
135:95?110.
[Ruo and Chen, 2007] Ruo, A. and Chen, F. (2007). Modified shallow water equations for inviscid
gravity currents. Physical Review E, 75.
[Schoklitsch, 1917] Schoklitsch, A. (1917). Uber dambruchwellen. Sitzungsber. Akad Wis-
senschaf, 126(1489).
[Shin, 2001] Shin, J. (2001). Colliding gravity currents. PhD thesis, University of Cambridge.
[Shin et al., 2004] Shin, J., Dalziel, S., and Linden, P. F. (2004). Gravity currents produced by
lock exchange. J. of Fluid Mech., 521:1?34.
[Simpson, 1997] Simpson, J. E. (1997). Gravity Currents In the Environment and the Laboratory.
Cambridge University Press.
[Simpson and Britter, 1979] Simpson, J. E. and Britter, R. E. (1979). The dynamics of the head of
a gravity current advancing over a horizontal surface. J. Fluid Mech., 94:477?495.
[Slim, 2006] Slim, A. C. (2006). High Reynolds number gravity currents. PhD thesis, University
of Cambridge.
103
[Stommel and Farmer, 1952] Stommel, H. and Farmer, H. G. (1952). Abrupt change in width in
two-layer open channel flow. J. Mar. Res, 11:205?214.
[Sturm, 2010] Sturm, T. W. (2010). Open Channel Hydraulics. McGraw-Hill, second edition.
[Timmermans et al., 2001] Timmermans, M. E., Lister, J. R., and Huppert, H. E. (2001). Com-
pressible particle-driven gravity currents. J. Fluid Mech., 445:305?325.
[Toro, 2001] Toro, E. F. (2001). Shock-Capturing Methods for Free-Surface Shallow Flows. John
Wiley and Sons.
[Turner, 1973] Turner, S. J. (1973). Buoyancy Effects in Fluids. Cambridge Unveristy Press.
[Ungarish, 2007] Ungarish, M. (2007). A shallow-water model for high-reynolds-number gravity
currents for a wide range of density differences and fractional depths. J. Fluid Mech., 579:373?
382.
[Ungarish, 2009] Ungarish, M. (2009). An Introduction to Gravity Currents and Intrusions. Chap-
man and Hall/CRC, first edition.
[Ungarish and Huppert, 2004] Ungarish, M. and Huppert, H. E. (2004). On gravity currents prop-
agating at the base of a stratified ambient: effects of geometrical constraints and rotation. J.
Fluid Mech., 521(69-104).
[Ungarish and Zemach, 2005] Ungarish, M. and Zemach, T. (2005). On the slumping of high
reynolds number gravity currents in two-dimensional and axisymmetric configurations. Euro-
pean J. Mech. B/Fluids, 24:71?90.
[Vasconcelos and Nunes, 2009] Vasconcelos, J. G. and Nunes, A. M. S. (2009). Evaluating nu-
merical schemes for the shallow water equations applied to dam break flows. In ABRH-Brazil,
editor, Proc. XVIII Brazilian Symposium in Water Resources, Campo Grande, Brazil.
[Vasconcelos et al., 2009] Vasconcelos, J. G., Wright, S. J., and Roe, P. (2009). Numerical oscil-
lations in pipe-filling bore predictions by shock-capturing models. J. Hydraulic Eng., 135(296).
[von Karman, 1940] von Karman, T. (1940). The engineer grapples with nonlinear problems. Bull.
M. Math Soc., 46:615?683.
[Wright et al., 1990] Wright, S. J., Kim, Y., and Buhler, J. (1990). Density current propagation
in flowing receiving fluid. Technical report, Int. Conf. on Physical Modeling of Transport and
Dispersion. Boston, MA. ASCE, New York, pp. 12A19-12A24.
[Wright and Paez-Rivadeneira, 1996] Wright, S. J. and Paez-Rivadeneira, D. (1996). Density in-
trusions with large relative thickness. Dynamics of Atmospheres and Oceans, 24:129?137.
[Zoppou and Roberts, 1999] Zoppou, C. and Roberts, S. (1999). Catastophic collapse of water
supply reservoirs in urban areas. J. Hydraulic Eng., 125(7):686?695.
[Zoppou and Roberts, 2003] Zoppou, C. and Roberts, S. (2003). Explicit schemes for dam-break
simulations. J. Hydraulic Eng., 129:11?34.
104