Dynamics of Saltwater Intrusion Processes in Saturated Porous Media Systems
by
Sun Woo Chang
A dissertation submitted to the Graduate Faculty of
Auburn University
in partial fulfillment of the
requirements for the Degree of
Doctor of Philosophy
Auburn, Alabama
December 8, 2012
Keywords: saltwater intrusion, density-dependent flow, climate change
Copyright 2012 by Sun Woo Chang
Approved by
T. Prabhakar Clement, Committee Chair, Arthur H. Feagin Professor of Civil Engineering
Ming-Kuo Lee, Committee member, Professor of Geology
Xing Fang, Committee member, Associate Professor of Civil Engineering
Jose G. Vasconcelos, Committee member, Assistant Professor of Civil Engineering
ii
Abstract
The overarching goal of this dissertation is to contribute towards a better understanding
of the dynamics of saltwater intrusion and the associated transport processes in coastal
groundwater systems. As a part of this effort, we have also investigated the impacts of various
climate-change impacted variables such as average sea level and regional recharge fluxes on
saltwater intrusion processes. Climate change effects are expected to substantially raise the
average sea level. It is widely assumed that sea-level rise would severely impact saltwater
intrusion processes in coastal aquifers. In the first phase of this study we hypothesize that a
natural mechanism, identified here as the ??lifting process,?? has the potential to mitigate, or in
some cases completely reverse, the adverse intrusion effects induced by sea-level rise. A
detailed numerical study using the MODFLOW-family computer code SEAWAT was completed
to test the validity of this hypothesis in both confined and unconfined systems. Our conceptual
simulation results show that if the ambient recharge remains constant, the sea-level rise will have
no long-term impact (i.e., it will not affect the steady-state salt wedge) on confined aquifers. Our
transient confined-flow simulations show that a self-reversal mechanism, which is driven by the
lifting of the regional water table, would naturally drive the intruded saltwater wedge back to the
original position. In unconfined systems, the lifting process would have a lesser influence due to
changes in the value of effective transmissivity. A detailed sensitivity analysis was also
completed to understand the sensitivity of this self-reversal effect to various aquifer parameters.
iii
The outcomes of the first phase of this research indicated that the changes in groundwater
fluxes due variations in rainfall patterns is one of the major climate-change-induced hydrological
variable that can impact saltwater intrusion in coastal aquifers. In the second phase of this study,
we use a combination of laboratory experiments and numerical simulations to study the impacts
of changes in various types of groundwater fluxes on saltwater intrusion dynamics. We have
completed experiments in a laboratory-scale sand tank model to study the changes in two types
of groundwater fluxes? areal-recharge flux and regional flux. The experimental results were
modeled using the numerical code SEAWAT. The transient datasets collected in this
experimental study are found to be useful benchmark data for testing numerical models that
employ flux-type boundary conditions. Also, based on the experimental observations we
hypothesized that when the fluxes are perturbed it would require relatively less time for a salt
wedge to recede from an aquifer when compared to the time required for the wedge to advance
into the aquifer. This rather counter-intuitive hypothesis implies that saltwater intrusion and
receding processes are asymmetric and the time scales associated with these processes will be
different. We use a combination of laboratory and numerical experiments to test this hypothesis
and use the resulting dataset to explain the reason for the difference in salt wedge intrusion and
recession time scales.
In coastal aquifers, presence of a salt wedge divides the groundwater flow system into
two distinct regions which includes a freshwater region above the wedge and saltwater region
below the wedge. Typically, the freshwater transport processes occurring above a wedge are
much faster than the saltwater transport processes occurring beneath the wedge. Recently, many
modeling and laboratory studies have investigated the movement of salt wedges and the
associated transport processes. Most of these transport studies, however, have focused on
iv
understanding the groundwater plume transport above a wedge. As per our knowledge, so far no
one has completed controlled laboratory experiments to study the transport processes occurring
beneath a saltwater wedge. In this study, we have completed contaminant transport experiments
to map the dynamics of saltwater flow patterns beneath a wedge and relate it to the freshwater
flow patterns present above the wedge. We used a novel experimental approach that employed
variety dyed neutral density tracers to map the mixing and transport processes occurring above
and below a salt wedge. The experimental datasets were simulated using the SEAWAT code.
The model was then used to investigate contaminant transport scenarios occurring beneath a
saltwater wedge in larger field-scale problems.
v
Acknowledgments
I would like to acknowledge my advisor, Dr. T. Prabhakar. Clement for his support and
mentorship. My achievements are evidence of his mentorship and only a small portion of why
he is a named outstanding mentor. I would like to thank to my committee members, Dr. Ming-
Kuo Lee, Dr. Xing Fang, Dr. Jose. G. Vasconcelos and also thank to Dr. Joel G. Melville. I
would like to thank to Dr. Mark Dougherty for giving useful suggestions and for careful edit of
my dissertation. I also thank to Dr. Kang-Kun Lee, Dr. Matthew Simpson, Dr. Rohit Goswami,
and the reviewers of Water Resources Research and Advances in Water Resources. This work
was, in part, supported by the Samuel Ginn College of Engineering Dean Fellowship and by the
Civil Engineering Teaching Assistantship. Finally, I would like to dedicate this work to my
family.
vi
Table of Contents
Abstract ...................................................................................................................................................... ii
Acknowledgments .................................................................................................................................... v
List of Tables ................................................................................................................................ x
List of Illustrations ....................................................................................................................... xi
Chapter 1
INTRODUCTION .................................................................................................................................... 1
1.1 Fundamentals of Saltwater Intrusion ........................................................................ 1
1.2 Impacts of Climate Change on Saltwater Intrusion Processes.................................... 4
1.3 Objectives ................................................................................................................... 5
Chapter 2
INVESTIGATION OF THE IMPACT OF SEA-LEVEL RISE ON SALTWATER INTRUSION
PROCESSES - A CONCEPTUAL MODELING STUDY ................................................... 8
2.1 Introduction ............................................................................................................... 8
2.2 Problem Formulation and Conceptual Modeling ..................................................... 11
2.3 Details of Numerical Experiment ............................................................................ 14
2.4 Results ...................................................................................................................... 17
2.4.1 Impacts of Sea-level Rise on Confined Flow Conditions ......................... 17
2.4.2 Sensitivity Analysis of the Self-reversal Process in Confined Systems .... 23
2.4.3 Impacts of Sea-level Rise on Unconfined Aquifers ................................... 29
vii
2.5 Conclusions .............................................................................................................. 35
Chapter 3
INVESTIGATION OF SALTWATER INTRUSION PROCESSES USING LABORATORY
EXPERIMENTS INVOLVING FLUX-TYPE BOUNDARY CONDITIONS .............. 37
3.1 Introduction ............................................................................................................. 37
3.2 Method .................................................................................................................... 40
3.2.1 Experimental Approach ............................................................................ 40
3.2.2 Numerical Modeling Approach ................................................................ 43
3.3 Results ...................................................................................................................... 45
3.3.1 Regional Flux Experiment ........................................................................ 45
3.3.2 Areal-recharge Flux Experiment ............................................................... 50
3.3.3 Analysis of Time Scales Associated with Intruding and Receding Salt
Wedge Transport Processes ............................................................................... 54
3.4 Conclusions ............................................................................................................... 64
Chapter 4
LABORATORY AND NUMERICAL INVESTIGATION OF TRANSPORT PROCESSES
OCCURING BENEATH A SALTWATER WEDGE. .................................................. 66
4.1 Introduction ............................................................................................................. 66
4.2 Methods .................................................................................................................... 72
4.2.1 Experimental Method................................................................................. 72
4.2.2 Details of Tracer Transport Experiments ................................................... 73
4.2.2 Numerical Modeling Method ..................................................................... 74
4.3 Results .................................................................................................................................78
4.3.1 Transport of Tracer Slugs Injected above the Wedge ................................ 78
4.3.2 Characterizing Tracer Transport Patterns beneath a Saltwater .................. 81
viii
4.3.3 Sensitivity of Recirculating Saltwater Flux to Dispersion and the Type of
Boundary Condition ................................................................................... 86
4.4 Conclusions ............................................................................................................. 91
Chapter 5
CONCLUSIONS AND RECOMMENDATIONS ..................................................................... 93
6.1 Summary and Conclusions ....................................................................................... 93
6.2 Recommendation for Future Work ........................................................................... 95
Appendix 1 Transient Variations in the Toe Position (XT) for 10%, 50% and 90% Isochlor of
Saltwater Wedge ............................................................................................................. 97
Appendix 2 Comparison of Salt-wedge between Central-difference and Upstream-weighting for
density-weighting scheme ............................................................................................... 98
Appendix 3 Comparison of Salt-wedge between FDM and TVD scheme ................................ 99
Appendix 4 Analytical Comparison of Salt-wedge Toe Position in Unconfined and Confined
Aquifers......................................................................................................................... 100
Appendix 5 Derivation of Dimensionless of Henry problem Parameters and Application to
Modified Henry Problem with Areal-recharge flux ...................................................... 102
Appendix 6 Development of Two-dimensional Numerical Code for Simulating Saltwater
Transport in Groundwater Systems .............................................................................. 106
Appendix 7 Code Formulation based on TVD scheme for Advection Package ...................... 128
References ................................................................................................................................ 133
ix
List of Tables
Table 2.1 Aquifer properties and hydrological properties used for the base-case problem ....... 16
Table 3.1 Summary of numerical model parameters .................................................................. 44
Table 3.2 Summary measured and modeled flows ..................................................................... 44
Table 4.1 Model parameters used for simulating the lab problem ............................................. 75
Table 4.2 Model parameters used for simulating the field problem (modified from Chang et al.
2011) ............................................................................................................................... 76
x
List of Illustrations
Figure 1.1. Groundwater flow patterns near a pumping well in a coastal aquifer (USGS, 2000) 3
Figure 2.1. Comparison of conceptual models used for visualizing the impacts of sea-level rise
on a saltwater wedge: (a) salt wedge profile after sea-level rise based on a traditional
conceptual model that ignores the lifting effect and (b) a new conceptual model that
includes the lifting effect. ......................................................................................... 13
Figure 2.2. Comparison of steady-state salt wedges predicted before and after the sea-level rise in
the confined system................................................................................................... 20
Figure 2.3. Transient salt wedge profiles for the confined system ............................................. 21
Figure 2.4. Transient variations in the toe position (XT) and freshwater head (hf) at the inland
boundary. .................................................................................................................. 22
Figure 2.5. Sensitivity of the toe position (XT) to (a) specific storage, (b) the magnitude of sea-
level rise, (c) the rate of sea-level rise velocity and (d) the dispersivity coefficients ..
................................................................................................................................... 26
Figure 2.6. Sensitivity of the toe position (XT) to: (a) hydraulic conductivity and (b) freshwater
flux. ........................................................................................................................... 28
Figure 2.7. Comparison of steady-state salt wedges predicted before the sea-level rise (sea level
at 30m) in the unconfined- and confined-flow systems ............................................ 32
Figure 2.8. Comparison of steady-state salt wedges predicted before the sea-level rise in the
unconfined system, after the sea-level rise in the unconfined system, and in an
equivalent 34-m thick confined-flow system ............................................................ 33
xi
Figure 2.9. Comparison of transient variations in toe positions (XT) predicted for the confined
and unconfined flow systems.. .................................................................................. 34
Figure 3.1. Schematic diagram of the (a) regional-flux and (b) areal-recharge flux experiments
and (c) photograph of laboratory setup. .................................................................... 42
Figure 3.2. Transient variations in the salt-wedge patterns due to changes in the regional flux at
the right boundary. The top figure shows intruding transport conditions when the flux
was reduced from 0.833 cm3/s to 0.467 cm3/s, and the bottom figure shows receding
transport conditions when the flux was increased from 0.467 cm3/s to 0.833 cm3/s .
................................................................................................................................... 47
Figure 3.3. Comparison of experimental data model predictions of transient salt wedges for the
regional-flux system: a) advancing-wedge experiment, and b) receding-wedge
experiment................................................................................................................. 48
Figure 3.4. Comparison of experimental and model-simulated toe positions estimated at 0.75 cm
above tank bottom in the regional-flux experiment. ................................................. 49
Figure 3.5. Transient variations in the salt-wedge patterns due to changes in the areal flux values.
The top figure shows intruding transport conditions when the flux was reduced from
1.133 cm3/s to 0.733 cm3/s, and the bottom figure shows receding transport
conditions when the flux is increased from 0.733 cm3/s to 1.133cm3/s ................... 51
Figure 3.6. Comparison of experimental data model predictions of transient salt wedges for the
areal-flux system: a) advancing-wedge experiment, and b) receding-wedge
experiment................................................................................................................. 52
Figure 3.7. Comparison of experimental and model-simulated toe positions estimated at 0.75 cm
above tank bottom in the areal-flux experiment ....................................................... 53
Figure 3.8. Comparison of toe migration rates for intruding and receding wedges for laboratory-
scale problem of (a) regional flux experimental problem and (b) field scale problem
(reported in Chang et al. 2011) ................................................................................. 57
xii
Figure 3.9. Comparison of toe migration rates for intruding and receding wedges for laboratory-
scale problem of areal-recharge flux experimental problem .................................... 58
Figure 3.10. Simulated steady-state velocity fields for the regional flux experiment (a) SS1 and
(b) SS2 ...................................................................................................................... 61
Figure 3.11. Simulated transient velocity fields for the regional flux experiment (a) 2 minutes
after SS1 and (b) 2 minutes after SS2 ....................................................................... 62
Figure 3.12. Visualization of the velocity fields near the salt wedge plotted in log scale at (a) SS1,
(b) 2 minutes after SS1, (c) SS2, and (d) 2 minutes after SS2 .................................. 63
Figure 4.1. Conceptual model illustrating various types of contaminant transport processes
occurring in a coastal groundwater system ............................................................. 67
Figure 4.2. Schematic diagram of the physical model used in this study ................................. 77
Figure 4.3. Comparison of observed data with numerical predictions for the neutral plume
transport experiment (Top figure: lab data; Bottom figure: SEAWAT predictions
with color filled contour -- red > 20 %; yellow 20-14 %; green 14-8 %; and turquoise
8-1 %; also, white line is the 50% saltwater contour of the wedge ......................... 80
Figure 4.4. Comparison of observed data with numerical predictions for the dense plume
transport experiment (Top figure: lab data; Bottom figure: SEAWAT predictions
with color filled contour -- red > 20 %; yellow 20-14 %; green 14-8 %; and turquoise
8-1 %; also, white line is the 50% saltwater contour of the wedge ......................... 80
Figure 4.5. Comparison of observed data with numerical predictions for the transport experiment
completed beneath the wedge (Top figure: lab data; Bottom figure: SEAWAT
predictions with color filled contour -- red > 20 %; yellow 20-14 %; green 14-8 %;
and turquoise 8-1 %; also, white line is the 50% saltwater contour of the wedge ... 84
Figure 4.6. Model simulated results illustrating the transport processes occurring above and
below the saltwater wedge for the field-scale test problem: (a) simulations for low
xiii
regional flow of 0.105 m3/day and (b) simulations for high regional flow of 0.210
m3/day. Initial contaminant concentration at the sources was 100 kg/m3 ................ 85
Figure 4.7. Relationship between freshwater flow and recirculate-saltwater flow for regional and
areal recharge problems with different dispersivity values ...................................... 89
Figure 4.8. Non-dimensional analysis of the relationship between freshwater and saltwater fluxes
[The symbols indicate the following literature data points: + from Chang et al. (2011)
study, square6 regional flux experimental data from Chang and Clement (2012), box3 areal
flux experiment data from Chang and Clement (2012), triangleup Cyprus case study from
Destouni and Prieto (2003), triangleopenup Israel case study from Destouni and Prieto (2003)]
................................................................................................................................... 90
1
Chapter 1
INTRODUCTION
1.1 Fundamentals of Saltwater Intrusion
Groundwater is an important freshwater resource used extensively for water supply.
Groundwater also supports industrial activities in many areas that have limited number of surface
water resources. Unlike a surface water reservoir, which is a concentrated source that can serve
only a narrow region, groundwater is a distributed source which can serve populations
distributed over a larger region. If the areas where groundwater is present coincided with areas
of demand, groundwater can be made available to the local public without any major investment
in water transmission infrastructure. Groundwater is also a relatively stable source since climate
fluctuations induce relatively small variations in groundwater levels in most cases; large volumes
of water stored in groundwater aquifers may serve as a buffer to supply water during lengthy
drought periods [Bear and Cheng, 2010] .
Currently, within the US, it has been estimated that groundwater aquifers located along
the Atlantic Coast supply drinking water to over 30 million residents living in coastal towns from
Maine to Florida [USGS, 2000]. These precious coastal aquifers are constantly being
contaminated by saline water through a natural process known as saltwater intrusion. Figure 1.1
illustrates the dynamics of groundwater flow and mixing of freshwater and saltwater in a coastal
aquifer. As shown in the figure, the fresh groundwater region is in direct contact with the saline
seawater region in coastal aquifers. Under natural conditions, a dynamic equilibrium exists
between these two regions. There are numerous human and environmental factors that can
adversely impact this equilibrium and that lead to severe saltwater contamination of the
2
freshwater region. Once contaminated, it is difficult and expensive to clean up saltwater
contaminated aquifers. In some cases the contamination could even be irreversible [Bear and
Cheng, 2006].
Saltwater intrusion into freshwater aquifers primarily occurs due to the difference in the
density of seawater (which has a density value of 1.025 g/cm3) and fresh groundwater (which has
a density value of 1.000 g/cm3). This small density difference can play a significant role in
controlling two types of saltwater intrusion processes: lateral intrusion of seawater beneath the
regional aquifer and up-coning of seawater near pumping wells. The up-coning of seawater
normally occurs near a pumping well when a large amount of groundwater is withdrawn from
the aquifer. Anthropogenic pumping activities induce a cone of depression around the well,
which can lead to upward migration of seawater directly into pumping wells. Lateral intrusion of
seawater, on the other hand, would occur naturally when denser saltwater seeps inland at the
bottom of the freshwater aquifer (see Figure 1.1). Lateral intrusion of seawater would result in a
distinct curved interface that separate the freshwater and saltwater regions; this interface is
known as the regional ?saltwater wedge.? The shape and extent of the interface is determined by
various factors such as the geology of the aquifer, climate patterns, variations in natural
groundwater flow and the sea level, etc.
3
Figure 1.1. Groundwater flow patterns near a pumping well in a coastal aquifer (USGS, 2000)
4
1.2. Impacts of Climate Change on Saltwater Intrusion Processes
In recent years, population growth and land use changes in coastal areas have increased
the demand for freshwater; in addition, it is believed that climate change effects have further
reduced the water availability in coastal areas [Ranjan et al., 2006]. Climate change could lead
to rise in sea levels, and climate change could also alter groundwater recharge fluxes due to
changes in rainfall patterns. Droughts induced by climate change would result in reduction of
regional freshwater recharge and this would enhance saltwater intrusion. In the published
literature, climate change models have forecasted significant sea-level rise ranging from 0.18 to
0.51 m by the end of this century [IPCC, 2007]. According to Titus et al. (2009), during the
twentieth century the majority of Atlantic Coast and Gulf of Mexico Coast regions experienced
rates of sea-level rise ranging from 2 to 10 mm/year. These estimates are considerably higher
than the current global average of 1.7 mm/year, and have led to the concern that rising oceans
levels would impact the quality of both surface and ground water supplies and reduce the world?s
existing fresh water reserves. Impact of climate change on saltwater intrusion in coastal
groundwater systems is a serious environmental concern since eighty percent of the world
population resides along the coastal regions and depends on coastal aquifers for their water
supply. Most of these regions have water-scarce coastal towns that have limited freshwater
supply, and are also undergoing rapid population growth. Saltwater intrusion is an important
concern for these towns.
Researchers have used modeling studies to understand the coupling of climate change
effects with groundwater problems. Dausman and Langevin (2005) evaluated the relative
importance of drought, well-field pumping, sea-level rise, and canal management on saltwater
intrusion for the Biscayne aquifer in Southern Florida. The simulation results for a coastal
5
aquifer in Broward County, Florida, demonstrated that if the sea-level rise becomes greater than
48 cm over the next 100 years local well fields would be vulnerable to chloride contamination.
Aquifers in small barrier islands are particularly more vulnerable to climate change effects than
other aquifers because extreme events such as hurricanes can deposit abnormal amount of salt by
directly inundating these islands. Gillett et al. [2000] reported that the occasional hurricanes that
strike the Gulf Coast can cause tidal surges up to 5 to 25 feet above the normal level. These
storm surges could dump large amount of saltwater over sandy soil and contaminate the shallow
groundwater. In addition to storm surges, catastrophic events such as tsunamis can also dump
large volume of saltwater above freshwater coastal aquifers. Ilangasekare et al. (2006) discuss a
detailed field survey completed in Sri Lanka that documented saltwater intrusion caused by the
tsunami wave generated by the 2004 Indian Ocean earthquake on Sunday, 26 December 2004,
with an epicenter off the west coast of Sumatra, Indonesia. The tsunami surge not only caused
the salt wedge to temporarily move inward, but it also dumped a dense seawater directly over
sandy unconfined aquifer along the entire Sri Lankan coast. Ilangasekare et al. (2006) used a
combination laboratory analysis together with simple analytical calculations to estimate the time
taken to naturally restore these aquifers. These field studies have demonstrated the importance
of understanding impacts of climate change and other extreme hydrological events on saltwater
intrusion processes.
1.3 Objectives
The overall goal of this research is to develop a better understanding of the dynamics of
saltwater intrusion processes in coastal groundwater aquifers. This study has four research
objectives, and this dissertation has four distinct chapters that document the results of
6
investigations related to each of these objectives. Each of the four chapters is organized in a
journal format and has an introduction section and also reviews all relevant literature.
The first objective of this study is to investigate the transient impacts of the sea-level
rise on saltwater intrusion processes in confined and unconfined coastal aquifers that are
driven by natural recharge fluxes. The results of this numerical study are used to develop an
intuitive understanding for saltwater intrusion dynamics, which can help better manage the
potential long term impacts of sea-level rise on coastal aquifers.
The result of the first phase of this study demonstrated the importance of flux-type
boundary conditions. Our review indicated that past experimental efforts related to saltwater
intrusion processes have primarily focused on conducting experiments involving constant head
boundary conditions. Therefore, the second objective of this work is to conduct
experimental studies to investigate the dynamics of saltwater intrusion processes in an
unconfined aquifer using different types of flux boundary conditions. We investigated the
influence of variations in recharge patterns (resulting from changes in rainfall conditions) on a
transient movement salt wedge using a laboratory-scale aquifer model and modeled the data
using a numerical code. Both the laboratory data and the numerical simulations were used to
develop a better understanding of the dynamics of salt wedges under transient flow conditions.
Many of the shallow unconfined aquifers in coastal areas are contaminated by various
type of contaminants released from a variety of anthropogenic sources [Westbrook et al., 2005].
Furthermore, regional groundwater flow can also transport large amount of nutrients, which then
could be exchanged across the salt wedge. Therefore, the third objective of this work is to
conduct laboratory experimental studies to investigate the contaminant transport processes
near a salt wedge. We completed laboratory experiments to study tracer transport in both
7
freshwater and salt water flow zones in the vicinity of the salt wedge. The results were used to
analyze contaminant transport scenarios above and below a saltwater wedge.
This dissertation has five chapters. The first chapter (the current chapter) provides a
basic introduction to this research, lists the objectives, and provides a brief summary of each
chapter. The second chapter focuses on the first objective, with the outcome of this effort
already published in Advances in Water Resources [Chang et al., 2011]. The third chapter
focuses on the second objective and the outcome of this effort has been accepted in Water
Resources Research [Chang and Clement, 2012b]. The fourth chapter focuses on the third
objective and the outcome of this effort has been communicated to Journal of Contaminant
Hydrology [Chang and Clement, 2012a] The fifth chapter provides a brief summary of the key
outcomes of this research and offers recommendations for future work.
8
Chapter 2
INVESTIGATION OF THE IMPACT OF SEA-LEVEL RISE ON SALTWATER INTRUSION
PROCESSES - A CONCEPTUAL MODELING STUDY
2.1 Introduction
Saltwater intrusion is a serious environmental issue since eighty percent of the world?s
population live along the coast and utilizes local aquifers for their water supply. In the US alone,
it is estimated that freshwater aquifers along the Atlantic coast supply drinking water to 30
million residents living in coastal towns located from Maine to Florida [USGS, 2000]. Under
natural conditions, these coastal aquifers are recharged by rainfall events, and the recharged
water flowing towards the ocean would prevent saltwater from encroaching into the freshwater
region. However, over exploitation of coastal aquifers has resulted in reducing groundwater
levels (hence reduced natural flow) and this has led to severe saltwater intrusion. Cases of
saltwater intrusion, with varying degrees of severity and complexity, have been documented
throughout the Atlantic coastal zone. For example, in May County, New Jersey, more than 120
water supply wells have been abandoned because of saltwater contamination [Lacombe and
Carleton, 1992]. A recent USGS [2000] study provides a summary of saltwater intrusion
problems and various mitigation techniques. International organizations have identified
saltwater intrusion as one of the major environmental issues faced by several coastal cities in
India, China, and Mexico [WWD, 1998]. Researchers have also reported that variations in the
sea level and the associated wedge movement can influence the near-shore and/or large-scale
submarine discharge patterns and impact nutrient loading levels across the aquifer-ocean
9
interface [Li and Jiao, 2003b; Li et al., 2008; Li et al., 1999; Michael et al., 2005; Robinson et al.,
2007]. Therefore, understanding the dynamics of saltwater intrusion in coastal aquifers and its
interconnection to anthropogenic activities is an important environmental challenge.
While anthropogenic activities such as over pumping and excess paving in urbanized
areas are the major causes of saltwater intrusion, it is anticipated that increases in sea level due to
climate change would aggravate the problem. Nevertheless, only a few studies have focused on
understanding the combined effects of climate change and anthropogenic impacts [Li and Jiao,
2003a; b]. Feseker [2007] completed a numerical modeling study to assess the impacts of
climate change and changes in land use patterns on the salt distribution in a coastal aquifer. The
model used parameters to reflect the conditions similar to those observed at the CAT-field site
located in the northern coast of Germany. The study concluded that rising sea level could induce
rapid progression of saltwater intrusion. Furthermore, the time scale of changes resulting from
the altered boundary conditions could take decades or even centuries to impact groundwater
flows and hence the present day salt distribution might not reflect the long term equilibrium
conditions. Leatherman [1984] investigated the effects of rising sea-level on salinization in an
aquifer in Texas. Meisler et al.[1985] used a finite-difference computer model to analyze the
effect of sea-level changes on the development of the transition zone between fresh groundwater
and saltwater in the northern Atlantic Coastal Plain (from New Jersey to North Carolina). Navoy
[1991] studied aquifer-estuary interactions to assess the vulnerability of groundwater supplies to
sea level rise-driven saltwater intrusion in a coastal aquifer in New Jersey. Oude Essink [2001]
used a three-dimensional transient density-driven groundwater flow model to simulate saltwater
intrusion in a coastal aquifer in the Netherlands for three types of sea-level rise scenarios: no
rise, a sea-level rise of 0.5 m per century, and a sea-level fall of 0.5 m per century. They
10
concluded that sea-level rise of 0.5 m per century would increase the salinity in all low-lying
regions closer to the sea. Dausman and Langevin [2005] completed SEAWAT simulations for a
coastal aquifer in Broward County, Florida, and demonstrated that if the sea-level rise becomes
greater than 48 cm over the next 100 years then several local well fields would be vulnerable to
chloride contamination. Melloul and Collin [2006] evaluated the potential of sea-level rise to
cause permanent freshwater reserve losses in a coastal aquifer in Israel. They quantified the
saltwater intrusion effect due the lateral movement of seawater and due changes in the
groundwater head. For the assumed sea-level rise of 0.5 m, about 77% of the loss was due to the
lateral movement and about 23% was due to the head change. Ranjan et al. [2006] used the
sharp interface assumption to analyze the effects of climate change and land use on coastal
groundwater resources in Sri Lanka. Giambastiani et al. [2007] conducted a numerical study to
investigate saltwater intrusion in an unconfined coastal aquifer of Ravenna, Italy. Their result
showed that the mixing zone between fresh and saline groundwater will be shifted by about 800
m farther inland for a 0.475 m per century of sea-level rise. Loaiciga et al. [2011] employed
hydrogeological data and the finite-element numerical model FEFLOW to assess the likely
impacts of sea-level rise and groundwater extraction on seawater intrusion in the Seaside Area
aquifer of Monterrey County, California, USA. Sea-level rise scenarios were consistent with
current estimates made for the California coast, and varied between 0.5 and 1.0 m over the 21st
century. These authors concluded that sea-level rise would have a minor contribution to seawater
intrusion in the study area compared to the contribution expected from groundwater extraction.
The primary focus of most of the field-scale modeling studies discussed above was to
understand saltwater intrusion problems related to a specific field site. More recently, Werner
and Simmons [2009] completed a conceptual modeling study using a steady-state, sharp-
11
interface analytical model focused on developing a general understanding of the impacts of sea-
level rise on groundwater aquifers that might have different types of boundary conditions. They
used a relatively simple analytical model to provide a first-order assessment of the impacts of sea
level changes on saltwater intrusion on aquifer systems with two types of boundary conditions.
Their results showed that the level of intrusion would depend on the type of inland boundary
condition assumed in the model. The steady-state, sharp-interface, analytical expression used in
the study did not consider saltwater mixing effects and transient effects. The transient effects in
unconfined aquifers were later investigated by Webb and Howard [2010], using a numerical
model that employed constant head boundary conditions to study the changes in the rates of
intrusion for a range of hydro-geological parameters.
The objective of this study is to complete a comprehensive investigation of the transient
impacts of sea-level rise on saltwater intrusion processes in both confined and unconfined coastal
aquifer systems that are driven by natural recharge fluxes. The results are then used to develop
an intuitive understanding of saltwater intrusion dynamics, which can help better assesses and
manage the potential long term impacts of sea-level rise on coastal aquifers.
2.2 Problem Formulation and Conceptual Modeling
Under natural flow conditions, dense sea water would have the tendency to intrude
beneath fresh groundwater. The spatial extent of the intruded saltwater wedge (designated as the
?toe position XT?) would depend on several aquifer parameters including recharge rate, regional
aquifer discharge rate, hydraulic properties, and sea level. Recent climate change studies have
shown that the global sea-level, on average, is expected to rise between 18 and 59 cm this
century [IPCC, 2007]. Worst-case projections show it could be as high as 180 cm [Vermeer and
12
Rahmstorf, 2009]. Therefore, environmental planners worldwide are concerned about the
impacts of sea-level rise on saltwater intrusion processes, especially in over-utilized, urbanized
coastal aquifers that already have low groundwater levels. Currently, it is expected that the
rising sea level will enhance saltwater intrusion and potentially contaminate many freshwater
reserves. Figure 2.1a depicts a commonly assumed conceptual model [Cheng et al., 2000;
Custodio and Bruggeman, 1987; Fetter, 2001; Mantoglou, 2003; Strack, 1976; Werner and
Simmons, 2009] that illustrates how the rising sea level would impact the groundwater quality by
forcing the wedge to migrate inland. This conceptual model, however, ignores the fact that when
the seawater rises at the sea-side boundary the system would pressurize and the water-table level
might be ?lifted? throughout the aquifer. Figure 2.1b illustrates a revised conceptual model.
This revised model accounts for the lift in the groundwater level over the entire system due to the
changes in the sea-side boundary condition. It is expected that after a long period (i.e., at or near
steady state) this lifting effect would approximately raise the entire fresh water body (measured
from the bottom of the aquifer) by an extent similar to the sea-level rise (i.e., a similar order of
magnitude). One could intuitively expect this lifting mechanism to counteract and reduce the
impacts due to the sea-level rise. However, it is unclear to what extent this lifting process could
reduce the overall impacts due to sea-level rise. In this study, we employ this revised conceptual
model to hypothesize that changes in the sea-level might have little or no impact on saltwater
intrusion when the net flux through the system is unchanged. The overall goal of this conceptual
study is to test the validity of this hypothesis, and understand its ramifications under transient
conditions in idealistic confined and unconfined systems.
13
Figure 2.1. Comparison of conceptual models used for visualizing the impacts of sea-level rise
on a saltwater wedge: (a) salt wedge profile after sea-level rise based on a traditional conceptual
model that ignores the lifting effect, and (b) a new conceptual model that includes the lifting
effect.
Regional flow, Q
Expected intrusion level
Recharge, W
( a )
New Sea Level
Old Sea Level
Wedge
before rise
Wedge
after rise
Regional flow, Q
New Groundwater Level
Old Groundwater LevelNew Sea Level
Old Sea Level
Expected intrusion level
after accounting for groundwater table rise
Wedge
before rise
( b )
Recharge, W
Wedge
after rise
14
2.3 Details of Numerical Experiment
The base-case problem considered in this study was adapted from Werner and Simmons?
conceptual model [Werner and Simmons, 2009] of a seawater intrusion field study completed in
Pioneer Valley, Australia. The original problem only considered unconfined flow conditions; in
this work, the problem was modified to investigate both confined and unconfined conditions.
Our numerical study considered a two-dimensional aquifer system which is 1000 m long and 30
m thick. A rectangular numerical grid, with ?x = 4 m, ?y = 1 m and ?z = 0.4 m, was used. The
initial sea level (prior to the rise) was assumed to be at 30 m and the sea level was then allowed
to rise instantaneously to 34 m. The net sea-level rise assumed was 4 m, a theoretical worst-case
scenario which is approximately double the extreme value predicted by Vermeer and Rahmstorf
[2009]. The assumed sea-level rise was chosen to approximately mimic the last interglacial
period?s rapid sea-level rise that reached up to 4 to 6 m due to rapid loss of ice-sheet [Blanchon
et al., 2009]. It is important to note that this is not a ?true? scenario simulation exercise; our
objective is not to forecast the future location of saltwater wedge for a specific system, rather it is
to develop a generic conceptual understanding to assess the potential impacts of sea-level rise on
saltwater intrusion processes. Therefore, some extreme sea-level rise scenarios were initially
simulated (as the base-case) to better illustrate certain subtle transport mechanisms. Later, as a
part of the sensitivity analysis, we explore some realistic sea-level rise scenarios that range from
0.2 m (minimum rise predicted by IPCC) to 4 m (double the maximum rise predicted by
Vermeer and Rahmstorf [2009] ). Also, in all base-case simulations we assumed instantaneous
rise (a worst-case scenario) and later in the sensitivity section we have presented the results for
other finite rates of sea-level rise.
15
The seawater density was assumed to be 1,025 kg/m3 with salt concentration of 35 kg/m3.
The regional groundwater flux supplied through the inland boundary using a set of artificial
injection well nodes. The flux (q, flow per unit depth) supplied from the right boundary was
0.005 m2/day (or a total flow rate (Q) of 0.15 m3/day over the 30 m thick aquifer). Recharge (W)
flow was delivered from the top boundary at the rate of 5?10-5 m/day. These regional and
recharge flows were simply selected to locate the initial location of the saltwater toe in between
300 to 500 m. The hydraulic conductivity was set to 10 m/day, specific storage was set to 0.008
m-1 in each confined layer, and porosity to 0.35. Longitudinal and transversal dispersivity values
were assigned to be 1 m and 0.1 m, respectively. The initial sea level was set at 30 m to simulate
a base-case steady-state wedge. This steady-state wedge was later used as the initial condition in
all subsequent sea-level rise simulations. Impacts of density-weighting schemes (upstream and
central) and various advection schemes (TVD and upstream finite-difference) were investigated
in using tests simulations (results summarized in Appendix 2 and 3). In the final simulations we
used upstream density weighting and the finite difference scheme advection package. Other
model parameters used are summarized in Table 2.1.
The MODFLOW-family variable density flow code SEAWAT [Langevin C et al., 2003]
was used in this study with the central difference weighting option. The modeling approach was
validated by solving Henry?s steady-state solution [Henry, 1964; Simpson and Clement, 2004]
and also laboratory data provided by Goswami and Clement [2007], and Abarca and Clement
[Abarca and Clement, 2009]. Several sets of numerical experiments were completed to explore
various types of flow conditions and parameter values. The first set of experiments only
considered confined-flow conditions and the results were used to explore certain novel salt-
16
wedge reversal mechanisms that have not been reported in the published literature. Later
simulations consider both confined and unconfined flow conditions.
Table 2.1. Aquifer properties and hydrological properties used for the base-case problem.
Property Symbol Value
Horizontal aquifer length L 1000 m
Vertical aquifer thickness B 30 m
Inland groundwater flow rate Q 0.15 m3/day
Uniform recharge rate from top W 5?10-5 m/day
Hydraulic conductivity K 10 m/day
Specific storage, Ss 0.008 m-1
Longitudinal dispersivity ?L 1 m
Transverse dispersivity ?T 0.1 m
Saltwater concentration Cs 35 kg/m3
Saltwater density ?s 1,025 kg/m3
Freshwater density ?f 1,000 kg/m3
Porosity n 0.35
17
2.4 Results
2.4.1. Impacts of Sea-level Rise on Confined Flow Conditions
The goal of the first steady-state simulation is to estimate the initial saltwater wedge
profile that would exist in the system prior to sea-level rise. Figure 2.2 (continuous line) shows
the 50% isochlor of the initial steady-state concentration profile for the base-case, confined-flow
system. This 50% concentration profile is defined as the base saltwater wedge. The figure
shows that under the assumed groundwater flow conditions, the toe of the salt wedge advanced
to 382 m into the aquifer before sea-level rise. This steady-state condition was used as the initial
condition in all subsequent simulations.
The second steady-state simulation aimed to predict the long-term salt-wedge profile
after the sea-level rise. The water level at the seaside boundary was abruptly increased from 30
to 34 m to simulate an instantaneous sea-level rise. The system was allowed to evolve for
80,000 days to reach steady-state conditions. The new steady-state solution for the salt wedge
simulated by the model after raising the sea level is also shown in Figure 2.2 (using circle data
points). Interestingly, and rather surprisingly, the model-predicted saltwater profiles for both
conditions (pre and post sea-level rise) were identical, indicating that the sea-level rise will have
no impact on the location of the steady-state wedge when the freshwater flux transmitted through
the system remained constant (i.e., if the rainfall/recharge pattern was not changed). This non-
intuitive result has several important practical implications. The result indicates that the lifting
effect postulated in the conceptual model described in Figure 1c will fully offset the negative
impacts of sea-level rise in constant flux confined flow systems.
18
To further explore this result, we present the transient evolution of the location of the
saltwater wedge in Figure 2.3. These profiles show that when the sea-level was instantaneously
raised, the wedge initially started to move inward; however, after about 8,000 days, the direction
of wedge movement reversed and the wedge started to move backward until it reached the initial
steady-state profile. As far as we are aware, no one has predicted or postulated this self-reversal
mechanism. Understanding this self-reversal process has an enormous implication on how water
resources managers would perceive and manage the impacts of sea-level rise on saltwater
intrusion. 10% isochlor and 90% isochlor shows the same transient movement (See appendix 1)
Figure 2.4 shows the transient variations in the toe position of the saltwater wedge (XT).
The data shows that, under assumed flow conditions, the salt wedge first advanced into the
aquifer for about 8,000 days and reached a maximum distance of 432 m (from the sea boundary)
and then started to recede. The figure also shows the simulated transient freshwater head levels
at the right side boundary. This freshwater head data is an excellent surrogate for quantifying the
progression of the aquifer ?lifting? process, postulated in Figure 2.1b. This dataset shows that
the initial head at the right boundary, as predicted by the model prior to the sea-level rise, was 31
m. When the sea-level was raised instantaneously from 30 to 34 m, it took about 2,000 days for
the right boundary to fully respond to this change. After about 2,000 days, the fresh-water level
in the right boundary reached a constant value of about 35 m. The net change in the freshwater
level was about 4 m, almost identical to the sea-level rise forced at the left boundary. This
implies that the raising seaward-boundary head lifted the entire fresh groundwater system by a
similar order of magnitude. This lifting process was able to fully reverse the salt wedge location.
The rate of the reversal process would depend on the rate at which the groundwater system was
lifted (or how quickly the groundwater heads would respond), which would in turn depend on the
19
storage properties of the aquifer and the rate of sea-level rise. In the following section, we
provide a detailed analysis that quantifies the sensitivity of the self-reversal process to the value
of storage coefficient, the rate of sea-level rise and the magnitude of rise. All these simulations
used the base-case steady-state solution as the initial condition.
20
Figure 2.2. Comparison of steady-state salt wedges predicted before and after the sea-level rise in
the confined system.
0
10
20
30
0 200 400 600 800 1000
z (
m)
x (m)
Initial condition
80,000 days, steady state
21
Figure 2.3. Transient salt wedge profiles for the confined system.
0
10
20
30
0 200 400 600 800 1000
z (
m)
x (m)
Initial condition
3,000 days
8,000 days
15,000 days
80,000 days, steady state
22
Figure 2.4 Transient variations in the toe position (XT) and freshwater head (hf) at the inland
boundary.
30
32
34
36
380
400
420
440
460
0 10,000 20,000 30,000 40,000
h f
(m
)
X T
(m
)
Time (day)
Toe position
Freshwater head
23
2.4.2. Sensitivity Analysis of the Self-reversal Process in Confined Systems
2. 4.2.1. Sensitivity to Specific Storage
We first explored the sensitivity of the intrusion mechanism to various values of specific
storage (Ss). In this analysis, the value of Ss was varied by an order of magnitude with the range
0.0008 to 0.008 m-1. In all sensitivity simulations, the other parameters were fixed at base-case
levels shown in Table 2.1. Figure 2.5a shows the temporal variations in the toe position (XT) for
different values of Ss. The data shows that when Ss was small the system responded rapidly and
the intrusion effect was reversed quickly. Also, the maximum intrusion length was small for
smaller values of Ss. The maximum values of the predicted saltwater toe position (XT) were 432,
412, 394, and 386 m for the Ss values 0.008, 0.005, 0.002 and 0.0008 m-1, respectively. The
figure shows that when Ss is small (at 0.0008 m-1) the change in the salt wedge location was
relatively small even for an extreme sea-level rise of 4 m. The total time required to reach the
maximum intrusion level (defined as the duration of intrusion) was 6,000 to 8,000 days. These
values appear to be relatively insensitive to Ss. Overall, the extent of saltwater intrusion as
indicated by the maximum value of XT was more sensitive to Ss than the duration of intrusion.
2.4.2.2. Sensitivity to the Magnitude of Sea-level Rise
The sensitivity of the intrusion length to the magnitude of the sea-level rise was explored
by varying the sea-level rise within the range of 0.2 to 4 m. Figure 2.5b shows the transient
variations in XT values for various values of sea-level rise. The maximum values of XT were
approximately proportional to the magnitude of sea-level rise. Also, similar to the previous
24
sensitivity experiment, the time taken to reach the maximum values was about 6,000 to 8,000
days, and this time was relatively insensitive to the magnitude of the sea-level rise.
2.4.2.3. Sensitivity to the Rate of Sea-level Rise
Simulations were completed to test the sensitivity of the intrusion process to changes in
the rate of sea-level rise. Literature data indicated that the current rate of sea-level rise observed
worldwide has ranged from 2 to 4 mm/year [McCarthy, 2009]. Climate change model
projections, however, show that the global rate could at least double by the end of this century
[Anderson et al., 2010]. A total rise of 4 m was simulated using six different rate scenarios:
instantaneous, 1 mm/day for 4,000 days, 0.1 mm/day for 40,000 days, 0.05 mm/day for 80,000
days, and infinite rise at a rate of 0.04 mm/day with Ss = 0.008 and Ss = 0.0008. Figure 2.5c
shows the temporal variations in the simulated toe position (XT) for all five scenarios. The data
demonstrates that the rate of the self-reversal process would depend on the sea-level rise rate.
When the rate of rise was low the reversal cycle had a longer duration. The maximum value of
the intrusion length, XT, decreased with decrease in the rate of sea-level intrusion. The time
taken to reach the maximum level of intrusion was 7,500, 10,000, 40,000 and 80,000 days for
instantaneous, 1 mm/day, 0.1 mm/day, 0.05 mm/day sea-level rise rates, respectively. When sea-
level was allowed to rise infinitely at a fixed rate, the results reached an irreversible, quasi
steady-state level. It is important to note this irreversible intrusion level was primarily an artifact
due the relatively large storage value (of 0.008 m-1) assumed as the base-case parameter. When
we reduced the specific storage value by an order of magnitude (to 0.0008 m-1, a more typical
value for confined flow) the head information propagated quickly and the lifting effect became
continuously active. Therefore, the system with a higher Ss value experienced immediate
25
reversal and experienced very little intrusion under the continuous-rise scenario. Overall,
changes in the rate of sea-level rise influenced the maximum level of intrusion as well as the
time required to reach the maximum level. The time required to reach the maximum would, to a
large extent, depended on how long the sea rise occurred.
2.4.2.4. Sensitivity to Dispersivity Values
We explored the sensitivity of the intrusion mechanism to dispersivity coefficients by
varying the value within the range of 0.5 m to 2 m for longitudinal dispersivity, ?L, and 0.05 m to
0.2 m for transverse dispersivity, ?T. As Abarca et al [2007b] pointed out, changing the values of
dispersivity would impact the value of XT (initial XT decreased when the dispersivity coefficients
were increased). Therefore, we defined a parameter ?net change in salt wedge location,? which
was computed as: maximum value of XT - initial position of XT, to quantify the changes. Figure
2.5d shows the temporal variations in the toe position (XT) for different values of dispersivity.
The net change in the salt wedge location was in between 44 and 53 m. The time taken to reach
the maximum intrusion level and the peak of XT values were relatively insensitive to changes in
the value of dispersivity coefficient.
26
Figure 2.5. Sensitivity of the toe position (XT) to (a) specific storage, (b) the magnitude of sea-
level rise, (c) the rate of sea-level rise velocity and (d) the dispersivity coefficients.
380
400
420
440
0 10,000 20,000 30,000 40,000
X T
(m
)
Time (day)
Ss = 0.008
Ss = 0.005
Ss = 0.002
Ss = 0.0008
(a)
380
400
420
440
460
0 20,000 40,000 60,000 80,000 100,000
X T
(m
)
Time (day)
instantaneous rise
1 mm/day for 4,000 days
0.1 mm/day for 40,000 days
0.05 mm/day for 80,000 days
0.04 mm/day, infinite rising (Ss = 0.008)
0.04 mm/day, infinite rising (Ss = 0.0008)
(c)
380
400
420
440
0 10,000 20,000 30,000 40,000
X T
(m
)
Time (day)
4 m sea level rise
2 m sea level rise
1 m sea level rise
0.2 m sea level rise
(b)
300
360
420
480
540
0 10,000 20,000 30,000 40,000
X T
(m
)
Time (day)
half dispersivity
alphaL = 1 m
alphaL = 2 m
?L= 0.5 m, ?T= 0.05 m
?L= 1.0 m, ?T= 0.1 m
?L= 2.0 m, ?T= 0.2 m
(d)
27
2.4.2.5. Sensitivity to Other Hydrological Parameters
It is important to note that the transient reversal patterns would depend on the values of
hydraulic conductivity, recharge rate and ambient groundwater flow. Sensitivity to variations in
all these model parameters was also explored in this study. Simpler analytical solutions can be
used to intuitively infer the sensitivity to individual variations in K, W and Q values under steady
conditions. Werner and Simmons [Werner and Simmons, 2009] followed this approach and
completed a detailed sensitivity assessment for a steady-state unconfined problem. In this study,
we present the results of a selected number of sensitivity tests completed for a transient confined
aquifer system using SEAWAT. Figure 2.6a shows the transient variations in XT values for
different hydraulic conductivity values; K values used are: 5 m/day, 10 m/day and 15 m/day. It
should be noted that the simulated profiles have different initial XT since individually changing
any one of these parameters (K, W or Q) would alter the steady-state solution of the base
problem. The results show that the peak value of XT (relative to the initial XT value) would
increase, and the responding time would decrease (as expected) when the K value was decreased.
Figure 3.6b shows the transient variations in XT values for various values of Q and W (0.67?base
values: Q = 0.1 m3/day and W = 3.3?10-5 m/day; base values: Q = 0.15 m3/day and W = 5?10-5
m/day; doubled the base values: Q = 0.3 m3/day and W = 1?10-4 m/day). The data show that
when the amount of freshwater flow was increased the peak value of XT increased and system
also had a shorter responding time. It should be noted that if the transport parameters were
scaled consistently (for example, increase K and reduce the Q by a similar factor) then there will
be very little variation in the peak value of XT.
28
Figure 2.6. Sensitivity of the toe position (XT) to: (a) hydraulic conductivity and (b) freshwater
flux.
200
400
600
800
0 10,000 20,000 30,000 40,000
X T
(m
)
Time (day)
K = 5 m/day
K = 10 m/day
K = 15 m/day
(a)
200
400
600
800
0 10,000 20,000 30,000 40,000
X T
(m
)
Time (day)
base freshwater flux
2 base flux
0.67 base flux
(b)
?
?
29
2.4.3. Impacts of Sea-level Rise on Unconfined Aquifers
To examine the impacts of sea-level rise on unconfined flow conditions, we completed
numerical simulations of an unconfined aquifer with dimensions identical to those used in the
confined simulation, except the top layer was modified to simulate unconfined flow. The length
of the unconfined aquifer was 1,000 m and the total thickness was 35 m. A numerical grid with
?x = 4 m, ?y = 1 m and ?z = 0.4 m (75 confined layers), and a top unconfined layer of ?z = 5 m
was used. In all unconfined flow simulations, the value of average hydraulic conductivity was set
to 10 m/day, total regional freshwater flow (Q) from the right boundary was set to 0.15 m3/day,
and the areal recharge flux (W) was set to 5?10-5 m/day. Initially, the regional flow rate was
assigned to 0.0019 m3/day for all confined layers with 0.0048 m3/day for top unconfined layer
(note the top unconfined layer was 1 m as opposed confined nodes which are 0.4 m). After the
sea-level rise, the flow rate was again redistributed over the unconfined zone of 4 m. In order to
redistribute the same flow of 0.15 m3/day over an increased saturated aquifer depth, all the
confined nodes were assigned a flow of 0.0017 m3/day and the top 4-m unconfined layer was
assigned 0.021 m3/day. In order to compare unconfined flow simulation results against confined
flow results, an identical instantaneous sea-level rise (rise from 30 to 34 m) scenario was
assumed. The specific storage, Ss, was set to 0.008 m-1 for all confined layers, and specific yield,
Sy, was set to 0.1 for the top unconfined layer.
We first completed a base-case, steady-state simulation for the unconfined system to generate
the initial conditions that existed prior to the sea-level rise. Figure 2.7 compares the steady-state
salt wedge predicted for the unconfined flow system with the wedge predicted for a similar
confined flow system (data from Figure 2.2). The figure shows the both wedges are almost
30
identical, indicating that both confined and unconfined systems would behave in a similar
manner at steady-state conditions. The similarity between the two steady-state solutions can also
be explained mathematically, as shown in Appendix 5.
To understand the impacts of sea-level rise on the unconfined aquifer, we instantaneously
raised the sea-level by 4 m and let the system reach steady state. Figure 2.8 compares the initial
and final steady-state saltwater wedge profiles in the unconfined system (dotted and continuous
lines). These profiles show that, unlike the confined system (compare with Figure 3 results), the
saltwater intrusion process is not reversible for the unconfined system. In unconfined aquifers,
the initial location of the toe XTi and the final location of the toe position XTf are distinctly
different. This is because, under unconfined conditions, the sea-level rise increases in the
saturated thickness (or transmissivity) of the aquifer. This increased transmissivity allows the
wedge to penetrate further into the system, resulting in a new steady-state condition. Comparing
the initial and final salt-wedge profiles indicates the initial salt wedge was approximately raised
by 4 m, similar to the level of sea-level rise. The groundwater level also rose over the entire
aquifer by about 4 m (as illustrated in Figure 2.1b), all the way to the inland boundary (the
predicted groundwater raise at the inland boundary was similar to the data shown in Figure 2.4).
These data indicate that the aquifer lifting process is active in the unconfined system. The lifting
process, however, was not able to fully reverse the wedge location since the system evolved from
a steady-state solution for the 30-m thick aquifer (with initial toe, XTi, at about 382 m) to a new
steady-state solution for the 34-m thick aquifer (with final toe, XTf, at about 510 m). The data
points (marked with circles) shown in the figure are steady-state wedge data predicted for a 34-m
thick ?equivalent confined aquifer.? As expected (see Appendix A for an analytical analysis),
31
the final unconfined steady-state solution matched with an equivalent confined aquifer solution
with an appropriate value of aquifer thickness.
Figure 2.9 compares the results of transient changes in toe length predicted for the
following three systems: Case-1) standard base-case confined system of 30 m thickness (same as
the data shown in Figure 2.2); Case-2) unconfined flow system with base-case parameters; and
Case-3) unconfined system with very high Ss value (0.02 m-1). It is important to note the Ss value
used for Case-3 is unrealistically high and this conceptual simulation was completed to
demonstrate the existence of certain subtle self-reversal mechanisms. Several observations can
be made from the results presented in Figure 2.9. All three solutions started at the initial toe
location XTi = 382 m and for about 1,000 days the unconfined solutions were approximately
equal to the confined flow solution. After this time, the unconfined model solutions started to
diverge. The Case-2 unconfined flow simulation reached the final steady-state toe position, XTf,
after about 50,000 days, and the system did not show any reversal effect. However, when we
repeated the unconfined simulation using high Ss values (Case-3) we could clearly observe the
reversal effect. Also, as expected, Case-3 required more time (about 100,000 days) to reach
steady-state final toe position XTf. In summary, the aquifer lifting effect is active in unconfined
systems but it is difficult to observe any reversal effects in transient unconfined systems due to
the changes in the aquifer transmissivity. When an appropriate value of aquifer thickness was
used, the steady-solution for an unconfined flow system would be almost identical to an
?equivalent confined flow system.? The storage available within fully-saturated layers (which
act as confined layers) are typically very low and this would allow rapid propagation of head
perturbations; therefore, it is difficult to observe reversal effects in any realistic unconfined
aquifers.
32
Figure 2.7. Comparison of steady-state salt wedges predicted before the sea-level rise (sea level
at 30m) in the unconfined- and confined-flow systems.
0
10
20
30
0 200 400 600 800 1000
z (
m)
x (m)
Confined flow system
Unconfined flow system
33
Figure 2.8. Comparison of steady-state salt wedges predicted before the sea-level rise in the
unconfined system, after the sea-level rise in the unconfined system, and in an equivalent 34-m
thick confined-flow system.
0
10
20
30
0.00 200.00 400.00 600.00 800.00 1000.00
z (
m)
x (m)
Unconfined flow system before sea-level rise
Unconfined flow system after sea-level rise
34 m - equivalent confined flow system
?4 m
?4 m
XTi XTf
34
Figure 2.9. Comparison of transient variations in toe positions (XT) predicted for the confined
and unconfined flow systems.
380
420
460
500
540
580
0 20,000 40,000 60,000
X T
(m
)
Time (day)
Case 1 : 30 m confined flow ( base case )
Case 2 : Unconfined flow system (Ss = 0.008)
Case 3 : Unconfined flow system (Ss = 0.02)
XTi
XTf
35
2.5 Conclusions
Detailed numerical experiments were completed using the MODFLOW-family computer
code SEAWAT to study the transient effects of sea-level rise on saltwater wedge in confined and
unconfined aquifers with vertical sea-land interface. The simulation results show that if the
ambient recharge remains constant, the sea-level rise will have no impact on the steady-state salt
wedge in confined aquifers. The transient confined-flow simulations help identify an interesting
self-reversal mechanism where the wedge, which initially intrudes into the formation due to sea-
level rise, would be naturally driven back to the original position. However, in unconfined-flow
systems this self-reversal mechanism would have lesser effect due to changes in the value of the
effective transmissivity (or average aquifer thickness). The increased trasmissivity resulted in a
redistribution of freshwater flux at the inland boundary condition. Both confined and unconfined
simulation experiments show that rising seas would lift the entire aquifer and this lifting process
would help alleviate the overall long-term impacts of saltwater intrusion.
The sensitivity simulations show that the rate and the extent of the self-reversal process
would depend on the value of specific storage (Ss) and the rate of sea-level rise. When the rate
of rise was low the reversal cycle had a larger duration. Overall, the changes in the rate of sea
rise had relatively less influence on the maximum level of intrusion and more influence on the
time taken to reach the maximum level. On the other hand, variation in Ss values was more
sensitive to the maximum level of intrusion than the duration of intrusion (the time taken to reach
the maximum value).
It is important to note that the results presented in the conceptual modeling study are
based on simulations completed for idealized rectangular aquifers with homogeneous aquifer
36
properties. While the results are useful for developing a large-scale conceptual understanding of
the impacts of sea-level rise, evaluation of true impacts would require detailed site-specific
modeling efforts [Abarca et al., 2007a; Lo?iciga et al., 2011]. A better understanding the self-
reversal mechanism (both its spatial and time scales) identified in this study would have an
enormous implication on managing the impacts of sea-level rise in coastal groundwater aquifers.
The results, however, do not imply that one could simply ignore climate change effects on
saltwater intrusion process. Rather it implies that we can minimize its risks based on a sound
scientific understanding of the transport processes and by developing pro-active management
strategies that are appropriate to unconfined aquifers and confined aquifers. It is important to
note that this study assumes the fluxes in the system would remain constant. However, site-
specific climate change effects could greatly alter the recharge and regional fluxes (these natural
hydrological fluxes can change due to variations in rainfall patterns), therefore the overall
problem should be managed in the context of large-scale variations in hydrological fluxes
expected to be induced by climate change effects. In addition, Lo?iciga et al. [Lo?iciga et al.,
2011] pointed out that variations in groundwater extraction (anthropogenic fluxes) was the
predominant driver of sea water intrusion in a model that simulated sea-level rise scenarios for
the City of Monterrey, California. Finally, it is very likely that rising heads might increase
evapo-transpiration fluxes. This could impact the overall hydrological budget resulting in less
recharge reaching the coastal area and this will hinder the self-reversal process. Therefore, site-
specific management models for coastal areas should carefully integrate changes in both natural
and anthropogenic fluxes with various sea-level rise scenarios.
37
Chapter 3
INVESTIGATION OF SALTWATER INTRUSION PROCESSES USING LABORATORY
EXPERIMENTS INVOLVING FLUX-TYPE BOUNDARY CONDITIONS
3.1.Introduction
Saltwater intrusion is a natural process where sea water would advance into coastal
groundwater aquifers due to the density difference between saline and fresh waters creating a
wedge that evolves landward. The horizontal extent of saltwater intrusion could range from a
few meters to many kilometers [Barlow and Reichard, 2010]. Several natural and anthropogenic
processes can enhance the effects of saltwater intrusion. For example, the water supply
operations in coastal regions may depend on pumping freshwater from local aquifers, and these
pumping activities can exacerbate both vertical and horizontal salt intrusion processes. Droughts
induced by climate change effects could result in the loss of freshwater resources and enhance
saltwater intrusion. Several modeling and field studies have provided evidence that climate
change could decrease the net freshwater input to groundwater resources [Feseker, 2007;
Lo?iciga et al., 2011; Masterson and Garabedian, 2007; Oude Essink, 2001; Oude Essink et al.,
2010; Ranjan et al., 2006; Rozell and Wong, 2010; Yu et al., 2010], and this could aggravate the
impacts of saltwater intrusion. Climate change studies estimate that the global sea level could
increase between 18 cm and 59 cm this century [IPCC, 2007]; other worst-case scenario
predictions have forecast higher sea-level rises up to 180 cm [Vermeer and Rahmstorf, 2009],
which could result in more severe saltwater intrusion.
The number of scientific investigations related to climate change effects have increased
rapidly in the last decade [Green et al., 2011]. Most of these investigations are based on the
38
premise that climate change will have severe adverse impacts on almost all natural systems
including groundwater systems. As pointed out by Kundzewics and D?lle [2009], flat areas such
as deltas and coral islands are expected to be inundated by rising seas leading to reduced
freshwater availability. According to an extreme case reported by Yu et al. [2010], about a meter
sea-level rise would inundate and contaminate over 18% of the total land in Bangladesh,
threatening to physically displace about 11% of their population.
Chang et al. [2011] recently presented some interesting counter-intuitive results
indicating that sea-level rise would have no impact on saltwater intrusion in certain idealized
confined systems that ignore marine transgression effects, such as a tidal effect and seasonal
variation. Their results show that rising seas would lift the entire water table, and this lifting
process has the potential to alleviate the overall long-term impacts of saltwater intrusion via a
natural self-reversal process. White and Falkland [2009] reviewed previously published field
datasets and concluded that the freshwater lenses in several atoll islands might not be severely
affected by sea-level rise up to 1 m, as long as the land mass is not lost due to inundation of the
edges of these islands. On the other hand, potential variations in the recharge patterns, resulting
from changes in the rainfall levels, are expected to have much larger impact on the overall water
budget of these islands. Rozell and Wong (2010) used climate change scenarios available in the
Intergovernmental Panel on Climate Change 2007 (IPCC) report and found that only a small
portion of Shelter Island, located New York, USA, would experience inundation and salt-water
intrusion even when the sea level was assumed to increase more than three times higher than
those values reported in the IPCC report. Some of these conflicting results imply that the
impacts of climate change on saltwater intrusion should be carefully predicted and managed at a
39
local scale, on a case-by-case basis, since the overall impacts would depend on various
hydrological factors including spatial and temporal variations in rainfall/recharge patterns.
Werner and Simmons [2009] categorized saltwater intrusion problems into two types
depending on how the freshwater would enter the aquifer. The first type of system is known as a
head-controlled system, where the groundwater level fixed at the inland boundary would
transmit various rates of freshwater flow depending on aquifer properties. The second type is a
flux-controlled system where the net freshwater flow entering at the aquifer inland boundary is
set to a constant flow rate. From a fundamental point of view, even in the head-controlled
system the changes in saltwater intrusion patterns would be essentially due to the net decrease or
increase in freshwater flow transmitted through the system. Therefore, the amount of freshwater
flux moving through the system is an important fundamental parameter that would control
saltwater intrusion effects. Chang et al. [2011] demonstrated how the impacts of sea-level rise
on both unconfined and confined systems could be strongly influenced by the changes in fresh-
water fluxes. They also found that when the fresh water flux is held constant the sea-level rise
would have no impact on the steady-state saltwater wedge profiles of confined aquifers. Kuan
[2012] reported that the saltwater wedge and the resulting upper saline plume in unconfined
coastal systems will be affected by the changes in the regional freshwater influx induced by
various tidal conditions. These studies have demonstrated the importance of understanding the
effects of changes in different types of groundwater fluxes on saltwater intrusion mechanisms.
The objective of this research is to conduct experimental studies to investigate the
impacts of variations in regional and areal-recharge fluxes on saltwater intrusion in unconfined
aquifers. In this effort, we have explored two types of flux-controlled systems: one driven by
regional flow delivered from the inland freshwater boundary and another driven by a fixed
40
amount of recharge flow distributed from the top. In the following sections, we refer to these
two systems as regional-flux system and areal-recharge flux system, respectively. The results
from the experiments were used to develop a better understanding of the influence of different
types of flux-boundary conditions on the saltwater intrusion process, especially under transient
flow conditions.
3.2 Methods
3.2.1 Experimental Approach
Figure 3.1 shows the experimental methods employed in this study. The laboratory setup
shown in Figure 3.1a was used to conduct regional-flux aquifer (RFA) experiments, and Figure
3.1b was used to conduct areal-recharge flux aquifer (AFA) experiments. All the experiments
were completed in a sand tank with dimensions: 50.5 cm ? 28 cm ? 2.2 cm. We used a relatively
narrow tank to simulate a two-dimensional system that represents cross-sectional flow in an
unconfined aquifer. As shown in the figure, the sand tank has three distinct chambers: an inlet
chamber, porous media chamber, and an outlet chamber. In the regional flux experiment, the
right inlet chamber was used to deliver freshwater flow to the system. The sea level variations
were simulated by adjusting the saltwater head at the left outlet chamber. The central porous
media chamber was filled with silica beads of diameter 1.1 mm. The beads were packed under
saturated conditions to prevent the trapping of air bubbles inside the tank. The average hydraulic
conductivity value of the system was estimated to be 1.46 cm/s. Other experimental methods
used here were similar to those used in our previous experiments reported in Goswami and
Clement [2007] and Abarca and Clement [2009]. It is important to note that our previous
experimental studies and those published by others [Karasaki et al., 2006] have used constant-
41
head boundary conditions; here we have employed two different types of flux boundary
conditions. To simulate the areal recharge flux boundary condition, we used a multi-head
peristaltic pump that delivered flow at various injection points. The flow was divided equally
and was distributed from four outlet tubes. The tubes were placed on the top of a porous media
chamber (well above the water table) at distances of 10 cm, 20 cm, 30 cm, and 40 cm away from
the inlet chamber. To simulate the constant regional flux condition, we used the same multi-
head pump that injected a constant amount of flow via a single outlet which was placed into the
right chamber. All the saltwater intrusion experiments were recorded using a Panasonic video
camera (HDC-HS250) and the photographs were later post-processed using methods described in
Goswami and Clement [2007]. A commercial food dye was added at a ratio of 50 ml per 20 L of
saltwater to help differentiate the saltwater from ambient fresh water. The measured density of
the colored saltwater was 1.027 g/cm3.
42
(c)
Figure 3.1. Schematic diagram of the (a) regional-flux and (b) areal-recharge flux experiments
and (c) photograph of laboratory setup.
Freshwater PumpSaltwater Pump
50.5 cm
Hs
= 2
1.0
cm
(a)
Freshwater PumpSaltwater Pump
50.5 cm
H s
= 2
1.0
cm
(b)
43
3.2.2 Numerical Modeling Approach
The MODFLOW-family variable density flow code SEAWAT [Guo and Langevin, 2002]
was used to model the experimental data. The fully implicit option and the total variation
diminishing (TVD) package were used in all our simulations. The code has been widely tested
by solving various benchmark problems including Simpson and Clement [2003] and Simpson
and Clement [2004]. In our model, the origin of an x-z coordinate system was set at the lower
left-hand corner of the plane. A uniform finite-differences grid of size 0.5 cm (total of 102 cells
in the x-direction) was used to discretize the horizontal direction. Forty-two uniform confined
layers of width 0.5 cm were used to discretize the vertical direction; in addition, a single
unconfined layer of width 1 cm was used to represent the top region. The depth of the top layer
was selected such that the layer would remain unconfined (i.e., it did not fully drain or overfill)
in all our simulations (note that the freshwater head in the tank varied from 21.6 to 22.0 cm,
measured from the bottom, in our transient experiments and the saltwater head was set to 21.0
cm). The lower boundary nodes were set to no-flow conditions. The longitudinal (?L) and
transversal (?T) dispersivity coefficients were assumed to be 0.05 cm and 0.005 cm, respectively;
these are the values previously used by Abarca and Clement [2009]. The value of specific
storage (Ss) was set to 10-6 cm-1 for the confined layers and the value of specific yield (Sy) was
set to 0.1 for the top unconfined layer. The measured values of freshwater density and the
saltwater density were 1.000 and 1.027 g/cm3, respectively. Using the standard density-
concentration slope factor of 0.714, the value of salt concentration was estimated to be 0.378
g/cm3. The concentration value was also verified using data for the amount of salt used to
prepare the saltwater solution. The entire flow domain was assumed to be initially filled with
44
saltwater, and a constant rate of freshwater flow was injected at the boundary to simulate the
desired initial steady-state condition.
Table 3.1. Summary of numerical model parameters
Property Symbol values
Horizontal length L 50.5 cm
Saltwater elevation hs 21 cm
Hydraulic conductivity K 1.46 cm/s
Specific yield, Sy 0.1
Specific storage, Ss 0.00001cm-1
Longitudinal dispersivity ?L 0.05 cm
Transverse dispersivity ?T 0.005 cm
Saltwater concentration Cs 0.0378 g/cm3
Saltwater density ?s 1.027 g/cm3
Freshwater density ?f 1.000 g/cm3
Table 3.2. Summary measured and modeled flows
Experiment
Experimentally
measured flows in the
sand tank (2.2 cm wide)
Numerical flows for unit width model
RFA-test
(SS1, SS3) 0.832 cm
3/s 0.378 cm3/s (0.0086 cm3/s for 42 confined cells,
0.017 cm3/s for one unconfined cell)
RFA-test (SS2) 0.463 cm3/s 0.211 cm
3/s (0.0048 cm3/s for 42 confined cells,
0.009 cm3/s for one unconfined cell)
AFA-test
(SS1, SS3) 1.111 cm
3/s 0.505 cm3/s (recharge rate is 0.0100 cm/s)
AFA-test
(SS2) 0.755 cm
3/s 0.343 cm3/s (recharge rate is 0.0068 cm/s)
45
3.3 Results
3.3.1 . Regional-Flux Experiments
We first completed a regional flux experiment to study the transient changes in saltwater
intrusion patterns when the flux was varied at the right boundary. The saltwater level at the left
boundary was fixed at 21 cm. Figure 3.2 shows the digital data collected for both advancing-
wedge and receding-wedge experiments. The time levels required for conducting these
experiments were estimated a priori from preliminary numerical simulations. To simulate the
initial condition, a constant rate of freshwater (0.832 cm3/s, see Table 4.2) was injected into the
right chamber and the system was allowed to reach a steady-state pattern (designated as SS1 in
Figure 3.2). To start the transient intruding conditions, the flow rate was instantaneously
reduced from 0.832 to 0.463 cm3/s (about 56% of the initial flow). The transient data were
recorded for 90 min, after which the system reached the second steady-state condition
(designated as SS2 in Figure 3.2). The receding-wedge experiment was then initiated by
increasing the amount of freshwater flow back to the initial rate of 0.832 cm3/s. This allowed the
wedge to recede back towards the saltwater boundary and the system reached the third steady-
state condition (designated as SS3 in Figure 3.2) after 60 mins.
Numerical simulations were completed using SEAWAT. The numerical model is
assumed to be a unit-width aquifer and hence the total flow rate used was set to 0.378 cm3/s; the
real sand tank was, however, 2.2 cm wide and the measured experimental flow rate was 0.832
cm3/s (Table 3.2 provides a summary of these conversions). In the model, the initial total flow of
0.378 cm3/s was distributed into 0.0086 cm3/s over 42 confined cells of thickness 0.5 cm; in
addition, a 1-cm thick unconfined cell was used and it received a flow rate of 0.017 cm3/s. To
simulate the reduced freshwater flow condition, the flow rates for the confined cells were set to
46
0.0048 cm3/s and the top unconfined layer flow was set to 0.009 cm3/s. The results of the
numerical model are compared against the experimental data in Figure 3.3. In the figures, -dot-
shaped symbols represent experimental observations and the continuous lines represent
numerical prediction of 50% concentration contour. The figures show excellent match between
the experimental data and model predictions. The toe length, XT, measured at the bottom of the
aquifer is a standard metric used for delineating the location of a wedge [Chang et al., 2011;
Watson et al., 2010; Werner and Simmons, 2009; Werner et al., 2011]. However, in this study,
due to some physical irregularities at the bottom joint of the tank (occurred due to smearing of
the excess welding glue outside the tank) prevented us from accurately viewing the wedge at the
bottom of the sand tank. This has introduced some uncertainties in measuring the length of the
saltwater toe XT (wedge extent at the bottom). To avoid any mismatch due to these uncertainties,
we compared the values of XT@0.75 (i.e., the toe length estimated at 0.75 cm above the tank
bottom) against numerical predictions for XT@0.75. Note we selected 0.75 cm since this height
conveniently corresponded to the center of a second numerical model layer. Using this upper
location allowed us to avoid the observational problem at the bottom of the physical model.
Figure 3.4 compares XT@0.75 predicted by SEAWAT against those values estimated for the
transient experimental data for both advancing and receding experiments. The figure shows that
the numerical model was able to predict the saltwater toe lengths of both intruding and receding
wedges.
Transient advancing wedge
0 min (SS1)
Transient receding wedge
0 min (SS2)
Figure 3.2. Transient variations in the
the right boundary. The top figure shows intruding transport conditions when the flux was
reduced from 0.833 cm3/s to 0.467 cm
conditions when the flux was increased from
47
10 min 90 min (SS2)
10 min 60 min (SS3)
salt-wedge patterns due to changes in the regional
3/s, and the bottom figure shows receding transport
0.467 cm3/s to 0.833 cm3/s.
flux at
48
Figure 3.3. Comparison of experimental data model predictions of transient salt wedges for the
regional-flux system: a) advancing-wedge experiment, and b) receding-wedge experiment.
0
5
10
15
20
0 10 20 30 40 50
z (
cm
)
x (cm)
time = 0 minute
time = 10 minutes
time = 30 minutes
time = 90 minutes
SEAWAT
(a)
0
5
10
15
20
0 10 20 30 40 50
z (
cm
)
x (cm)
time = 0 minute
time = 10 minutes
time = 60 minutes
SEAWAT
(b)
49
Figure 3.4. Comparison of experimental and model-simulated toe positions estimated at 0.75 cm
above tank bottom in the regional-flux experiment
0
10
20
30
40
50
0 30 60 90 120 150 180
To
e p
os
iti
on
o
f s
alt
-w
ed
ge
, X
T 0
.75
Time (minutes)
Experiment
SEAWAT
50
3.3.2 Areal-Recharge Flux Experiments
The areal-recharge flux experiments were completed to study the extent of saltwater
intrusion when the recharge-flux varied at the top of the model. In all these experiments, the
saltwater level was fixed at 21 cm at the left boundary. Figure 3.5 shows the digital data for both
advancing- and receding-wedge experiments. The time levels and flow rates required for
conducting these areal-recharge flux experiments were estimated from preliminary test
simulations. Note that the employed inflow freshwater rate in the areal-recharge test was larger
than the regional-flux test; this increase was required to create a wedge profile similar to the
previous test (see Table 3.2 for flow data). The wedge in the areal-recharge experiment would
have advanced less into the flow tank if we would have set the volumetric inflow rates for the
area-recharge flux experiment to be identical to the regional- flux experiment. To begin the
areal-flux experiment, a constant total freshwater flow rate of 0.505 cm3/s was injected at the top
of the tank and the system was allowed to reach a steady-state pattern (designated as SS1 in
Figure 3.5). To initiate transient conditions, flow rate was instantaneously reduced to 0.343 cm3/s
(about 68% of the initial flow rate), which then allowed the wedge to advance to the right,
towards the flux boundary, by approximately 20 cm. This transient data was recorded for 90 min,
after which the system reached the second steady-state condition (designated as SS2 in Figure
3.5). Then, the receding-wedge experiment was initiated by increasing the amount of freshwater
flow back to the initial rate of 0.505 cm3/s. The wedge started to recede backward towards the
saltwater boundary and reached the third steady state (designated as SS3 in Figure 3.5) and the
transport data was recorded for 60 mins.
Numerical simulations were completed using SEAWAT. In the numerical model, the
initial areal-recharge rate of 0.010 cm/s was distributed using the recharge package over the top
unconfined layer. To simulate the
recharge rate was reduced to 0.00
experimental data in Figure 3.6.
and numerical predictions. Similar to previous experiments, in order to avoid the observational
errors associated with our tank bottom, we compared the experimental values of
0.75 cm against the corresponding
shown in Figure 3.7 compare
SEAWAT against data for both intruding and receding areal
indicate that the model was able to predict the
Transient advancing wedge
0 min (SS1)
Transient receding wedge
0 min (SS2)
Figure 3.5. Transient variations in the
The top figure shows intruding transport conditions when the flux was reduced from
to 0.733 cm3/s, and the bottom figure shows receding transport conditions when the flux is
increased from 0.733 cm3/s to 1.133
51
instantaneous reduction in the freshwater flow rate
68 cm/s. The numerical predictions are compared against the
The figure shows good match between the experimental data
XT values computed by the numerical model. The results
the transient variations in toe length (XT@0.75
-recharge flux experiment
experimental data well.
10 min 90 min (SS2)
10 min 60 min (SS3)
salt-wedge patterns due to changes in the areal flux values.
cm3/s.
, the
XT estimated at
) predicted by
s. Results
1.133 cm3/s
52
Figure 3.6. Comparison of experimental data model predictions of transient salt wedges for the
areal-flux system: a) advancing-wedge experiment, and b) receding-wedge experiment.
0
5
10
15
20
0 10 20 30 40 50
z (
cm
)
x (cm)
time = 0 minute
time = 10 minutes
time = 30 minutes
time = 90 minutes
SEAWAT
(a)
0
5
10
15
20
0 10 20 30 40 50
z (
cm
)
x (cm)
time = 0 minute
time = 10 minutes
time = 60 minutes
SEAWAT
(b)
53
Figure 3.7. Comparison of experimental and model-simulated toe positions estimated at 0.75 cm
above tank bottom in the areal-flux experiment
0
10
20
30
40
50
0 30 60 90 120 150 180
To
e p
os
iti
on
o
f s
alt
-w
ed
ge
, X
T 0
.75
Time (minutes)
Experiment
SEAWAT
54
3.3.3 Analysis of Time Scales Associated with Intruding and Receding Salt Wedge
Transport Processes
The transient data for the regional flux experiment, presented in Figure 3.4, show that the
intruding wedge requires more time to reach steady-state when compared to the receding wedge.
This trend is also observed in the recharge flux experiment shown in Figure 3.5. A previously
published experimental dataset presented by Goswami and Clement [2007] for a constant head
system also indicated a similar trend. This is an interesting, counter-intuitive result which
implies that it will take relatively less time to remediate (or reverse) saltwater intrusion effects
when compared to the time taken to contaminate (or advance the wedge into) the system. This is
a subtle yet an important result that can play a significant role in managing water quality of
coastal aquifers. However, so far no one has either quantified or analyzed the disparities in the
time scales associated with saltwater intrusion and recession processes.
In order to quantify the difference between intruding- and receding-wedge migration
rates, we first define a term ?XT(t), which is the distance that remains to be traversed by a
migrating saltwaterwedge toe to reach its final steady-state condition (note, we will be using
numerical data in this analysis and report the estimates of XT at the bottom of the aquifer). Using
this definition, for an intruding wedge ?XT(t) = ABS(XT(t)-XT,for SS2) and for a receding wedge
?XT(t) = ABS(XT(t)-XT for SS3). Note in our experiments the final steady state SS3 is identical to
the initial steady state SS1. For the regional-flux test, the values of XT for SS1 (or XT for SS3) and XT
for SS2 were estimated using SEAWAT as 19.3 and 36.8 cm, respectively. The computed values
of ?XT(t) for both intruding and receding wedges are presented in Figure 3.8(a).
The figure clearly shows that the receding wedge moves considerably faster than the
intruding wedge as seen in Figure 3.8(a). The function ?XT(t) resembles a first-order decay
55
process, and hence one could fit an exponential decay curve of the form ?XT(t) = ?XT0 exp(-t/?)
to these data; where, the parameter ?XT0 is the maximum distance traversed by a moving wedge,
and ? is a characteristic time constant. It should be noted that this approximate fitting analysis
was completed just to estimate a characteristic time for this transient problem. For the regional-
flux test, the value of ?XT0 was calculated as 17.5 cm [ABS (XT of SS1 - XT of SS2)], and the fitted
values of ? for intruding and receding wedges were, 15.0 and 7.9 minutes, respectively. The time
constant ? is a measure of the characteristic time associated with the underlying transport process
(note, here the value of ?XT(t) starts from a maximum and decays with time). The value of ?
computed for the intruding system is larger than the value computed for the receding system
indicating that the intrusion process is inherently slower than the recession process. Based on
this result one could hypothesize that it will take relatively less time for a salt wedge to recede
from an aquifer when compared to the time taken for the wedge to intrude into an aquifer. To
test the validity of this hypothesis at a larger scale, we completed a simulation study for a field
scale problem reported in Chang et al. [2011]. The problem considered a two-dimensional
unconfined aquifer system which is 1000 m long and 31 m thick with unit width. The sea level
was set at 30 m along the left boundary. A regional flux of 0.186 m3/day was injected at the
right boundary. The initial steady-state (SS1) location of the salt wedge toe, XT for SS1, was 380 m
away from the coastal boundary. To initiate transient intrusion conditions, the boundary flux
was instantaneously reduced to 0.140 m3/day (about 75% of the initial flux) which allowed the
wedge to advance towards the landward direction. It took approximately 60,000 days for the
system to reach the second steady state (SS2), and value of XT for SS2 was 485 m. The receding-
wedge experiment was then initiated by restoring the freshwater flux back to the initial rate of
0.186 m3/day. The wedge started to recede back and reached the third steady state (SS3) after
56
about 40,000 days. In Figure 3.8b, we provide the ?XT(t) profiles for both intruding and
receding simulations. The value of ?XT0 for the system is 105 m (485 m ? 380 m) and the fitted
values of the time constant ? for intruding and receding wedges are 13,600 and 9,400 days,
respectively.
To compare the differences in intruding and receding transport times scales under areal-
recharge-flux boundary conditions, we completed a new set of simulations and compiled the
?XT(t) data for the recharge condition. The initial amount of areal recharge used in the
simulation was 0.0102 cm/sec and the corresponding value of XT for SS1 was 19.7 cm. The flux
was then reduced to 0.0071 cm/sec the corresponding value of XT for SS2 was 37.3 cm. Note that
the recharge rates were selected to match the values of XT close to the regional flux scenario.
Figure 3.9 shows the SEAWAT simulated data and the fitted exponential profiles for this
problem. Similar to previous problems, the intruding wedge moved slower than the receding
wedge even under areal-recharge flux conditions. The simulated data for the receding wedge
however, did not show a smooth exponential trend. We evaluated the parameters of best fit
exponential functions to both datasets and used it to estimate the values of the time constants for
the intrusion and receding wedges as 22.7 and 10.1 mins, respectively. These results once again
confirm that the transport time required for the intruding wedge is more than the time required
for the receding wedge.
57
Figure 3.8. Comparison of toe migration rates for intruding and receding wedges for laboratory-
scale problem of (a) regional flux experimental problem and (b) field scale problem (reported in
Chang et al. 2011).
0
5
10
15
20
0 50 100 150 200
? X
T
(cm
)
Time (minutes)
SEAWAT data - Intruding wedge
SEAWAT data - Receding wedge
Fitted curve for intruding data
Fitted curve for receding data
(a)
0
30
60
90
120
0 20,000 40,000 60,000 80,000 100,000
?X
T
(m
)
Time (days)
SEAWAT data - Intruding wedge
SEAWAT data - Receding wedge
Fitted curve for intruding data
Fitted curve for receding data
(b)
58
Figure 3.9. Comparison of toe migration rates for intruding and receding wedges for laboratory-
scale problem of areal-recharge flux experimental problem.
0
5
10
15
20
0 50 100 150 200
? X
T
(cm
)
Time (minutes)
SEAWAT data - Intruding wedge
SEAWAT data - Receding wedge
Fitted curve for intruding data
Fitted curve for receding data
59
In order to understand the fundamental reason for the difference in the transport time
scales, we analyzed the underlying velocity distribution (or flow field) near the intruding and
receding wedges. Figure 3.10 shows the model-predicted vector field for the system at the two
steady state conditions SS1 and SS2 (for the regional flux system). The length of each velocity
vector presented in Figure 3.10 is directly proportional to the magnitude of the velocity. The
figures also show the 50% concentration contour to delineate the wedge location. Several
interesting observations can be made from these steady-state data. First, the vector field shows
that the velocities in the freshwater region above the wedge are much higher than the velocities
in the saltwater region (velocity vectors of salt water are small and we will use log scale to better
visualize these vectors in a later figure). It can be observed that the freshwater flow from the
inland boundary converges into the narrow zone above the salt wedge and discharges within a
high velocity region as it approaches the sea-side boundary. On the contrary, saltwater flow
velocities below the wedge are small forming a relatively stagnant region. Upon closer
observation one could observe a convection cell circulating saltwater beneath the wedge. Also,
there is a stagnation point (with opposing flow vectors) present near the toe of the wedge at the
bottom of the aquifer.
To better understand the transport dynamics under transient conditions, we plotted the
velocity fields for both intruding and receding wedge after two minutes of transport (i.e., two
minutes after the steady state condition was perturbed). This data is shown in Figure 3.11. In
order to better visualize the low velocity patterns present beneath the salt wedge, we transformed
the velocity vector field in log scale and Figure 3.12 shows the vector field computed at various
times. The first picture, Figure 3.12a, simply presents the steady-state (SS1) data shown in
Figure 3.10a in logarithmic scale. One could now clearly see the circulation cell and the
60
stagnation point (present at around 20 cm) in this log-transformed velocity vector plot. Note,
under the wedge, saltwater continuously enters the aquifer from the bottom and exists slightly
beneath the wedge, forming a convection cell.
Figure 3.11a and Figure 3.12b show the velocity field of the transient intruding salt
wedge (2 mins of transport after SS1). The figures show that the transport process distorts the
convection cell when wedge is migrating into the aquifer. All the vectors below the wedge are
now pointing towards the landward direction. Within 2 minutes, the stagnation point has been
pushed into the system by about 26 cm. The relative magnitude of saltwater velocity is still
small (see Figure 3.11a) and the saltwater flow vectors are diametrically opposite to the
freshwater flow vectors forming a stagnation region where the velocities reverse; also a distinct
stagnation point continue to persist at the aquifer bottom (see Figure 3.12b). In the intruding
system, the saltwater has to advance slowly against an opposing fresh water flow field. Once the
steady state condition is reached (Figure 3.12c) the convection cell was restored.
Figures 3.11b and 3.12d show the velocity field of the transient receding wedge (2 mins
of transport after SS2). A direct comparison of the flow fields intruding and receding flow fields
presented in Figures 3.11a and 3.11b shows that the receding system is characterized by much
higher velocities since the net freshwater flux transmitted through the system is higher in this
case. Figure 3.12d clearly shows that once the fresh water flux was increased to initiate the
receding condition, the direction of saltwater water flow is reversed towards seaward direction.
This flow reversal is a transient phenomenon and was due to the sudden increase in freshwater
flux and the associated changes in groundwater heads. There is no stagnation point in this
system. Unlike the intruding system, which supports opposing flow fields, here the flow field is
unidirectional with all the vectors pointing towards the saltwater boundary. Therefore, the
61
transport processes associated with a receding wedge can flush the saltwater rapidly out of the
system, facilitating faster transport of the salt wedge.
Figure 3.10. Simulated steady-state velocity fields for the regional flux experiment (a) SS1 and
(b) SS2.
0
5
10
15
20
0 10 20 30 40 50
z (
cm
)
x (cm)
(a)
0
5
10
15
20
0 10 20 30 40 50
z (
cm
)
x (cm)
(b)
62
Figure 3.11. Simulated transient velocity fields for the regional flux experiment (a) 2 minutes
after SS1 and (b) 2 minutes after SS2.
0
5
10
15
20
0 10 20 30 40 50
z (
cm
)
x (cm)
(a)
0
5
10
15
20
0 10 20 30 40 50
z (
cm
)
x (cm)
(b)
63
Figure 3.12. Visualization of the velocity fields near the salt wedge plotted in log scale at (a) SS1,
(b) 2 minutes after SS1, (c) SS2, and (d) 2 minutes after SS2.
0
5
10
15
0 10 20 30 40
z (
cm
)
x (cm)
12(a) ? SS1
0
5
10
15
0 10 20 30 40
z (
cm
)
x (cm)
12(b) ? 2 minutes after SS1
0
5
10
15
0 10 20 30 40
z (
cm
)
x (cm)
12(c) ? SS2
0
5
10
15
0 10 20 30 40
z (
cm
)
x (cm)
12(d) ? 2 minutes after SS2
64
3.4 Conclusions
Sea-level rise and reduction of recharge due to changes in rainfall patterns and the
associated changes in groundwater fluxes are the two major climate-change-induced
hydrological processes that can affect saltwater intrusion processes in coastal aquifers. The
objective of this research was to conduct laboratory-scale experiments to investigate the impacts
of variations in regional and areal-recharge fluxes on saltwater intrusion in unconfined aquifers.
We have successfully completed laboratory experiments to study two types of flux-controlled
saltwater intrusion problems: one controlled by the regional flow delivered from the inland
freshwater boundary, and the other controlled by the recharge flow delivered from the top
boundary. The results of the experiments are used to develop a better understanding for the
influence of changes in boundary fluxes on transient saltwater intrusion processes. The results
were also modeled using the numerical code SEAWAT. The experimental data presented in this
study are useful benchmarks for testing the validity of density-coupled flow and transport models
involving flux-type boundary conditions.
Based on our experimental data we hypothesized that the time scales associated with an
intruding wedge would be higher than a receding wedge. This hypothesis was confirmed using
numerical simulations completed for problems involving two different scales (1-m and 1000-m
scales). Our analysis also shows that the intruding system supports two opposing flow fields
with a distinct stagnation point at the bottom of the aquifer, and a diffused stagnation zone along
the wedge where the flow velocities reverse direction. On the other hand, receding systems
support a unidirectional flow field with all the velocity vectors pointing towards the saltwater
boundary. Due to the presence opposing flow fields and the stagnation zone, the time taken for
saltwater to intrude into an aquifer will be relatively high when compared to the time taken for
65
the saltwater to recede out of an aquifer (which has a well-aligned unidirectional flow field).
The insights gained from this study help us better understand the transport processes occurring
within transient saltwater wedges, and these processes can be coupled to multi-species reactive
transport formulations [Clement et al., 1998; Sun and Clement, 1999] to predict the exchange of
nutrients and other contaminants between terrestrial and marine waters. Overall, this study has
the following two major contributions: 1) the experimental efforts report transient saltwater
intrusion data for systems involving two types of flux boundary conditions; and 2) the
experimental observations together with the numerical results point out the time-scale
differences between intruding and receding saltwater wedges and offer a mechanistic explanation
for their difference.
66
Chapter 4
LABORATORY AND NUMERICAL INVESTIGATION OF TRANSPORT PROCESSES
OCCURRING BENEATH A SALTWATER WEDGE
4.1 Introduction
Due to rapid urban development, shallow aquifers in coastal areas have become highly
vulnerable to contaminants dissolved in urban runoff water and/or contaminants inadvertently
discharged from leaking underground storage tanks or landfills that are located close to the
shoreline. These contaminants not only pollute the local unconfined aquifer, but they also have
the potential to contaminate the entire shoreline beach environment. This is due to the presence
of a saltwater wedge, a common feature in coastal aquifers, which can influence transport
pathways and force groundwater plumes to discharge along the beach face. Westbrook et al.
[2005] completed a field study to document the transport of a dissolved hydrocarbon
groundwater plume flowing towards a tidally and seasonally forced saline water body in Perth,
Western Australia. Their field data indicate shore-focused discharge where the groundwater
seepage occurred along beach-face sediments within a zone less than 10 m from the high tide
mark. The plume discharge was impacted by the transient movement of the saltwater wedge.
This study illustrated the need for understanding the groundwater flow dynamics in the vicinity
of the wedge while predicting the exchange of contaminants between terrestrial to marine waters.
Figure 4.1 is a conceptual diagram that illustrates how contaminant discharged in coastal
aquifers will be impacted by the presence of the salt wedge. The anthropogenic contaminants
67
discharged within the shallow freshwater zone (well above the wedge) would be transported
towards the beach face. In addition, as shown in the figure, coastal groundwater systems can
also circulate nutrients or other contaminants of marine origin within the saline water present
beneath the wedge. The density contrast between the fresh groundwater and saline ocean water
induces a circulation cell below the wedge which can enable exchange solutes between marine
and terrestrial waters. Saltwater flow induced by this circulator motion is an important
component of the submarine groundwater discharge (SGD) [Destouni and Prieto, 2003; Li et al.,
1999; Smith, 2004]. It has been well established that chemicals entering coastal water via SGD
have considerable consequence on the marine ecosystems [Johannes, 1980; Moore, 1996; Moore
and Church, 1996; Robinson et al., 2007; Simmons, 1992]. Typical groundwater flow velocities
within this circulatory cell would be relatively small (compared to the freshwater velocities
above the wedge) and hence the transport within the wedge is expected to be dispersion
dominated.
Figure 4.1. Conceptual transport model in unconfined coastal aquifer.
Sea Level
Salt wedge toe, XT
Contaminant
Contaminant
Freshwater level
?
68
In the literature, only a few studies have focused on conducting physical experiments to
investigate the contaminant transport processes occurring near a saltwater wedge. Zhang et al.
[2002] conducted laboratory experiments to visualize contaminant transport above a wedge.
They completed a set of tracer transport experiments that tracked the movement of a dense
contaminant plume discharged above the salt wedge. Their observations showed that a dense
contaminant plume became more diffusive compared to the less dense plumes as it approached
the narrow discharge zone near the coastline. Non-dimensional parameters were developed to
qualitatively characterize the plume; no attempts were made to reproduce the observations using
any numerical models. Volker et al. [2002] completed a numerical modeling study to investigate
the transport of a dense contaminant plume in an unconfined aquifer. Their simulation results
were also compared against experimental data. The study modeled the sea-side boundary using
two types of boundary conditions: a) as constant head freshwater boundary that neglected density
effects, and b) as a constant head saltwater boundary that included density effects. Their results
indicated that when the plume is away from the boundary, predictions based on both approaches
matched reasonably well with experimental results. However, as the plume approached the coast,
the experimental data indicated that the plume will be uplifted along the saltwater interface and
finally exited along the shoreline. This observed uplifting and shore-focused discharge patterns
were recreated only by the model that considered density effects. The model that ignored density
effects predicted an incorrect travel path where the plume travelled further into the sea and exited
under the seabed. The study concluded that ignoring seawater intrusion effects would not only
underestimate the amount of contaminant mass exiting along the beach face, but it would also
result in simulating unrealistic transport pathways.
69
A few laboratory studies have focused on understanding the effects of tides on
controlling transport patterns in coastal aquifers near the groundwater-ocean boundary. For
example, Boufadel [2000] studied the transport pathways of nutrients used for remediating
beach-face contaminants resulting from the past oil spills. Laboratory experiments using two
types of density gradients were completed to investigate the effects of tides and buoyancy on
beach hydraulics. In both experiments, the authors observed that the groundwater flow from the
aquifer moved seaward via the intertidal zone and pinched out near the low tide mark. The
experimental data was modeled using a numerical simulator, which was then used to design
nutrient delivery strategies for bioremediating oil contaminated beaches. Kuan et al. [2012]
completed laboratory experiments to study tide-induced circulation processes that would lead to
the formation of an upper saline plume within the intertidal zone. They found that the presence
of the upper saline plume moved the fresh groundwater discharge zone towards the low tide
mark and restricted the movement of the salt wedge. Illangasekare et al. [2006] presented a
laboratory dataset that characterized the movement of dense saline plumes in a tsunami impacted
aquifer. Their laboratory study first simulated the pre-tsunami initial conditions by setting up a
constant regional hydraulic gradient. Their experiments used a variety of dyes to study the
mixing of different types of salt waters and their interactions with regional freshwater. Initially,
an uncolored fresh groundwater (identified as ??freshwater??) plume was allowed to flow from
right to left toward the dense seawater boundary, which was marked green. The green seawater
was allowed to intrude the aquifer from the left to form a steady-state saltwater wedge. The
inundation caused by the tsunami wave was then simulated by evenly discharging a fixed amount
of red-colored saltwater, identified as ??tsunami?? water. In addition, a blob of red-colored
tsunami water was also injected directly into the aquifer to mimic saltwater flow via open wells.
70
Qualitative observations were made to understand the fate and transport of different saltwater
plumes. The study used non-dimensional analysis to characterize the transport times, however
no quantitative numerical assessments were made. Use of different types of colored waters
allowed researchers to visualize various types of mixing effects; in the current effort, we employ
a similar approach to study contaminant plume migration processes occurring near a saltwater
interface.
Above studies have demonstrated the importance of understanding the role of a salt
wedge when predicting the movement of a contaminant plume in coastal aquifers. While several
experimental efforts have focused on studying contaminant transport above a wedge, as per our
knowledge, so far no one has studied the fate and transport of contaminants beneath a saltwater
wedge. Kohout [1960] was the first to complete a field study in Florida and estimated that the
saltwater flux circulated through the system is about 10% of the fresh groundwater flux flowing
towards the sea. Destouni and Prieto [2003] completed a modeling study to develop empirical
correlations between groundwater fluxes present above and below a saltwater wedge. Their
correlations also show that the saltwater flux is roughly an order of magnitude smaller than the
freshwater flux. Smith [2004] completed a comprehensive modeling study and concluded that
the amount of saltwater flux circulated within the system can vary widely and would depend on
the level of dispersion in the system. However, most of these studies are based on theoretical
modeling data and so far no one has provided experimental data to quantify the circulating
saltwater flux and the associated transport processes.
The objective of this effort is to conduct well-controlled contaminant transport
experiments to study the dynamics of groundwater flow transport patterns beneath a saltwater
wedge and relate the transport to the flow patterns present above the wedge. We used a novel
71
experimental approach that employed a variety of dyed waters to map and compare the mixing
and transport patterns occurring above and below a saltwater wedge. The experiments were
simulated using a numerical model. The model was then used to investigate these transport
problems in a larger scale problem.
72
4.2 Methods
4.2.1 Experimental Method
Figure 4.2 shows the experimental setup used in this study. The dimension of the flow
tank was 50.5 cm ? 28 cm ? 2.2 cm. The origin of the Cartesian coordination system used to
record experimental observations was set at the left bottom corner of the tank. The tank was
fitted with three distinct flow chambers. The right inlet chamber was used to supply freshwater
flow delivered using a peristaltic pump. The sea level was set to 25.1 cm by adjusting the
saltwater head at the left outlet chamber. The central porous media chamber was filled with
silica glass beads of diameter 1.1 mm. The beads were packed under saturated conditions to
prevent the trapping of air bubbles inside the tank. Thin stainless-steel meshes were used to
separate the side chambers (boundary conditions) from the porous media chamber. Freshwater
level at the right boundary was set at 26.03 cm. The hydraulic conductivity (K) of the tank was
measured using the in-situ method described in Simpson et al. [2003]. The average value of
conductivity was 1.03 cm/s. The measured values of saltwater and freshwater densities were
1.025 g/cm3 and 1.000 g/cm3, respectively. All the experimental parameters used in this study
are summarized in Table 1. The transport experiments were recorded using a video camera
(Panasonic HDC-HS250). Other experimental procedures used were similar to the ones reported
in our previous studies [Abarca and Clement, 2009; Goswami and Clement, 2007; Goswami et
al., 2009]. However, in this study, as shown in Figure 4.2, the flow tank was slightly modified to
accommodate a long needle and a short needle of lengths 300 mm and 50 mm, respectively, and
thickness of 1.6 mm. These needles were buried in the tank to allow tracer injection above and
below the wedge. Both needles were set in the tank when the porous medium was wet packed.
In addition to this physical modification, in this study we employed a multi-colored tracer
73
scheme to map various transport processes. The scheme used three different types of waters to
discern the colored freshwater (dyed red) and dense tracer water (dyed green) colorless the
seawater (colorless). Note that in all our previous experiments, the saltwater wedge was colored
red. In this effort, we reversed the coloring scheme and created a colorless salt wedge, which
allowed us to better visualize the tracer transport patterns occurring under the wedge. The
freshwater was marked using 1 ml of red food color for every liter of water. The dense tracer
solution was prepared by mixing 1 ml of green dye for every 250 ml of saltwater.
4.2.2 Details of Tracer Transport Experiments
Three types of tracer solutions with density values of 1.000, 1.0125 and 1.025 g/cm3 were
prepared to conduct three different transport experiments (identified as Test-1, Test-2 and Test-
3). Test-1 and Test-2 are baseline reference transport experiments that involved injection of 22
ml of tracer slug with density values of 1.000 (neutral tracer) and 1.0125 (dense tracer),
respectively, above the saltwater wedge very close to the freshwater boundary (see Figure 4.2 for
details). The volume of the tracer slug was 22 ml which was injected in about 8 seconds,
yielding an injection rate of 2.75 ml/s. The x-z coordinates of the injection point are: x = 45 cm
and z = 22.5 cm. The tracer slugs were allowed to transport with the freshwater flow and the
images were recorded at different times until the plume reached the exit boundary. The result of
Test-1 was used to calibrate a numerical model by adjusting the dispersivity values, and Test-2
data was used to test the model performance.
Test-3 was the most important tracer study that recorded the transport patterns occurring
beneath the wedge. The coordinates of the injection point are: x = 4.0 cm and z = 0.5 cm. About
11 ml of saltwater solution with a density value of 1.025 g/cm3 (which is identical to the
74
seawater in the left tank) was injected beneath the wedge. The injection time was 40 s, yielding
an average injection rate of 0.275 ml/s. Note since the ambient density of water present beneath
the wedge was 1.025 g/cm3, and hence the tracer slug acted like a neutral tracer within the salt
wedge.
4.2.3 Numerical Modeling Method
The MODFLOW-family variable density flow code SEAWAT [Guo and Langevin, 2002]
was used to model the experimental data. A uniform two-dimensional grid of dimension 0.5 cm
(102 cells in the x-direction and 50 cells in the z-direction) was used to discretize the flow tank.
Thickness of numerical cells was set to 2.2 cm, which is identical to the thickness of the tank.
The saltwater head at the left boundary was set to 25.1 cm for the entire simulation. Freshwater
was allowed flow from the constant-head boundary set in the right hand side tank. The value of
specific storage (Ss) was set to 10-6 cm-1. Similar to the Goswami and Clement [2007] study, we
used the confined flow approximation to simulate this steady-state experiment. All fifty model
layers (of thickness 0.5 cm) were set as confined layers, and the boundary heads (i.e., the water
levels at both saltwater and freshwater boundaries) always remained slightly above the top layer.
The bottom of the model was set to a no-flow boundary condition. The values of longitudinal (?L)
and transversal (?T) dispersivity were initially selected based on the values suggested by Abarca
and Clement [2009], and were adjusted to fit the neutral tracer transport (Test-1) data and were
later verified using dense tracer transport (Test-2) data. The total variance diminishing (TVD)
scheme solver was selected for simulating all the lab experiments. A Courant number stability
constraint of Cr ? 0.5 was used in all the simulations.
75
Table 4.1: Model parameters used for simulating the lab problem
Parameter Symbol Value
Horizontal length L 50.5 cm
Saltwater elevation Hs 25.1 cm
Saltwater elevation Hf 26.03 cm
Hydraulic conductivity K 1.03 cm/s
Specific storage, Ss 1?10-6 cm-1
Longitudinal dispersivity ?L 0.05 cm
Transverse dispersivity ?T 0.005 cm
Saltwater concentration Cs 0.035 g/cm3
Saltwater density ?s 1.025 g/cm3
Freshwater density ?f 1.000 g/cm3
Density slope E 0.714
Porosity n 0.385
76
Table 4.2: Model parameters used for simulating the field problem (modified from Chang et al.
2011)
Parameter Symbol Value
Horizontal aquifer length L 1000 m
Vertical aquifer thickness B 30 m
Inland groundwater flow rates Q1 (low flux) 0.105 m3/day
Q2 (High flux) 0.210 m3/day
Hydraulic conductivity K 10 m/day
Specific storage, Ss 1?10-5 m-1
Longitudinal dispersivity ?L 0.4 m
Transverse dispersivity ?T 0.04 m
Saltwater concentration Cs 35 kg/m3
Saltwater density ?s 1,025 kg/m3
Freshwater density ?f 1,000 kg/m3
Density Slope E 0.714
Porosity n 0.35
77
Figure 4.2. Schematic diagram of the flow tank equipped with solute injection syringes.
78
4.3 Results
4.3.1 Transport of Tracer Slugs Injected above the Wedge
Figure 4.3 shows time series of pictures of the neutral tracer slug as it travelled from the
injection point, located near the freshwater boundary, to the saltwater boundary. The
photographs were taken at 0, 5, 10 and 15 minutes after injecting the slug. The experimental
data indicated that the plume required about 20 minutes for traversing the entire flow domain.
The figure also shows the model simulated plume at different times. In the figure, the salt wedge
is identified by white solid line which represents the 50% isochlor. The color filled plumes were
generated by the contouring software that used the following color coding scheme: red > 0.8 %;
yellow 0.8?0.7 %; green 0.7?0.3 %; turquoise 0.3?0.05 %; and blue < 0.05 %. As a part of this
simulation study, the dispersion parameters ?L and ?T were adjusted within an expected range to
reproduce the observed plume spreading levels. Based on qualitative comparison, dispersivity
values of ?L = 0.05 cm and ?T = 0.005 cm were able to recreate the observed plume spreading
patterns. These values are within the range of values reported in a previous study [Abarca and
Clement, 2009] that used a similar type of flow tank. All the model parameters used for
simulating this baseline experiment are summarized in Table 1. Both lab and model-simulated
results (see Figure 4.3) indicate that the shape and spatial extent of the plume lifted upward as it
approached towards the discharge boundary. Overall, the model simulated tracer transport
patterns agreed well with the experimental data.
The second tracer experiment was completed to primarily assemble a dataset to test the
performance of the model. In addition, the second tracer test also explored how a dense plume
would behave as it sunk and interacted with the wedge. All the model parameters including the
dispersion coefficients used were identical to the values used for simulating the Test-1 dataset.
79
Figure 4.4 shows the experimental data for the dense tracer plume as it travelled from the
injection point to the saltwater boundary. The photographs were taken at 0, 5, 10 and 15 minutes
after slug injection. The dense plume also took about 20 minutes to traverse the flow tank. The
figure shows that the model-simulated plumes at different times
It is interesting to note that the dense plume initially started to sink and the sinking
process continued for about 10 minutes; during this time, the core of the plume sank by about 10
cm in the vertical direction. The model was able to reproduce this sinking pattern. As the dense
plume approached the salt wedge, the flow direction started to change and the tracer plume
started to rise upwards due to the presence of saltwater wedge that almost acted as a flow barrier.
Since the freshwater flow regime was characterized by strong advection the dense plume
remained stable throughout the experiment. Kanel et al. [2008] studied the transport of a dense
plume of nano-particles under strong advection-dominated flow conditions, and their data also
showed such a stable dense plume transport pattern. Goswami et al. [2012] recently published a
comprehensive study to investigate the dynamics of sinking and buoyant plumes and
demonstrated the role of advection in inhibiting plume instabilities.
80
Time = 0min Time = 5 min Time = 10 min Time = 15 min
Fig.4.3. Comparison of observed data with numerical predictions for the neutral plume transport
experiment (Top figure: lab data; Bottom figure: SEAWAT predictions with color filled
contour -- red > 20 %; yellow 20-14 %; green 14-8 %; and turquoise 8-1 %; also, white line is
the 50% saltwater contour of the wedge)
Time = 0min Time = 5 min Time = 10 min Time = 15 min
Fig. 4.4. Comparison of observed data with numerical predictions for the dense plume transport
experiment (Top figure: lab data; Bottom figure: SEAWAT predictions with color filled
contour -- red > 20 %; yellow 20-14 %; green 14-8 %; and turquoise 8-1 %; also, white line is
the 50% saltwater contour of the wedge)
81
4.3.2 Characterizing Tracer Transport Patterns beneath a Saltwater
The above two tests allowed us to estimate the transport parameters for the experimental
system, and also allowed to visualize the transport patterns occurring within the freshwater zone,
which served as baseline information. Test-3 was the most important experiment that directly
mapped the transport patterns occurring beneath the wedge. Our preliminary experiments
indicated that transport under the wedge would have very little advection, and hence even a small
fraction of density difference between the ambient water and the tracer slug would lead to
buoyancy effects. Therefore, we carefully matched the density of the green-dyed dense tracer
solution as identical to the ambient saltwater solution. Figure 4.5 shows the migration patterns
of this neutral (neural since the tracer density was matched to the ambient saltwater density)
plume release beneath the saltwater wedge. The digital pictures were taken at 0, 60, 120 and 180
minutes after injecting the tracer slug. The figure also shows the model predicted results, which
are in good agreement with observed data. Unlike the transport in the freshwater domain, the
transport beneath the wedge was extremely slow. The tracer slug was circular in shape at the
time of injection; however, the circular slug started to degenerate rapidly and the evolved into an
elongated plume of irregular geometry within a relatively short period. The bulk of the plume
mass remained at the bottom of the tank with a ?finger like? structure evolving upward; the
fingers rapidly extended and were flushed out of the system by the strong advective flow active
near the interface. Due to dilution of the plume mass, it was difficult to delineate a distinct
plume in the digital photographs taken after about 4 hours of transport; however, some residual
mass was still visible in the tank. It took approximately 6 hours (estimated based on visual
82
observation of the remnant green residues in the tank) for the tracer slug to flush out of the
saltwater region.
Several interesting observations can be made from this laboratory dataset. First, the
plume residence time within the saltwater region was about 6 hours, which is much higher than
the 20 minute residence time estimated for the freshwater region. The numerical model also
required about 6 hours to fully flush the tracer slug out of the system. Second, the plume was not
tracking the idealized circular flow path commonly depicted in textbooks [Fetter, 2001;
Goswami et al., 2009]. Due to dispersion effects, the plume attained an elongated shape as it
approached the interface. Once the plume reached the interface, a diffused fraction of the plume
transported rather rapidly along the mixing zone and discharged out as a narrow finger. The
model was able to reproduce these complex transport patterns well.
To test the scalability of the observed transport patterns, we completed numerical
simulations for a field scale problem, which based on a modified version a problem discussed in
Chang et al. [2011]. The problem considered a two-dimensional confined aquifer which is 1000
m long, 30 m deep and 1 m wide. Dispersivity coefficients are set to: ?L = 0.4 m and ?T = 0.04
m. Porosity was assumed to be 0.35. The sea level was fixed at 30 m along the left side of the
model boundary. A constant amount of regional freshwater flow was allowed to flow from the
inland boundary, and areal-recharge was set to zero. Numerical simulations were completed for
two conditions that modeled low and high regional groundwater flow scenarios which used the
flow rates of 0.105 m3/day and 0.210 m3/day, respectively. These relatively low levels of fluxes
were selected to allow simulation of relatively long salt wedges. In order to study the tracer
transport patterns, two tracer slugs (one above the wedge and one beneath the wedge) were
released simultaneously by setting the initial condition of four model grid cells (total grid volume
83
of 64 m3) to 100 kg/m3. To be consistent with Chang et al. (2011) study, we used the fully
implicit method option available in SEAWAT for all our field scale simulations. Other model
parameters used were identical to the ones used in Chang et al. [2011] study (details are
summarized in Table 4.2). Steady-state location of the salt wedge toe for the high flow rate
problem was found to be around 390 m away from the coastal boundary and for the low flow rate
problem it was around 730 m away from the boundary. Figure 4.6 summarizes the tracer
transport scenarios predicted under both high and low flow conditions. In the figure, the salt
wedge is identified by a black solid line which represents the 50% isochlor. The color filled
plumes were generated by the contouring software that used the following color coding scheme:
red > 0.8 kg/m3; yellow 0.8?0.7 kg/m3; green 0.7?0.3 kg/m3; turquoise 0.3?0.05 kg/m3; and blue
< 0.05 kg/m3. The figures show the observed plume patterns after 360, 20,000, 40,000, and
60,000 days of transport. The data show that after 40,000 days of transport, the center of the
freshwater slug that was released above the wedge travelled about 400 meters and 800 meters
under low and high flux conditions, respectively. Based on these transport times the average
pore-scale transport velocity calculated for these systems are 1 and 2 cm/day, respectively, which
agrees with the expected values. The presence of the salt wedge lifted the freshwater plumes,
and at later times, due to converging flow conditions, the center mass moved slightly faster than
the average freshwater velocity. During the same transport period (of 40,000 days), the saltwater
slug, which was released beneath the wedge, travelled only about 170 m and 220 m under low
and high flux conditions, respectively. Once the saltwater plume reached the plume wedge, it
fingered along the wedge and flushed out rapidly from the system. The overall transport patterns
predicted by the field scale model were quite similar to those observed in our experimental tank.
Due to dispersion effects, the circular transport pathway could not be established within the
84
system; once the slug reached the wedge the bulk of the contaminant mass transported rather
rapidly along the wedge. The flow budget analysis indicated that the amount of saltwater flow
circulated within the low and high flux systems were 0.069 and 0.096 m3/day, respectively.
These fluxes primarily control the transport processes occurring beneath the wedge. Therefore,
we completed further simulation studies to quantify the sensitivity these fluxes to changes in
model parameters (dispersivity values) and boundary forcings.
Time = 0 min Time = 60 min Time = 120 min Time = 180 min
Fig.4.5. Comparison of observed data with numerical predictions for the transport experiment
completed beneath the wedge (Top figure: lab data; Bottom figure: SEAWAT predictions with
color filled contour -- red > 20 %; yellow 20-14 %; green 14-8 %; and turquoise 8-1 %; also,
white line is the 50% saltwater contour of the wedge)
85
Fig. 4.6. Model simulated results illustrating the transport processes occurring above and below
the saltwater wedge for the field-scale test problem: (a) simulations for low regional flow of
0.105 m3/day and (b) simulations for high regional flow of 0.210 m3/day. Initial contaminant
concentration at the sources was 100 kg/m3.
360 day
20,000 days
40,000 days
60,000 days
(a)
0
10
20
30
0 200 400 600 800 1000
0
10
20
30
0 200 400 600 800 1000
0
10
20
30
0 200 400 600 800 1000
0
10
20
30
0 200 400 600 800 1000
x (m)
(b)
0
10
20
30
0 200 400 600 800 1000
0
10
20
30
0 200 400 600 800 1000
0
10
20
30
0 200 400 600 800 1000
0
10
20
30
0 200 400 600 800 1000
x (m)
86
4.3. 3 Sensitivity of Recirculating Saltwater Flux to Dispersion and the Type of Boundary
Condition
Smith [2004] pointed out that the saltwater flux circulated within the system under steady
state conditions is strong function of the level of mixing (value of dispersivity) in the system.
Therefore, we first completed numerical experiments to characterize impacts of dispersion on
our field-scale problem at different boundary flux values. Simulations were completed by
varying the regional flux from 0.15 to 1.8 m3/day; the model was run for a sufficiently long
period until steady state conditions were reached. Two sets of simulations were completed using
longitudinal dispersivity values of 0.4 m and 4 m; the ratio of the transverse to longitudinal
dispersivity value was set to 0.1. Furthermore, additional simulations were also completed using
identical model parameters and freshwater flux values, but the freshwater deliver point was
modified; instead of injecting the flow from the right regional boundary, we distributed the same
amount of flow from the top of the model to simulate a system involving areal recharge flux
boundary condition. The regional flux was set to zero in these simulations. Chang and Clement
[2012b] completed laboratory experiments to compare saltwater intrusion problems involving
regional and areal recharge boundary conditions. Their data indicate that while the extent of
saltwater intrusion is highly sensitive to the magnitude of the freshwater flux, the extent is almost
insensitive to the type of boundary (recharge or areal) condition. Their results also revealed that
the wedge in the regional-flux experiment was thick and almost linear, whereas the wedge in the
areal-flux experiment was relatively thin and shaped like an exponential decay curve. However,
the impact of these differences on saltwater recirculation patterns was not investigated in Chang
and Clement [2012b] study.
87
Figure 4.7 summarizes the results of this sensitivity study. Firstly, as expected, the
amount of saltwater flux entering system was proportional to the amount of freshwater flux
entering the system. Furthermore, when the dispersivity value was increased, the amount of
saltwater flux circulated within the system also increased. This is because, higher dispersion
resulted in more mixing which facilitated saltwater flushing; this resulted in driving a cyclic
process that allowed more saltwater to enter the system. Interestingly, the relationship between
saltwater and freshwater fluxes was not sensitive to the type of boundary condition used for
modeling the inlet boundary; i.e., the location where the freshwater flux is allowed to enter the
system. The results show that both regional and area flux boundary conditions resulted in
inducing almost identical amount of saltwater flux. We completed an analytical analysis to
generalize some of these results using the non-dimensional parameters proposed by Smith [2004].
In order to normalize the relationship between freshwater inflow Qf and saltwater inflow Qs, we
used a modified version of a dimensionless freshwater flux parameter proposed by Smith (2004),
which is defined as:
uni0056g1499g3404uni03B2uni004Buni0042uni0057uni0051
g2916
(1)
Where uni03B2 is the fluid density ratio, uni004Buni0020g4670uni004Cuni0054g2879g2869g4671 is the hydraulic conductivity, B [L] is the saturated
aquifer depth, uni0051g2916uni0020g4670uni004Cg2871uni0054g2879g2869g4671 is the rate of freshwater inflow, and W is the aquifer width of the
aquifer [L]. The amount of saltwater circulated within the system was quantified using the
following dimensionless parameter (defined as the percent saltwater circulation (PSC) in Smith
[2004] study):
uni0050uni0053uni0043g3404uni0051g2929uni002Funi0051g2916 (2)
88
Where uni0051g2929uni0020g4670uni004Cg2871uni0054g2879g2869g4671 is the rate of saltwater flow into the system. As pointed out by Smith (2004),
the parameter uni0056g1499 is the ratio of free convection velocity to forced convection in the system. Also,
uni0056g1499 is simply the inverse of the non-dimensional parameter ?a? commonly used in Henry
problems [Henry, 1964; Simpson and Clement, 2004]. The derivation for Henry problem?s
parameters are shown in Appendix 6. Figure 4.8 presents the data reported in Figure 7 in a non-
dimensionless form. Visualizing the data in this non-dimensional format allowed us to compare
the data for the current problem with other literature-derived data (which are shown marked
using various symbols in Figure 4.8). The figure shows that the percentage of saltwater
circulated within the system scales well with the non-dimensional freshwater flux parameter V*.
Furthermore, the values V* and PSC predicted for the current problem are within the range of
values observe in other published problems of different scales. Similar to Smith (2004) study,
the relationship between PSC and V* shows a non-linear trend, and it is a strong function of the
mixing level (or value of dispersivity). The figure shows that the non-dimensional format
proposed by Smith [2004] provides a robust framework for correlating the relationship between
freshwater and saltwater fluxes for a variety of problems of different scales.
89
.
Fig.4.7. Relationship between freshwater flow and recirculate-saltwater flow for regional and
areal recharge problems with different dispersivity values
0
0.2
0.4
0.6
0.000 0.500 1.000 1.500 2.000
Sa
ltw
at
er
Fl
ow
Q
s
(m
3 d
-1
m
-1
)
Inland Freshwater Flow Qf (m3d-1m-1)
Regional flux - Dispersivity 0.4
Regional flux - Dispersivity 4
Areal flux - Dispersivity 0.4
Areal flux - Dispersivity 4
90
Fig. 4.8 Non-dimensional analysis of the relationship between freshwater and saltwater fluxes
[The symbols indicate the following literature data points: + from Chang et al. (2011) study, square6
regional flux experimental data from Chang and Clement (2012), box3 areal flux experiment data
from Chang and Clement (2012), triangleup Cyprus case study from Destouni and Prieto (2003), triangleopenup
Israel case study from Destouni and Prieto (2003)]
0%
20%
40%
60%
80%
100%
120%
0 10 20 30 40
Pe
rc
en
t S
alt
wa
te
r C
irc
ul
at
io
n
PS
C
(%
)
Dimensionless flux parameter V *
Regional flux - Dispersivity 0.4
Regional flux - Dispersivity 4
Areal flux - Dispersivity 0.4
Areal flux - Dispersivity 4
91
4.4 Conclusions
In this effort we have completed both laboratory and numerical experiments to study the
contaminant transport patterns occurring beneath a saltwater wedge. For comparison purposes,
we have also completed transport experiments to map contaminant transport patterns occurring
above the wedge. Results indicate that the time scales associated with the transport processes
occurring below the wedge are distinctly different from the transport processes occurring above
the wedge. The typical solute residence times associated with steady state saltwater circulation
cell beneath a saltwater wedge are at least an order of magnitude higher than the residence times
associated with the freshwater flow occurring above the wedge. The experimental dataset
showed that it took approximately 6 hours to fully flush the tracer slug when it was released
beneath the wedge, and on the other hand, it took about 30 minutes to flush the tracer slug when
it was released above the wedge. The experimental data also indicated that the tracer slugs
released below the wedge did not follow an idealized circular flow path commonly depicted in
textbooks [Fetter, 2001; Goswami et al., 2009]. Due to dispersion effects, the plume attained an
elongated shape as it approached the interface. Once the interface was reached, a diffused
fraction of the plume transported rather rapidly along the mixing zone and discharged out as a
narrow finger. The modeling results for both small and large scale systems showed similar
transport patterns. The modeling data also confirmed that dispersion plays an important role in
controlling the amount of saltwater flux recirculated beneath a wedge. Model simulations show
that when mixing across the wedge was allowed to increase (by using high dispersivity values)
the amount saltwater flux circulated beneath the wedge also increased. The percentage of
recirculated saltwater flow in a system correlates well with a dimensionless flux parameter,
proposed by Smith [2004], for a variety of problems of different scales. While the experimental
92
results showed that the saltwater flux was about an order magnitude less than the freshwater flux,
the theoretical modeling results show that when the value of dispersion is high the saltwater flux
could be comparable (or even be higher) than the freshwater flux. These theoretical results were
consistent with the results reported in Smith (2004).
The insights gained from this study help us better understand the solute exchange
processes occurring between terrestrial and marine waters. The current and Smith (2004) studies,
however, have focused only on understanding steady-state systems with a constant sea level
boundary condition. Further studies are needed to study dynamic systems which can be
influenced by transient boundary conditions which are impacted by tidal cycles and other long
term changes such as sea level rise. It will be also interesting to quantify the impacts of mixing
due to fluctuations transport velocities induced by changes in boundary conditions, changes in
recharge fluxes, and perturbations arising from spatial heterogeneities.
93
Chapter 5
CONCLUSIONS AND RECOMMENDATIONS
6.1 Summary and Conclusions
The overall goal of this research is to develop a better scientific understanding of the
dynamics of a saltwater intrusion and the associated transport processes in coastal groundwater
aquifers. Our first objective of this investigation is to understand the long term impacts of sea-
level rise on saltwater intrusion processes in confined and unconfined coastal aquifers. We
completed numerical simulations for a field-scale problem with a sea level boundary rising with
time. Numerical results provided conceptual insights into the transient effects of sea-level rise
on saltwater intrusion in coastal aquifers. The simulation results demonstrated that the rising sea
would initially force the wedge to migrate inland; however, the rise would also result in lifting of
the regional groundwater table which in-turn would naturally drive the wedge back towards the
sea. Due to this self-reversal process, the impact of sea level rise will be negligible on the
steady-state salt wedge if the total freshwater discharge in the system remained constant. The
self-reversal process will also be active in unconfined aquifers; however, in unconfined systems
there would be some level of intrusion caused due to the increase in the value of aquifer
transmissivity (or saturated aquifer thickness). However, it is important to note that this
numerical investigation ignored inundation of low lying coastal areas due to sea-level rise and
the associated movement of the coastline. Future studies should investigate these effects.
The result of the first phase of this study demonstrated the amount of freshwater flux
transmitted through the system is most important fundamental parameter that would have the
94
greatest impact on the dynamics of saltwater wedges in coastal aquifers. Therefore, the second
objective of this study is to conduct experimental studies to investigate the dynamics of saltwater
intrusion processes in an unconfined aquifer involving different types of flux-type boundary
conditions. We investigated the influence of variations in regional flow and recharge flows
(resulting from changes in local and regional rainfall patterns, respectively) on a salt wedge
using a laboratory-scale aquifer model. Both transient and steady-state data were collected from
two types of experiments namely the areal-recharge experiment and regional-flow experiment.
The experimental results were simulated using the SEAWAT numerical code. The experimental
data themselves are unique and can be used as benchmarks for testing numerical codes that use
different types of flux-type boundary conditions. In addition, based on the results of the transient
experiments, we hypothesized that the times scales associated with the salt wedge intrusion step
(where the wedge moved in towards the land due to reduction in groundwater recharge) is larger
than the times scales associated with salt wedge recession step (where the wedge moved back
towards the sea due to increase in recharge). This implies that if natural fluxes are restored, it
will require relatively less time to drive the wedge backward and restore a saltwater
contaminated aquifer. Clearly, this is an important result that has an enormous practical
implication for managing coastal aquifers.
Many of the shallow coastal unconfined aquifers are routinely contaminated by various
type of chemical discharges released from a variety of anthropogenic sources (Westbrook et al.,
2005). Furthermore, regional groundwater flow can also transport large amount of nutrients,
which then could be exchanged across the salt wedge. Therefore, the third objective of this study
is to conduct laboratory studies to investigate the contaminant transport processes near a salt
wedge. We completed laboratory experiments to study solute transport patterns in both fresh and
95
salt water regions present above and below a saltwater wedge, respectively. The experimental
data showed that it took approximately 6 hours to fully flush a tracer slug released beneath the
wedge, whereas it took only about 30 minutes to flush a tracer slug released above the wedge.
This result indicates the rate of transport beneath a saltwater wedge will be an order of
magnitude slower than the transport rate above the wedge. Based on the experimental and
numerical results we estimated that the transport times are expected to differ approximately. The
data also indicated that the tracer slugs released below the wedge did not follow the idealized
circular flow path commonly depicted in textbooks. Due dispersion effects, the circular slug
released beneath the wedge attained an elongated shape as it approached the interface. In the
vicinity of the interface, a fraction of the diffused plume transported rather rapidly along the
mixing zone and moved towards the boundary as a narrow finger. The modeling data simulated
for a large scale system also showed similar solute transport patterns. The insights gained from
this study help us better understand the solute transport processes occurring near a saltwater
wedge.
6.2 Recommendations for Future Work
Several additional investigations could be conducted to further extend our understanding
of saltwater intrusion mechanisms in coastal aquifers. First, the findings of this study which
were primarily based on laboratory data and numerical simulation results could be verified in the
field study impacted by climate change effects (i.e., sea-level rise, rainfall variations). The
laboratory efforts focused on understanding an idealized sea level boundary condition that
ignored marine transgression effects, such as a tidal condition and seasonal variation. Future
studies are needed to study dynamic systems which can be influenced by transient boundary
forces that are impacted by tidal cycles and other long term changes such as inundation.
96
Secondly, the developed two-dimensional code has several computational limitations and should
be improved to reduce the simulation time. More robust advection tracking techniques and/or
more efficient numerical solvers can be used to speed-up code. The code could be combined
with an optimization method to analyze water resource management options in coastal aquifers.
Finally, future modeling efforts should focus on coupling density-coupled transport within a
multi-species reactive transport formulation to predict the exchange of nutrients and other
reactive contaminants between terrestrial to marine waters.
97
Appendix 1: Transient Variations in the Toe Position (XT) for 10%, 50% and 90% Isochlor of
saltwater wedge
Figure A1. The temporal variation of XT 10%, 50% and 90% of saltwater wedge.
The behavior of XT 10% and 90% isochlor of saltwater wedge follow similar pattern to 50%
isochlor. Therefore, XT 50% islcholor is useful to represent the temporal change of saltwater
wedge.
380
400
420
440
460
0 20,000 40,000
Ti
tle
Title
10% isochlor
50% isochlor
90% isochlor
98
Appendix 2 Comparison of Salt-wedge between Central-difference and Upstream-weighting for
density-weighting scheme
Figure A2. The comparison of saltwater wedge influenced by central difference and upstream
weighting for density-weighting scheme.
Based on ?x = 4 m and ?L = 1 m (Table 2.1), the grid Peclet number in the x-direction is Pex =
?x/ ?L = 4.0. Also, based on ?z = 0.4 m and ?T = 0.1 m (Table 2.1), the grid Peclet number in
the z-direction is Pez = ?z/ ?T = 4.0. These are relatively large values that may have led to
instabilities in the solutions if the central-difference weighting method was used to obtain the
solutions. However. we have not observed any instability issues in our solution. We also
completed some preliminary simulation to test various schemes and model convergence and
found that the wedge location was insensitive to the differencing scheme used.
0
10
20
30
0.00 200.00 400.00
z (
m)
x (m)
Central difference
upstream weighting
99
Appendix 3 Comparison of Salt-wedge between FDM and TVD scheme
Figure A3. The comparison of saltwater wedge influenced by dispersivity and by advection
schemes (Implicit FDM and TVD scheme).
The analytical toe position that computed using by equations of Appendix 1 is 625.5 m. Figure
shows that the length of XT is increased by decreasing dispersivity. For more accurate prediction
of toe position, TVD scheme is applied and compared to implicit finite difference scheme. In
results, the XT, 565 m, located close to analytical solution when dispersivity value was 0 for
TVD scheme.
0
10
20
30
0 200 400 600 800 1000
z (
m
)
x (m)
base case FDM
base case TVD
no dispersivity FDM
no dispersivity TVD
Analytical Toe
100
Appendix 4 Analytical Comparison of Salt-wedge Toe Positions in Unconfined and Confined
Aquifers
The basic concepts used for deriving analytical solutions for salt wedge locations using
the sharp-interface approximation are presented by Strack and others [Cheng et al., 2000;
Custodio and Bruggeman, 1987; Mantoglou, 2003; Strack, 1976]. Using the sharp-interface
approach, the toe position (XT) in a steady-state, unconfined flow system can predicted using the
following expression [Custodio and Bruggeman, 1987]:
g1850g3021 g3404uni0020g3436g1869
g1499
g1849g3397uni0020g1838g3440g3398g3496g3436
g1869g1499
g1849g3397uni0020g1838g3440
g2870
g3398g2012g1828g3042g2870g4666uni0031g3397g2012g4667g3436g1837g1849g3440uni0020 (1A)
Where, q*is is the depth-averaged flow through the boundary per unit width of the aquifer
[L2T-1], uni0042g2925 is the depth of aquifer bottom measured from the mean sea level [L]. W is the
uniform recharge rate [LT-1], uni03B4 is equal to g2977g3177g2879g2977g3164g2977
g3164
, where uni03C1g2916 is the density of fresh water [ML-3]
and uni03C1g2929 is the density of saltwater [ML-3]. Note the above equation is purely a function of
groundwater flows (or fluxes) and it does not depend on boundary head levels. Hence, the
location of the saltwater toe position will be insensitive to changes in head levels at the sea-side
boundary.
A similar expression for estimating the toe position in a confined flow system can be
derived from the expressions presented by Cheng et al [Cheng et al., 2000] as:
g1850g3021 g3404uni0020g3436g1869
g1499
g1849g3397uni0020g1838g3440g3398g3496g3436
g1869g1499
g1849g3397uni0020g1838g3440
g2870
g3398g2012g1828g2870g3436g1837g1849g3440uni0020 (2A)
uni0020
101
Where, B is the thickness of the confined aquifer [L]. Equation (2A) is similar to (1A) and the
only difference is the (uni0031g3397uni03B4g4667 term. For seawater, the value of the dimensionless parameter uni03B4 =
0.025, and hence the term (uni0031g3397uni03B4g4667uni0020uni007Euni0020uni0031. Therefore, when uni0042g2870uni004Buni007Euni0020uni0042g2925g2870uni004Buni0020(or when aquifer
thicknesses are matched approximately) the XT values predicted for confined and unconfined
flow systems would almost be the same.
102
Appendix 5 Derivation of Dimensionless Parameters for the Henry Problem and Analysis of
Henry-type Problem with Areal-recharge flux and zero regional flux
Notations
a = uni0020 g2901uni0020g2977g3116g2921g2914g4666g2977
g3177g2879uni0020g2977g3116g4667
b = D/Q
? = l/d, aspect latio.
c = concentration of salt, in mass per unit volume of solution [ML-3].
cs = concentration of salt in sea water [ML-3].
Q = net fresh-water discharge per unit length of beach [L2].
u = horizontal velocity [LT-1].
v = vertical velocity [LT-1].
? = density of solution [ML-3].
?o= density of freshwater [ML-3].
?s= density of saltwater [ML-3].
uni03C6?= dimensionless stream function
k = permeability of sand [L2].
d = thickness of aquifer [L].
l = length of aquifer [L].
P = fluid pore pressure [ML-1T-2].
g= acceleration of gravity [LT-2]
uni00B5 = viscosity of freshwater [ML-1T -1].
g1869g1499g3042 = flux at the coastline or left boundary [L2T-1].
N = recharge rate [LT-1].
103
1. Dimensionless parameter of Henry Problem with regional flux
Scalar form of Darcy?s law can be written as:
g1869g3051 g3404g1873g3404uni0020g3398g1863g2020uni0020g4666g2034g1842g2034g1876g4667 1)
g1869g3052 g3404g1874g3404uni0020g3398g1863g2020uni0020g4666g2034g1842g2034g1877g3397uni0020g2025g1859g4667 2)
Stream function definition in dimensionless form
g1873? g3404uni0020uni0020uni0020g2034g2030
?
g2034g1877?
g1874? g3404uni0020g3398g2034g2030
?
g2034g1876?
g1487g2870g2030? g3404uni0020 g2034g2034g1876?g4678g2034g2030
?
g2034g1876?g4679g3397
g2034
g2034g1877?g4678
g2034g2030?
g2034g1877?g4679g3404uni0020uni0020uni0020
g2034g1873?
g2034g1876? g3398uni0020
g2034g1874?
g2034g1876?
3)
From the dimensionless form
g1873? g3404g3048g3031g3044g1499 g1874? g3404g3049g3031g3044g1499 g1876? g3404g3051g3031,uni0020g1877? g3404g3052g3031, g1855? g3404 g3004g3004
g3294
, g2025? g3404 g3096g2879g3096g3290g3096
g3294g2879g3096g3290
Dimensionless derivation
g2034g1876
g2034g1876? g3404g1856
g2034g1877
g2034g1877? g3404g1856
g2034g1873?
g2034g1877? g3404uni0020
g2034uni0020
g2034g1877?g3436
g1873g1856
g1869g1499g3440g3404uni0020
g1856
g1869g1499
g2034g1873
g2034g1877
g2034g1877
g2034g1877? g3404uni0020
g1856g2870
g1869g1499
g2034g1873
g2034g1877
g2034g1874?
g2034g1876? g3404uni0020
g2034uni0020
g2034g1876?g3436
g1874g1856
g1869g1499g3440g3404uni0020
g1856
g1869g1499
g2034g1874
g2034g1876
g2034g1876
g2034g1876? g3404uni0020
g1856g2870
g1869g1499
g2034g1874
g2034g1876
g2034g1829?
g2034g1876? g3404uni0020
g2034uni0020
g2034g1876?g3436
g1829
g1829g3046g3440g3404uni0020
uni0031
g1829g3046
g2034g1829
g2034g1876
g2034g1876
g2034g1876? g3404uni0020
g1856
g1829g3046
g2034g1829
g2034g1876
104
g2034g2025?
g2034g1876? g3404uni0020
g2034uni0020
g2034g1876?g3436
g2025g3398g2025g3042
g2025g3046g3398g2025g3042g3440g3404uni0020
uni0031
g2025g3046g3398g2025g3042
g2034g2025
g2034g1876
g2034g1876
g2034g1876? g3404
g1856
g2025g3046g3398g2025g3042
g2034g2025
g2034g1876
Using these dimensionless form
Equation (3) is converted into dimensioned form
g1487g2870g2030? g3404uni0020uni0020uni0020g2034g1873
?
g2034g1876? g3398uni0020
g2034g1874?
g2034g1876? g3404
g1856g2870
g1869g1499
g2034g1873
g2034g1877g3398
g1856g2870
g1869g1499
g2034g1874
g2034g1876uni0020uni0020g3404
g1856g2870
g1869g1499g4666uni0020
g2034g1873
g2034g1877g3398
g2034g1874
g2034g1876g4667uni0020
g3105g3048
g3105g3052,
g3105g3049
g3105g3051 terms are transformed in terms of pressure,
g2034g1873
g2034g1877g3404uni0020uni0020
g2034uni0020
g2034g1877g3428g3398
g1863
g2020uni0020g3436
g2034g1842
g2034g1876g3440g3432g3404uni0020
g1863
g2020
g2034g1842
g2034g1876g2034g1877
5)
g2034g1874
g2034g1876g3404
g2034
g2034g1876g3428g3398
g1863
g2020uni0020g3436
g2034g1842
g2034g1877g3397uni0020g2025g1859g3440g3432g3404uni0020
g1863
g2020
g2034g1842
g2034g1876g2034g1877g3398
g1863g1859
g2020
g2034g2025
g2034g1876
6)
When Subtract between equation (6) from (5) is
g2034g1873
g2034g1877g3398
g2034g1874
g2034g1876g3404uni0020g1859
g2034g2025
g2034g1876
Therefore
g1487g2870g2030? g3404g1856
g2870
g1869g1499g3436uni0020
g2034g1873
g2034g1877g3398
g2034g1874
g2034g1876g3440g3404uni0020
g1856g2870
g1869g1499g3436uni0020
g1863g1859
g2020
g2034g2025
g2034g1876g3440g3404
g1856g2870
g1869g1499uni0020g4678
g2020g3033g1837
g2020g2025g3042
g2025g3046g3398g2025g3042
g1856
g2034g2025?
g2034g1876?g4679g3404uni0020
g1837g1856g4666g2025g3046g3398g2025g3042g4667
g1869g1499
g2034g1855?
g2034g1876?
using
g1837g3404uni0020g1863g2025g3042g1859uni0020g2020
g3033
g2034g2025
g2034g1876g3404uni0020
g2025g3046g3398g2025g3042
g1856
g2034g2025?
g2034g1876?
g2034g2025?
g2034g1876? g3404
g2034g1855?
g2034g1876?
105
2. Dimensionless parameter of Henry Problem with areal-recharge flux
g1869g1499g3404uni0020g1869g1499g3042g3398g1840g1876g3404g1858g4666g1876g4667
g2034g1869g1499
g2034g1876 g3404g3398g1840
From the dimensionless form
g1873? g3404 g3048g3031g3044g1499
g3290g2879g3015g3051
g1874? g3404 g3049g3031g3044g1499
g3290g2879g3015g3051
g1876? g3404g3051g3031,uni0020g1877? g3404g3052g3031, g1855? g3404 g3004g3004
g3294
, g2025? g3404 g3096g2879g3096g3290g3096
g3294g2879g3096g3290
Dimensionless derivation
g2034g1876
g2034g1876? g3404g1856
g2034g1877
g2034g1877? g3404g1856
g2034g1873?
g2034g1877? g3404uni0020
g2034uni0020
g2034g1877?g3436
g1873g1856
g1869g1499g3440g3404uni0020
g1856
g1869g1499
g2034g1873
g2034g1877
g2034g1877
g2034g1877? g3404uni0020
g1856g2870
g1869g1499
g2034g1873
g2034g1877
g2034g1874?
g2034g1876? g3404uni0020
g2034uni0020
g2034g1876?g3436
g1874g1856
g1869g1499g3440g3404uni0020g1856g4678g3398
uni0020g1874
g1869g1499g2870
g2034g1869g1499
g2034g1876
g2034g1876
g2034g1876? g3397
uni0031
g1869g1499
g2034g1874
g2034g1876
g2034g1876
g2034g1876?g4679g3404uni0020g1856g4678g3398
uni0020g1874
g1869g1499g2870g4666g3398g1840g4667g1856g3397
uni0031
g1869g1499
g2034g1874
g2034g1876g1856g4679uni0020uni0020
g3404 g1856
g2870
g1869g1499g2870g4666g1840g1874g3397g1869
g1499g2034g1874
g2034g1876g4667
g2034g1829?
g2034g1876? g3404uni0020
g2034uni0020
g2034g1876?g3436
g1829
g1829g3046g3440g3404uni0020
uni0031
g1829g3046
g2034g1829
g2034g1876
g2034g1876
g2034g1876? g3404uni0020
g1856
g1829g3046
g2034g1829
g2034g1876
g2034g2025?
g2034g1876? g3404uni0020
g2034uni0020
g2034g1876?g3436
g2025g3398g2025g3042
g2025g3046g3398g2025g3042g3440g3404uni0020
uni0031
g2025g3046g3398g2025g3042
g2034g2025
g2034g1876
g2034g1876
g2034g1876? g3404
g1856
g2025g3046g3398g2025g3042
g2034g2025
g2034g1876
Equation (3) is converted into dimensioned form
g1487g2870g2030? g3404uni0020uni0020uni0020g2034g1873
?
g2034g1876? g3398uni0020
g2034g1874?
g2034g1876? g3404
g1856g2870
g1869g1499
g2034g1873
g2034g1877g3398
g1856g2870
g1869g1499g2870g3436g1840g1874g3397g1869
g1499g2034g1874
g2034g1876g3440uni0020g3404
g1856g2870
g1869g1499g3436uni0020
g2034g1873
g2034g1877g3398
g2034g1874
g2034g1876g3440g3398
g1856g2870
g1869g1499g2870g1840g1874
Due to the presence of the last term, which has dimensions, we were unable to cast this
problem in a fully non-dimensional form. Further analysis is needed to study this problem.
106
Appendix 6 Development of Two-dimensional Numerical Code for Simulating Saltwater
Transport in Groundwater Systems
A7.1.Introduction
In the published literature, several numerical codes are available for simulating density-
coupled flow problems. The most commonly used codes are: SEAWAT [Guo and Langevin,
2002] which is based on a finite-difference formulation and FEFLOW [Diersch, 2002] and
SUTRA [Voss and Provost, 2002] which are based on a finite-element formulation. The source
code for SEAWAT uses a modified version of the USGS groundwater model MODFLOW to
simulate equivalent freshwater head as the principal dependent variable. In addition, it uses the
popular solution transport code MT3DMS to simulate salt and other solute transport processes.
Despite the availability of these public-domain/commercial codes, many researchers have
pointed out the need for more robust density coupled flow and transport models to study more
complex problems (e.g. metal leachate migration, transport of zero-valent iron, and redox
controlled dense metal plume transport, to name a few) [Goswami et al., 2012; Oostrom et al.,
1992; Schincariol and Schwartz, 1990]. Goswami et al. [2012], for example, has shown that
there are still several inaccuracies in simulating certain variable-density systems and hence more
research is needed to formulate computationally efficient codes to solve these systems.
Henry problem and Elder problem are commonly used benchmark problems for verifying
density-dependent groundwater codes. Simpson and Clement [2003] studied these two
benchmark problems using a coupled and uncoupled approach. The test showed that the Elder
problem is better than the Henry problem for testing variable-density codes.
107
The Henry problem, however, is the only variable density problem that has a semi-
analytical solution [Henry, 1964]. Henry used an analytical framework to define the location
and shape of the interface under a constant seaward freshwater flow moving toward an oceanic
boundary. Since Henry?s initial effort, several other analytical and numerical solutions have
been developed for the original Henry problem by Pinder and Cooper [1970], Lee and Cheng
[1974], Segol et al. [1975], Frind [1982], Huyakorn et al.[1987], and Voss and Souza [1987].
These studies have tested and modified some of the original parameters to fit Henry?s solution;
Croucher and O?sullivan [1995], reviewed all these studies and concluded that Henry?s original
analytical solution might have a minor error. Segol [1994] reinvestigated the errors associated
with the Henry solution and suggested some corrections. Simpson and Clement [2004]
developed an improved Henry problem. Goswami and Clement [2007] completed experimental
studies to develop a new benchmark dataset and used the coupled-uncoupled analysis to test the
worthiness of the new experimental problem for testing density-coupling effects. The results of
the analysis showed that the proposed experimental benchmark problem is a more robust
alternative to the traditional Henry problem.
The objective of the present work is to develop a computationally efficient finite-
difference algorithm that can solve a wide variety of two-dimensional, density-dependent flow
problems in saturated groundwater systems. The code was then tested using several benchmark
problems discussed in the previous section. The numerical code is built in Microsoft Excel
based Visual Basic framework. The finite-difference MODFLOW-family variable density flow
code SEAWAT [Guo and Langevin, 2002] was also used to simulate the benchmark datasets.
The model uses the pressure-based, finite-difference approach for formulating the governing
flow equations. The pressure based approach is relatively simple and direct compared to other
108
head-based formulations. For solving the solute transport, the advection process was solved
using the third-order TVD scheme and the dispersion/source terms are solved using implicit
finite-difference scheme.
A7.2 Governing Equation
The governing equations for density-dependent saturated flow through a two-dimensional
homogeneous and isotropic porous medium can be expressed as:
uni2202
uni2202uni0078g3428
uni004B
uni0067uni0020uni0020g3436
uni2202uni0050
uni2202uni0078uni0020g3440g3432g3397uni0020
uni2202
uni2202uni007Ag3428
uni004B
uni0067uni0020uni0020g3436
uni2202uni0050
uni2202uni007Ag3398uni03C1uni0067uni0020uni0020g3440g3432g3404uni0020uni03C1uni0053g2926
uni2202uni0050
uni2202uni0074g3397uni03B8uni0045
uni2202uni0043
uni2202uni0074 (A7.1)
where, Sp [M-1LT2]is the specific storage in terms of pressure g4672g3404 g2903g3177g2977
g3164g2917
g4673 , P [ML-1T-2]is
pressure, uni03C1 [ML-3]is the fluid density, uni03B8 is the water content, K [LT-1]is the hydraulic
conductivity, t [T] is time and x and z are the Cartesian coordinates, E is the density
concentration slope. The first term in the right-hand side of equation (A7.1) represents the rate
of fluid mass accumulation due to ground-water storage effects (for example, due to the
compressibility of the bulk porous material and fluid compressibility). The second term in the
right-hand side of the equation represents the rate of fluid mass accumulation due to the change
of solute concentration.
The two-dimensional transport is modeled using the advection-dispersion equation,
which is solved using the operator split strategy. The governing equation and the operator split
version of the transport equations can be written as:
uni2202uni0043
uni2202uni0074 g3404g3398uni0056g2934
uni2202uni0043
uni2202uni0078g3398uni0056g2935
uni2202uni0043
uni2202uni0079g3397
uni2202
uni2202uni0078g4666uni0044g2934
uni2202uni0043
uni2202uni0078g4667uni0020g3397uni0020
uni2202
uni2202uni0079g4666uni0044g2935
uni2202uni0043
uni2202uni0079g4667uni0020
(A7.2)
109
uni2202uni0043uni0020
uni2202uni0074 g3404g3398uni0056g2934
uni2202uni0043uni0020
uni2202uni0078g3398uni0056g2935
uni2202uni0043uni0020
uni2202uni0079
(A7.2a)
uni2202uni0043g3364
uni2202uni0074 g3404
uni2202
uni2202uni0078g4666uni0044g2934
uni2202uni0043g3364
uni2202uni0078g4667uni0020g3397uni0020
uni2202
uni2202uni0079g4666uni0044g2935
uni2202uni0043g3364
uni2202uni0079g4667uni0020
(A75.2b)
where, uni0056g2934uni0020uni0061uni006Euni0064uni0020uni0056g2935g4670uni004Cuni002Funi0054g4671 are velocity, C [ML-3] is the solute concentration, uni0044g2934uni0020uni0061uni006Euni0064uni0020uni0044g2935 [L2/T] are
are the hydrodynamic dispersion coefficient.
uni0020
Iuni006Euni0020uni0074heuni0020operuni0061uni0074iouni006Euni0020spliuni0074uni0020meuni0074houni0064,uni0020uni0074heuni0020uni0061uni0064vecuni0074iouni006Euni0020puni0061runi0074uni0020ofuni0020uni0074heuni0020equuni0061uni0074iouni006Euni0020isuni0020solveuni0064uni0020uni0061suni0020showuni006Euni0020iuni006Euni0020
equuni0061uni0074iouni006Euni0020g4666A7.2uni0061g4667uni0020usiuni006Euni0067uni0020uni0054uni0056uni0044uni0020puni0061cg141uni0061uni0067e.uni0020uni0020uni0053ug132sequeuni006Euni0074luni0079,uni0020uni0074heuni0020couni006Eceuni006Euni0074runi0061uni0074iouni006Euni0020og132uni0074uni0061iuni006Eeuni0064uni0020fromuni0020
uni0061uni0064vecuni0074iouni006Euni0020g4666uni0043g3364g4667uni0020isuni0020useuni0064uni0020uni0074ouni0020solveuni0020uni0074heuni0020uni0064ispersiouni006Euni0020puni0061runi0074uni0020ofuni0020uni0074heuni0020equuni0061uni0074iouni006Euni0020uni0061suni0020showuni006Euni0020iuni006Euni0020equuni0061uni0074iouni006Euni0020
g4666A7.2g132g4667uni0020usiuni006Euni0067uni0020Impliciuni0074uni0020g9uni0044g16.uni0020uni0020To calculate velocity at nodes, the specific discharge is divided by
porosity. The general expression for specific discharge (Darcy flux) can be written as: uni0020
qg2934g3404g3398g141g2934uni00B5uni0020uni0020uni2202uni0050uni2202uni0078 (
qg2936g3404g3398g141g2936uni00B5uni0020g4666uni2202uni0050uni2202uni007Ag3397uni03C1uni0067uni0020uni0020g4667 (A7.3)
where qg2934 and qg2934uni0020are the individual component of specific discharge, P is pressure [ML-1T-2], uni03C1 is
the fluid density [ML-3], g is acceleration due to gravity [LT-2].
uni0045quuni0061uni0074iouni006Euni0020g4666A7.uni0031g4667uni0020uni0061uni006Euni0064uni0020g4666A7.2g4667uni0020uni0061reuni0020coupleuni0064uni0020uni0074hrouuni0067huni0020equuni0061uni0074iouni006Euni0020g4666A7.g886g4667.uni0020
uni03C1g3404uni0020uni03C1g2916g3397uni0045uni0043 (A7.4)
110
Based on SI unit, the value of typical seawater concentration is about 35 kg/m3 and the
freshwater density is set to 1,000 kg/m3 and saltwater density is set to 1,025 kg/m3 with the
density slope of 0.714 (g3404 g2914g2915g2924g2929g2919g2930g2935uni0020g2914g2919g2916g2916g2915g2928g2915g2924g2913g2915,uni0020uni0020uni0020g2870g2873uni0020g2921g2917uni002Fg2923g3119g2903g2911g2922g2930g2933g2911g2930g2915g2928uni0020g2913g2925g2924g2913g2915g2924g2930g2928g2911g2930g2919g2925g2924,uni0020uni0020uni0020g2871g2873uni0020g2921g2917uni002Fg2923g3119).
A7.3 Model Development
Flow Equation
uni2202
uni2202uni0078g3428
uni004B
uni0067uni0020uni0020g3436
uni2202uni0050
uni2202uni0078uni0020g3440g3432g3397uni0020
uni2202
uni2202uni007Ag3428
uni004B
uni0067uni0020uni0020g3436
uni2202uni0050
uni2202uni007Ag3398uni03C1uni0067uni0020uni0020g3440g3432 (A7.5a)
g3406 uni004Buni0067uni2206uni0078g4680uni0020uni0020g4678uni0050g2919g2878g2869,g2920
g2923,g2924g3398uni0050
g2919,g2920
g2923,g2924
uni2206uni0078 g4679g3398uni0020g4678
uni0050g2919,g2920g2923,g2924g3398uni0050g2919g2879g2869,g2920g2923,g2924
uni2206uni0078 g4679g4681
g3397uni0020 uni0020uni004Buni0067uni2206uni0079g4680uni0020g4678uni0050g2919,g2920g2878g2869
g2923,g2924g3398uni0050
g2919,g2920
g2923,g2924
uni2206uni0079 g4679g3398uni0020g4678
uni0050g2919,g2920g2923,g2924g3398uni0050g2919,g2920g2879g2869g2923,g2924
uni2206uni0079 g4679g4681
g3398uni0020 uni0020uni004Buni0020uni0067uni2206uni0079g3428uni03C1g2919g2878g2869
g2870,g2920
g2923g2879g2869,g2924uni0067g3398uni03C1
g2919g2879g2869g2870,g2920
g2923g2879g2869,g2924uni0067g3432
where n denotes the nth discrete time level, tn, when the solution is known, and m is the Picard
iteration level, uni0050g2919,g2920 is the pressure head at node ij. The subscripts i+1/2, i-1/2 refer to the value
between two neighboring cells. The density between two cells expressed as average value. The
nodal order matches the numbering order of VBA language for all formulation in this study.
The current and previous iteration levels are denoted as m and m-1, respectively.
uni03C1uni0053g2926uni2202uni0050uni2202uni0074g3397uni03B8uni0045uni2202uni0043uni2202uni0074 g3406uni03C1g2919,g2920uni0053g2926uni0050g2919,g2920
g2923,g2924g3398uni0050
g2919,g2920
uni0020g2924g2879g2869
uni2206uni0074 g3397uni03B8uni0045
uni0043g2919,g2920g2923g2879g2869,g2924g3398uni0043g2919,g2920g2924g2879g2869
uni2206uni0074
(A7.5b)
where, uni2206uni0074g3404uni0074g2924g3398uni0074g2924g2879g2869is the time step, which in general is set to maintain acceptably small
temporal variation in pressure head.
uni0054heuni0020fiuni006Eiuni0074euni0020uni0064iffereuni006Eceuni0020formuni0020reuni0064uceuni0064uni0020uni0074ouni0020
111
uni0061uni0031uni0020uni0050g2919g2879g2869,g2920g2923,g2924g3397uni0020g132uni0031uni0020uni0050g2919,g2920g2879g2869g2923,g2924g3397cuni0031uni0020uni0050g2919,g2920g2923,g2924g3397uni0064uni0031uni0020uni0050g2919,g2920g2878g2869g2923,g2924g3397euni0031uni0020uni0050g2919g2878g2869,g2920g2923,g2924g3404uni0052uni0048uni0053uni0020 g4666A7.g888g4667uni0020
where,uni0020uni0020
uni0061uni0031g3404uni0020g3398 uni004Bg2934uni0067uni2206uni0078g2870uni0020
g1322g3404g3398 uni004Bg2936uni0067uni2206uni007Ag2870uni0020
cuni0031g34042g3436 uni004Bg2936uni0067uni2206uni0078g2870g3397 uni004Bg2936uni0067uni2206uni007Ag2870g3440g3397uni0020uni03C1g2919,g2920uni0053g2926uni2206uni0074uni0020
uni0064uni0031g3404uni0020g3398uni0020 uni004Bg2936uni0067uni2206uni007Ag2870uni0020
euni0031g3404uni0020g3398 uni004Bg2934uni0067uni2206uni0078g2870uni0020
uni0052uni0048uni0053g3404uni0020 uni004Bg29362uni2206uni007Ag3427uni03C1g2919g2878g2869,g2920g2923g2879g2869,g2924uni0020g3398uni03C1g2919g2879g2869,g2920g2923g2879g2869,g2924g3431g3397uni03C1g2919,g2920uni0020uni0053g2926uni2206uni0074uni0050g2919,g2920uni0020g2924g2879g2869g3398uni03B8uni0045uni2206uni0074g4666uni0043g2919,g2920g2923g2879g2869,g2924g3398uni0043g2919,g2920g2924g2879g2869g4667uni0020
uni0020
Eq. (A7.6) applies to all interior nodes: at boundary nodes this equation is modified to
reflect the appropriate boundary conditions. This would lead to a system of linear equations of
the form:
uni0041uni00B7uni0050g3553g3404g132uni0020 g4666A7.g889g4667uni0020
Where, A is a banded square matrix, uni0050g3553 is the vector of unknown pressure at the current time level
n and iteration level, m, and b is the forcing vector. The set of linear algebraic equations are
solved using a banded version of a Gauss?elimination solver.
Transport Equation
The linear formulation for TVuni0044uni0020schemeuni0020foruni0020uni0061uni0064vecuni0074iouni006Euni0020puni0061cg141uni0061uni0067euni0020isuni0020uni0064escrig132euni0064uni0020g132elow.uni0020uni0020uni0045quuni0061uni0074iouni006Euni0020
g4666A7.g889uni0061g4667uni0020isuni0020spuni0061uni0074iuni0061luni0020formuluni0061uni0074iouni006Euni0020uni0061uni006Euni0064uni0020equuni0061uni0074iouni006Euni0020g4666A7.g889g132g4667uni0020isuni0020uni0074emporuni0061luni0020formuluni0061uni0074iouni006E.uni0020uni0020
112
g3398uni0056g2934uni2202uni0043uni2202uni0078g3398uni0056g2935uni2202uni0043uni2202uni0079
g3406uni0020g3398uni0056g2919,g2920g2878g2869uni002Fg2870uni0043g2919,g2920g2878g2869uni002Fg2870
g2924 g3398uni0056g2919,g2920g2879g2869uni002Fg2870uni0043
g2919,g2920g2879g2869uni002Fg2870
g2924
uni0394uni0079 g3398
uni0056g2919g2878g2869uni002Fg2870,g2920uni0043g2919g2878g2869uni002Fg2870,g2920g2924 g3398uni0056g2919g2879g2869uni002Fg2870,g2920uni0043g2919g2879g2869uni002Fg2870,g2920g2924
uni0394uni0078
g4666A7.g889uni0061g4667uni0020
uni2202uni0043uni0020
uni2202uni0074 g3406uni0020
uni0043g2919,g2920g2924g2878g2869g3398uni0043g2919,g2920g2924
uni0394uni0074 uni0020
g4666A7.g889g132g4667uni0020
The TVD scheme adopts a universal flux limiting procedure to minimize numerical oscillations
which may occur if sharp concentration fronts are involved. This scheme is mass conservative,
without excessive numerical dispersion, and essentially oscillation-free. Appendix 8 illustrates
the detailed TVD formulation procedure and the application of universal limiter.
The implicit FDM formulation is described below for dispersion package.uni0020uni0045quuni0061uni0074iouni006Euni0020g4666A7.g890uni0061g4667uni0020isuni0020
spuni0061uni0074iuni0061luni0020formuluni0061uni0074iouni006Euni0020uni0061uni006Euni0064uni0020equuni0061uni0074iouni006Euni0020g4666A7.g890g132g4667uni0020isuni0020uni0074emporuni0061luni0020formuluni0061uni0074iouni006E.uni0020uni0020
g3398 uni2202uni2202uni0078g4666uni0044g2934uni2202uni0043uni2202uni0078g4667uni0020g3397uni0020uni2202uni2202uni0079g4666uni0044g2935uni2202uni0043uni2202uni0079g4667uni0020
g3406 uni0031uni0020uni2206uni0078g3429uni0020uni0044g2919g2878g2869
g2870,g2920uni0020
g2923g2879g2869,g2924uni0020g4678uni0043g2919g2878g2869,g2920
g2923,g2924g3398uni0043
g2919,g2920
g2923,g2924
uni2206uni0078 g4679g3398uni0044g2919g2879g2869g2870,g2920uni0020
g2923g2879g2869,g2924uni0020uni0020g4678uni0043g2919,g2920
g2923,g2924g3398uni0043
g2919g2879g2869,uni006Ag2923,g2924
uni2206uni0078 g4679g3433
g3397uni0020 uni0031uni0020uni2206uni0078g3429uni0020uni0044g2919,g2920g2878g2869
g2870uni0020
g2923g2879g2869,g2924uni0020g4678uni0043g2919,g2920
g2923,g2924g3398uni0043
g2919,g2920g2879g2869
g2923,g2924
uni2206uni0078 g4679g3398uni0044g2919,g2920g2879g2869g2870uni0020
g2923g2879g2869,g2924uni0020uni0020g4678uni0043g2919,g2920
g2923,g2924g3398uni0043
g2919,g2920g2879g2869
g2923,g2924
uni2206uni0078 g4679g3433uni0020
g4666A7.g890uni0061g4667uni0020
uni2202uni0043uni0020
uni2202uni0074 g3406uni0020
uni0043g2919,g2920g2924g2878g2869g3398uni0043g2919,g2920g2924
uni0394uni0074
g4666A7.g890g132g4667uni0020
Note that the average dispersion tensor used in the spatial terms is the value that was updated by
the previous iteration step, m-1. The value is calculated explicitly so that the equation has only
concentration as an independent variable to solve in linear matrix.
The general expression for finite difference form is
uni0061g2870uni0043g2919g2879g2869,g2920g2923,g2924g3397uni0020g132g2870uni0020uni0043g2919,g2920g2879g2869g2923,g2924g3397cg2870uni0020uni0043g2919,g2920g2923,g2924g3397uni0064g2870uni0020uni0043g2919,g2920g2878g2869g2923,g2924g3397eg2870uni0020uni0043g2919g2878g2869,g2920g2923,g2924g3404uni0052uni0048uni0053 (A7.9)
113
where,
uni0061g2870g3404uni0020uni0020uni0044g2934g2919g2879g2869
g2870,g2920uni0020
g2923g2879g2869,g2924 uni0031
uni0020uni2206uni0078g2870
g132g2870g3404uni0044g2935g2919,g2920g2879g2869
g2870uni0020
g2923g2879g2869,g2924 uni0031
uni0020uni2206uni0079g2870uni0020
cg2870g3404g3398 uni0031uni0020uni2206uni0078g2870uni0020g4678uni0044g2934g2919g2879g2869
g2870,g2920uni0020
g2923g2879g2869,g2924g3397uni0044
g2934g2919g2879g2869
g2870,g2920uni0020
g2923g2879g2869,g2924g4679g3398 uni0031
uni0020uni2206uni0079g2870g4678uni0044g2935g2919,g2920g2879g2869g2870uni0020
g2923g2879g2869,g2924g3397uni0044
g2935g2919,g2920g2878g2869
g2870uni0020
g2923g2879g2869,g2924g4679g3398uni0031uni0020
uni0064g2870g3404uni0020uni0044g2935g2919,g2920g2878g2869
g2870uni0020
g2923g2879g2869,g2924 uni0031
uni0020uni2206uni0079g2870uni0020
eg2870g3404uni0020uni0044g2934g2919g2878g2869
g2870,g2920uni0020
g2923g2879g2869,g2924 uni0031
uni0020uni2206uni0078g2870uni0020
uni0052uni0048uni0053g3404uni0020uni0043
g1499g2919,g2920g2923,g2924
uni2206uni0074 uni0020uni0020
Eq. (A7.9) applies to all interior nodes: at boundary nodes this equation is modified to reflect the
appropriate boundary conditions. The resulting system of linear equations isuni0020uni0061lsouni0020solveuni0064uni0020g132uni0079uni0020the
banded Gauss-elimination solver.
Boundary Conditions
Test problems-1, 2 and 4 use Dirichlet boundary condition for freshwater boundary. Test
problem-3 uses Numann boundary condition for freshwater inflow. Also, the tank bottom is set
to no-flow boundary for all test problems. Saltwater intrusion modifies the concentration at the
boundary. The detailed modification is expressed in the model setup for each test problem.
Time Step Determination
Since advection is solved by explicitly, Courant condition is required to limit transport
time step, ?t. Courant condition is applied to all simulations to maintain sufficiently small time
steps. The length of time step is calculated by Courant condition to satisfy the stability
114
constraints and accuracy requirements of the transport equation. In this study, saltwater
simulation generally set the Courant number to 0.5. The Elder problem decreases the Courant
number to 0.1.
115
A7.3 Results
Test Problem 1: Sinking plume model
The first test problem is an example of dense solute transport. Goswami et al. [2012]
provided an experimental and numerical solution conducted in a two-dimensional flow tank
using three different numerical techniques available in SEAWAT/MT3DMS. The experiments
were designed to represent a sinking groundwater plume and a rising groundwater plume. This
problem is also helpful when applied to practical environmental issues such as contaminant
transport and tsunami issues in coastal aquifers. The domain and boundary conditions describe a
confined aquifer; constant freshwater pressure heads are maintained on the left and the right
boundary of the main, respectively (see Figure A7.1). The two-dimensional modeling domain is
225 mm long and 180 mm high. The longitudinal dispersivity value was set to 0.1 mm and
transverse dispersivity is set to 0.01 mm. Three stress periods were employed to simulate the
following three experimental stages: (1) steady-state conditions prior to injection, (2) injection of
freshwater/saltwater slug, and (3) migration of slug through the tank. The time step size was set
to an initial minimum value of 0.1 s and a maximum value of 1 s. The Courant number constraint
was set to 0.5. Values for this model are presented in Table A7.1.
The simulated results of third state (migration stage) were compared to SEAWAT results
and shown in Figure A7.2. The simulation shows good agreement with SEAWAT (See Figure
A7.4). The density-dependent sinking study can be also applied to a different composition in
fluid from typical seawater or the salt concentration. For example, Kanel et al. [2008] simulated
the fate and transport of Zero-valent iron nanoparticles (INP) and stabilized zero-valent iron
nanoparticles using poly acrylic acid (S-INP) in porous media under saturated, steady-state flow
conditions.
116
Figure A7.1 Conceptual description of the numerical model used for simulating the sinking-
plume experiments (Goswami et al, 2012).
Table A7.1 Summary of numerical model parameters
Property Symbol values
Horizontal length L 23 cm
Vertical Length hs 18 cm
Hydraulic conductivity K 1.16 cm/s
Porosity n 0.38
Specific storage, Ss 0.000001cm-1
Longitudinal dispersivity ?L 0.01 cm
Transverse dispersivity ?T 0.001 cm
Saltwater concentration Cs 0.035 g/cm3
Saltwater density ?s 1.025 g/cm3
Freshwater density ?f 1.000 g/cm3
Density slope E 0.714
117
(a) VDFT code (b) SEAWAT
Figure A7.2 Comparison of VDFT code simulation results against to SEAWAT for the sinking
plume experiment at various times at 0 min (left plume) , 2 min (middle) and 4 min (right),
Red color > 10% and Blue color 1 ? 10 % .
0 5 10 15 200
2
4
6
8
10
12
14
16
18
0 5 10 15 200
2
4
6
8
10
12
14
16
18
118
Test Problem 2: Saltwater Intrusion Model
A conceptual model based on Goswami and Clement [2007] was simulated. The
saltwater behavior in the saltwater intrusion problem is relatively stable. Goswami and Clement
[2007] proposed the experimental data set that can be used as a benchmark problem. The
domain and boundary conditions describe a confined aquifer; constant pressure heads are
maintained on the left and the right boundary of the main, respectively (see Figure A7.3). This
Henry-type problem is useful to understand the role of several hydrodynamic and geologic
factors affecting the flow and transport process in the simplified flow system. In order to
construct the numerical model, the origin of an x-z coordinate system was set on the lower left-
hand corner of the plane. A uniform grid of size 1.0 cm (54 cells in the x-direction and 26 cells
in the z-direction) was used to discretize the 53.0 cm ? 26.0 cm dimensional model. The
saltwater head is imposed on left boundary where the head is set to 25.5 cm for the entire
simulation time. The value of specific storage (Ss) was set to 10-6 cm-1 for all nodes. All layers
are set to confined layer and bottom of the model is set to no-flow condition. Freshwater flux is
allowed to penetrate from the right constant-head boundary. The value of 0.05 cm and 0.005 cm
were selected for the longitudinal and transversal dispersivity coefficient, respectively. Values
for this model are presented in Table A7.2.
Three steady state conditions were selected and simulated using VDFT code and
compared to SEAWAT results. The simulation shows good agreement with SEAWAT for the
steady-state condition (See Figure A7.4). Transient simulations for an intruding condition and a
receding condition also show good agreements (See Figure A7.5).
119
Figure A7.3 Computation domain and boundary conditions used in the numerical model
(Goswami and Clement, 2007)
Table A7.2 Summary of numerical model parameters
Property Symbol values
Horizontal length L 53 cm
Saltwater elevation hs 26 cm
Hydraulic conductivity K 1.215 cm/s
Porosity n 0.385
Specific storage, Ss 0.000001cm-1
Longitudinal dispersivity ?L 0.01 cm
Transverse dispersivity ?T 0.001 cm
Saltwater concentration Cs 0.0371 g/cm3
Saltwater density ?s 1.026 g/cm3
Freshwater density ?f 1.000 g/cm3
Density Slope E 0.714
120
Figure A7.4 Comparison of VDFT code simulation results against to SEAWAT for steady state
in saltwater problem.
0
5
10
15
20
25
0 10 20 30 40 50
z (
cm
)
x (cm)
VDFT SS1
VDFT SS2
VDFT SS3
SEAWAT
121
(a)
(b)
Figure A7.5 Comparison of VDFT code simulation results against to SEAWAT for (a) intruding
and (b) receding condition in saltwater problem.
0
5
10
15
20
25
0 10 20 30 40 50
z (
cm
)
x (cm)
VDFT 5 minutes
VDFT 15 minutes
VDFT 55 minutes
SEAWAT
0
5
10
15
20
25
0 10 20 30 40 50
z (
cm
)
x (cm)
VDFT10 minutes
VDFT 15 minutes
VDFT 25 minutes
SEAWAT
122
Test Problem 3: Flux-type Henry Problem
Conceptual models based on original Henry problem [1964] and modified Henry
problem [Simpson and Clement, 2004] were simulated. The domain and boundary conditions for
the Henry and Modified Henry problem describe a confined aquifer. A constant pressure head is
maintained on the right boundary and constant flux is maintained on the left boundary of the
main (see Figure A7.6). In order to construct a numerical model, the origin of an x-z coordinate
system was set on the lower left-hand corner of the plane. A uniform grid of size 5.0 cm (41
cells in the x-direction and 20 cells in the z-direction) was used to discretize the 200 cm ? 100
cm dimensional model. The saltwater head is imposed on the right boundary where the head is
set to 100 cm for the entire simulation time. The value of specific storage (Ss) was set to 10-6 cm-
1 for all nodes. All layers are set to confined layer and bottom of the model is set to no-flow
condition. The longitudinal and transversal dispersivity coefficients are set to zero. Values for
this model are presented in Table A7.3.
The simulated results for (a) original Henry problem and (b) modified Henry problem
were compared to SEAWAT results and shown in Figure A7.7. This difference indicates that the
relative importance of the density effects is greater for the modified case where the influence of
the boundary forcing is reduced [Simpson and Clement, 2004]. In this figure, bigger discrepancy
is shown at the bottom of the model when the modified Henry problem was simulated. The
VDFT code is still under development to predict more accurate density effects.
123
Figure A7. 6 Boundary conditions applied for solving Henry saltwater intrusion Problem (Henry,
1964; Simpson and Clement, 2003).
Table A7.3 Summary of numerical model parameters
Property Symbol values
Horizontal length L 200 cm
Saltwater elevation hs 100 m
Freshwater darcy flow q 6.6 ? 10-3 cm/s
Hydraulic conductivity K 1.0 cm/s
Porosity n 0.35
Specific storage, Ss 0 cm-1
Longitudinal dispersivity ?L 0 cm
Transverse dispersivity ?T 0 cm
Saltwater concentration Cs 35 kg/m3
Saltwater density ?s 1,025 kg/m3
Freshwater density ?f 1,000 kg/m3
Density slope E 0.714
Diffusion coefficient Dm 1.886 ? 10-5 cm2/s
0
20
40
60
80
100
0 20 40 60 80 100 120 140 160 180 200
Ele
va
tio
n (
cm
)
Length (cm)
g1817g1824 g3404g2783.g2783uni0020g3400g2778g2777uni0020g3398g2783uni0020g2195uni0020g2201g3398g2778uni0020g4666g2164g2187g2196g2200g2207uni0020g2198g2200g2197g2184g2194g2187g2195g4667
uni0020uni0020g3404g2780.g2780uni0020g3400g2778g2777uni0020g3398g2783uni0020g2195uni0020g2201g3398g2778uni0020g4666g2169g2197g2186g2191g2188g2191g2187g2186uni0020g2164g2187g2196g2200g2207uni0020g2198g2200g2197g2184g2194g2187g2195g4667
g1777g3404g2777
g1817g1826g3404g2777, g2770g1777g2770g1826g3404uni0020g2777
g1817g1826g3404g2777, g2770g1777g2770g1826g3404uni0020g2777 g1790uni0020g3404g2761
g1819g1807uni0020g1808g1819
uni0020uni0020uni0020g1777g3404g1777g1819uni0020g1809g1806uni0020g1796g1824g3407g2777
g2770g1777
g2770g1824g3404g2777uni0020g1809g1806uni0020g1796g1824g3410g2777
124
Figure A7.7. Steady-state 50% isochlor distribution for the Henry Problem and the modified
Henry problem.
0
20
40
60
80
100
0 20 40 60 80 100 120 140 160 180 200
z (
cm
)
x (cm)
VDFT - Modified Henry problem (2004)
VDFT - Henry problem (1964)
SEAWAT
125
Test Problem 4: Elder Problem
A conceptual model based on original Elder Problem [Elder, 1967] was simulated. This
study follow the aquifer property and parameter values suggested by Simpson and Clement
[2003]. The domain and boundary conditions for the Elder problem describe a confined aquifer,
a zero pressure head is maintained on the left and right upper corner of the domain (see Figure
A7.8). In order to construct a numerical model, the origin of an x-z coordinate system was set on
the lower left-hand corner of the plane. A grid of size of 16.8 m for x-direction and 6 m of z-
direction (44 cells in the x-direction and 25 cells in the z-direction) was used to discretize the
600 m ? 150 m dimensional model. The value of specific storage (Ss) was set to zero for all
nodes. All layers are set to confined layer and bottom of the model is set to no-flow condition.
The longitudinal and transversal dispersivity coefficients are set to zero. Cells in layer 1 from
columns 12 to 33 are imposed by constant concentration with a value of 285.7 kg/m3. The
Courant number is set to 0.1. The model summary is presented in Table A7.3.
The simulated results were compared to SEAWAT results and shown in Figure A7.9.
VDFT results show relatively smooth migration downward where SEAWAT shows clear
vortices in the model. The VDFT code is still under development to predict more accurate
density effects.
126
Figure A7.8. Domain and boundary conditions for the Elder problem (Simpson and Clement,
2003).
Table A7.4 Summary of numerical model parameters
Property Symbol values
Horizontal length L 600 m
Saltwater elevation hs 150 m
Hydraulic conductivity K 4.76 ? 10-6 cm/s
Porosity n 0.1
Specific storage, Ss 0 cm-1
Longitudinal dispersivity ?L 0 cm
Transverse dispersivity ?T 0 cm
Saltwater concentration Cs 286 kg/m3
Saltwater density ?s 1,200 kg/m3
Freshwater density ?f 1,000 kg/m3
Density Slope E 0.7
Diffusion coefficient Dm 3.57 ? 10-6 cm2/s
127
(a) 2 years at SEAWAT (b) 2 years at VDFT
(c) 4 years at SEAWAT (d) 4 years at VDFT
Figure A7.9. 20% (blue line) and 60% (red line) of the flow pattern of the dense fluid into the
aquifer for the Elder Problem.
A7.4 Discussion
The VDFT code currently reproduces the experiment data well; however, some discrepancies
where observed for certain benchmark problems (Henry problem and Elder problem). The
implicit transport package of VDFT code is still under development, and further work is needed
to improve this package.
0 100 200 300 400 500 6000
50
100
150
0 100 200 300 400 500 6000
50
100
150
0 100 200 300 400 500 6000
50
100
150
0 100 200 300 400 500 6000
50
100
150
128
Appendix 7 Code Formulation based on TVD Scheme for Advection Package
1. General Formulation (assume constant g2722g1824)
uni03B8uni0043g2919,g2920
g2924g2878g2869g3398uni0043
g2919,g2920
g2924
uni0394uni0074 g3404uni0020g3398
qg2919,g2920g2878g2869uni002Fg2870uni0043g2919,g2920g2878g2869uni002Fg2870g2924 g3398qg2919,g2920g2879g2869uni002Fg2870uni0043g2919,g2920g2879g2869uni002Fg2870g2924
uni0394uni0079 g3398
qg2919g2878g2869uni002Fg2870,g2920uni0043g2919g2878g2869uni002Fg2870,g2920g2924 g3398qg2919g2879g2869uni002Fg2870,g2920uni0043g2919g2879g2869uni002Fg2870,g2920g2924
uni0394uni0078
Rearranging terms
uni0043g2919,g2920g2924g2878g2869g3404uni0043g2919,g2920g2924 g3398g4666uni03B3g4667g2928g2919g2917g2918g2930uni0043g2919,g2920g2878g2869uni002Fg2870g2924 g3397g4666uni03B3g4667g2922g2915g2916g2930uni0043g2919,g2920g2879g2869uni002Fg2870g2924 g3398g4666uni03B3g4667g2912g2925g2930g2930g2925g2923uni0043g2919g2878g2869uni002Fg2870,g2920g2924 g3397g4666uni03B3g4667g2930g2925g2926uni0043g2919g2879g2869uni002Fg2870,g2920g2924
Where the Courant number, ?, is defined as
g4666uni03B3g4667g2922g2915g2916g2930g3404qg2919,g2920g2879g2869uni002Fg2870uni0394uni0074uni03B8uni0394uni0078 g3404vg2934g2919,g2920g2879g2869uni002Fg2870uni0394uni0074uni0020uni0394uni0078 uni0020
g4666uni03B3g4667g2928g2919g2917g2918g2930g3404qg2919,g2920g2878g2869uni002Fg2870uni0394uni0074uni03B8uni0394uni0078 g3404vg2934g2919,g2920g2878g2869uni002Fg2870uni0394uni0074uni0020uni0394uni0078 uni0020
g4666uni03B3g4667g2930g2925g2926g3404qg2919g2879g2869uni002Fg2870,g2920uni0394uni0074uni03B8uni0394uni0079 g3404vg2935g2919g2879g2869uni002Fg2870,g2920uni0394uni0074uni0020uni0394uni0079 uni0020
g4666uni03B3g4667g2912g2925g2930g2930g2925g2923g3404qg2919g2878g2869uni002Fg2870,g2920uni0394uni0074uni03B8uni0394uni0079 g3404vg2935g2919g2878g2869uni002Fg2870,g2920uni0394uni0074uni0020uni0394uni0079 uni0020
Applying the same procedure for third-order polynomial interpolation of nodel
concentrations as described in chapter 3 (MT3DMS manual), the concentration values at the left
interface can be derived as
uni0043g2919,g2920g2879g2869uni002Fg2870g2924 g3404uni0043g3364g2919,g2920g2879g2869uni002Fg2870g2924 g3398uni0020uni0064g29342 fg2934uni0020g3398uni0020uni0064g29352 fg2935g3398g3435uni0394uni0079
g2870g3398uni0064g2935g2870g3439
g888 fg2935g2935g3397g4678
uni0064g2934g2870
g888 g3398
uni0394uni0078uni0064g2934
g886 g4679fg2934g2934g3397g4678
uni0064g2934uni0064g2935
uni0033 g3398
uni0394uni0079uni0064g2934
g886 g4679fg2935g2934
Where,
uni0064g2935 g3404uni0020vg2934g2919,g2920g2879g2869uni002Fg2870uni0394uni0074
uni0064g2934 g3404uni0020vg2935g2919,g2920g2879g2869uni002Fg2870uni0394uni0074
129
uni0043g3364g2919,g2920g2879g2869uni002Fg2870g2924 g3404uni0043g2919,g2920g2879g2869
g2924 g3397uni0043
g2919,g2920
g2924uni0020uni0020
2 ,iuni006Euni0020uni0074hisuni0020suni0074uuni0064uni0079
The normal gradient
fg2935g3404uni0043g2919,g2920
g2924 g3398uni0043
g2919,g2920g2879g2869
g2924 uni0020uni0020
uni0394uni0079
The upwind gradient along the y-direction depends on the directions of the velocities
fg2934g3404g3435uni0043g2919,g2920g2879g2869g2924 g3398uni0043g2919g2879g2869,g2920g2879g2869g2924 g3439uni002Funi0394uni0078uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919,g2920g2879g2869uni002Fg2870g3408uni0030,vg2934uni0020g2919,g2920g2879g2869uni002Fg2870g3408uni0030uni0020uni0020uni0020uni0020uni0020uni0020uni0020
uni0020uni0020uni0020uni0020g3404g3435uni0043g2919g2878g2869,g2920g2879g2869g2924 g3398uni0043g2919,g2920g2879g2869g2924 g3439uni002Funi0394uni0078uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919,g2920g2879g2869uni002Fg2870g3408uni0030,vg2934uni0020g2919,g2920g2879g2869uni002Fg2870g3407uni0030uni0020uni0020uni0020uni0020uni0020uni0020uni0020
uni0020uni0020uni0020uni0020g3404g3435uni0043g2919,g2920g2924 g3398uni0043g2919g2879g2869,g2920g2924 g3439uni002Funi0394uni0078uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919,g2920g2879g2869uni002Fg2870g3407uni0030,vg2934uni0020g2919,g2920g2879g2869uni002Fg2870g3408uni0030uni0020uni0020uni0020uni0020uni0020uni0020uni0020
uni0020uni0020uni0020uni0020g3404g3435uni0043g2919g2878g2869,g2920g2924 g3398uni0043g2919,g2920g2924g3439uni002Funi0394uni0078uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919,g2920g2879g2869uni002Fg2870g3407uni0030,vg2934uni0020g2919,g2920g2879g2869uni002Fg2870g3407uni0030uni0020uni0020uni0020uni0020uni0020uni0020uni0020
Based on the first-order derivatives,
Second order derivatives can be calculated below:
fg2935g2935g3404uni0020uni0020g3428g3435fg2935g3439g2919,g2920g2879g2869
g2870
g3398g3435fg2935g3439g2919,g2920g2879g2869g2879g2869
g2870
g3432uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919,g2920g2879g2869uni002Fg2870g3408uni0030
uni0020uni0020uni0020uni0020uni0020g3404uni0020g3428g3435fg2935g3439g2919,g2920g2878g2869
g2870
g3398g3435fg2935g3439g2919,g2920g2879g2869
g2870
g3432uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919,g2920g2879g2869uni002Fg2870g3407uni0030
fg2934g2934g3404uni0020uni0020g3428g4666fg2934g4667g2919g2878g2869
g2870,g2920g2879g2869
g3398g4666fg2934g4667g2919g2879g2869
g2870,g2920g2879g2869
g3432uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919,g2920g2879g2869uni002Fg2870g3408uni0030
uni0020uni0020uni0020uni0020uni0020g3404uni0020g3428g4666fg2934g4667g2919g2878g2869
g2870,g2920
g3398g4666fg2934g4667g2919g2879g2869
g2870,g2920
g3432uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919,g2920g2879g2869uni002Fg2870g3407uni0030
Mixed second-order derivatives (twist terms)
fg2935g2934g3404uni0020uni0020g3428g4666fg2934g4667g2919g2879g2869
g2870,g2920
g3398g4666fg2934g4667g2919g2879g2869
g2870,g2920g2879g2869
g3432uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2934g2919,g2920g2879g2869uni002Fg2870g3408uni0030
uni0020uni0020uni0020uni0020uni0020g3404uni0020g3428g4666fg2934g4667g2919g2878g2869
g2870,g2920
g3398g4666fg2934g4667g2919g2878g2869
g2870,g2920g2879g2869
g3432uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2934g2919,g2920g2879g2869uni002Fg2870g3407uni0030
uni0043g2919g2879g2869uni002Fg2870,g2920g2924
130
uni0043g2919g2879g2869uni002Fg2870,g2920g2924 g3404uni0043g3364g2919g2879g2869uni002Fg2870,g2920g2924 g3398uni0020uni0064g29342 fg2934uni0020g3398uni0020uni0064g29352 fg2935g3398g4666uni0394uni0078
g2870g3398uni0064g2934g2870g4667
g888 fg2934g2934g3397g4678
uni0064g2935g2870
g888 g3398
uni0394uni0079uni0064g2935
g886 g4679fg2935g2935g3397g4678
uni0064g2934uni0064g2935
uni0033 g3398
uni0394uni0078uni0064g2935
g886 g4679fg2934g2935
Where,
uni0064g2934 g3404uni0020vg2934g2919g2879g2869uni002Fg2870,g2920uni0394uni0074
uni0064g2935 g3404uni0020vg2935g2919g2879g2869uni002Fg2870,g2920uni0394uni0074
uni0043g3364g2919g2879g2869uni002Fg2870,g2920g2924 g3404uni0043g2919g2879g2869,g2920
g2924 g3397uni0043
g2919,g2920
g2924uni0020uni0020
2 ,iuni006Euni0020uni0074hisuni0020suni0074uuni0064uni0079
The normal gradient
fg2934g3404uni0043g2919,g2920
g2924 g3398uni0043
g2919g2879g2869,g2920
g2924 uni0020uni0020
uni0394uni0078
The upwinding gradient along the y-direction depends on the directions of the velocities
fg2935g3404g3435uni0043g2919g2879g2869,g2920g2924 g3398uni0043g2919g2879g2869,g2920g2879g2869g2924 g3439uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919g2879g2869uni002Fg2870,g2920g3408uni0030,vg2934g2919g2879g2869uni002Fg2870,g2920g3408uni0030uni0020uni0020uni0020uni0020uni0020uni0020uni0020
uni0020uni0020uni0020uni0020g3404g3435uni0043g2919,g2920g2924 g3398uni0043g2919,g2920g2879g2869g2924 g3439uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919g2879g2869uni002Fg2870,g2920g3408uni0030,vg2934g2919g2879g2869uni002Fg2870,g2920g3407uni0030uni0020uni0020uni0020uni0020uni0020uni0020uni0020
uni0020uni0020uni0020uni0020g3404g3435uni0043g2919g2879g2869,g2920g2878g2869g2924 g3398uni0043g2919g2879g2869,g2920g2924 g3439uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919g2879g2869uni002Fg2870,g2920g3407uni0030,vg2934g2919g2879g2869uni002Fg2870,g2920g3408uni0030uni0020uni0020uni0020uni0020uni0020uni0020uni0020
uni0020uni0020uni0020uni0020g3404g3435uni0043g2919,g2920g2878g2869g2924 g3398uni0043g2919,g2920g2924g3439uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919g2879g2869uni002Fg2870,g2920g3407uni0030,vg2934g2919g2879g2869uni002Fg2870,g2920g3407uni0030uni0020uni0020uni0020uni0020uni0020uni0020uni0020
Based on the first-order derivatives,
Second order derivatives can be calculated below:
fg2934g2934g3404uni0020uni0020g3428g4666fg2934g4667g2919g2879g2869
g2870,g2920
g3398g4666fg2934g4667g2919g2879g2869g2879g2869
g2870,g2920
g3432uni002Funi0394uni0078uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2934g2919,g2920g2879g2869uni002Fg2870g3408uni0030
uni0020uni0020uni0020uni0020uni0020g3404uni0020g3428g4666fg2934g4667g2919g2878g2869
g2870,g2920
g3398g4666fg2934g4667g2919g2879g2869
g2870,g2920
g3432uni002Funi0394uni0078uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2934g2919,g2920g2879g2869uni002Fg2870g3407uni0030
fg2935g2935uni0020uni0020uni0020uni0020g3404uni0020g3428g3435fg2935g3439g2919g2879g2869,g2920g2878g2869
g2870
g3398g3435fg2935g3439g2919g2879g2869,g2920g2879g2869
g2870
g3432uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2934g2919,g2920g2879g2869uni002Fg2870g3408uni0030
uni0020uni0020uni0020uni0020uni0020g3404uni0020g3428g3435fg2935g3439g2919,g2920g2878g2869
g2870
g3398g3435fg2935g3439g2919,g2920g2879g2869
g2870
g3432uni002Funi0394uni0079uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2934g2919,g2920g2879g2869uni002Fg2870g3407uni0030
131
Mixed second-order derivatives (twist terms)
fg2934g2935g3404uni0020uni0020g3428g3435fg2935g3439g2919,g2920g2879g2869
g2870
g3398g3435fg2935g3439g2919g2879g2869,g2920g2879g2869
g2870
g3432uni002Funi0394uni0078uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919g2879g2869uni002Fg2870,g2920g3408uni0030
uni0020uni0020uni0020uni0020uni0020g3404uni0020g3428g3435fg2935g3439g2919uni0020,g2920g2878g2869
g2870
g3398g3435fg2935g3439g2919g2879g2869,g2920g2878g2869
g2870
g3432uni002Funi0394uni0078uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020vg2935g2919g2879g2869uni002Fg2870,g2920g3407uni0030
2. Limiter Application
uni0043g3560g3404uni0020uni0043g3398uni0043g2920g2879g2870uni0043
g2920g3398uni0043g2920g2879g2870
Note that uni0020
uni0043g3560g2920g2879g2870g3404uni0020uni0030uni0020uni0061uni006Euni0064uni0020uni0020uni0043g3560g2920uni0020g3404uni0020uni0031uni0020uni0020
If the concentration profile across the three nodes are monotonic,
uni0043g3560g2920g2879g2869uni0020g3409uni0043g3560g2920g2879g2869
g2870uni0020
g3409uni0031
Which equal to
uni0043g2920g2879g2869g3398uni0043g2920g2879g2870
uni0043g2920g3398uni0043g2920g2879g2870 g3409
uni0043g2920g2879g2869
g2870
g3398uni0043g2920g2879g2870
uni0043g2920g3398uni0043g2920g2879g2870 g3409uni0031
uni0043g2920g2879g2869g3398uni0043g2920g2879g2870g3409uni0043g2920g2879g2869
g2870
g3398uni0043g2920g2879g2870g3409uni0043g2920g3398uni0043g2920g2879g2870
The right-hand inequality in above equation is insured
uni0043g2920g2879g2869
g2870
uni0020g3409uni0043g2920uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020couni006Esuni0074runi0061iuni006Euni0074uni0020g4666uni0031g4667uni0020
The left-hand inequality in above equation is insured
uni0043g2920g2879g2869
g2870
uni0020uni0020g3410uni0043g2920g2879g2869uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020couni006Esuni0074runi0061iuni006Euni0074uni0020g46662g4667
Also, Explicit finite difference solution at note j-1 is
132
uni0043g2920g2879g2869g2924g2878g2869g3404uni0043g2920g2879g2869g2924 g3398uni0043runi0020g4666uni0043g2920g2879g2869
g2870
g2924uni0020 g3398uni0043
g2920g2879g2869g2879g2869g2870
g2924uni0020 g4667uni0020
Is transformed in terms of normalized variables,
uni0043g3560g2920g2879g2869g2924g2878g2869g3404uni0043g3560g2920g2879g2869g2924 g3398uni0043runi0020g4666uni0043g3560g2920g2879g2869
g2870
g2924uni0020 g3398uni0043g3560
g2920g2879g2869g2879g2869g2870
g2924uni0020 g4667uni0020
Since
uni0043g3560g2920g2879g2869
g2870
g2924uni0020 uni0020g3410uni0043g3560
g2920g2879g2869g2879g2869g2870
g2924uni0020 uni0020uni0020uni0020uni0020uni0064ueuni0020uni0074ouni0020mouni006Eouni0074ouni006Eic
uni0030g3409uni0043g3560g2920g2879g2869g2924g2878g2869uni0020g3409uni0043g3560g2920g2879g2869g2924 g3409uni0031
Left-hand inequality implies
uni0030g3409uni0043g3560g2920g2879g2869g2924g2878g2869
uni0030g3409uni0043g3560g2920g2879g2869g2924 g3398uni0043g2928uni0020g4666uni0043g3560g2920g2879g2869
g2870
g2924uni0020 g3398uni0043g3560
g2920g2879g2869g2879g2869g2870
g2924uni0020 g4667uni0020
uni0043g3560g2920g2879g2869
g2870
g2924uni0020 g3409uni0043g3560
g2920g2879g2869g2879g2869g2870
g2924uni0020 g3397uni0043g3560
g2920g2879g2869
g2924 uni0020uni002Funi0043g2928
The worst ?case estimate for uni0043g3560g2920g2879g2869g2879g3117
g3118
g2924uni0020 g3404uni0030 results in additional constraints on uni0043g3560
g2920g2879g3117g3118
g2924uni0020 :
uni0043g3560g2920g2879g2869
g2870
g2924uni0020 g3409uni0020uni0043g3560
g2920g2879g2869
g2924 uni0020uni002Funi0043g2928uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020uni0020g1855g1867g1866g1871g1872g1870g1853g1872g1861g1866g1872uni0020g4666uni0033g4667
If uni0043g3560g2920g2879g3117
g3118
g2924uni0020 does not lie with above boundary,
uni0043g3560g2920g2879g2869
g2870
g2924uni0020 g3404uni0020uni0043g3560
g2920g2879g2869
g2924uni0020
Therefore
uni0043g2920g2879g2869
g2870
g2924uni0020 g3404uni0020uni0043
g2920g2879g2869
g2924uni0020
133
References
Abarca, E., and T. P. Clement (2009), A novel approach for characterizing the mixing zone of a
saltwater wedge, Geophysical Research Letters, 36(6).
Abarca, E., J. Carrera, X. S?nchez-Vila, and C. I. Voss (2007a), Quasi-horizontal circulation
cells in 3D seawater intrusion, Journal of Hydrology, 339(3-4), 118-129.
Abarca, E., J. Carrera, X. S?nchez-Vila, and M. Dentz (2007b), Anisotropic dispersive Henry
problem, Advances in Water Resources, 30(4), 913-926.
Anderson, J., K. Milliken, D. Wallace, A. Rodriguez, and A. Simms (2010), Coastal Impact
Underestimated From Rapid Sea Level Rise, Eos Trans. AGU, 91(23), 205-206, doi:
10.1029/2010EO230001
Barlow, P. M., and E. G. Reichard (2010), Saltwater intrusion in coastal regions of North
America, Hydrogeology Journal, 18(1), 247-260.
Bear, J., and A. H. D. Cheng (2010), Modeling groundwater flow and contaminant transport, 1st
edition ed., Springer.
Blanchon, P., A. Eisenhauer, J. Fietzke, and V. Liebetrau (2009), Rapid sea-level rise and reef
back-stepping at the close of the last interglacial highstand, Nature, 458(7240), 881-884.
Boufadel, M. C. (2000), A mechanistic study of nonlinear solute transport in a groundwater-
surface water system under steady state and transient hydraulic conditions, Water Resour. Res.,
36(9), 2549-2565.
Chang, S. W., and T. P. Clement (2012a), Laboratory and numerical investigation of transport
processes occurring beneath a saltwater wedge Journal of Contaminant Hydrology, sumbitted.
Chang, S. W., and T. P. Clement (2012b), Experimental and numerical investigation of saltwater
intrusion dynamics in flux controlled groundwater systems, Water Resources Research,
doi:10.1029/2012WR012134, In press.
Chang, S. W., T. P. Clement, M. J. Simpson, and K.-K. Lee (2011), Does sea-level rise have an
impact on saltwater intrusion?, Advances in Water Resources, 34(10), 1283-1291.
Cheng, A. H. D., D. Halhal, A. Naji, and D. Ouazar (2000), Pumping optimization in saltwater-
intruded coastal aquifers, Water Resources Research, 36(8), 2155-2165.
134
Clement, T. P., Y. Sun, B. S. Hooker, and J. N. Petersen (1998), Modeling Multispecies Reactive
Transport in Ground Water, Ground Water Monitoring & Remediation, 18(2), 79-92.
Croucher, A. E., and M. J. O'Sullivan (1995), The Henry Problem for Saltwater Intrusion, Water
Resour. Res., 31(7), 1809-1814.
Custodio, E., and G. A. Bruggeman (1987), Groundwater problems in coastal areas. Studies and
reports in hydrology 45, edited, UNESCO, Paris.
Dausman, A., and C. Langevin (2005), Movement of the saltwater interface in the Surficial
Aquifer System in response to hydrologic stresses and water-management practices, Broward
County, Florida: USGS Scientific Investigations Report : SIR 2004-5256Rep.
Destouni, G., and C. Prieto (2003), On the possibility for generic modeling of submarine
groundwater discharge, Biogeochemistry, 66(1), 171-186.
Diersch, H. J. G. (2002), FEFLOW finite element subsurface flow and transport simulation
system -user's manual/reference manual/whilte papers. Release 5.0Rep., Berlin.
Elder, J. W. (1967), Steady free convection in a porous medium heated from below, Journal of
Fluid Mechanics, 27(01), 29-48.
Feseker, T. (2007), Numerical studies on saltwater intrusion in a coastal aquifer in northwestern
Germany, Hydrogeology Journal, 15(2), 267-279.
Fetter, C. W. (2001), Applied Hydrogeology, Prentice Hall, N.J.
Frind, E. O. (1982), Simulation of long-term transient density-dependent transport in
groundwater, Advances in Water Resources, 5(2), 73-88.
Giambastiani, B. M. S., M. Antonellini, G. H. P. Oude Essink, and R. J. Stuurman (2007),
Saltwater intrusion in the unconfined coastal aquifer of Ravenna (Italy): A numerical model,
Journal of Hydrology, 340(1-2), 91-104.
Goswami, R. R., and T. P. Clement (2007), Laboratory-scale investigation of saltwater intrusion
dynamics, Water Resources Research, 43(4), doi: 10.1029/2006WR005151.
Goswami, R. R., B. Ambale, and T. P. Clement (2009), Estimating errors in concentration
measurements obtained from image analysis, Vadose Zone Journal, 8(1), 108-118.
Goswami, R. R., T. P. Clement, and J. H. Hayworth (2012), Comparison of Numerical
Techniques Used for Simulating Variable-Density Flow and Transport Experiments, Journal of
Hydrologic Engineering, 17(2), 272-282.
135
Green, T. R., M. Taniguchi, H. Kooi, J. J. Gurdak, D. M. Allen, K. M. Hiscock, H. Treidel, and
A. Aureli (2011), Beneath the surface of global change: Impacts of climate change on
groundwater, Journal of Hydrology, 405(3-4), 532-560.
Guo, W., and C. D. Langevin (2002), User's Guide to SEWAT: A Computer Program for
Simulation of Three-Dimensional Variable-Density Ground-Water Flow, United States
Geological Survey.
Henry, H. R. (1964), Effects of dispersion on salt encroachment in coastal aquifers, U.S.
Geological Survey Water-Supply Paper Rep. 1613-C, C71-C84. pp.
Huyakorn, P. S., P. F. Andersen, J. W. Mercer, and H. O. White, Jr. (1987), Saltwater intrusion
in aquifers: Development and testing of a three-dimensional finite element model, Water Resour.
Res., 23(2), 293-312.
Illangasekare, T., et al. (2006), Impacts of the 2004 tsunami on groundwater resources in Sri
Lanka, Water Resour. Res., 42(5), W05201.
IPCC (2007), Climate change 2007: Impacts, adaptation and vulnerability, Contribution of
Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate
Change, edited by M. L. Parry, O. F. Canziani, J. P. Palutikof, P. J. van der Linden and C. E.
Hanson, Cambridge University Press Cambridge.
Johannes, R. E. (1980), The ecological significance of the submarine discharge of groundwater,
Marine Ecology Progress Series 3, 365?-373.
Kanel, S. R., R. R. Goswami, T. P. Clement, M. O. Barnett, and D. Zhao (2008), Two
dimensional transport characteristics of surface stabilized zero-valent iron nanoparticles in
porous media, Environ Sci Technol, 42(3), 896-900.
Karasaki, K., K. Ito, and K. Maekawa (2006), Simulation of salt water intruison, paper presented
at TOUGH Symposium, Lawrence Berkeley National Laboratory, Berkely, California.
Kohout, F. A. (1960), Cyclic Flow of Salt Water in the Biscayne Aquifer of Southeastern Florida,
J. Geophys. Res., 65(7), 2133-2141.
Kuan, W. K., G. Jin, P. Xin, C. Robinson, B. Gibbes, and L. Li (2012), Tidal influence on
seawater intrusion in unconfined coastal aquifers, Water Resour. Res., 48(2), W02502.
Kundzewicz, Z. W., and P. D?lle (2009), Will groundwater ease freshwater stress under climate
change?, Hydrological Sciences Journal, 54(4), 665-675.
Lacombe, P. J., and G. B. Carleton (1992), Saltwater intrusion into fresh ground-water supplies,
southern Cape May County, New Jersey, 1890-1991.
136
Langevin C, Shoemaker WB, and G. W (2003), Modflow-2000, The USGS Geological Survey
modular ground-water model ? documentation of the Seawat-2000 version with the variable
density flow process (VDF) and the integrated MT3DMS transport process (IMT)Rep., US
Geological Survey.
Leatherman, S. P. (1984), Coastal geomorphic response to sea level rise : Galveston Bay, Texas,
in Barth M.C., Titus J.G. (eds) Greenhouse effect and sea level rise : a challenge for this
generation.Rep., 151-178 pp, Reinhold, New York.
Lee, C. H., and R. T. S. Cheng (1974), On seawater encroachment in coastal aquifers, Water
Resources Research, 10, 1039-1043.
Li, H., and J. J. Jiao (2003a), Influence of the tide on the mean watertable in an unconfined,
anisotropic, inhomogeneous coastal aquifer, Advances in Water Resources, 26(1), 9-16.
Li, H., and J. J. Jiao (2003b), Tide-induced seawater-groundwater circulation in a multi-layered
coastal leaky aquifer system, Journal of Hydrology, 274(1-4), 211-224.
Li, H., M. C. Boufadel, and J. W. Weaver (2008), Tide-induced seawater-groundwater
circulation in shallow beach aquifers, Journal of Hydrology, 352(1-2), 211-224.
Li, L., D. A. Barry, F. Stagnitti, and J. Y. Parlange (1999), Submarine groundwater discharge
and associated chemical input to a coastal sea, Water Resour. Res., 35(11), 3253-3259.
Lo?iciga, H. A., T. J. Pingel, and E. S. Garcia (2011), Sea Water Intrusion by Sea-Level Rise:
Scenarios for the 21st Century, Ground Water, no-no.
Mantoglou, A. (2003), Pumping management of coastal aquifers using analytical models of
saltwater intrusion, Water Resources Research, 39(12), 1335.
Masterson, J. P., and S. P. Garabedian (2007), Effects of Sea-Level Rise on Ground Water Flow
in a Coastal Aquifer System, Ground Water, 45(2), 209-217.
McCarthy, J. J. (2009), Reflections On: Our Planet and Its Life, Origins, and Futures, Science,
326(5960), 1646, doi: 10.1126/science.1184937.
Meisler, H., Leahy PP and Knobel LL. (1985), The effect of eustatic sea-level changes on
saltwaterfreshwater relations in the northern Atlantic Coastal Plain: U.S. Geological Survey
Water-Supply Paper 2255, 28 p.Rep.
Melloul, A., and M. Collin (2006), Hydrogeological changes in coastal aquifers due to sea level
rise, Ocean & Coastal Management, 49(5-6), 281-297.
Michael, H. A., A. E. Mulligan, and C. F. Harvey (2005), Seasonal oscillations in water
exchange between aquifers and the coastal ocean, Nature, 436(7054), 1145-1148.
137
Moore, W. S. (1996), Large groundwater inputs to coastal waters revealed by 226Ra enrichments,
Nature, 380(6575), 612-614.
Moore, W. S., and T. M. Church (1996), Submarine Groundwater Discharge, reply to Younger.,
Nature, 382, 122.
Navoy, A. S. (1991), Aquifer-estuary interaction and vulnerability of groundwater supplies to sea
level rise-driven saltwater intrusion, Pennsylvania State University, USA.
Oostrom, M., J. S. Hayworth, J. H. Dane, and O. G?en (1992), Behavior of dense aqueous phase
leachate plumes in homogeneous porous media, Water Resour. Res., 28(8), 2123-2134.
Oude Essink, G. H. P. (2001), Salt Water Intrusion in a Three-dimensional Groundwater System
in The Netherlands: A Numerical Study, Transport in Porous Media, 43(1), 137-158.
Oude Essink, G. H. P., E. S. van Baaren, and P. G. B. de Louw (2010), Effects of climate change
on coastal groundwater systems: A modeling study in the Netherlands, Water Resour. Res., 46,
W00F04.
Pinder, G. F., and H. H. Cooper, Jr. (1970), A Numerical Technique for Calculating the
Transient Position of the Saltwater Front, Water Resour. Res., 6(3), 875-882.
Ranjan, P., S. Kazama, and M. Sawamoto (2006), Effects of climate change on coastal fresh
groundwater resources, Global Environmental Change, 16(4), 388-399.
Robinson, C., L. Li, and H. Prommer (2007), Tide-induced recirculation across the aquifer-ocean
interface, Water Resour. Res., 43(7), W07428.
Rozell, D. J., and T.-F. Wong (2010), Effects of climate change on groundwater resources at
Shelter Island, New York State, USA, Hydrogeology Journal, 18, 1657-1665.
Schincariol, R. A., and F. W. Schwartz (1990), An experimental investigation of variable density
flow and mixing in homogeneous and heterogeneous media, Water Resour. Res., 26(10), 2317-
2329.
Segol, G. (1994), Classic Groundwater Simulations Proving and Improving
Numerical Models, Old Tappan, N. J.
Segol, G., G. F. Pinder, and W. G. Gray (1975), A Galerkin-finite element technique for
calculating the transient position of the saltwater front, Water Resour. Res., 11(2), 343-347.
Simmons, G. M. (1992), Importance of submarine groundwater discharge (Sgwd) and seawater
cycling to material flux across sediment water interfaces in marine environments., Marine
Ecology Progress Series 84, 173-184.
138
Simpson, M. J., and T. P. Clement (2003), Theoretical analysis of the worthiness of Henry and
Elder problems as benchmarks of density-dependent groundwater flow models, Advances in
Water Resources, 26(1), 17-31.
Simpson, M. J., and T. P. Clement (2004), Improving the worthiness of the Henry problem as a
benchmark for density-dependent groundwater flow models, Water Resour. Res., 40(1), W01504.
Simpson, M. J., T. P. Clement, and T. A. Gallop (2003), Laboratory and Numerical Investigation
of Flow and Transport Near a Seepage-Face Boundary, Ground Water, 41(5), 690-700.
Smith, A. J. (2004), Mixed convection and density-dependent seawater circulation in coastal
aquifers, Water Resour. Res., 40(8), W08309.
Strack, O. D. L. (1976), A single-potential solution for regional interface problems in coastal
aquifers, Water Resources Research, 12(6), 1165-1174.
Sun, Y., and T. P. Clement (1999), A Decomposition Method for Solving Coupled
Multi-Species Reactive Transport Problems, Transport in Porous Media, 37, 327-346.
USGS (2000), Ground-Water Resources for the Future - Atlantic Coastal ZoneRep.
Vermeer, M., and S. Rahmstorf (2009), Global sea level linked to global temperature, Proc. Natl.
Acad. Sci., 106(51), 21527-21532, doi: 10.1073/pnas.0907765106
Volker, R. E., Q. Zhang, and D. A. Lockington (2002), Numerical modelling of contaminant
transport in coastal aquifers, Mathematics and Computers in Simulation, 59(1??), 35-44.
Voss, C. I., and W. R. Souza (1987), Variable density flow and solute transport simulation of
regional aquifers containing a narrow freshwater-saltwater transition zone, Water Resour. Res.,
23(10), 1851-1866.
Voss, C. I., and A. M. Provost (2002), A model for saturated-unsaturated, variable-density
ground-water flow with solute or energy transportRep., U.S. Geol. Surv., Reston, Va.
Watson, T. A., A. D. Werner, and C. T. Simmons (2010), Transience of seawater intrusion in
response to sea level rise, Water Resour. Res., 46(12), W12533.
Webb, M. D., and K. W. F. Howard (2010), Modeling the Transient Response of Saline Intrusion
to Rising Sea-Levels, Ground Water.
Werner, A. D., and C. T. Simmons (2009), Impact of Sea-Level Rise on Sea Water Intrusion in
Coastal Aquifers, Ground Water, 47(2), 197-204, doi: 10.1111/j.1745-6584.2008.00535.x.
Werner, A. D., J. D. Ward, L. K. Morgan, C. T. Simmons, N. I. Robinson, and M. D. Teubner
(2011), Vulnerability Indicators of Sea Water Intrusion, Ground Water, no-no.
139
Westbrook, S. J., J. L. Rayner, G. B. Davis, T. P. Clement, P. L. Bjerg, and S. J. Fisher (2005),
Interaction between shallow groundwater, saline surface water and contaminant discharge at a
seasonally and tidally forced estuarine boundary, Journal of Hydrology, 302, 255-269.
White, I., and T. Falkland (2009), Management of freshwater lenses on small Pacific islands,
Hydrogeology Journal, 18(1), 227-246.
WWD (1998), Groundwater: the Invisible Resource.
Yu, W., C. Voss, H. Michael, K. M. Ahmed, L. Feinson, M. R. Khan, and A. Tuinhof (2010),
Implications of climate change for fresh groundwater resources in coastal aquifers in
BangladeshRep., The world bank
US Geological Survey.
Zhang, Q., R. E. Volker, and D. A. Lockington (2002), Experimental investigation of
contaminant transport in coastal groundwater, Advances in Environmental Research, 6(3), 229-
237.