PHYSICAL AND NUMERICAL MODELING OF
BUOYANT GROUNDWATER PLUMES
Except where reference is made to the work of others, the work described in this thesis is
my own or was done in collaboration with my advisory committee. This thesis does not
include proprietary or classified information.
Linzy Kay Brakefield
Certificate of Approval:
Joel Melville
Professor
Civil Engineering
T. Prabhakar Clement, Chair
Professor
Civil Engineering
Elena Abarca
Fulbright Fellow
Civil Engineering
Jacob Dane
Professor
Agronomy and Soils
Joe F. Pittman
Interim Dean
Graduate School
PHYSICAL AND NUMERICAL MODELING OF
BUOYANT GROUNDWATER PLUMES
Linzy Kay Brakefield
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
May 10, 2008
PHYSICAL AND NUMERICAL MODELING OF
BUOYANT GROUNDWATER PLUMES
Linzy Kay Brakefield
Permission is granted to Auburn University to make copies of this thesis at its discretion,
upon request of individuals or institutions and at their expense. The author reserves all
publication rights.
Signature of Author
Date of Graduation
iii
iv
VITA
Linzy Kay Brakefield, daughter of Brenda Jensen and James Brakefield, was born
July 21
st
, 1979, in Birmingham, Alabama. She graduated from Jefferson County
International Baccalaureate School, formerly known as Shades Valley Resource Learning
Center, with an I.B. diploma in 1997. She attended Maryville College in Maryville,
Tennessee and graduated cum laude with a Bachelor of Arts degree in Mathematics in
May, 2000. After completing some coursework towards a Master of Science degree in
Mathematics at Montana State University, she moved to Auburn, Alabama where she
worked for three years for a local construction company, Shannon, Strobel & Weaver,
Constructors & Engineers, Inc. In fall of 2005, she began coursework to pursue her M.S.
degree in Civil Engineering under the guidance of T. Prabhakar Clement. In June of
2006, she entered Graduate School at Auburn University. She married Rohit Raj
Goswami, son of Chanchal and the late Shuchindra Goswami, on September 12
th
, 2007.
v
THESIS ABSTRACT
PHYSICAL AND NUMERICAL MODELING OF
BUOYANT GROUNDWATER PLUMES
Linzy Kay Brakefield
Master of Science, May 10, 2007
(B.A., Maryville College, 2000)
107 Typed Pages
Directed by T. Prabhakar Clement
In coastal states, the injection of treated, relatively fresh, wastewater into deep
saline aquifers offers a disposal alternative to ocean outfalls and discharge directly into
local waterways. The density contrast causes upward buoyant movement of the
wastewater plume during and after injection. Since some wastewater treatment plants
inject more than 100 MGD of this treated wastewater, it is important to be able to
determine the fate and transport rates of the plume. In this study, both physical and
numerical modeling were undertaken to investigate and understand buoyant plume
behavior and transport.
vi
Physical models using a 2D tank filled with glass beads were carried out under
different ambient density scenarios. The experiments consisted of injection of a
freshwater pulse-source bubble in an initially static system with no ambient flow.
Using the finite-difference numerical code SEAWAT v.4 to simulate variable
density flow, the experiments were numerically modeled and compared with the physical
model results. Due to the sensitivity of this problem to spatial resolution, results from
three different grid sizes were compared to determine a reasonable compromise between
computing times and numerical accuracy. Furthermore, a comparison of advection
solvers was undertaken to identify the best solver to use for this specific problem. From
these scenarios, the Method of Characteristics (MOC) advection solver with the fine
resolution grid (0.1 cm x 0.1 cm x 2.7 cm cells) resulted in a simulation that was in good
agreement with the physical experiments. This model was determined to be the base-case
problem for further sensitivity analysis. The finite element based numerical code
SUTRA_MS was also used for intercode comparison with SEAWAT v.4.
Dimensionless analysis of the flow and transport governing equations was
undertaken to determine important physical problem parameters and a characteristic
plume travel time. From the derived dimensionless numbers, it was hypothesized that
density, hydraulic conductivity and dispersivity should all play important roles in this
problem. A parameter sensitivity analysis on these parameters was performed. The
parameter sensitivity investigation involved a quantitative comparison based on moment
analysis. It was determined that the problem was most sensitive to density contrast and
hydraulic conductivity in regards to vertical velocity rates. Dispersivity changes played
an important role in affecting fingering.
vii
ACKNOWLEDGEMENTS
The author would like to thank her mother, Brenda Jensen, for her continued
support throughout her education as well as her brother, Wesley Brakefield, and the rest
of her family. Thanks are also due to Mr. C. Hadley Weaver, Jr. and his wife, Kristin, for
their support and encouragement early in the author?s graduate career. The author would
also like to thank her committee members for their ongoing guidance and support,
especially Dr. T. Prabhakar Clement who gave the author several opportunities to
develop professionally and encouragement to succeed. Thanks are extended to Dr.
Christian Langevin and Ms. Alyssa Dausman at the USGS for their professional and
personal encouragement. Thanks are due to the author?s colleagues: Rohit Goswami,
Matthew Hogan, Justin McDonald, Gautham Jeppu, Venkatraman Srinivasan, Anjani
Kumar, Vijay Loganathan, Che An Kuo, Julia Bower, Douglas Kilgour, Jagadish and
Shyamsundar. Thanks are also extended to the author?s friends for their encouraging
words through the years: Kristy Sharpton, Laura Graham, Chasity Simpson, Tonya
Hollis, Wendy Cantrell, Noemi Caridad Mendez-Sanchez, Elena Abarca, Matt Hogan and
Alyssa Dausman. Thanks are also due to the author?s ?best friend? Bridger and his
unending loyalty and patience. The author would also like to thank her husband Rohit R.
Goswami for his guidance, patience and love. This thesis is dedicated to the author and
husband?s unborn daughter.
viii
Style manual or journal used Auburn University Graduate School Guide to Preparation of
Master?s Thesis, Water Resources Research Journal
Computer Software used Microsoft Windows XP, Microsoft Office XP Professional
(Microsoft Word, Microsoft Excel), SEAWAT v. 4, SUTRA_MS, Groundwater Vistas
5.0, Endnote 9.0, MATLAB, Intel Visual FORTRAN, Surfer 8, Modelviewer 1.1,
UltraEdit 32
ix
TABLE OF CONTENTS
Chapter
1. Introduction 1
1.1 Background 1
1.2 Research Objectives 2
1.3 Organization of Thesis 4
2. Literature Review 5
2.1 Introduction 5
2.2 Review of Analytical Modeling of Variable-density Plumes 6
2.3 Review of Physical Modeling of Variable-density Plumes 7
2.4 Review of Numerical Modeling of Variable-density Plumes 12
3. Physical Experiments 15
3.1 Introduction 15
3.2 General Experimental Set-up 16
3.3 Description of the Experiments 17
3.3.1 Experiment A ? Single Injection in Density-stratified System 18
3.3.2 Experiment B ? Single Injection in Confined System 19
3.3.3 Experiment C ? Multi-injection in Confined System 19
3.4 Experimental Results and Discussion 20
3.5 Need for Numerical Modeling 25
4. Numerical Modeling of Experimental Results 27
4.1 Objective 27
4.2 Numerical Modeling of Density-dependent Flow and Transport 27
4.2.1 Governing Equations 27
4.2.2 Introduction to SEAWAT 29
4.2.3 Introduction to SUTRA 30
4.3 Numerical Model Description 30
x
4.3.1 Initial Conditions 31
4.3.2 Boundary Conditions 32
4.3.3 Injection Well 33
4.3.4 SUTRA Model Comparison 33
4.4 Numerical Model Results 34
4.5 Grid Discretization and Solver Sensitivity 37
4.6 Method of Characteristics (MOC) 41
4.7 An Investigation of Mass Balance Closure in the Numerical Model 41
4.8 Need for Detailed Parameter Sensitivity Analysis 44
5. Parameter Sensitivity Analysis 46
5.1 Introduction 46
5.2 Dimensional Analysis 46
5.3 Methodology for Qualitative Analysis 51
5.4 Methodology for Quantitative Analysis 52
5.4.1 Quantifying the Plume Using Vertical Bulk Velocity and Spatial 53
Variances
5.4.2 Quantitative Method Testing ? Advection Solver Comparison 55
5.5 Parameter Sensitivity Results 63
5.5.1 Application of Moment Analysis to Base-case Simulation Results 63
5.5.2 Density Contrast Sensitivity Results 66
5.5.3 Hydraulic Conductivity Sensitivity Results 71
5.5.4 Dispersivity Sensitivity Results 79
5.5.4.1 Dispersivity Sensitivity ? Investigation I 80
5.5.4.2 Dispersivity Sensitivity ? Investigation II 84
6. Summary and Conclusions 89
7. Future Work 92
References 95
1
CHAPTER 1
INTRODUCTION
1.1 Background
According to the U. S. Geological Survey, in the year 2000, approximately 26
percent of the freshwater used in the U.S. was derived from groundwater sources [2006].
Also, more than 98 percent of the world?s potable water supply comes from groundwater
[Fetter, 1988]. Therefore, it is of utmost importance that underground sources of
drinking water (USDW) be protected from possible contamination, especially in urban
areas where there is higher wastewater production. In urban coastal areas, for example in
Southeast Florida, the population produces enormous amounts of wastewater ? over 500
mgd ? which must be properly disposed of [Muniz, et al., 2005]. Historically, this
wastewater, after completing proper treatment in a waste water treatment facility, was
discharged to surface water bodies or to the ocean. However, due to escalating concerns
over possible unknown effects to local ecologies surrounding discharge areas and the
increased cost associated with regulation and monitoring of ocean outfalls, a growing
number of areas are taking advantage of the option of deep-well injection as a means of
wastewater disposal.
As part of the Underground Injection Control Program (UIC), the EPA has
classified wells into five major types. The first type ? Class I Injection Wells ? is used
2
for injection of hazardous and non-hazardous wastes, liquids and municipal wastewater
beneath the lowermost USDW. As of December 12
th
, 2007, an inventory of 549 Class I
Injection Wells exist in the United States [EPA, 2007]. The popularity of these types of
injection wells can be measured by the fact that they are now used in 19 states. In
Florida, these wells have been used for over 30 years with the typical injection rate for a
well being approximately 18 mgd. In Southeast Florida, the total daily injection rate for
all wells is about 265 mgd [Muniz, et al., 2005]. Due to the number of these wells and
their associated high injection volumes, it is important to be able to understand the fate
and transport of the waste plumes originating from these wells.
1.2 Research Objectives
Since the injected water undergoes primary and secondary treatment at a
wastewater treatment facility, its density is close to that of freshwater (~1.000 g/cc). This
is in contrast to the ambient density of approximately 1.025 g/cc, typically found in saline
Floridian aquifers which are used for waste disposal. This density contrast causes
buoyant upward movement of the plume during and after injection. Also, wastewater is
typically injected under continuous injection conditions in the presence of an ambient
flow field in the receiving aquifer. However, in this study, we only focus on scenarios
involving pulse-source plumes under static regional flow conditions. Future work in this
area should involve analysis of buoyant plumes under ambient flow conditions and/or
with continuous injection.
The objective of this study was to investigate the fate and transport of pulse-
source buoyant freshwater contaminant plumes in saltwater aquifers under hydrostatic
conditions. This was accomplished using both physical and numerical modeling
techniques. The interest was to study not only the upwards movement of the plume but
also the plume behavior near the overlying confining unit. Two types of confining units
were investigated: first, a density stratification set-up with freshwater overlying a saline
layer, and second, a physical impermeable confining unit overlying a fully saline system.
It was also of interest to explore the stratification of subsequently injected buoyant plume
slugs underneath the confining unit.
The numerical code SEAWAT was used to model the experimental data. A
detailed sensitivity analysis was completed to study the sensitivity of the system to
various parameters. In the sensitivity analysis, the first set of simulations considered grid
and advection solver sensitivity. Secondly, a dimensional analysis was carried out to
determine the parameters important in the problem. A chosen set of these parameters
were then perturbed by ? 50% to determine their effects on vertical advection, fingering,
and the fate of the plume. This comparison was initially carried out qualitatively by
contrasting model outputs against one another. Later, a quantitative analysis, using the
method of moments, was completed to better understand the effects of parameter changes
on concentration variations in the plume.
3
4
1.3 Organization of Thesis
The goal of this thesis was to gain a better understanding of the fate and transport
dynamics of a buoyant groundwater contaminant plume. Chapter 2 presents a literature
review covering physical, numerical and analytical studies of buoyant plumes. Chapter 3
describes the physical model experiments of injection into a density-stratified aquifer and
a confined aquifer. Results from a multi-injection experiment are also discussed.
Chapter 4 introduces the governing equations and discusses the numerical modeling
investigation. Chapter 5 includes a dimensional analysis where the important problem
parameters are identified. Furthermore, both a qualitative and quantitative analysis of the
parameter sensitivity results are summarized in this chapter. The overall discussion and
summary of results are presented in Chapter 6. Finally, Chapter 7 discusses future
research that should be extended from this body of work.
5
CHAPTER 2
LITERATURE REVIEW
2.1 Introduction
According to Anderson and Woessner [1992] a model ?is any device that
represents an approximation of a field situation.? Models can be divided into two types:
mathematical models and physical models. Physical models, in general, consist of
laboratory experiments that allow visualization and quantification of both groundwater
flow and solute transport by directly simulating these phenomena. Mathematical models,
on the other hand, replicate these physical processes through mathematical governing
equations, boundary conditions, and initial conditions. This may be done in one of two
ways ? by analytical solutions or numerical approximations of the governing equations
[Anderson and Woessner, 1992].
Analytical models yield a true solution to a given problem, although this usually
comes at the expense of one or more simplifying assumptions. Numerical models, on the
other hand, can generally solve more complex problems, albeit at the expense of
computational errors [Chapra and Canale, 2002]. The popularity and usefulness of
numerical modeling have increased over the past 30 years. This is because model
6
simulation results can ?lead to scientifically defensible answers to specific questions
about ground water systems [McDonald and Reilly, 2007].?
2.2 Review of Analytical Modeling of Variable-density Plumes
Bear and Jacobs [1965] published a paper involving analytical solutions for the
movement of water bodies injected into confined aquifers through wells. They
investigated two scenarios: one of steady injection and one of non-steady injection.
However, they assumed density and viscosity effects were negligible and that the two
fluids were immiscible.
Paschke and Hoopes [1984] investigated analytical solutions of negatively
buoyant (or sinking) plumes in groundwater. These solutions described the movement,
concentration and boundaries of plumes in an ambient flow field. The aquifer conditions
were assumed to be homogeneous and isotropic. Hydrodynamic dispersion was
considered in this problem, along with buoyancy induced flow and ambient flow
[Paschke and Hoopes, 1984]. However, their analysis was limited to dense sinking
plumes..
Buoyant plumes and jets have been investigated in the atmosphere and in surface
water bodies [Fischer, et al., 1979]. However, groundwater scenarios are not the same as
those in air or open water since in these cases dispersion through porous media is
involved. Also, flow is generally assumed to be laminar in groundwater, whereas
turbulence controls atmospheric and open water flows.
7
More recently, Degan et al. [2003] published analytical results related to buoyant
plumes above a heat source, which is an analogous problem to buoyant contaminant
plumes. However, their goal was deriving analytical solutions for an anisotropic
medium. This paper focused on effects of anisotropic parameters and the Rayleigh
number on plume transport.
2.3 Review of Physical Modeling of Variable-density Plumes
Physical modeling was conducted by Paschke and Hoopes as verification for their
analytical solutions (discussed in the previous section). They used a sand packed
laboratory model to test the analytical model for the plume boundary, transport and
concentration. Results from the physical model were shown to correspond well to those
of the analytical model [Paschke and Hoopes, 1984]. However, as stated before, their
investigation involved the study of dense sinking plumes in groundwater instead of
buoyant plumes.
Physical model studies of the behavior of dense sinking plumes were also carried
out by Schincariol and Schwartz [1990], Traylor [1991] and Oostrom et al. [1992a;
1992b]. A large portion of the investigation by Schincariol and Schwartz involved study
of sinking dense plume dynamics in heterogeneous porous media, which is beyond the
scope of this work. We were only interested in buoyant plume dynamics in isotropic,
homogeneous systems. Also, their work involved effects from an ambient flow field and
with a continuous injection of the contaminant where as our work involved static tank
8
conditions and a pulse-source plume. Schincariol and Schwartz investigated effects of
density contrast changes on plume instabilities and the resulting extent of vertical mixing.
Traylor [Traylor, 1991] conducted laboratory experiments of dense sinking pulse-
source plumes in a two-dimensional tank. He also investigated stability of the plume
interface using a simplified analytical approach. The laboratory experiments were
conducted under hydrostatic conditions. He investigated plume transport and stability
sensitivity with respect to changes in density contrast and injection volume. His results
show that vertical velocities of the plume were sensitive to density contrast changes, but
not to changes in volume. However, for unchanging density contrast, he states that there
is a plume volume below which instabilities do not occur.
Oostrom et al. investigated the fate, instabilities, and transport of dense sinking
plumes in a laboratory setting [1992a; 1992b]. They investigated the effects of saturated
hydraulic conductivity, horizontal Darcy velocity, density contrast between ambient and
injected waters and the leakage rate on the stability of plumes. Their experiments were
carried out in homogeneous unconfined aquifer models. A dimensionless analysis was
also undertaken and two important dimensionless numbers were formulated. The
determination of instabilities was seen to be important because unstable plumes can
contaminate much larger regions of the domain than stable plumes. Oostrom et al. also
estimated dispersivity values and indicated that little mixing occurred in the transverse
direction. Although these results were important for understanding dense sinking plumes,
such results may not be applicable to buoyant rising plumes.
9
Hogan [2006] undertook physical modeling of different dense sinking plume
scenarios. This work involved physical modeling using glass bead porous media under
different types of scenarios typical of saltwater inundation due to a tsunami-type event.
The fate and movement of dense saltwater plumes were investigated under three different
set-ups. The first was to mimic a tsunami injection into a well and along the beach face
of a shallow coastal aquifer and was reported by Illangasekare et al. [2006]. The second
was representative of only beach face infiltration (without a shallow well). The third was
illustrative of the instance of saltwater being deposited further inland by a tsunami. The
physical modeling of these different set-ups helped in understanding the dynamics of
plumes injected into coastal aquifers during tsunami-type events.
Other than in Florida, deep well injection has also been used in the state of Hawaii
for handling excess wastewater. In 1977, Peterson and his coworkers from the University
of Hawaii Water Resources Research Center, published a paper investigating buoyant
wastewater plume migration in stratified groundwater bodies [Peterson, et al., 1977].
The groundwater system they investigated consisted of a saline layer overlain by a
freshwater lens with a transition zone in between. The experiment was carried out in a
sand-packed hydraulic model. Blue dye was mixed with the freshwater in the injection
well and the saltwater layer was dyed with a green dye to help visualize the different
layers. The important problem parameters they investigated were the injection rate,
injection depth, length of injection well screening, the density contrast between ambient
and injected waters, and the ambient groundwater flux. None of these parameters were
determined to play a major role in controlling the fate of the contaminant plume. The
10
injected plumes were observed to always break through the transition zone and move into
the overlying freshwater layer. They also found that little mixing between the ambient
saline water and injected water was occurring since very little saltwater entrained into the
plume. With injection in the saltwater zone, however, more saltwater entrainment
occurred than for injection in the transition zone, although it was determined to occur
along the plume?s outer margins. Concentrations of the plume were measured through
sampling ports along one wall of the tank. However, the tank dimensions (1.8 m ? 0.9 m
? 1.2 m) cannot be considered two dimensional although the assumption was made that
what was visualized along the tank boundary was representative of what was occurring
for the plume over the entire thickness. The experiments were carried out under both
static and ambient groundwater flow conditions.
Under static conditions, it was determined that the only parameter exhibiting any
significant control over the plume shape and transport was that of the depth of injection
relative to the transition zone. This parameter however showed little effect over the
ultimate fate of the plume. Injection occurred both in the saltwater zone and in the
transition zone. Plumes injected in the saltwater zone tended to form ?roughly
hemicylindrical column[s]? rising until reaching the freshwater zone where the plume
began to expand outward laterally [Peterson, et al., 1977]. Plumes injected in the
transition zone, however, were markedly different in that they developed no initial
column, but began to immediately spread outward and upward.
Under dynamic conditions, once the plumes penetrated the freshwater layer, their
vertical growth slowed as the plumes moved with the ambient flow field. This was
11
different than the results seen under static conditions where the vertical movement of the
plume was almost unrestricted. It was determined with the experiments under dynamic
conditions that a critical (though unidentified) relationship exists between the ambient
flow rate and the waste injection rate [Peterson, et al., 1977]. Peterson et al. also carried
out a physical model investigation of the same scenarios and parameter sensitivities with
a Hele-shaw cell for comparison to their sand-box model [Peterson, et al., 1978]. Similar
results were found.
Due to the dimensions of the sand box used in the study by Peterson et al, the
flow cannot be assumed to be two-dimensional. Also, their study focused on continuous
injection into a layered system of freshwater overlying saltwater with a transition zone in
between. We were interested in investigating plume development and transport not only
with an overlying freshwater and transition zone, but also with an overlying confining
layer, which is more typical of actual field scenarios at deep-well injection sites in the
state of Florida. Furthermore, rigorous numerical modeling has not been completed for
any buoyant plume study.
A more recent example of physical modeling of buoyant groundwater plumes was
presented by Richardson et al [2004]. This work involved one and two dimensional
physical models and field tests. They investigated the effects of ambient water salinity
on plume migration. According to the results, they suggested that increased background
salinities caused an increase in pore water velocities and dispersion in the system. Also,
as the distance from the injection point increased, the pore water velocities decreased.
However, the focus of this paper was to check the important parameters related to the
12
validity of Rhodamine dye as a tracer in saline systems. Rhodamine dye was found to be
a useful dye for qualitative purposes of determining preferential pathways and general
flow directions. However, due to its complex adsorbtion/desorption relationship in
organically rich waste waters and waters with background salinity, it was not
recommended for quantitative analysis.
2.4 Review of Numerical Modeling of Variable-density Plumes
Numerical modeling of dense sinking plumes has been undertaken by various
researchers. Schincariol et al. [1994] investigated the development of instabilities due to
density effects. Their physical modeling work was discussed in the previous section.
They determined that numerical errors can cause instabilities to form. However, these
errors are difficult to control or predict. Therefore, they used perturbations in the
concentration boundary to generate instabilities. They concluded that increasing
dispersion may cause a decrease in instabilities.
Liu and Dane also undertook numerical modeling of dense sinking plumes [Liu
and Dane, 1997]. Their investigation focused on instabilities in a three dimensional
system. They discuss three possibilities for causes of instabilities: numerical errors (as
also reported by Schincariol et al. [1994]), random small scale perturbations in the
permeability field, and width of the porous media. In a three dimensional flow domain,
plumes can exhibit greater instability in the transverse cross-sectional area than in the
longitudinal direction. Their research confirms conclusions from Schincariol regarding
effect of dispersion on instabilities.
13
Hogan [2006] also completed a parameter sensitivity and stability analysis of the
inland tsunami infiltration physical model set-up described in the previous section. This
sensitivity analysis was undertaken to investigate effects on the plume behavior and
stability caused by changes in the ambient flow field, injection density, injection loading
rate, and dispersivities. The plume behavior and stability was found to be sensitive to all
these parameters. However, the plumes were least sensitive to changes in the injection
loading rate. The parameter sensitivity and stability analysis was completed using
SEAWAT.
Dorgarten and Tsang [1991] numerically studied density-driven transport in
sloping aquifers. This involved numerical modeling of both dense sinking and buoyant
rising plumes. They used a two-dimensional finite element code for their numerical
modeling. Their analysis attempted to quantify the risk of upward movement of waste
water injected into saline systems. It was determined that the slope of the aquifer could
play a significant role in the transport of these injected waters. They suggest a thorough
analysis should be undertaken before and after injection of these wastes including
determining the aquifer characteristics and density contrast.
Numerical modeling of deep-well injection of industrial wastewater at a field site
in Louisiana was conducted by Jin et al [1996]. They used the finite difference based
numerical model HST3D for simulating their scenarios. They reported only minor
differences between effects of deep and shallow wells with regards to the fate and
transport of the injected plume in the aquifer. Due to the large scale of their problem, the
numerical grid used was very coarse in order to minimize the computational time.
14
Furthermore, their results are relevant only to their particular problem parameters since
no general conclusions were made.
A more recent investigation of numerical modeling of buoyant groundwater
plumes was carried out by Maliva et al. [2007] in South Florida. Using the variable-
density code SEAWAT, they investigated the impact of vertical hydraulic conductivity
on buoyant plume migration rates. These plumes were derived from treated wastewater
injected into saline formations. In the presence or absence of fracturing, which can create
vertical conduits for faster flow, the mean residence times still remain large enough to
allow natural degradation to reduce effluent contaminant levels. This work only
investigated the effects of variations in vertical hydraulic conductivity.
15
CHAPTER 3
PHYSICAL EXPERIMENTS
3.1 Introduction
We designed physical models to replicate problems related to buoyant
groundwater plumes under controlled laboratory conditions. This allowed for easy
visualization of real-world scenarios. Furthermore, the physical experimental results can
be used later for comparison against numerical modeling results for code validation. In
this study, three physical experiments were carried out and they will be referred to as
experiments A, B, and C. Experiment A involved injection of a freshwater slug into a
static, density-stratified system of saline water overlain by freshwater. Experiment B was
a simpler set-up involving a freshwater slug injection into a static, fully saline, confined
aquifer. Experiment C involved multiple injections of freshwater slugs into a fully saline,
confined aquifer to compare the effects of previous injections on plume dynamics. The
same two-dimensional tank was used in all three injection experiments.
Experiment A was conducted to mimic a real-world scenario where deep well
injection is used as a means of disposal of treated wastewater. In Florida, for example,
these injection wells can reach to depths of over 3000 m to dispose of the wastewater.
Since the wastewater has a density similar to that of freshwater and is being injected into
an aquifer with density close to seawater density, a buoyancy force is created which plays
a large role in controlling the ultimate fate of the plume. The saltwater aquifers used for
disposal can generally be overlain by freshwater layers which may or may not be used as
a potable drinking water source. This leads to understandable concerns over the safety of
this source of drinking water. Generally, one or several confining layers may separate the
injected aquifer from the overlying freshwater aquifer. Using the experiment A set-up,
we investigated how the plume would act in the absence of a confining layer.
Experiment B was designed to investigate the plume behavior in a confined
system. The tank was filled with porous medium to the top. The open air boundary was
used to conceptually model a non-leaky, impermeable confining unit. An assumption
was made in this experiment that the air-water boundary acts as a confining unit.
The goal of experiment C was to model the effects of previous injections on any
subsequent freshwater plume transport behavior. At most wastewater treatment plants
that utilize deep-well injection, the injection is carried out on and off over long spans of
time. Therefore, the purpose of experiment C was to investigate scenarios involving
multiple plumes. It was our interest to study the effects of previous injections on the
freshwater plume movement and to visualize how multiple plumes stratify beneath the
confining layer.
3.2 General Experimental Set-up
A two-dimensional flow tank (53 cm ? 30.5 cm ? 2.7 cm) constructed of 6 mm
thick Plexiglas walls was used as the physical model. The tank consisted of two lateral
16
17
head chambers and a central porous media chamber, as can be seen in Figure 1. A
careful packing of 1 mm diameter glass beads was performed under saturated conditions
to avoid air entrapment. The beads in the tank were tamped during packing to minimize
layering. Separating the lateral chambers from the porous media were two fine screens of
mesh size US #16 to prevent loss of porous media into the head chambers. Two holes
(approximate diameter of 1 cm), drilled at different heights, were used to inject the
freshwater (see Figures 1 and 2). The holes were fitted with a rubber stopper with an 18
gauge needle in which the injection syringe could be fastened.
The saline water was prepared in a clean bucket using NaCl and de-ionized water.
The density was measured using a ASTM 11H hydrometer. Three food coloring dyes
were used: green (FD&C Yellow # 5 and Blue # 1), blue (FD&C Blue # 1) and red
(FD&C Red # 40) for plume visualization. Food coloring, which is highly soluble in
water and causes negligible effects on density, has been previously used as a tracer in
laboratory experiments [Goswami and Clement, 2007; Pringle, et al., 2002].
Digital images of the experiments were captured using a standard 6 mega-pixel
CCD digital camera and the times of these pictures were recorded.
3.3 Description of the Experiments
Three different experiments were performed as discussed previously. The relevant
parameters used in these physical experiments are summarized in Table 1.
3.3.1 Experiment A ? Single Injection in Density-stratified System
The first experiment (experiment A) represents the injection of freshwater into a
density stratified aquifer. A volume of 24 mL of blue-dyed freshwater was injected by a
syringe over a period of 14 seconds, into a layered system of green-dyed saltwater
overlain by colorless freshwater (Figure 1). Much care was taken in ensuring the
injection was at a continuous rate. The length of injection time was measured using a
stopwatch.
30.5 cm
Saltwater Layer
Head Chambers
Freshwater Layer
53 cm
26.5 cm
15 cm
12.5 cm
Injection Needle
Figure 1 ? Conceptual model of layered experiment A
Table 1 ? Parameters for physical model experiments
Experiment A Experiment B Experiment C
Static head 25.5 cm 29 cm 28.4 cm
?
fresh
1.000 g/mL N/A g/mL N/A g/mL
?
salt
1.025 g/mL 1.020 g/mL 1.030 g/mL
?
injectate
1.000 g/mL 1.000 g/mL 1.000 g/mL
Volume of Injection 24 mL 60 mL 60 (3) mL
Injection Rate 1.71 mL/sec 1.46 mL/sec 2.00 mL/sec
18
3.3.2 Experiment B ? Single Injection in Confined System
The second experiment (experiment B) represents buoyant plume movement in a
confined saltwater aquifer. A 60 mL volume of red-dyed freshwater was injected by
syringe into colorless saltwater over 41 seconds. In this experiment, the second injection
point ? located 5 cm from the bottom of the tank ? was utilized. This allowed us to better
visualize the entire circular plume as it moved through the system (Figure 2).
Saltwater
53 cm
Porous Media Chamber
27 cm
5 cm
Head
29
cm
Injection Needle
Figure 2? Conceptual model of single injection confined experiment B
3.3.3 Experiment C ? Multi-injection in Confined System
Experiment C involved three consecutive injections of freshwater. All three
injections occurred in the same tank as used in previous experiments under fully saline
confined conditions. The ambient water density was higher in this experiment than in
experiment B. The first part of this experiment involved injection of 60 mL of red-dyed
19
20
freshwater into colorless saltwater. The plume was allowed to migrate vertically to the
top of the tank. One hour after the previous injection, a second injection of 60 mL of
green-dyed freshwater was injected into the same port. One hour after the second
injection, a third and final injection of freshwater occurred. Similar to the first injection,
this freshwater was dyed red and a 60 mL volume was injected. Each injection occurred
over a 30 second time-frame.
In all three experiments, the injection rates given in Table 1 are estimated values
as each injectate volume was injected over a measured amount of time (~14 seconds (A),
~41 seconds (B) and ~30 seconds (C)). These injection rates are assumed to be
continuous over these given time frames.
3.4 Experimental Results and Discussion
Results from experiment A are shown in Figure 3. As shown in the figure, the
freshwater plume developed three major fingers as it migrated upwards. These fingers
were presumably caused by small scale porous media heterogeneities and instabilities
caused by dense fluid sinking into the rising light fluid. The three main fingers appeared
to rise at similar rates until reaching the overlying freshwater layer. However, the middle
finger advected slightly faster than the outer two fingers initially. In general, faster
advection of the middle finger was observed in most of the physical experiments. As the
plume rose, processes of mechanical dispersion and diffusion caused mixing between the
plume and the ambient saltwater. This mixing occurred along the boundary between the
plume and the saline water. As fingering became more complex, the surface area of the
21
plume in contact with the ambient water increased, thereby increasing mixing and salt
entrainment into the plume. Saltwater entrainment into the plume can be visualized by a
dilute concentration of the blue dye tracer being left behind in a ?tail? as the plume rose
(see Figure 3). As mixing increased, the density contrast between the plume and the
ambient water decreased which lessened the uplifting buoyant force. With decreasing
buoyancy, vertical velocities diminished and the ratio between advective and dispersive
forces lessened. This cycle continued as the vertical velocity decreased, experiencing a
rapid decline when the plume came into contact with the overlying freshwater layer.
When the boundary was reached, as seen in Figure 3, some of the freshwater injectate
entered the bottom of the freshwater layer (by displacement and/or mixing) while the rest
of the injectate spread out laterally along the boundary. The entrainment of saltwater into
the plume caused the average density of the plume to increase as it rose. Therefore, when
the plume reached the freshwater interface, its average density was slightly more than
that of the freshwater layer. Consequently, the behavior of the freshwater layer was
similar to that of a confining unit in that vertical transport was considerably reduced.
a b c
d e f
g
Figure 3 ? Experiment A results: [times are from end of injection (seconds)] (a) 0,
(b) 90, (c) 210, (d) 330, (e) 510, (f) 750 and (g) 1410
The results from experiment B are shown in Figure 4. Similar to results from
experiment A, the plume in experiment B also rose vertically and fingered until it reached
the boundary. In this case, the zero pressure boundary was the top of the tank open to the
atmosphere which simulates a ?confining unit?. The fingering between plumes in
Figures 3 and 4 is similar early on in the experiments as slight instabilities form due to
small-scale heterogeneities and the density difference between the injected and ambient
water. Also, at later times, the initial instabilities develop into distinct fingers, or
pathways for flow. Upon reaching the upper boundary both plumes exhibited similar
behavior in that vertical advection decreased and the plume began to spread out laterally
along the boundary. It can also be observed that the complex fingering pattern evolved
into a cone-shaped formation at the upper boundary.
22
Figure 4 ? Experiment B results: [times are from end of injection (seconds)] (a) 27,
(b) 369, (c) 685, (d) 1385, (e) 2631
a b c
d e
One main difference that can be noted by comparing the results from these two
experiments is the variations that occur when the plume reaches the overlying boundary.
Even though both plumes begin an outward lateral migration after reaching the boundary,
the plume underneath the confining unit (experiment B) spreads laterally considerably
more as time elapses. This is due to the inherent difference in the overlying boundaries
in each experiment. In experiment A, when the freshwater plume comes into contact with
the overlying boundary, a small part of the freshest part of the plume pushes up into the
freshwater layer displacing some of the overlying water and mixing with it. Therefore,
the plume?s lateral migration rate is slow. In contrast, in experiment B, the freshwater
plume cannot continue its vertical path once in contact with the confining layer, thereby
causing immediate outward expansion. This is the major difference observed between
the two plumes (Figures 3 and 4).
The results from experiment C are shown in Figure 5. It can be seen in this figure
that when a plume mixes with the ambient saltwater as it moves upwards, it leaves a
small tail of slightly fresher water behind. This local change in the ambient density may
23
affect the upward movement of subsequent plumes. Pictures in the left-most column
show results for the first injected plume. These results look very similar to results from
experiment B, the main difference being the travel time of the plume. Experiment C
traveled vertically faster than the plume in experiment B. Three fingers developed with
similar vertical velocity rates. A tail was observed as the plume moved upward. After
reaching the overlying confining interface, the plume spread laterally to form a classic
conical shape.
Figure 5 ? Results of multi-injection confined experiment C
24
The second and third injections are shown in the central and right columns of
Figure 5, respectively. In general, in all of the injections, similar upward behavior is
25
observed, exhibiting similar mean travel times. However, small changes are observed in
the shape of the plume. In the second plume, the middle finger is advanced with respect
to the lateral fingers and, also, with respect to the fingers in the previous injection. This is
caused by the tail of the previous injection. Density is slightly lower in the tail than in the
surroundings as saltwater entrainment as occurred here. Therefore, less mixing of the
next injected plume with the ambient water occurred in the tail zone. Also, the new
plume should have met less resistance to flow in this area as the water is slightly lighter.
Therefore, in the tail area, the buoyant force is reduced more slowly than in other areas of
the plume. As a result, the central part of the plume moved faster following the trail of
the previous plume. When the second plume reached the confining unit, less lateral
spreading of the plume is observed. The third injection shows similar patterns: the central
part of the plume moved faster and less fingering is observed although the general shape
of the plume is the same. This effect will be enhanced with a higher number of sequential
injections. The final effect is a preferential path for the freshwater to reach the confining
unit rather rapidly, following the trail of previous plumes.
3.5 Need for Numerical Modeling
In general, similar fate and transport behavior was observed for all the injected
plumes in the physical experiments. However, two major differences were observed
between injected plumes. First, the amount of lateral spreading was significantly
different between plumes that came in contact with the overlying freshwater layer and
plumes that came in contact with the confining layer. Also, there were changes in plume
travel times (or plume velocities) caused by different density contrasts between ambient
and injected water. It was noticed that the plume in the first injection from experiment C
(where
s
? = 1.03 g/mL) moved slightly faster than the plume in experiment B (where
s
? =1.02 g/mL).
A major disadvantage to physical modeling is the constraint on varying problem
parameters that may not be easy to alter to investigate the changes in plume behavior.
Changing the density contrast is relatively easy; however, altering the permeability value
requires a change in the porous media selection and repacking the tank. Perturbing
parameters such as dispersivity, while maintaining a constant permeability, is unrealistic
in a laboratory experiment. This is where numerical modeling can be helpful. Numerical
modeling allows flexibility in selecting problem parameters making it better suited for
predicting results under different types of scenarios. The following sections discuss
numerical modeling that was undertaken for comparison to physical model results, as
well as for parameter sensitivity analysis.
26
27
CHAPTER 4
NUMERICAL MODELING OF EXPERIMENTAL RESULTS
4.1 Objective
The objective of this section was to develop a numerical model for the buoyant
plume experiments and compare the results against laboratory data. Experiment B was
chosen for this numerical study because out of all three experiments it allowed the largest
travel time for the plume.
4.2 Numerical Modeling of Density-dependent Flow and Transport
For the numerical portion of this study, we used the finite-difference numerical
code SEAWAT, published by the U. S. Geological Survey (USGS). In addition, we also
completed a simulation using the USGS finite element code SUTRA for intercode
comparison.
4.2.1 Governing Equations
Fluid flow in porous media is governed by Darcy's law and for density-dependent
systems it can be written in the form:
?
?
?
?
?
?
?
? ?
+??=
zf
ehKq
0
0
?
??
(4.1)
where q is specific discharge [LT
-1
], K is the freshwater hydraulic conductivity [LT
-1
], h
f
is the equivalent freshwater head [L], ? is the density of the water [ML
-3
],
0
? is the
reference water (freshwater) density [ML
-3
] and is the unit vector in the z-direction.
z
e
Mass continuity of the fluid in transient state, in the absence of sources and sinks,
can be written as:
)( q
tt
h
S
f
f
?
?
?? ???=
?
?
+
?
?
(4.2)
where is specific storage in terms of freshwater head [L
-1
] and
f
S ? is porosity.
Salt transport is described by the advection-dispersion equation written as:
))(( cIDDcq
t
c
m
?+??+??=
?
?
?? (4.3)
where c is the concentration of salt in the water [ML
-3
], D is the mechanical dispersion
coefficient [L
2
T
-1
], D
m
is the coefficient of molecular diffusivity of salt [L
2
T
-1
] and I is
the identity matrix.
A linear density state equation is assumed which couples the flow (4.2) and
transport (4.3) equations written as:
)1(
0
s
c
c
??? += (4.4)
28
where
0
0
)(
?
??
?
?
=
s
will be referred to as the relative density contrast and c
s
is the
maximum concentration of salt corresponding to saltwater density
s
? [ML
-3
].
4.2.2 Introduction to SEAWAT
SEAWAT couples a modified version of the MODFLOW code with the
MT3DMS transport code through the density term in order to solve variable density flow
problems [Langevin and Guo, 2006]. SEAWAT was used to numerically model
experiment B data and for the parameter sensitivity analysis.
The SEAWAT version of equation (4.2) used for modeling two-dimensional flow,
with the assumption that the x and z axes are aligned with the principal permeability
directions, can be written as:
s
f
fz
f
fz
f
fx
q
t
c
E
t
h
Se
z
h
K
zx
h
K
x
?????? ?
?
?
+
?
?
=?+
?
?
?
?
+
?
?
?
?
))(()( (4.5)
where:
fzfx
KK , are the freshwater hydraulic conductivities in the x and z directions,
respectively, [LT
-1
],
E is a slope constant relating the change in density to the change in salt
concentration( ) and is equivalent to 0.7143,
? is the density of source water [ML
-3
], and
s
q is the source/sink water flux [LT
-1
] [Guo and Langevin, 2002].
29
The SEAWAT version of the transport equation given by (4.3) for our two-
dimensional system without sorption or reactions can be written as:
() ()
30
?
ss
zxzzzxxzxx
cq
cv
z
cv
xz
c
D
x
c
D
zz
c
D
x
c
D
xt
c
?
?
?
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+
?
?
?
?
+
?
?
?
?
?
?
?
?
+
?
?
?
?
=
?
?
(4.6)
where:
zxxzzzxx
DDDD ,,, are the longitudinal and transverse hydrodynamic dispersion
coefficients in the x and z directions [L
2
T
-1
],
zx
vv , are the fluid velocities in the x and z directions [LT
-1
], and
is the external source concentration [ML
-3
] [Guo and Langevin, 2002].
4.2.3 Introduction to SUTRA
A second numerical code, SUTRA, was used for the purpose of intercode
comparison. SUTRA is a finite-element based numerical code also developed by the
USGS [Voss and Provost, 2003]. SUTRA includes unsaturated flow and transport
processes but is primarily intended for use in modeling saturated systems. SUTRA?s
equations for flow and transport are similar to those of SEAWAT with the main
difference being that the SUTRA equations are expressed in terms of pressure, whereas
the SEAWAT equations are expressed in terms of equivalent freshwater heads.
4.3 Numerical Model Description
As stated previously, experiment B was chosen to be numerically modeled. In
order to simulate advection with minimal numerical dispersion, SEAWAT offers several
advection solvers. The MOC (Method of Characteristics) was used as the advection
solver and the two-dimensional grid resolution was 0.1 cm ? 0.1 cm ? 2.7 cm. Time
stepping was determined by the Courant number. The initial and boundary conditions are
outlined below. This SEAWAT model of experiment B will be referred to as the base-
case model.
4.3.1 Initial Conditions
The initial condition for this problem was that of a fully saline aquifer with an
ambient density of 02.1=? g/mL. The top and bottom boundaries were no-flow
boundaries. The initial head in the problem was 29 cm for all cells. The assumed
problem parameters are summarized in Table 2.
The injection rate was measured by injecting a known volume of freshwater
(60mL) into the tank over a fixed time period (41 sec). Porosity as well as injection and
ambient water density of the system are given in Table 2. The porosity of the porous
medium had been previously measured in our laboratory using both volumetric and
gravimetric methods [Goswami and Clement, 2007]. The density of the injected water, as
well as the ambient saline water, were measured using a standard ASTM 11H hydrometer
which has an accuracy of ? 0.001 g/mL. The value for molecular diffusivity of salt was
taken from published literature [Hughes et al., 2005]. The other parameter values listed
in Table 2 were found by calibration of the numerical model, or by assuming a previously
measured estimate or range and using it as a basis for calibration. For example, from
previous in situ measurements for the tank, the known range of hydraulic conductivity
values is from 950 to 1200 m/day. The in situ conductivity can vary because of the
31
changes in packing conditions. A value from this range was used as a starting point for
numerical model calibration. The value of K = 950 m/day (1.0995 cm/sec) was shown to
give the best matching quantitative results for experiment B data. The value of
longitudinal dispersivity was calibrated using the numerical model and is on the order of
a tenth of the grain size. This low value was required in the numerical model to replicate
fingering patterns seen in the physical experiments. This can be seen in the results from
the parameter sensitivity analysis as discussed later in Chapter 5. A 1/10 ratio for
transverse to longitudinal dispersivity is used, as suggested in the literature [Johannsen,
et al., 2002]. The value for the Courant number was set at 0.25 and comes from a
suggestion that ?each particle needs at least four time steps in order to pass one cell [Spitz
and Moreno, 1996].?
Table 2 ? Details of the numerical model for the base-case scenario (experiment B)
Model Type Confined Cr (Courant) 0.25
Injection Rate 60 mL/41 s
Injection Density 1.000 g/mL Hydraulic
Conductivity
1.0995 cm/s
Longitudinal
Dispersivity
0.01 cm
Ambient Density 1.019 g/mL
Transverse to
Longitudinal
Ratio
0.1
Porosity 0.39
Molecular
Diffusivity of Salt
1.477E-5 cm
2
/s
4.3.2 Boundary Conditions
32
No-flow conditions surround the outside of the domain in the numerical grid by
default. Therefore, only two zero-pressure boundary conditions were added, one at the
leftmost top corner and one at the rightmost top corner of the grids. These constant head
boundary modifications ensured that flow and concentration could be allowed to enter
33
and exit the system and to keep appropriate constant heads at these boundaries so that
static ambient flow conditions could be maintained. Similar boundary conditions have
been used previously in the literature, such as with the Elder problem [Voss and Souza,
1987]. Furthermore, these conditions force close to zero pressure head along the top
layer of model cells, which is similar to the air-water interface of the physical
experiment.
4.3.3 Injection Well
The injection point for all of the numerical simulations was simulated by placing
a well boundary condition at the corresponding coordinates in SEAWAT. Using the
multi-species capabilities of SEAWAT, both freshwater and a tracer representing the dye
were injected in the well in the numerical models so that dual molecular diffusivities of
NaCl and the dye could be incorporated. It is important to note that only the NaCl salt
concentration controls the density in the density state equation. However, based on our
past experience the dye is known to be highly soluble and has negligible effects on
density; also, the dye is nonsorbing [Goswami and Clement, 2007].
4.3.4 SUTRA Model Comparison
A numerical simulation of experiment B was also carried out with the finite
element code SUTRA [Voss and Provost, 2003] using quadrilateral elements. Implicit
finite differences are used for the time integration. The iterative methods chosen to solve
the linear system of equations are the conjugate gradient method for the flow equation
and ORTHOMIN for the transport equation. The Picard method is used to solve the non-
34
linear system of equations in the numerical problem. The grid resolution used (0.2 cm ?
0.2 cm ? 2.7 cm) is coarser than the grid used with SEAWAT due to limitations on the
number of elements that can be used in SUTRA.
4.4 Numerical Model Results
Figure 6 displays results for our numerical modeling effort for experiment B. In
Figure 6, we can see the comparison at specific times between the physical model (left-
most column) and the SEAWAT model (central column), also known as the base-case
numerical model. The shape of the plume appears to be matched well by the SEAWAT
results. Similar to the physical model, the plume in the SEAWAT model contains initial
instabilities that develop into preferential pathways and elongated fingers as time
progresses. Also, similar to the physical model, the fingers reach the confining unit at
slightly different times. Once the confining unit is reached, the fingers join to form a
conic shape, leaving behind a tail due to saltwater entrainment from mixing of the plume
with the ambient saline water.
The numerical model reveals similar vertical advection rates as the experimental
results, as well as similar lateral spreading of the plume once the overlying boundary is
encountered. Since actual concentration contours could not be obtained from the
physical experiment, the numerical model reveals much more qualitative information in
regards to concentration stratification at the end of the run and in reference to the amount
of mixing that is occurring between the plume and the ambient saline water. As can be
seen from the numerical model, freshwater begins mixing with the saline water once
35
injection begins. This mixing appears to increase as the plume moves vertically upward
causing an increase in plume density which leads to a subsequent decrease in the buoyant
force. This decrease in buoyancy causes an ensuing decrease in the vertical advection
rate. The amount of mixing with the ambient water increases as it vertically rises until it
reaches the confining layer. Mixing is what causes salt entrainment into the plume and
the subsequent tail of brackish water to be left behind. However, relatively freshwater
does appear to reach the confining layer, although it cannot be determined how fresh this
water is, except to say from Figure 6 that C/C
0
is in the range of 0 ? 0.25 which
represents ?relatively fresh? water.
(e)
(d)
(c)
(b)
(a)
Physical Experiment SEAWAT SUTRA
C/C
0
0
0.25 0.50 0.75 0.999
Figure 6 ? Experiment B physical and numerical model results (SEAWAT and
SUTRA). All times are from end of injection (seconds) (a) 27, (b) 369, (c) 685, (d)
1385, (e) 2631
The results of the simulation with SUTRA are shown in the right most column of
Figure 6. The results show a good agreement with the experimental results as well as with
the MOC simulations. Three fingers are developed at the early stages of the upwards
movement with the middle finger moving slightly faster than the lateral ones. Overall,
the vertical velocity is very similar to that of MOC and the physical model. Therefore,
SEAWAT (with MOC as the advection solver) when compared to the widely-used code
36
SUTRA and to the physical experiment, is shown to qualitatively model this system quite
well, although it exhibits more complex fingering.
It should be noted that at earlier times, as can be seen in Figure 6, the plume in
both numerical models does seem to have a slightly faster vertical velocity than the
plume in the laboratory experiment. This can be attributed to the uncertainty in
laboratory hydraulic conductivity measurements.
Furthermore, errors may also be present due to uncertainties in both ambient
water density and the well location. The error range for the density measurement is
equivalent to the precision of the hydrometer (? 0.001 g/mL) which is very small. Also,
measurement errors for the location of the well are equivalent to the precision of the scale
on the tank, which is
1 mm, equal to one cell width in the base-case numerical model.
?
4.5 Grid Discretization and Solver Sensitivity
While modeling the base-case scenario, we found that the simulated plume was
highly sensitive to the choice of advection solver and grid resolution. Therefore, an
investigation of plume sensitivity to the advection solver and grid resolution was
completed.
The grid resolutions used are summarized in Table 3. The finest grid resolution
had cells the same diameter as the glass bead grain size (~1mm).
37
Table 3 ? Numerical models grid discretization table
Discretization ?x ?y ?z
Coarse 0.4 cm 2.7 cm 0.4 cm
Medium 0.2 cm 2.7 cm 0.2 cm
Fine 0.1 cm 2.7 cm 0.1 cm
In addition to comparing effects of discretization levels on plume fingering and
transport, results from four advective transport solvers in SEAWAT were also compared.
The four solvers used were: the Total-Variation-Diminishing scheme (TVD?
ULTIMATE), standard Method of Characteristics (MOC), finite-difference with
upstream weighting (FDU) and finite-difference with central-in-space weighting (FDC).
The sensitivity analysis involved comparing the model outcomes qualitatively for
each of the twelve scenarios (three grid meshes ? four advection solvers). These results
are reported in Figure 7. This figure contains outputs at 369 seconds after the end of
injection for the base-case model described in Table 2. As can be seen from the results,
the outcome from the models using MOC as the advection solver appear to give the best
qualitative match of the physical model results. This can be attributed to MOC
encompassing less numerical dispersion than the other solvers. The results from the runs
using upstream weighting with finite-difference (FDU) (Figure 7) include too much
numerical dispersion prohibiting the plume from even fingering. The central-in-space
finite difference (FDC) results (Figure 7) show fingering, but appear to contain numerical
oscillations, which is common to results using this technique [Zheng and Bennett, 1995].
The results from the Total-Variation-Diminishing (TVD) method (Figure 7) also
demonstrate fingering but result in longer run times than the central-in-space finite
difference method and comprise more numerical dispersion than MOC. Even though
38
fingering is observed with both FDC and TVD methods, the shape of the plume in these
numerical models is not representative of the outline of the plume observed in the
physical model.
TVD
MOC
FDU
FDC
COARSE MEDIUM FINE
Figure 7 ? Grid resolution and solver sensitivity comparison: time 369 seconds after
end of injection
Also, from Figure 7, we can observe qualitatively how the buoyant plume is
sensitive to the grid resolution. The FDU method seems to be the least sensitive and can
again be contributed to too much numerical dispersion inherent in the model. Results
39
40
from FDC show that this method is very sensitive to the grid discretization as results
from all three grids show significantly different results. Oscillations observed in FDC for
the coarse and medium grids seem to stabilize in the fine resolution grid. MOC and TVD
also appear to be sensitive to mesh refinement as the fingering patterns change between
each grid size. However, in terms of overall shape of the plume, MOC exhibits a better
qualitative match to the physical experiment when compared to the other solvers and the
general shape appears to be less sensitive to grid refinement. Also, results from MOC are
known to yield less numerical dispersion than those from other solvers. Therefore, the
MOC solver on the fine mesh was chosen for modeling experiment B. The base-case
model with MOC on the fine mesh will also be used in the subsequent parameter
sensitivity analysis discussed in Chapter 5.
The results from the grid and solver sensitivity analysis are consistent with
guidelines set aside by Al-Maktoumi et al. [2007] in regards to grid refinement and
advection solver for unstable flow problems. They investigated sensitivity of the Elder-
Voss-Souza problem to spatial and temporal discretization and advection solver choice.
They suggested using ?a grid cell size of about 0.38% (dx) and about 0.6% (dz) of the
total length and depth of the domain, respectively.? This is only recommended with use
of a small enough Courant number. They recommend a Cr = 0.1 but we went with a
recommendation specific for MOC made by Spitz and Moreno [1996]. They recommend
that a Cr = 0.25 is sufficient, which was the value used in this numerical modeling study.
41
4.6 Method of Characteristics (MOC)
The method of characteristics (MOC) is a standard Lagrangian particle-tracking
method first used for solving transport in porous media by Garder et al. in 1964 and later
used by Konikow and Bredehoeft in 1978. SEAWAT?s use of MOC comes from the
implementation of MOC into the MT3D code by Zheng in 1990 [Zheng and Bennett,
1995]. The particles in MOC represent the corresponding concentration field, in which
each particle?s concentration corresponds to the concentration of the cell in the finite-
difference grid in which it is located. For these model runs, a dynamic pattern of
particles was used for the initial particle placement. One possible drawback to MOC, in
comparison to the other advection solvers, is that MOC is not mass conservative. Also
MOC ? along with TVD ? is the most computationally expensive advection solver. In
order to investigate more quantitatively how mass conservation for MOC compared to
results from other advection solvers on all grid sizes, the total mass of salt in the aquifer
over time was calculated and can be seen in Figure 8.
4.7 An Investigation of Mass Balance Closure in the Numerical Model
The total mass was calculated for each solver on each grid and compared to the
theoretical mass of salt in the system. The theoretical mass of salt was calculated to be
43.05 grams before injection and 41.45 grams after injection assuming 60 mL of saline
water, with a salt concentration of 0.0266 g/mL, was expelled out through the boundary
cells during injection. The measure of salt concentration for density of 1.019 g/mL
comes from the state equation given below:
(4.6)
EC
f
+= ??
where:
E is a slope constant relating the change in density to the change in salt
concentration( ) and is equivalent to 0.7143, and
C is the concentration of salt (g/mL) [Langevin and Guo, 2006].
E ?NOTE: The
shown above is not the same as the buoyancy term
given in equation
(4.4) although it is related via the equation
s
c
E
0
??
= .
In Figure 8, the values along the ordinate are based on relative concentration and
can be converted into actual mass by multiplying them by the salt concentration of
0.0266 g/mL. This theoretical line is the solid black line in the figure.
42
Figure 8 ? Graph of total mass of salt in model domain versus time after the
end of injection for grid resolution and solver sensitivity analysis
It should first be mentioned that for all model runs, total mass as shown in Figure
8 remained within 1 % error of the theoretical value. Therefore, although differences
between solvers were observed and will be discussed, these differences are minor. It can
be observed in Figure 8, that of the four scenarios on the coarse grid, TVD and the finite-
difference methods yield no oscillations in total mass of salt, while MOC loses mass.
The results from FDU and FDC are indistinguishable, as can be seen they actually
overlap in the figure, with the only decrease in mass occurring just a few seconds after
injection ceases. The total mass with the MOC solver steadily decreases from injection
until eventually stabilizing around 1000 seconds on each grid. This decrease is attributed
to the lack of mass conservation in MOC. For the coarse grid, all four advection solvers
yield overestimates of the theoretical total mass while for the medium grid, all four
43
44
advection solver scenarios yield underestimates. For the medium grid, the results for
TVD, FDU and FDC are identical on the graph and all three are closer to the theoretical
value than the results from MOC. The results from the fine grid yield better results for
MOC and the same results for TVD when compared to the medium grid. The outcomes
for FDU and FDC, however, are even less than the theoretical total mass value when
compared to the same solver results on the medium grid. In this case, TVD on the fine
grid yields the overall best matching and most stable results out of all solvers on all grids.
However, MOC on the fine grid performed similarly except for its characteristic loss of
mass after injection until stabilizing at around 1000 seconds. It should also be mentioned
that MOC on the coarse grid performed well according to this quantitative measure;
however, qualitatively, it was not as good of a match. Although MOC was shown from
this analysis to be less mass conservative than the other three advection solvers, its values
still only fluctuated by at most 0.3% throughout the model runs after injection ceased (42
seconds). MOC as the advection solver on the fine resolution grid appears to yield the
best overall results for this buoyant plume scenario (Figure 7), even though its
computational time was greater, in some cases more than four times the length of run
times using finite-difference methods.
4.8 Need for Detailed Parameter Sensitivity Analysis
As discussed from the physical modeling results (Chapter 3), changes in the
density contrast between ambient water and injected water caused a change in the travel
time for the plume to reach the overlying confining layer. The greater the density
contrast, the less the travel time. Therefore, it appears the time needed for the plume to
45
travel to the confining layer is inversely proportional to the density contrast. Because of
this, it is difficult to directly compare outputs of these experimental runs to one another.
It is also an understandable question to ask what kind of effect changes in hydraulic
conductivity would have on this travel time. Is the travel time also inversely proportional
to K ? Also, what other parameters are critical in this problem? For this, a dimensional
analysis of the governing equations is necessary. From this dimensional analysis, the
important problem parameters can be identified and they can be used for further
sensitivity analysis. This analysis is the focus of the next chapter.
46
CHAPTER 5
PARAMETER SENSITIVITY ANALYSIS
5.1 Introduction
In order to determine ?important problem parameters?, a dimensional analysis
was undertaken of the governing equations for flow and transport. According to
Wikipedia.com, dimensional analysis is ?a conceptual tool ? used to form reasonable
hypotheses about complex physical situations that can be tested by experiment or by
more developed theories of the phenomena [Wikipedia, 2008].?
5.2 Dimensional Analysis
NOTE: This dimensional analysis was completed and contributed by Dr. Elena Abarca as
part of her continuing work with this problem.
A dimensional analysis of the flow and transport equations given in Section 4.2.1
was completed to analyze the upwards movement of the plume, determine a characteristic
plume travel time, and determine which physical parameters are integral to the problem.
The general set-up used in the analysis is described in Figure 9. The physical parameters
derived through this analysis will also be used for the parameter sensitivity analysis.
47
L
H
d
x
z
Figure 9: Problem Dimensions
We first write the governing equations in a dimensionless form using the set of
dimensionless variables summarized in Table 4.
Table 4 ? Dimensionless Variables
d
x
x
D
=
d
z
z
D
=
'
,,
1
yxyx
d
?=?
d
h
h
D
=
)1(
0
??
?
?
+
=
D
?K
q
q
q
q
c
D
== ?K
d
t
vd
t
t
t
t
cc
D
===
Dc
ttt ?
=
?
? 1 ?
where:
d is the distance from the center of mass of the plume to the overlying confining unit (see
Figure 9),
h is the equivalent freshwater head, and
K is the freshwater hydraulic conductivity.
NOTE: The subscript D in the variables and parameters summarized in Tables 4 and 5
stand for dimensionless while the subscript C stands for characteristic.
Using the variables presented in Table 4, equations (4.1), (4.2) and (4.3) given in
Chapter 4 can be written as:
aahaq
DDD
?++??= )1(
'
?
(5.1)
DDDD
D
D
D
D
D
qq
tt
h
'' ????=
?
?
+
?
?
??
?
??? (5.2)
(5.3)
The dimensionless expression for Darcy?s flow is shown in (5.1) while the new
form of the flow equation is listed in (5.2). The boundary conditions of the salt water?s
equivalent freshwater head, zh ?= are specified in the lateral boundaries
( and ) and is imposed in both the upper and bottom impermeable
boundaries. The dimensionless forms of these boundary conditions are:
=
z
q0=x Lx = 0
for (5.4)
48
for (5.5)
The dimensionless form of the transport equation is shown in (5.3). Salt mass
flux across the boundary is zero at the bottom and upper boundaries. At the lateral
boundaries, we prescribe the concentration equal to the ambient saltwater (c
s
). The
dimensionless form of this boundary condition is:
for (5.6)
From these dimensionless equations, a number of dimensionless numbers are
obtained and are summarized in Table 5. Also, from dimensionalizing the time, a
characteristic time ? )/( ?Kdt
c
= ? was derived, as seen in Table 4. This characteristic
time is an average time required for the plume to traverse the distance from injection
point to the overlying confining unit based on the maximum buoyant force in the system
after injection.
Table 5 ? Dimensionless Problem Parameters
49
?
1
=a
dK
D
b
m
D
?
?
=
dS
s
= ?
?
d
b
L
L
?
=
d
b
T
T
?
=
where:
? is the porosity,
m
D is the molecular diffusivity of the salt,
s
S is the specific storage in terms of freshwater head,
L
? is the longitudinal dispersivity, and
T
? is the transverse dispersivity.
As can be seen in Table 5, several key problem parameters arose out of the
dimensionless analysis. These are
Lms
DSdK ??? ,,,,, and
T
? . For our parameter
sensitivity study, four important parameters were selected from this list. Sensitivity of
the plume to changes in the other parameters should be carried out in future work related
to this problem. The parameters chosen were the relative density contrast between
injected and ambient water ?, the saturated hydraulic conductivity K, and the longitudinal
and transverse dispersivity values
L
? and
T
? . The buoyancy, or relative density
contrast? , was chosen to be investigated because as we saw in our physical modeling
results, an increase in ambient water density caused a decrease in travel time of the plume
from the injection point to the overlying confining layer. It is also noted that this density
contrast appears in the formula for characteristic travel time of the plume. Hydraulic
conductivity is also found in this formulation of characteristic time. Therefore, it was
also of interest to investigate the sensitivity of the problem to and to determine if the K
?perturbed cases of and K scale well to the dimensionless time, as it is hypothesized
they should. The longitudinal dispersivity value was selected for the parameter
sensitivity analysis because the
L
? used for the numerical modeling in the last chapter
was on the order of a tenth of the grain size, which was unusually small compared to
what is normally used. If
L
? values are to be perturbed while maintaining the ratio of
T
?
50
to
L
? of 0.1, then
T
? values must also be perturbed. Therefore,
T
? was also
selected for the parameter sensitivity analysis although results for it are shown in the
L
? results.
To investigate the sensitivity of the plume to the selected parameters, each
parameter was perturbed by ? 50% about the base-case value (Table 2) while keeping all
other model parameters constant. Therefore, runs of the base-case model with each
changed parameter were undertaken to compare effects of changes in these parameters on
plume behavior. The values of these perturbations are summarized in Table 6.
Furthermore, it was also of interest to investigate the effect of changing the longitudinal
dispersivity
L
? to 0.1 cm (and
T
? to 0.01 cm). This value of
L
? is equivalent to the
actual average grain diameter of the porous media, a typical used value for longitudinal
dispersivity in numerical modeling studies. To further compare the changes due to this
parameter, an even larger value of
L
? = 0.2 cm (and
T
? = 0.02 cm) was also chosen to be
modeled.
Table 6 ? Parameter sensitivity values
-50% Base-case +50%
Density Contrast 1.0095
g/mL
1.0285
g/mL
1.019 g/mL
?
Hydraulic
Conductivity
0.5498
cm/s
1.6493
cm/s
1.0995 cm/s
K
Dispersivity
51
L
?
0.005 cm 0.01 cm 0.015 cm
0.0005 cm 0.001 cm 0.0015 cm
T
?
An investigation of the sensitivity of the plume behavior to changes in these
problem parameters was completed using two methods. In the first approach, a
qualitative analysis involved simple comparison of outputs of perturbed model runs to
one another and the base-case. In the second approach, a more rigorous technique that
consisted of a quantitative investigation of perturbed parameter model runs using moment
analysis was used. Results for the sensitivity analysis will first be qualitatively reported
in each section, followed by a discussion of quantitative results.
5.3 Methodology for Qualitative Sensitivity Analysis
Numerical model runs of the perturbed parameters were completed and pictorial
outputs from these runs were organized into figures for easy comparison.
A characteristic time, based on the formula given in Table 4, was calculated for
each of the model runs and, using this time, the model run times were non-
dimensionalized for easier comparison. These characteristic times needed to only be
calculated for the? and K sensitivity scenarios as dispersivity values were not present in
the formula for . Therefore, by scaling to dimensionless time, differences due to bulk
c
t
52
advection were omitted from the figures so that the results could be compared for
differences in plume fingering, spreading and mixing. A qualitative comparison of
mixing in this study involved comparing pictorially the amount of red water, representing
the freshest water, that reaches the overlying confining unit by T
D
= 0.67 in the figures.
It is hypothesized that the results from the parameter sensitivity will show that
when the values of hydraulic conductivity and density contrast increase, vertical
velocities of the plume increase. Furthermore, for lower dispersivity values, it is
hypothesized that the plume should have a larger vertical velocity because less mixing
with ambient water should be occurring.
5.4 Methodology for Quantitative Sensitivity Analysis
It was inferred from the qualitative results, which will be reported in sections
5.5.2 ? 5.5.4, the fingering of the plumes can be quite complex. This can be attributed to
the more random nature of MOC when compared to other advection solvers, its low
numerical dispersion, and the small values which were required for longitudinal and
transverse dispersivities to better calibrate the model to the physical experiment.
Therefore, qualitatively comparing the model outcomes proved to be difficult, and it did
not help us fully understand the effects of parameter changes on physical processes such
as spatial spreading, mixing of the plume with the ambient water, plume fingering and
vertical velocity. For this reason, there was a need to analyze the sensitivity analysis
results using a more rigorous quantitative method, which will be discussed in the
following section. The quantitative method was based on moment analysis.
5.4.1 Quantifying the Plume Using Vertical Bulk Velocity and Spatial Variances
In order to better quantify the plume movement, a description for the center of the
plume was needed. This would allow calculation of the vertical velocity of the plume as
a whole, which we will refer to as the bulk velocity. The center of mass of each plume
was used as an indicator of the center of the plume. Center of mass is equivalent to the
centroid of the plume when the plume has uniform density. However, our plume
increases in density as it moves vertically and mixes with the ambient water. Therefore,
we will use the center of mass to characterize the movement of our plume. The bulk
velocity will be calculated by measuring the rate of change of the z coordinate of the
center of mass with time. In the quantitative results, which will be reported throughout
the remainder of this chapter, the z coordinate of the center of mass will be plotted at
several dimensionless times. These data points will then be fitted with a linear trend line.
The slope of this trend line is representative of the bulk velocity of the plume.
The center of mass or the mean of the concentration distribution, can be
calculated from the spatial moments of the concentration distribution of the plume. The
lower spatial moments in the x-direction are as follows [Fischer, et al., 1979]:
zeroth moment = M
0
= (5.7)
?
?
?dtzyxC ),,,(
first moment = M
1
= (5.8)
?
?
?dtzyxxC ),,,(
second moment = M
2
=
(5.9)
?
?
?dtzyxCx ),,,(
2
where ? is the total volume of the domain and is the concentration of solute. ),,,( tzyxC
53
The center of mass is equivalent to M
1
/M
0.
The variance of the concentration
distribution about the center of the mass (or what is referred to in the ground water
literature as the dispersion from the center of mass in x and z directions) can be
calculated from the second central spatial moment normalized by the total solute mass
[Kitanidis, 1994]. For this problem, the x and z coordinates of the center of mass of the
plume, along with the corresponding variance in each direction were calculated using the
formulas:
54
?
?
?
?
?
?
=
dtzyxC
dtzyxxC
tx
),,,(
),,,(
)( (5.10)
?
?
?
?
?
?
=
dtzyxC
dtzyxzC
tz
),,,(
),,,(
)( (5.11)
?
?
?
?
?
??
=
dtzyxC
dtzyxCxx
t
Lx
),,,(
),,,()(
)(
2
2
? (5.12)
?
?
?
?
?
??
=
dtzyxC
dtzyxCzz
t
Lz
),,,(
),,,()(
)(
2
2
? (5.13)
In the above equations, C is the concentration of solute, Error! Bookmark not
defined. is the x coordinate of the center of mass of the plume, ?is the total volume, x
z is the z coordinate of the center of mass of the plume and and are variances or
dispersion measures in both the x and z directions respectively. The x present in
equations (5.10) and (5.12) is the distance from the z-axis to the center of the
2
Lz
?
2
Lx
?
corresponding cell while the z in equations (5.11) and (5.13) is the distance from the x-
axis to the center of the corresponding cell.
These values were calculated from the concentrations present in each cell in the
model domain at a specific point in time. This is based on the coordinate system with the
origin (0, 0) corresponding to the leftmost bottom corner of the grid as shown in Figure
10.
Figure 10 ? Coordinate system of numerical model for moment calculation
Using the model outputs from the previously completed advection solver
sensitivity analysis, discussed in section 4.5, the quantitative analysis methodology will
be tested and compared against the qualitative results.
5.4.2 Quantitative Method Testing ? Advection Solver Comparison
In section 4.5, qualitative results were reported (Figure 7) for the advection solver
sensitivity analysis. It was concluded, for results on the fine grid, that MOC gave the
55
56
overall best matching results when compared to the physical experiment. Also, results
from the finite difference methods exhibited significant numerical dispersion even on a
fine grid. As a way of testing our quantitative analysis methodology before applying it to
results from the parameter sensitivity analysis, we applied the quantitative methodology
to the results for the advection solver sensitivity on the fine grid. It should be noted here
that the results for the last time step of T
D
= 0.67 were ignored in some of the analysis
since the plume had interference from the boundary.
Figure 11 shows the z coordinates of the center of mass for the advection solver
comparison graphed against the dimensionless times of T
D
= 0, 0.11, 0.20, 0.36 and 0.67.
The data points have been fitted with linear regression trend lines and the equations of
these lines are shown in the colored boxes (the color corresponding to the color of the
data points). The slopes of these lines are equivalent to the rate of change of the vertical
location of the center of mass with respect to time. These are estimates of the vertical
bulk velocities of the plumes.
Figure 11 ? Z-coordinate of center of mass versus dimensionless time ? advection
solver comparison
It can be observed in Figure 11 that the centers of mass of all four plumes move at
approximately the same rates. However, there are slight differences. The bulk velocity
rates of results from TVD and MOC are the largest while the velocities of the finite-
difference methods are smaller. The FDU results exhibited the smallest bulk velocity.
This can be attributed to a higher amount of numerical dispersion inherent with finite-
difference, resulting in a slightly impeded velocity. The actual numerical values for the x
and z coordinates of the center of mass for these results are listed in Table 7. Overall, it
appears that the choice of advection solver does make a slight difference in terms of bulk
velocity rates, especially at later times.
57
Table 7 ? Center of mass coordinates and horizontal and vertical variance values ?
advection solver comparison on fine grid ? blue is low variance for time, red is high
variance for time
MOC (BASE-CASE) TVD
58
T
D
x
(cm)
z
2
Lx
?
2
Lz
?
(cm)
(cm
2
)
(cm
2
)
x
(cm)
z
2
Lx
?
2
Lz
?
(cm)
(cm
2
) (cm
2
)
0 27.044 6.177 4.715 4.498 27.0444 6.187 4.765 4.548
0.11 27.053 9.101 4.545 4.938 27.046 9.129 4.529 4.87
0.20 27.066 11.663 5.286 5.729 27.05 11.668 4.548 5.261
0.36 27.138 16.09 8.948 7.222 7.867 27.066 16.152 6.325
0.67 27.49 23.265 15.843 10.91914.863 27.128 22.799 9.646
FDU FDC
T
D
x
(cm)
z
2
Lx
?
2
Lz
?
(cm)
(cm
2
)
(cm
2
)
x
(cm)
z
2
Lx
?
2
Lz
?
(cm)
(cm
2
) (cm
2
)
0 27.044 6.142 5.095 4.734 27.044 6.165 4.945 4.627
0.11 27.045 8.713 5.782 6.152 27.045 8.853 5.661 5.603
0.20 7.294 27.047 10.858 5.407 27.048 11.103 5.414 6.233
0.36 27.05 14.587 4.635 10.701 27.051 15.022 4.87 8.62
0.67 3.415 23.50727.051 22.184 27.075 22.773 4.708 18.847
To further compare results from different SEAWAT advection solvers, the
horizontal and vertical spatial variances were calculated from results of each advection
solver model run on the fine grid. These values are summarized in Table 7. In the
horizontal direction, results from MOC displayed maximum changes in variance over
time, followed by TVD, FDC and FDU. This can be seen in Figure 12. The results from
finite-difference (upstream weighting) actually show a decrease in horizontal variance
over time, while variances for finite-difference (central) seem to remain relatively
constant around 5 cm
2
. The solver exhibiting the most lateral variance was MOC with its
largest variance being 15.8 cm
2
at T
D
= 0.67. However, at this last time, effects from the
upper model boundary may be interfering. However, even at T
D
= 0.36, MOC had the
largest lateral variance (8.9 cm
2
).
Figure 12 ? Horizontal variance versus dimensionless time ? advection solver
comparison
Results for calculated vertical variances are plotted in Figure 13. Unlike
horizontal direction, all four advection solvers displayed an increase in variance around
the center of mass as time progressed, with the largest variation occuring for FDU.
Variances observed for FDC were less than that of upstream weighting, but were more
than the results of MOC and TVD. TVD indicated the least amount of vertical variance
about the center of mass.
59
Figure 13 ? Vertical variance versus dimensionless time ? advection solver
comparison
It is of our interest to understand whether there is a direct relationship between the
horizontal and vertical variance estimates and spreading and/or mixing of the plume. In
order to understand this relationship, outputs from the model runs of advection solver
sensitivity scenarios were compared to numerical values of horizontal and vertical
variances from Table 7. At early times, little difference can be observed between results
from the different advection solvers, which can also be seen in Figures 12 and 13. At
time T
D
= 0.67, interference may exist from the upper boundary of the model, as
mentioned previously. Therefore, we analyzed the outputs from T
D
= 0.36, which are
shown in Figures 14 and 15.
60
Figure 14 ? Qualitative outputs for T
D
= 0.36 ? advection solver comparison ?
investigating horizontal variance about the center of mass
As observed in Table 7, the results from MOC exhibit the greatest horizontal
variance about the center of mass at time T
D
= 0.36. If spread is defined as the extent, or
distance of the edge of the plume about the center of mass in a particular direction, then
these four plumes have similar lateral spread as can be seen in Figure 14. However, if we
look at the actual values for horizontal variance, given in Table 7, we can see that the
numbers vary considerably. For example, horizontal variance for FDU (4.6 cm
2
) is
approximately half that of MOC (8.9 cm
2
). This cannot be physically attributed to just
spreading as the extent of both plumes in the horizontal direction appears to be the same
(Figure 13). However, it is clear that the variance about the center of mass is greatest
61
with MOC and least with FDU in the horizontal direction. The MOC plume contains
high amounts of freshwater at greater distances from the center of mass. The FDU
plume, on the other hand, appears to be more smoothly distributed with concentrations
decreasing as we move away from the center of mass. In other words, the concentration
of the FDU plume is more ?centered? about the center of mass.
Figure 15 ? Qualitative outputs for T
D
= 0.36 ? advection solver comparison ?
investigating vertical variance about the center of mass
For vertical variance, Table 7 shows that FDU has the greatest variance about the
center of mass (10.7 cm
2
) and TVD has the least (7.2 cm
2
) for time T
D
= 0.36. This can
also be observed in Figure 15 as FDU exhibits a significant amount of vertical variance
while TVD seems to show much less variance. The vertical variance of TVD is
approximately 67% of that of FDU. Again, it should be noted that this does not directly
relate to the amount of spreading of the plume nor to the amount of mixing, which cannot
be quantified from the pictures due to the complexities of the plume. It appears that the
best description of the spatial variance is that it is an estimate that quantifies the
distribution of the concentrations about the center of mass in a particular direction.
The quantification of the center of mass used to calculate the bulk velocity
appears to agree with qualitative results. Also, we have learned that the spatial variance
gives an estimate of the spatial concentration distribution of the plume about the center of
mass in both the horizontal and vertical directions. Overall, our quantitative
62
63
methodology appears to agree with the qualitative results, but also appears to yield more
information than can be extracted from the qualitative analysis alone. This analysis will
further be used to analyze the base-case simulation results.
5.5 Parameter Sensitivity Results
5.5.1 Application of Moment Analysis Method to Base-case Simulation Results
From testing the quantitative methodology with the advection solver sensitivity
outputs, the base-case model was analyzed (MOC advection solver). The purpose of this
section is to analyze the simulation results for the base-case using the quantitative
method. These results will be later used for comparison against different sensitivity
results.
First, as reported in Table 7, the center of mass was calculated at five
dimensionless times. The z coordinate of this center of mass was graphed versus the x
coordinate and the vertical movement of the plume over time can be visualized in Figure
16. The initial placement of the plume should theoretically be 5 cm from the bottom and
27 cm from the left (see Figure 10). This placement is based on instantaneous injection
with no mixing area around the plume. The initial placement of the base-case plume
from the numerical model, as can be seen in Table 7, is 27.04 cm from the left and 5.86
cm from the bottom. The results from the numerical model show a higher initial starting
point because in the actual model we do not have instantaneous injection; instead we
have buoyancy occurring while injection is occurring. This explains why the z-
coordinate of the center of mass is slightly higher than the theoretical value.
Figure 16 ? Center of mass movement in x-z plane over time for base-case
In Figure 17, the z coordinate of the center of mass of the plume in the base-case
scenario was plotted versus the elapsed time from the end of injection. The five data
points were fitted with a linear regression trend line. From the slope of this line, the
actual vertical bulk velocity of the plume was determined to be 0.0227 cm/sec. From the
characteristic time formulation in Table 4,
?K
d
v
d
t
c
c
== , the characteristic velocity ( )
is the product of the buoyancy and the hydraulic conductivity. For the base-case, this
velocity is calculated to be 0.021 cm/sec which corresponds well to the slope of the trend
line in Figure 17. Therefore, this characteristic time ( ) appears to correctly quantify the
time required for the plume to traverse the distance from injection point to confining
layer and should allow for proper scaling of the parameter sensitivity results to the
dimensionless time.
c
v
c
t
64
Figure 17 ? Z coordinate of center of mass versus real time for base-case
Variances in the x and z directions about the center of mass calculated for the
base-case scenario were graphed against dimensionless time in Figure 18. It was
observed that variances in both directions increased with time. This was a trend
generally observed in all scenarios simulated during the parameter sensitivity
investigation, which will be discussed later. Also, with the base-case model, the lateral
variance increased slightly more than the vertical variance at later times. However, this
difference was not statistically significant. Therefore, it appears that the variances in the
horizontal direction are approximately the same as the variances in the vertical direction
for the base-case model.
65
Figure 18 ? Variances in horizontal and vertical directions about the center of mass
versus T
D
? BASE-CASE
5.5.2 Density Contrast Sensitivity Results
Figure 19 shows results from the density sensitivity model runs. Differences
between plumes in Figure 19 can be seen as early as T
D
= 0.11 and become more
pronounced as time elapses. The times reported in Figure 19 are all dimensionless and
based on the maximum buoyancy in each scenario. The results show that the vertical
velocity of each plume scales well to this dimensionless time. However, the plume with
the larger density contrast appears to have a slightly faster vertical velocity. This can be
easily seen at T
D
= 0.36 in Figure 19. With increase in ambient water density (or density
contrast), the vertical variance of the plume increases while the horizontal variance
appears to remain the same. From Figure 18, it is difficult to qualitatively determine the
effects on mixing because of the complexities of fingering in each plume. Therefore, no
conclusions can be drawn, as far as mixing, except to say that it appears from these
pictorial outputs that density has little effect.
66
67
T
D
= 0.11
T
D
= 0.20
T
D
= 0.36
T
D
= 0.67
? = 1.0095 g/mL ? = 1.019 g/mL ? = 1.0285 g/mL
C/C
0
0
0.25 0.50 0.75
0.999
Figure 19 ? Density sensitivity results (time is dimensionless time after the end of
injection)
Table 8 provides a summary of x and z coordinates of the center of mass of each
plume for the density sensitivity scenarios. The first observation is that the horizontal
components of the center of mass remain close to 27 cm even though we do see some
deviations, as much as ? 0.5 cm, for the last time step for the base-case scenario. The
explanation for the deviation of the center of mass from 27 cm can be explained by the
asymmetry present in the plumes in the MOC simulations. Due to the random nature of
MOC, when compared to other advection solvers, symmetry about the centerline is lost.
This is caused by preferential pathways or fingers being formed and the mass of the
plume moving vertically through these fingers. All movement of the center of mass in
the lateral direction remains within 2 % of the theoretical value of 27 cm.
Table 8 ? Center of mass coordinates for base-case and density sensitivity scenarios
68
=
s
1.0285 g/mL ? 1.0095 g/mL
BASE-CASE
=
s
?
T
D
x (cm) z (cm) x (cm) z (cm) x (cm) z (cm)
0 27.0445 5.86292 27.0444 6.17697 27.0448 6.49743
0.11 27.0473 8.73511 27.0525 9.10062 27.0508 9.48668
0.2 27.0591 11.2981 27.0662 11.6627 27.073 12.1593
0.36 27.065 15.7201 27.1379 16.0901 27.0391 16.8992
0.67 26.9255 22.5067 27.4894 23.2646 26.8569 23.6433
Figure 20 shows results from linear fitting of data points from temporal changes
in the z coordinate of the center of mass of each plume resulting from the density
sensitivity scenarios. It can be seen from this figure that a perturbation of ?50% about
the base-case value corresponds to an increase or decrease in the velocity of the plume by
approximately ?50%. In other words, it is observed that the velocities for the three
density contrast scenarios scale well to the dimensionless time.
Figure 20 ? Z-coordinate of center of mass versus dimensionless time for density
scenarios
In both directions for density sensitivity, we can see that at early times little
difference can be noted between spatial variances while at later times, there are more
differences. With the x direction variance, results from low and high ambient water
density both resulted in less variance at later times when compared to the base-case
(Figure 21). However, the small differences exist within the error produced by the
randomness of MOC. With variance in the z direction (Figure 22), we can see that higher
density results in a larger variance in this direction while the lower density results have
approximately the same variance as the base-case. The highest variance was observed
with the high density results at T
D
= 0.36 in the vertical direction. It should be noted that
the differences between variances in the lateral direction were not statistically significant.
However, in the vertical direction, a larger density contrast appears to cause more
variance with increasing time. This is nonetheless caused by the faster advection of the
69
middle finger as can be seen in Figure 19. When the results for TD = 0.67 are
considered, the difference between vertical variances decreases. Therefore, it is possible
that the results for high density contrast at time TD = 0.36 could be an outlier and within
the range of errors attributed to the randomness of MOC. Therefore, it can be stated that
from both the qualitative and quantitative results, density contrast (?) does not play a
significant role in affecting the spatial variance of concentrations in the plume. It does,
however, play a significant role in affecting vertical bulk velocity rates. The bulk
velocity increases, or decreases, proportional to an increase, or decrease, respectively of
?. However, when scaled using the characteristic time, the effective vertical transport of
the plumes are approximately the same.
Figure 21 ? Horizontal variance versus dimensionless time for perturbed density
scenarios
70
Figure 22 ? Vertical variance versus dimensionless time for perturbed density
scenarios
5.5.3 Hydraulic Conductivity Sensitivity Results
Figure 23 shows results from hydraulic conductivity sensitivity runs. Similar to
the density contrast results, the differences in plumes can be seen as early as T
D
= 0.11.
The plume velocities for all three scenarios appear to be the same when scaled to the
dimensionless time. The spatial variance in both the vertical and horizontal directions do
not appear to be affected. Again, the complexity of the plume fingering prohibits any
conclusions to be reached regarding mixing of the plume with the ambient saline water.
71
72
T
D
= 0.11
T
D
= 0.20
T
D
= 0.36
T
D
= 0.67
K = 0.5498 cm/sec K = 1.0995 cm/sec K = 1.6493 cm/sec
C/C
0
0
0.25 0.50 0.75
0.999
Figure 23 ? Conductivity sensitivity results (time is dimensionless time after the
end of injection)
Table 9 shows the numerical results for computation of the coordinates of the
center of mass for each of the hydraulic conductivity scenarios. The results were similar
to those for the density contrast sensitivity scenarios. The horizontal component of the
center of mass was located approximately 27 cm from the left, with a deviation of
approximately 2% at most.
Table 9 ? Center of mass locations for base-case and hydraulic conductivity
sensitivity scenarios
K = 0.5498 cm/sec K = 1.6493 cm/sec BASE-CASE
73
T
D
x (cm) z (cm) x (cm) z (cm) x (cm) z (cm)
0 27.0444 5.86944 27.0444 6.17697 27.0447 6.49173
0.11 27.0466 8.74278 27.0525 9.10062 27.055 9.48356
0.2 27.06 11.3154 27.0662 11.6627 27.0849 12.1251
0.36 27.011 15.9282 27.1379 16.0901 27.1035 16.6756
0.67 26.6874 23.4084 27.4894 23.2646 26.9566 22.7375
The vertical bulk velocities for the conductivity scenario plumes can be seen in
Figure 24. Similar to the velocities for the different density scenarios, the velocity
increased with an increase in hydraulic conductivity, and decreased with a decrease in
conductivity. But when scaled against the dimensionless time as in Figure 24, the net
vertical transport velocities are almost the same.
Figure 24 ? Z-coordinate of center of mass vs. dimensionless time for hydraulic
conductivity scenarios
Figures 25 and 26 show results for lateral and vertical variances for the
conductivity scenarios. Similar to the density scenarios, we see little difference at early
times in the x and z direction variances between perturbed cases and the base-case. Even
the differences that are observed at later times do not seem to be statistically significant.
We detect slightly greater variance in the x direction than in the z direction, especially at
later times. From Figures 24 and 25, it can be deduced that changing the hydraulic
conductivity does not appear to affect spatial variance of the concentrations.
Figure 25 ? Horizontal variance versus dimensionless time for perturbed hydraulic
conductivity scenarios
74
Figure 26 ? Vertical variance versus dimensionless time for perturbed hydraulic
conductivity scenarios
Since both density contrast and hydraulic conductivity are present in the
characteristic time formula used to non-dimensionalize the time, it was of interest to
compare results from these scenarios to one another. In other words, the question arose
as to whether or not changes in ? yielded similar results as changes in K.
It was first of interest to compare initial locations of the z coordinate of the center
of mass at the end of injection. It was found to correctly correspond to the perturbed
density or hydraulic conductivity as can be noted in Tables 8 and 9. With an increase in
ambient water density or hydraulic conductivity, a higher location of the center of mass
was observed when compared to the base-case, while a decrease in either of these two
parameters caused a subsequently lower initial location. This again is related to the
plume exhibiting vertical buoyant movement while injection was still occurring since it
occurred over a 41-second time period, which is not instantaneous.
75
Figure 27 compares results of vertical center of mass movement for all density
contrast and conductivity scenarios. It can be observed that the slopes of the lines are
approximately the same for the base-case as they are for the other perturbed cases. This
shows again that the characteristic time we computed is a good measure to non-
dimensionalize the time as equivalent results can be observed between equivalent
changes in K and ?. It can also be observed from this figure that results from the 50%
decrease in density contrast correspond to results from the 50% decrease in conductivity.
Similar results were observed for an increase of K and ? values.
Figure 27 ? Z-coordinate of center of mass versus dimensionless time ? comparing
results from density to conductivity
In Figures 28 and 29, we compare the density and conductivity variance results on
the same graphs. Figure 28 shows the results for lateral variance comparison. At early
times, all scenarios seem to correspond to one another. However, later, slight differences
are observed. In the x direction, low and high perturbations of conductivity yield greater
spatial variance while low and high values of density yield less horizontal spatial
76
variance. In summary, the lateral variance appears to be more sensitive to changes in
conductivity than density.
Figure 29 shows results for the vertical variance comparison. In the vertical
direction and at early times, the variances between base-case, density and conductivity
scenarios appear to correspond well. Later, the high density perturbation exhibited
greater variance in this direction. However, as discussed previously, this result for high ?
at time T
D
= 0.36 could be an outlier as differences between variances decrease at later
times (not shown). If this is considered as an outlier, it appears perturbations in ? or K
cause little change in vertical variances.
Figure 28 ? Horizontal variance versus dimensionless time for density and
conductivity comparison
77
Figure 29 ? Vertical variance versus dimensionless time for density and
conductivity comparison
The characteristic time used for non-dimensionalizing the elapsed times in these
figures can be seen in Table 10. Since the characteristic time formula is not dependent on
dispersivity, there is no need to compute new characteristic times for these scenarios as
they are the same as the base-case. It is expected that for the same dimensionless times,
we should observe similar bulk vertical movement of each plume. And, in fact, this is
what was observed in both the qualitative and quantitative results.
Table 10 ? Characteristic times in seconds for density and hydraulic conductivity
scenarios
Characteristic Time (t
c
) Scenario
?
s
= 1.0095 g/mL
2297.6 seconds
or
78
K = 0.5498 cm/s
BASE-CASE
?
s
= 1.019 g/mL 1148.8 seconds
or K = 1.0995 cm/s
?
s
= 1.0285 g/mL
765.9 seconds
or K = 1.6493 cm/s
5.5.4 Dispersivity Sensitivity Results
Dispersivity sensitivity analysis occurred in two parts. First, similar to the other
parameters, the dispersivity values were perturbed ?50% from the base-case value to
examine effects on plume transport and spatial variance (Investigation I). However, since
the base-case dispersivity was so small compared to what is normally used in a numerical
model (grain size longitudinal dispersivity), a second investigation (Investigation II) was
undertaken to compare the base-case dispersivity to that of grain-size dispersivity and an
even larger value. It should be noted that the ratio of transverse to longitudinal
dispersivity (
LT
?? ) was maintained at 0.1 in all scenarios. Model results from both
investigations are shown in Figure 30. Investigation I results are outlined in red in the
figure, while results from Investigation II are outlined in green.
79
?
L
=0.005 cm ?
L
=0.01 cm ?
L
=0.015 cm ?
L
=0.1 cm ?
L
=0.2 cm
C/C
0
Figure 30 ? Longitudinal dispersivity sensitivity results (time is dimensionless time
after the end of injection) RED outlines results for Investigation I and GREEN
outlines results for Investigation II
0
0.25 0.50 0.75 0.999
T
D
= 0.11
BASE-CASE
T
D
= 0.20
T
D
= 0.36
T
D
= 0.67
80
5.5.4.1 Dispersivity Sensitivity ? Investigation I
Results from Investigation I (left three columns of Figure 30 outlined in red) show
that, qualitatively, vertical bulk velocities of the plumes appear to be the same. At earlier
times, the amount of fingering appears to increase with decreasing dispersivity, but this
distinction is lost at later times due to the complexity of the plume geometry. Differences
in plumes can be noted as early as T
D
= 0.11. For consistency, the times in this figure
were dimensionalized as discussed in the previous sections on density and conductivity
sensitivity. However, since dispersivity was not determined to be an integral parameter
in estimating characteristic plume travel time, these equivalent dimensionless times are
also equal real elapsed times from the end of plume injection. In other words, the
characteristic time used to dimensionalize the times in Figure 30 is the same as the
characteristic time for the base-case given in Table 10. It appears from Figure 30 that the
dispersivity does play a slight role in affecting the lateral width and vertical length of the
plume. This can be seen by the decrease in horizontal width and increase in vertical
length as the dispersivity increases. This is most evident in the figure in times T
D
= 0.36
and T
D
= 0.67. Again, as reported with previous parameter sensitivity analysis, the
complexity of the plume fingering prevented any conclusions from being drawn in terms
of effects of dispersivity on the amount of mixing with ambient water.
The center of mass locations for the dispersivity sensitivity scenarios are reported
in Table 11. Again the base-case plume exhibited the most lateral movement of the
center of mass with all lateral movements staying within 2% of the initial location of 27
cm. The vertical coordinates of the center of mass for each plume at time T
D
= 0 were
the same, showing that dispersivity played no role in vertical movement during injection.
Table 11 ? Center of mass coordinates for base-case and dispersivity scenarios ?
Investigation I
?
L
= 0.005 cm ?
L
= 0.015 cm BASE-CASE
z z zT
D
(cm) (cm) (cm) (cm) (cm) (cm)
81
x x x
0 27.0446 6.18633 27.0444 6.17697 27.0442 6.16938
0.11 27.0392 9.19715 27.0525 9.10062 27.0482 9.05109
0.2 27.0152 11.7664 27.0662 11.6627 27.0775 11.5487
0.36 26.9117 16.0642 27.1379 16.0901 27.0951 15.985
0.67 26.7091 22.6328 27.4894 23.2646 27.1122 22.9753
The vertical bulk velocity rates for each of the perturbed dispersivity scenarios
were calculated and are reported in Figure 31. Note, as stated before, the true elapsed
times for all plumes are the same. Therefore, since the vertical bulk velocities are the
same, for the same elapsed times, it can be deduced that dispersivity does not play an
important role in determining the vertical movement of the plume. Furthermore, while
changes in the dispersivity values do not enhance or deter vertical advection, changes in
the density contrast and hydraulic conductivity do.
Figure 31 ? Z-coordinate of center of mass versus dimensionless time for
dispersivity scenarios (Investigation I)
Vertical and horizontal variances were also calculated for the dispersivity
scenarios. Overall, for both directions, values for variance increased with time. In
regards to variance in the x direction (Figure 32), little differences between plumes are
observed. However, at later times, an increase in dispersivity caused a slight decrease in
variance, while a decrease in dispersivity caused a slight increase in variance. This was
also detected qualitatively in Figure 30 by comparing horizontal and vertical spread of
each plume. However, this slight difference may not be statistically significant and may
lie within the error produced by the randomness inherent in MOC.
82
Figure 32 ? Horizontal variance over time for dispersivity scenarios
(Investigation I)
With respect to changes in vertical variance, less difference can be noted between
model runs in this direction than in the horizontal direction. Again, at early times, no
difference in plume variance is observed. At later times, minor differences are observed
although these may not be statistically significant.
Figure 33 ? Vertical variance over time for dispersivity scenarios (Investigation I)
83
Overall, it can be deduced that vertical bulk velocity and spatial variance of the
plume are not sensitive to changes in longitudinal and transverse dispersivities (while
maintaining a constant ratio for
84
L
T
?
?
of 0.1).
5.5.4.2 Dispersivity Sensitivity ? Investigation II
This second dispersivity sensitivity investigation was undertaken because the
value of dispersivity used in the base-case model was small. It was therefore of interest
to compare results from the base-case to plumes with
L
? equal to the grain-size and
larger. Again, the ratio of transverse to longitudinal dispersivity was maintained at 0.1 in
these model runs. The results from these model outputs are shown in Figure 30 (outlined
in green).
The qualitative results show stark differences between plumes as early as T
D
=
0.11. The plumes for the two larger dispersivity values are characterized by a mixing
zone ring surrounding the freshwater. Also, the development of fingers in the base-case
scenario at this time are absent in the larger dispersivity plumes. With the base-case
scenario, these fingers developed and elongated with time to create a complex geometry
of the plume by the time it reached the overlying confining unit. This fingering pattern
also caused preferential pathways of vertical plume transport, and with these pathways,
symmetry around the centerline was lost. With the plumes with the larger dispersivities,
this fingering is not present, but is replaced by a characteristic symmetry around the
vertical centerline of the plume. With the presence of the mixing zone bands surrounding
the plumes with larger dispersivities, and due to their more simple geometries, it is easy
to see that less ?fresh? water reaches the confining unit of the plume with
L
? = 0.2 cm
than with the plume with the grain-size dispersivity of
L
? = 0.1 cm. It should also be
noted that neither of these two plumes have reached the overlying confining layer by T
D
= 0.67 which is not the case with the plumes with smaller dispersivity values. Therefore,
it can be inferred from these results that larger dispersivity values play a role in affecting
both bulk velocity rates and mixing.
Table 12 displays the center of mass locations for the larger dispersivity and base-
case scenarios. The last time of TD = 0.67 for the base-case is the largest lateral
deviation of the center of mass of the plume but still remains within 2% of the initial
location of 27 cm. Deviations between vertical locations of the center of mass at the end
of injection (T
D
= 0) were minor. This coincides with results from Investigation I where
dispersivity was not shown to have an effect on the vertical movement of the plume
during injection.
Table 12 ? Center of mass coordinates for base-case and larger dispersivity
scenarios ? Investigation II
BASE-CASE ?
L
= 0.1 cm ?
L
= 0.2 cm
z z zT
D
(cm) (cm) (cm) (cm) (cm) (cm)
85
x x x
0 6.17697 27.044 6.12882 27.0437 6.11729 27.0444
0.11 27.0525 9.10062 27.046 8.65093 27.0458 8.4031
0.2 11.6627 27.0478 10.7367 27.0475 10.2665 27.0662
0.36 27.1379 16.0901 27.0494 14.324 27.049 13.3725
0.67 27.4894 23.2646 27.1746 20.495 27.0525 18.7089
The vertical bulk velocity rates for results from Investigation II are shown in
Figure 34. Again, since dispersivity is not present in the formula used for determining ,
the dimensionless times are equivalent to real elapsed times from the cessation of
c
t
injection. The results in Figure 34 show a marked difference from results displayed in
Investigation I in reference to bulk velocity (Figure 31). The results from the ? 50%
perturbations around the base-case value of
L
? = 0.01 cm showed no changes in bulk
velocity rates. Figure 34, however, demonstrates that bulk velocity is sensitive to
changes in dispersivity at larger dispersivity values. With an increase in local
dispersivity, a decrease in bulk velocity ? indicated by the slope of the lines ? is
observed. This trend was not identified in Investigation I.
Figure 34 ? Z-coordinate of center of mass versus dimensionless time for larger
dispersivities (Investigation II)
Vertical and horizontal variances about the center of mass were calculated for
Investigation II outputs and are shown in Figures 35 and 36. First, it should be noted that
results from time T
D
= 0.67 where kept in these figures because the larger dispersivity
plumes had not yet reached the boundary. With all horizontal variance results for other
parameter perturbations, an increase in variance was observed over time. However,
86
results in Figure 35 demonstrate that with larger dispersivity values, the horizontal
variance does not change with time. This can be attributed to a lack of fingering in the
model outputs when compared to the base-case. As seen previously, at early times, no
difference in x direction variance can be noted between model runs. However, by time
T
D
= 0.36, a noticeable difference can be seen between model runs. This difference
increased even more by the time of T
D
= 0.67 because the variance for the base-case
model continued to increase, while the variances for the larger dispersivity plumes
remained around 5 cm
2
.
Figure 35 ? Horizontal variance over time for larger dispersivity values
(Investigation II)
In the vertical direction (Figure 36), we see the familiar trend of an increase in
variances for all values of dispersivity over time. Results for most dispersivities seemed
to follow the same trend except for
L
? of 0.2 cm which exhibited more variance in this
direction at all times. This could be attributed to more mixing occuring in the plume with
L
? = 0.2 cm, and this mixing occuring in the vertical direction more than in the
87
horizontal one as noted by increased tailing in Figure 30. Overall, there was much less
difference observed between vertical variance values for these dispersivities than seen
with the lateral, or horizontal, variances. In fact, the differences observed between
variances in Figure 36 may not be statistically significant when taking into account the
randomness of MOC.
Figure 36 ? Vertical variance over time for larger dispersivity values
(Investigation II)
In general, for larger dispersivity values, vertical bulk velocities were found to be
sensitive to changes in longitudinal and transverse dispersivity. An increase in
L
? and
T
? caused a decrease in bulk velocity. With larger dispersivities, the plume became more
symmetric around the vertical centerline and all fingering disappeared. Also, the
horizontal spatial variance about the center of mass was found to remain the same as
from the end of injection for larger dispersivities, instead of increasing in time. Vertical
variance was not shown to be very sensitive to these large dispersivity changes.
88
89
CHAPTER 6
SUMMARY AND CONCLUSIONS
The objective of this study was to investigate the fate and transport of pulse-
source buoyant freshwater contaminate plumes in saltwater aquifers under static
conditions. This was completed by performing both physical experiments in two-
dimensional tanks and by numerical modeling using SEAWAT.
The physical models were completed using a two dimensional tank under static
ambient flow conditions packed with glass beads. The experiments consisted of two
single injection set-ups (one in a fully saline confined aquifer and another in a density-
stratified system with freshwater overlying saline water) and a multi-injection confined
saline aquifer scenario. From the physical models, it can been that a buoyant plume in a
fully confined system will experience more lateral spreading, once coming into contact
with the confining unit, than will a plume coming into contact with a freshwater layer of
the same permeability. From the multi-injection experiment, it was observed that small
ambient density changes from previous freshwater injections have the effect of
influencing the movement of the plume by creating a more preferential pathway. Also, it
was noticed from the physical modeling results that changes in density can change the
time taken for the plume to traverse from the injection point to the overlying confining
layer. Overall, the buoyant plumes in all scenarios exhibited similar behavior in terms of
fingering and transport dynamics.
The numerical modeling showed that SEAWAT is capable of modeling these
types of scenarios. Also, results from SEAWAT compared well to results from another
widely known numerical code, SUTRA. The plume behavior was shown to be sensitive
to both grid resolution and advection solver and MOC on the fine grid was chosen as the
best solution for a good qualitative match to the physical results. From dimensional
analysis important problem parameters were derived and a selected number of these
parameters were investigated for sensitivity. The parameter sensitivity analysis consisted
of both a qualitative and quantitative analysis. The quantitative analysis was based on
moment analysis and involved calculations of bulk velocities and spatial variances about
the center of mass of each plume. The results from the parameter sensitivity analysis
were compared to a dimensionless time ( ) derived through the dimensional analysis.
The time was dimensionalized based on a characteristic time
D
t
?Kdt
c
= for each plume
scenario. The vertical bulk velocity rates of the buoyant plumes were observed to be
most sensitive to hydraulic conductivity (K) and the density contrast (? ) between
ambient and freshwater. Results from both density contrast and conductivity sensitivity
analysis scaled well when viewed in terms of this dimensionless time. Spatial variances
about the center of mass in both the horizontal and vertical directions were not found to
be sensitive to changes in any of the parameters. However, for large dispersivity values
(as observed in dispersivity sensitivity Investigation II) spatial variance in the horizontal
direction was shown to be sensitive to changes in dispersivity. The bulk velocity was
90
91
also shown to be sensitive to changes in dispersivity ? when dispersivities are large ? as
increases in dispersivity caused retardation in vertical transport.
92
CHAPTER 7
FUTURE WORK
Due to the results obtained from the sensitivity results, it was determined that
perhaps the randomness of particles in MOC does not yield trustworthy results of
qualitative analysis. This more random nature of MOC could have influenced results in
the subsequent quantitative analysis. It is recommended that sensitivity analysis, both
qualitative and quantitative, be completed using another advection solver, such as TVD,
to gain a more thorough understanding and to compare to the MOC results.
In an experimental sense, it is of interest to carry out more robust experiments
than these general ?proof of principle? experiments in the lab. For instance, in the real
world, it is rare that injection occurs as a pulse source. To more realistically understand
these scenarios, it is further recommended that both laboratory and numerical
investigations of behavior of buoyant plumes under constant injection scenarios be
undertaken.
Another significant drawback to the experiments discussed in this thesis is that the
actual concentration isochlors could not be determined from physical model results.
Therefore, application of an image analysis technique to the data images provided in this
body of work, or to subsequent similar experiments, would greatly aid in calibration of
93
the numerical model and further understanding of plume transport dynamics. Another
possibility is to use gamma radiation to determine concentration values at different points
in the experimental domain.
From numerical modeling of this problem, it was discovered that to best calibrate
to the physical results, again with the use of actual concentration contours, a
homogeneous, isotropic model could not really be used. Better calibration results may be
obtained with a stochastically generated permeability field with small-scale
heterogeneities.
In regards to MOC, it would be also be of future interest to carry out a detailed
particle number sensitivity analysis as this was not found to be done in the literature.
From iterations of different numerical models being built and run until finally converging
on a set base-case model scenario, it was discovered that results for MOC seem to be
sensitive to NPMIN, NPMAX, NPHIGH, and NPLOW parameters chosen. These
parameters control how the particle allocation is handled in the presence of a boundary,
in the presence of sinks and sources, and under different concentration gradients.
With changes in these parameters, it was noted that different qualitative results
could be obtained, which is not desirable. However, a quantitative analysis could shed
light on how much of an effect these parameters play in results. This is highly
recommended before future use of MOC as an advection solver in quantitative analysis.
The quantitative analysis of spatial variances was shown to give an estimate of the
variance of the concentration distribution about the center of mass in a particular
94
direction. However, this did not yield definitive results in terms of the amount of
spreading and mixing of the plume as hoped. Therefore, future work should involve a
better quantification of the spreading and mixing of the plume in and with the ambient
water.
95
REFERENCES
Al-Maktoumi, A., et al. (2007), SEAWAT 2000: modelling unstable flow and sensitivity
to discretization levels and numerical schemes, Hydrogeology Journal, 15, 1119-1129.
Anderson, M. P., and W. W. Woessner (1992), Applied Groundwater Modeling:
Simulation of Flow and Advective Transport, Academic Press, Inc., San Diego,
California.
Bear, J., and M. Jacobs (1965), On the movement of water bodies injected into aquifers,
Journal of Hydrology, 3, 37-57.
Chapra, S. C., and R. P. Canale (2002), Numerical Methods for Engineers with Software
and Programming Applications, McGraw Hill, New York, NY.
Degan, G., et al. (2003), Buoyant plume above a line source of heat in an anisotropic
porous medium, Heat and Mass Transfer, 39, 209-213.
Dorgarten, H. W., and C. F. Tsang (1991), Modeling the Density-Driven Movement of
Liquid Wastes in Deep Sloping Aquifers, Ground Water, 29, 655-662.
EPA (2007), Underground Injection Control Program: Classes of Wells, edited, U. S.
Environmental Protection Agency.
Fetter, C. W. (1988), Applied Hydrogeology, 2nd Ed. ed., Merrill Publishing Company,
Columbus, OH.
Fischer, H. B., et al. (1979), Mixing in Inland and Coastal Waters, Academic Press, New
York, New York.
Goswami, R. R., and T. P. Clement (2007), Laboratory-scale investigation of saltwater
intrusion dynamics, Water Resources Research, 43.
Guo, W., and C. D. Langevin (2002), User's Guide to SEAWAT: A Computer Program
for Simulation of Three-Dimensional Variable-Density Ground-Water Flow, U. S.
Geological Survey, Tallahassee, Florida.
Hogan, M. B. (2006), Understanding the Flow and Mixing Dynamics of Saline Water
Discharged into Coastal Freshwater Aquifers, Auburn University, Auburn, AL.
96
Hughes, J. D., et al. (2005), Numerical simulation of double-diffusive finger convection,
Water Resources Research, 41.
Illangasekare, T., et al. (2006), Impacts of the 2004 tsunami on groundwater resources in
Sri Lanka, Water Resources Research, 42.
Jin, P. K., et al. (1996), Numerical simulation of deep-well injection in Geismar,
Louisiana, Ground Water, 34, 989-1000.
Johannsen, K., et al. (2002), The saltpool benchmark problem - numerical simulation of
saltwater upconing in a porous medium, Advances in Water Resources, 25, 335-348.
Kitanidis, P. K. (1994), The Concept of the Dilution Index, Water Resources Research,
30, 2011-2026.
Langevin, C. D., and W. X. Guo (2006), MODFLOW/MT3DMS-based simulation of
variable-density ground water flow and transport, Ground Water, 44, 339-351.
Liu, H. H., and J. H. Dane (1997), A numerical study on gravitational instabilities of
dense aqueous phase plumes in three-dimensional porous media, Journal of Hydrology,
194, 126-142.
Maliva, R. G., et al. (2007), Vertical migration of municipal wastewater in deep injection
well systems, South Florida, USA, Hydrogeology Journal, 15, 1387-1396.
McDonald, M. G., and T. E. Reilly (2007), Models of ground water systems - Not just
tools but components in a scientific approach, Ground Water, 45, 1-1.
Muniz, A., et al. (2005), Why current regulations protect Florida's subsurface
environment, in Underground Injection: Science and Technology, edited by C.-F. Tsang
and J. A. Apps, Elsevier,
Amsterdam, The Netherlands.
Oostrom, M., et al. (1992a), Experimental Investigation of Dense Solute Plumes in an
Unconfined Aquifer Model, Water Resources Research, 28, 2315-2326.
Oostrom, M., et al. (1992b), Behavior of Dense Aqueous Phase Leachate Plumes in
Homogeneous Porous-Media, Water Resources Research, 28, 2123-2134.
Paschke, N. W., and J. A. Hoopes (1984), Buoyant Contaminant Plumes in Groundwater,
Water Resources Research, 20, 1183-1192.
Perlman, H. (2006), Ground water use in the U. S. , edited, U. S. Geological Survey.
97
Peterson, F. L., et al. (1977), Waste Injection into Stratified Groundwater Bodies, Water
& Sewage Works, 124, 60-65.
Peterson, F. L., et al. (1978), Waste Injection into a 2-Phase Flow Field - Sand-Box and
Hele-Shaw Models Study, Ground Water, 16, 410-416.
Pringle, S. E., et al. (2002), Double-diffusive Finger Convection in a Hele?Shaw Cell: An
Experiment Exploring the Evolution of Concentration Fields, Length Scales and Mass
Transfer, Transport in Porous Media, 47, 195-214.
Richardson, S. D., et al. (2004), Use of rhodamine water tracer in the marshland
upwelling system, Ground Water, 42, 678-688.
Schincariol, R. A., and F. W. Schwartz (1990), An Experimental Investigation of
Variable Density Flow and Mixing in Homogeneous and Heterogeneous Media, Water
Resources Research, 26, 2317-2329.
Schincariol, R. A., et al. (1994), On the Generation of Instabilities in Variable-Density
Flow, Water Resources Research, 30, 913-927.
Spitz, K., and J. Moreno (1996), A Practical Guide to Groundwater and Solute Transport
Modeling, John Wiley & Sons, Inc.
Traylor, D. H. (1991), Stability of Fluid Interfaces in Ground Water: Physical Modeling
and Analysis, Auburn University, Auburn, AL.
Voss, C. I., and A. M. Provost (2003), SUTRA: A Model for Saturated-Unsaturated,
Variable-Density Ground-Water Flow with Solute or Energy Transport, U.S. Geological
Survery, Reston, VA.
Voss, C. I., and W. R. Souza (1987), Variable Density Flow and Solute Transport
Simulation of Regional Aquifers Containing a Narrow Fresh-Water-Saltwater Transition
Zone, Water Resources Research, 23, 1851-1866.
Wikipedia (2008), Dimensional Analysis, edited.
Zheng, C., and G. D. Bennett (1995), Applied Contaminant Transport Modeling: Theory
and Practice, Van Nostrand Reinhold.