UNDERSTANDING THE FLOW AND MIXING DYNAMICS OF SALINE WATER DISCHARGED INTO COASTAL FRESHWATER AQUIFERS Except where reference is made to the work of others, the work described in this thesis is my own or was done in collaboration with my advisory committee. This thesis does not include proprietary or classified information. ________________________ Matthew Brooks Hogan Certificate of Approval: ________________________ ________________________ Joel G. Melville T. Prabhakar Clement, Chair Professor Associate Professor Civil Engineering Civil Engineering ________________________ ________________________ Dongye Zhao Joe F. Pittman Associate Professor Interim Dean Civil Engineering Graduate School UNDERSTANDING THE FLOW AND MIXING DYNAMICS OF SALINE WATER DISCHARGED INTO COASTAL FRESHWATER AQUIFERS Matthew Brooks Hogan A Thesis Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Master of Science Auburn, Alabama December 15, 2006 iii UNDERSTANDING THE FLOW AND MIXING DYNAMICS OF SALINE WATER DISCHARGED INTO COASTAL FRESHWATER AQUIFERS Matthew Brooks Hogan Permission is granted to Auburn University to make copies of this thesis at its discretion, upon the request of individuals or institutions and at their expense. The author reserves all publication rights. _______________________ Signature of Author _______________________ Date of Graduation iv VITA Matthew B. Hogan, son of Mark B. Hogan and Marsha K. Hogan, was born on May 18, 1981, in Leesburg, Virginia. He graduated from Broad Run High School, in Ashburn, Virginia in the spring of 1999. In the spring of 2004, he graduated from Virginia Polytechnic Institute and State University (Virginia Tech) with a Bachelor?s Degree in the field of Civil Engineering. He entered the Graduate School at Auburn University in the fall of 2004 to begin his graduate studies at Auburn University to pursue a Master of Science degree in the field of Civil Engineering specializing in groundwater hydrology. v THESIS ABSTRACT UNDERSTANDING THE FLOW AND MIXING DYNAMICS OF SALINE WATER DISCHARGED INTO COASTAL FRESHWATER AQUIFERS Matthew B. Hogan Master of Science, December 15, 2006 (B.S., Virginia Polytechnic and State University, 2004) 114 Typed Pages Directed by T. Prabhakar Clement Typically the origin of saltwater intrusion into coastal aquifers has been attributed to the reduction of freshwater flux (or head) in the aquifer. This intrusion is the most common natural source of saltwater contamination in coastal aquifers; however, it is not the only source. Other sources of saltwater contamination that include saltwater deposition caused by tsunamis and hurricanes have been observed. The country of Sri Lanka was devastated by a catastrophic tsunami event in 2004. In this study, laboratory experiments were completed using a physical two-dimensional groundwater model to observe subsurface saltwater transport processes in a coastal aquifer after a tsunami event. vi Experiments were completed to investigate the transport from different types of saltwater sources. The first experiment was conducted to determine the fate of the overland saltwater that was discharged through the vadose zone. This experiment specifically considered the infiltration of saltwater along the beach face and the injection of saltwater through an open well. The next experiment investigated an isolated beach face infiltration source. The final experiment considered the fate of the saltwater deposited into a coastal lagoon and later leached into a coastal aquifer. The experimental datasets were then modeled using the USGS code SEAWAT. The numerical solutions generated by SEAWAT matched the physical model datasets well. An analysis was performed to further investigate the sensitivity of the plume to certain several flow parameters include freshwater flux, density variations, source loading rate and dispersivity. The results of the sensitivity analysis indicate that the source loading rate is the least sensitive model parameter. In this study, a conceptual basis was established for understanding the saltwater transport processes occurring from multiple sources after a large-scale inundation (tsunami-type) event. These conceptual experiments are prototypes for future field studies, and they provide insight into the density driven transport processes that would occur in a local unconfined aquifer after a tsunami-type event. vii ACKNOWLEDGMENTS The author would like to thank his advisor Dr. T. Prabhakar Clement who has helped the author immensely. The author sincerely thanks Dr. Clement for his support and for changing his life by introducing him to a new level of excellence. The author extends his thanks to his committee members, Dr. Joel Melville and Dr. Dongye Zhao, who have been supportive and have the author?s utmost respect. The next group of people the author would like to thank is his office mates. The author?s experiences over the past two years at Auburn University would not have been as memorable without them. The author cannot put into words what they have meant to him. The author thanks Anjani, Gautham, Justin, Linzy, Rohit, Sushil, Venkat, and Vijay. Special thanks to Mike Murdock, Xu Chen, and Jimmy Walker who will make very fine colleagues someday. The author would also like to extend his absolute gratitude towards his friends in Auburn, back home in Virginia, and else where. Finally, thanks to Auburn University and the entire staff of the Department of Civil Engineering. The author of this study would like to dedicate this thesis to his family, who has been there for him time and time again. Thank you. viii Style manual or journal used Auburn University Graduate School Guide to Preparation of Master?s Thesis Computer Software used Microsoft Office XP Professional (Microsoft Word, Microsoft Excel), SEAWAT 2000, Groundwater Vistas 4.0, Endnote 9.0, Solid Edge (CAD) ix TABLE OF CONTENTS List of Tables .................................................................................................................... xi List of Figures .................................................................................................................. xii Nomenclature .................................................................................................................. xiv Chapter 1. Introduction ................................................................................................................... 1 1.1 Background ................................................................................................................. 1 1.2 Sri Lanka Tsunami Scenario........................................................................................ 7 1.3 Research Goals and Objectives ................................................................................. 10 1.4 Organization of Thesis .............................................................................................. 11 2. Literature Review ....................................................................................................... 12 2.1 Introduction ............................................................................................................... 12 2.2 Groundwater Models Definition ............................................................................... 13 2.2.1 Analytical Models ........................................................................................... 13 2.2.2 Numerical Models ........................................................................................... 14 2.2.2.1 Review of Density Coupled Flow Codes ........................................... 19 2.2.3 Physical Models .............................................................................................. 21 2.2.3.1 Review of Physical Model Studies .................................................... 32 3. Governing Equations of Density Driven Flow ........................................................... 33 3.1 Governing Equations ................................................................................................. 33 3.1.1 Equivalent Freshwater Head ........................................................................... 34 3.1.2 Derivation of General Groundwater Flow Equation ....................................... 36 3.1.3 Equivalent Freshwater Head Form of Groundwater Flow Equation .............. 38 3.1.4 Solute Transport Equation .............................................................................. 41 3.2 SEAWAT Coupling Procedure ................................................................................. 43 4. Materials and Methods of Laboratory Study .............................................................. 44 4.1 Flow Container Design and Construction ................................................................. 44 x 4.2 Experimental Setup ................................................................................................... 48 4.3 Beach Face Infiltration with Well Experiment ......................................................... 50 4.4 Beach Face Infiltration Experiment .......................................................................... 50 4.5 Sinking Plume Experiment ....................................................................................... 52 4.6 Porous Medium ......................................................................................................... 53 4.7 Dye ............................................................................................................................ 54 4.8 Hydraulic Conductivity ............................................................................................. 54 4.8.1 Technique 1 ..................................................................................................... 54 4.8.2 Technique 2 ..................................................................................................... 56 4.8.3 Technique 3 ..................................................................................................... 57 4.9 Porosity and Dispersivity Values .............................................................................. 57 5. Results of Physical Model Experiments ..................................................................... 59 5.1 Problem Definition ................................................................................................... 59 5.2 Beach Face Infiltration with Well Experiment ......................................................... 62 5.3 Beach Face Infiltration Experiment .......................................................................... 65 5.4 Sinking Plume Experiment ....................................................................................... 67 5.4.1 Discussion of Plume Stability ......................................................................... 80 5.4.2 SEAWAT Simulations .................................................................................... 81 5.4.3 Sensitivity Analysis ........................................................................................ 84 6. Summary and Conclusions .......................................................................................... 92 References ........................................................................................................................ 97 xi LIST OF TABLES Table 1: Porosity Measurement ...................................................................................... 58 Table 2: Summary of Sinking Plume Experiment .......................................................... 68 Table 3: Plume Stability ................................................................................................. 81 Table 4: Numerical Simulation Base Case Parameters .................................................. 85 xii LIST OF FIGURES Figure 1.1: Stable Interface Problem ................................................................................ 2 Figure 1.2: Unstable Interface Problem ............................................................................ 3 Figure 1.3 Plume Stability ................................................................................................ 4 Figure 1.4: Cross Sectional View of Coastal Aquifer ...................................................... 6 Figure 1.5: Inundation Distance and Run-up Elevation ................................................... 7 Figure 1.6: Maruthumunai, Sri Lanka .............................................................................. 8 Figure 2.1: Henry Saltwater Intrusion Problem .............................................................. 16 Figure 2.2: Elder Problem ............................................................................................... 17 Figure 2.3: Hydrocoin Benchmark Problem ................................................................... 18 Figure 2.4: Diagram of Oostrom Flow Container ........................................................... 27 Figure 2.5: Shackelford Fingering Experiment .............................................................. 30 Figure 3.1: Conceptual Figure of Hydraulic Head ......................................................... 34 Figure 3.2: Representative Elementary Volume (REV) ................................................. 36 Figure 3.3: SEAWAT Coupling ..................................................................................... 43 Figure 4.1: Flow Container ............................................................................................. 46 Figure 4.2: Flow Container ............................................................................................. 47 Figure 4.3: Flow Container ............................................................................................. 47 Figure 4.4: Experiment Setup ......................................................................................... 48 Figure 4.5: Beach Face Infiltration with Open Well ...................................................... 50 Figure 4.6: Beach Face Infiltration ................................................................................. 51 Figure 4.7: Sinking Plume Experiment ........................................................................... 53 Figure 4.8: Hydraulic Conductivity Measurement (Column).......................................... 55 Figure 4.9: Photograph of Column Setup ....................................................................... 56 Figure 4.10: Hydraulic Conductivity Measurement (Tracer) ......................................... 57 Figure 4.11: Porosity Measurement................................................................................. 58 Figure 5.1: Conceptual Diagram of Contaminant Sources ............................................. 60 Figure 5.2: Contaminant Sources and Two-dimensional Approximation ...................... 61 xiii Figure 5.3: Beach Face Infiltration with Open Well ...................................................... 64 Figure 5.4: Beach Face Infiltration ................................................................................. 66 Figure 5.5: Sinking Plume Dataset A-1 (Stable Plume) .................................................. 70 Figure 5.6: Sinking Plume Dataset A-2 (Stable Plume) ................................................. 72 Figure 5.7: Sinking Plume Dataset A-3 (Stable Plume) ................................................. 74 Figure 5.8: Sinking Plume Dataset A-4 (Unstable Plume) ............................................. 75 Figure 5.9: Sinking Plume Dataset A-5 (Unstable Plume) ............................................. 77 Figure 5.10: Sinking Plume Dataset A-6 (Unstable Plume) ........................................... 79 Figure 5.11: SEAWAT Simulation of Stable Plume ...................................................... 83 Figure 5.12: SEAWAT Simulation of Unstable Plume .................................................. 84 Figure 5.13: Sensitivity Analysis: Head Gradient .......................................................... 87 Figure 5.14: Sensitivity Analysis: Source Density ......................................................... 88 Figure 5.15: Sensitivity Analysis: Loading Rate ............................................................ 89 Figure 5.16: Sensitivity Analysis: Dispersivity .............................................................. 90 NOMENCLATURE H f = Total Freshwater Hydraulic Head (L) H = Total Hydraulic Head (L) h f = Freshwater Pressure Head (L) h s = Saltwater Pressure Head (L) ? f = ? 0 = Freshwater Density (ML -3 ) ? s = Saltwater Density (ML -3 ) P = Fluid Pore Pressure (M -1 LT 2 ) C = Solute Concentration (ML -3 ) ? = Density ? Gravity (ML -1 T -1 ) Z = Reference Elevation (L) Q = Flowrate (L 3 T -1 ) Q L = Loading Rate (L 3 T -1 ) q = Groundwater Flux (L 3 L -2 T -1 ) V = Average Porewater Velocity K sat = Saturated Hydraulic Conductivity (LT -1 ) k = Intrinsic Permeability (L 2 ) ? = Porosity ? = Viscosity (FTL -2 ) ? L = Longitudinal Dispersivity (L) ? T = Transverse Dispersivity (L) D hl = Hydrodynamic Dispersion (L 2 T -1 ) D d = Effective Molecular Diffusion (L 2 T -1 ) ? = Compressibility of a Bulk Porous Material (M -1 LT 2 ) ? = Coefficient of Water Compressibility (M -1 LT 2 ) S p = Specific Storage in terms of Pressure (M -1 LT 2 ) S f = Specific Storage in terms of Freshwater Head (L -1 ) ? = Vector Operator (L -1 ) g = Acceleration Due to Gravity (LT -2 ) H sat = Height of the Saturated Zone (L) H dp = Penetration Depth (L) xiv 1 CHAPTER 1 INTRODUCTION 1.1 Background A large percentage of the world?s population is dependent upon groundwater resources and utilizes them via public water supply or private domestic wells. However, groundwater resources are presently insufficient for human needs in many regions of the world due to various contamination problems. This is particularly true in coastal regions where natural seawater intrusion processes and catastrophic storm surges may disperse large volumes of saline water into the unconfined coastal aquifers. Saltwater intrusion in the subsurface is the most prevalent source of contamination to coastal aquifers, but it is not the only source. Other sources of contamination include leachate from municipal landfills, leaky underground storage tanks that store saline water (Hanf et al. 2005), contamination from natural brine formations (Penny et al. 2003), and surface saltwater deposition caused by tsunamis and hurricanes (Illangasekare et al. 2006). Therefore, a thorough understanding of the interaction of saltwater with fresh groundwater under various source conditions is needed to effectively manage these contaminant problems. Saltwater contamination in the saturated zone often propagates from surface sources located above the vadose zone. Once the contaminants reach the saturated zone, the contaminants migrate with the ambient groundwater to form a plume with lower contaminant concentrations at the fringes and higher contaminant concentrations within the core. The shape and contaminant concentration distribution of the plume depends on many site specific soil and fluid properties such as the density differences between the contaminant solution and the ambient groundwater. Fluid flow problems in which density influences the flow pattern are commonly referred to as density-driven flow (DDF) problems (Holzbecher 1998). DDF problems can be categorized by the location of the denser fluid. A DDF problem in which a dense fluid is situated below a less dense fluid is described as having a stable density configuration. These problems are referred as stable-interface problems. Likewise, a problem in which the dense fluid is situated above the less dense fluid has an unstable density configuration and is referred to as an unstable-interface problem. Examples of a stable-interface problem are illustrated in Figure 1.1. The left-most example in the figure shows a laboratory simulation of a classic saltwater intrusion problem in which a dense saltwater wedge has intruded beneath the less dense ambient fresh groundwater. The second example in the figure (middle) gives another example of a stable-interface problem in which the less dense fluid (dark colored) is situated above a more dense fluid (light colored). FW SW 2 Figure 1.1: Examples of stable interface problems An example of an unstable-interface problem, in which the more dense fluid is introduced above the less dense fluid, is shown in Figure 1.2. In this example the interface of the dense fluid (dark colored) tends to sink rapidly and mix with the ambient fresh groundwater. FW SW Figure 1.2: Example of an unstable-interface problem Within the category of unstable interface problems, the behavior of the contaminant plume can be classified as stable or unstable. This plume classification refers to the transport behavior of the plume with respect to the mixing and sinking patterns of the solute with the ambient fresh groundwater. An example of a stable and an unstable plume is shown in Figure 1.3. A stable plume is formed under high groundwater flow conditions when the plume forms a clear stable interface as shown in the figure. An unstable plume will generally propagate under low flow conditions which can cause the plume to penetrate deeper into the aquifer. Unstable plumes have the potential to contaminate larger regions of an aquifer and can mix with the ambient water to a larger degree (Oostrom et al. 1992). Plume stability is influenced by several variables such as the horizontal darcy flux (q x ), loading rate (Q L ), longitudinal (? L ) and transverse (? T ) 3 dispersivities, and the density difference (??) between the contaminant source and the ambient water. Figure 1.3: Plume stability in unstable interface problems: Left: stable plume. Right: unstable plume (Freeze, Cherry 1979). Plume instability can also be influenced by the interfacial disturbances caused by porous medium heterogeneities. For density driven flow problems involving immiscible fluids, plume instabilities are manifested by lobe shaped perturbances which first form at the edge of the plume and then within the plume (Schincariol et al. 1994). In miscible plumes, dispersion tends to smooth out the fingers before they have a chance to grow. The mixing zone reduces the fingering effect, which changes the fluid properties and may contribute to the development of plume instability. The stratification of freshwater above seawater in coastal environments is a commonly encountered stable-interface DDF problem. The mixing of 3% to 4% of seawater is enough to render freshwater unfit for most uses and may have a serious environmental impact (Custodio 2002). The United States Environmental Protection Agency (USEPA) has set a secondary maximum contamination level (MCL) for many of the constituents existing in seawater. The MCL of chloride, the most abundant constituent in seawater, is set at 250 mg/L which is less than one fiftieth of the chloride concentration of seawater, which is 19000 mg/L. Ingestion of groundwater with high chloride concentrations has been reported to cause vomiting, diarrhea, prostration, and 4 5 dehydration. The higher chloride concentrations can also have adverse effects on agriculture as many plant species are intolerant to higher salinity levels. In addition, many developing countries are dependent on tourism and fishing industries in coastal regions to provide income for the local people. Contamination of fresh groundwater resources may force these industries to relocate which can result in periods of economic downturn for the region. Therefore, the quality of life in many coastal regions is extremely sensitive to the quality of the groundwater available in local aquifers. A typical cross section of a coastal aquifer is given in Figure 1.4. This conceptual figure illustrates a practical stable interface problem where the less dense freshwater is situated above the denser saltwater intruded from the ocean. As shown in the figure, coastal aquifers are connected to the sea through the subsurface forming a freshwater- saltwater interface, also known as the transition zone. A density difference between the fresh ground water and the seawater is the primary cause for the intrusion of a saltwater wedge beneath the fresh groundwater. The position of the saltwater wedge can move laterally in the subsurface depending on the total head gradient. A reduction of freshwater flux, due to groundwater pumping or low recharge reduces the gradient and induces the saltwater wedge to migrate laterally inland. Once the lateral movement of the wedge is counteracted by the freshwater discharge the saltwater wedge can attain a state of equilibrium, which is called a steady-state saltwater wedge. In a theoretical sense, the concentration profile of the steady-state salt wedge in the subsurface does not change with time. However, the transport within the wedge is dynamic and there will always be a dispersion zone (or transition zone) where the salt is carried away by the freshwater flowing along the interface. Figure 1.4: Cross-sectional view of a coastal aquifer illustrating a saltwater intrusion wedge (Barlow 2003). Large-scale saltwater inundation events, such as tsunamis and storm surges, represent another potential saltwater contaminant source to coastal aquifers. These types of problems involving the deposition of dense seawater on the land surfaces above the fresh groundwater are classified as unstable interface problems. Typically, tsunamis occur as a result of energy transfer due to tectonic movements of the plates located beneath the ocean which can suddenly alter the elevation of the sea floor. This sudden shift causes a seismic water wave to spread producing tsunamis. Storm surges are the result of wind stresses or atmospheric low pressure zones (hurricanes) near the coast. The effects of storm surges can be amplified during the high tide. Tsunamis and storm surges can cause elevated sea levels which if high enough can deposit large volumes of seawater into coastal lagoons and depressions overlying the shallow unconfined aquifers. The amount of saltwater deposited over the coastal aquifer is controlled by two main parameters, inundation distance and run-up level. The inundation distance is defined as the maximum horizontal penetration of the tsunami from the shoreline. The run-up level is defined as the difference between the elevation of maximum tsunami penetration (inundation line) and the sea level at the time of the tsunami as described in Figure 1.5. 6 Figure 1.5: Diagram illustrating the inundation distance and run-up elevation (UNESCO-IOC 2006) 1.2 Sri Lanka Tsunami Scenario Groundwater contamination due to inundation events has received considerable attention in recent months because of the 2004 Asian tsunami event. The Sumatra- Andaman earthquake, which occurred on December 26, 2004 in the Indian Ocean, measured 9.3 on the Richter scale. The earthquake produced an ocean-wide tsunami that devastated many countries including Indonesia, Sri Lanka, India, and Thailand. An estimated 280,000 people were killed and more than 1.5 million were displaced (Illangasekare et al. 2006). This disaster greatly impacted the infrastructure, sanitation, economy, morale, and overall quality of life in the region. It was considered the worst tsunami catastrophe in history according to a report published by UNESCO (UNESCO- IOC 2006). More than a year after the tsunami, the groundwater along the eastern coast of Sri Lanka was still inadequate for human consumption. Groundwater serves as the primary water source for many of the inhabitants in the affected regions who access fresh water through open wells which were hand dug in the shallow unconfined aquifer. Due to the high inundation levels, large quantities of saltwater were directly deposited into the open 7 dug wells. Saltwater was also directly injected into the aquifer due to the residence time of the tsunami over the highly permeable sandy soil. The devastation of the tsunami at the coastal town of Maruthumunai on the eastern coast of Sri Lanka and the high density of open dug wells (concrete cylinders and open excavations) is shown in Figure 1.6. Reports indicated that the water level at this site was approximately 7-9 meters (Liu et al. 2005) at this site which is higher than the top of the palm trees shown in the figure. 8 Indian Ocean Open Dug Wells Figure 1.6 ? Picture illustrating the high density of open wells in the coastal town of Maruthumunai, Sri Lanka Currently, there are no guidelines available for the decontamination of aquifers and the open dug wells that have been contaminated by saltwater. One strategy is to pump the contaminated water out of the wells and the aquifer. However, any increase in the drawdown of the aquifer due to pumping would increase the likelihood of contamination due to saltwater upconing. Another possible decontamination strategy is to allow the aquifer to remediate itself naturally. The latter approach appears to be a more pragmatic solution because the hydraulic conductivity is relatively high in coastal sandy aquifers (10-500 m/day) and Sri Lanka experiences an annual recharge of 1000 ? 1700 mm, which is concentrated during the monsoon season (Villholth et al. 2005). However, more 9 information is required to accurately predict the length of time for the natural remediation of the contaminated aquifers. Massive efforts were undertaken in the short period of time following the tsunami event to restore the water quality to pre-tsunami conditions. However, many of these efforts proved to be unsuccessful. For example, pumped wells were prone to collapse and pumped contaminant water was often dumped on the lands surface in the vicinity of the well where it re-infiltrated the open dug wells. For reasons such as these, a team of scientists headed by the International Water Management Institute (IWMI) completed a detailed impact assessment study to develop a better remediation strategy. Villholth et al. (2005) provided field data and professional insight into the problems incurred as a result of the tsunami and its aftermath. Prior to this study, there was limited information available about the specific site details of this large-scale contaminant problem. More recently, a field monitoring team was setup on the east coast of Sri Lanka in three towns named Kallady, Kaluthavalai, and Oluvil in the Batticaloa and Ampara districts. They measured salinity, groundwater level, turbidity, etc. on a regular basis (20 to 40 day interval). One hundred and fifty wells were monitored within a 2 km distance from the coastline. These wells included both affected and unaffected areas. The salinity level in the wells had decreased significantly by the start of the study which was expected because of dilution caused by the seasonal monsoon. However, in the months following the monsoon the salinity levels decreased slowly (Illangasekare et al. 2006). It was discovered that a large fraction of both un-flooded and flooded wells were unsuitable for drinking at the end of the study in mid-July of 2005. 10 The coastal aquifers of Sri Lanka consist of structureless sand, ranging from fine to moderately coarse. The aquifers stretch on average 2 ? 8 km inland (Villholth et al. 2005). There are also many lagoons in these coastal regions which are connected to the top of the coastal aquifer. The lagoons and wells become sources of saltwater contamination for a period of time following inundation. The vulnerability of these aquifer systems to contamination is high due to the highly permeable porous media which allows for the rapid infiltration of contaminants. Furthermore, the aquifers are shallow and unconfined with little retention capacity. A conceptual physical study of groundwater contamination due to a tsunami type event is lacking in the literature. 1.3 Research Goals and Objectives The main goal of this study was to develop a laboratory-scale physical model to simulate the migration of a dense miscible plume in a coastal aquifer. The study was performed in a thin rectangular tank which represented a two-dimensional strip of an unconfined coastal aquifer. The experiments were designed to give a conceptual basis to the behavior of dense miscible plumes with regard to sinking and mixing patterns within a sandy coastal aquifer. Furthermore, the interaction of a saltwater wedge and a solute plume were observed simultaneously within the model to better understand transport conditions under the field setting. This study began with the construction of a flow container with in-situ flow visualization capability. The experimental setup included saltwater and freshwater boundary conditions. The development of the experimental setup was an iterative procedure that evolved over the course of this two year study. 11 Upon completion of the physical model, many trial experiments were performed. The outcomes of three key experiments are presented in this thesis. The first study involved the combined effects of saltwater infiltration along the beach face, in addition to saltwater injection through an open well with a salt wedge present. The second study qualitatively modeled saltwater infiltration along the beach face in the presence of a saltwater wedge. Understanding this source of contamination is important since so many of the tsunami- affected residents had freshwater supply wells within a couple of hundred meters from the coast line. The final dataset qualitatively models the behavior of a dense solute plume propagating from a salt pond and infiltrating into a shallow coastal aquifer. The head gradient of the system was altered to observe the effects on plume behavior. The physical data was compared with the solutions of the numerical code SEAWAT 2000. While density-driven flow problems have received some attention in past years, a conceptual physical study of a tsunami type event has not been performed. 1.4 Organization of Thesis The goal of the thesis is to investigate the transport of saltwater in a coastal aquifer using a laboratory scale physical model. A thorough literature review was performed concerning physical modeling of density-driven flow and related topics. The literature review is presented in Chapter 2. Chapter 3 describes the experimental setup of the physical model. The design of the flow container and the measurement of certain flow parameters are described in detail. Chapter 4 presents a series of physical and numerical simulations of various types of saltwater inundation events in a laboratory model. In Chapter 5, a summary and suggestions for future research are discussed. 12 CHAPTER 2 LITERATURE REVIEW 2.1 Introduction Groundwater represents 98% of the freshwater resources available for human use (Fetter 1988). As the population and economic development in many regions of the world continues to increase the demand for groundwater resources is expected to increase. The dependence on groundwater is of special concern for the inhabitants of coastal regions where the demand for freshwater and the risks related to the saline contamination of coastal aquifers are high. These coastal aquifers can often serve as a sustainable source of fresh water if managed and used according to recharge and local hydrogeological characteristics. Otherwise, remediation of saltwater contaminated aquifers may prove to be a costly and an extremely difficult problem to correct as evidenced by several regions in Florida, USA (Barlow 2003). An important first step in coastal aquifer management is to obtain a representative model of the coastal aquifer system to model the flow and salinity transport (Custodio 2002). Scientists and engineers have developed several types of groundwater models for predicting the impact of groundwater pumping and the fate of groundwater contaminants in the subsurface. These models are necessary tools to help manage the resource for 13 future generations. The types of groundwater models are described in the following section. 2.2 Groundwater Models Definition A model is defined as a representation of an object, system, or idea in some form other than the entity itself (Holzbecher 1998). Groundwater models are the most widely used tools for understanding groundwater flow behavior and contaminant transport in the subsurface. Groundwater flow models are always an incomplete representation of a real system. The reliability of predictions from a groundwater model depends on how well the model approximates the field situation (Wang and Anderson 1982). The three types of groundwater models which will be discussed in this section are analytical, numerical, and physical models. 2.2.1 Analytical Models Analytical models provide the exact mathematical solution to the groundwater problem but in order to obtain this solution a series of significant simplifying assumptions must be made. In general, obtaining the exact analytical solution to the partial differential equations requires that the properties and boundaries of the flow system be highly and perhaps unrealistically idealized. For most field problems, the mathematical benefits of obtaining an exact analytical solution are probably outweighed by the errors introduced by the required simplifying assumptions of the complex field environment (Konikow 1996). 14 Analytical models cannot model complex contaminant transport in general groundwater flow problems governed by advection, diffusion, dispersion, decay, and sorption reactions. For this type of problem, which requires the flow and transport equations to be solved simultaneously, an iterative approach is needed to couple the governing equations. Density coupled flow problems are another example of a problem in which a scheme is required to couple the flow and transport equations. The traditional verification procedure by use of analytical solutions is generally not applicable due to the nonlinear nature of density driven flow problems (Diersch and Kolditz 2002). Hence, sophisticated numerical techniques are required to solve these coupled problems. 2.2.2 Numerical Models Numerical models are capable of solving large systems of equations, nonlinearities, and complex geometries that are often impossible to solve analytically (Chapra and Canale 2002). This is beneficial for solving groundwater problems which typically involve a set of partial differential equations (PDE) with many variables. However, a drawback of numerical models is that numerical errors are introduced during the solution procedure. These errors are classified into two categories: round-off error and truncation error. Round-off errors occur because computers can only use a finite number of digits to represent a number. Truncation errors occur because the numerical method may employ approximations to represent exact mathematical operations and quantities (Chapra and Canale 2002). Density driven groundwater flow problems are greatly affected by heterogeneities in the subsurface which introduces randomness. While numerical errors inherently 15 introduce some degree of randomness, it has been shown that these numerical errors are not representative of randomness caused by field heterogeneities (Schincariol et al. 1994). Finer meshes and smaller time steps can be incorporated into models to reduce the effects of the numerical errors. However, the question which arises is what level of error in the numerical scheme can be tolerated (Diersch and Kolditz 2002). A numerical solution that has not been calibrated against a known benchmark solution may be significantly different than what transpires in the actual field setting. This presents a challenging problem for scientists and engineers responsible for maintaining a required level of confidence in model predictions. To reduce the uncertainty associated with numerical solutions, several benchmark solutions have been developed for testing density-driven flow problems. These benchmark problems include the Henry Problem (Henry 1964), Elder Problem (Elder 1967), and Hydrocoin Problem (Konikow et al. 1997). However, each of these solutions has drawbacks. The Henry problem is a semi-analytical solution which has not been reproduced in a physical model (Simpson and Clement 2003). The Elder problem is an experimental heat transport problem which was performed in a Hele-Shaw cell. It has since been modified into a solute transport problem. It is extremely difficult to numerically reproduce the Elder solution due to the sensitivity of the problem to model parameters and mesh discretization level (Diersch and Kolditz 2002). The Elder problem for solute transport has yet to be reproduced in a laboratory setting. Finally, the Hydrocoin problem is a theoretical problem that is only applicable for inter-code testing. Details of these benchmark problems, as described in the SEAWAT 2000 User?s Guide (Guo and Langevin 2002), are summarized in Figures 2.1, 2.2, and 2.3. Henry Saltwater Intrusion Problem Figure 2.1 ? Solution to the Henry?s steady state saltwater intrusion problem (Guo and Langevin 2002) 16 Elder Heat Problem Figure 2.2 ? Solution to Elder?s transient heat problem (Guo and Langevin 2002) 17 Hydrocoin Sloping Pressure Boundary Problem Figure 2.3 ? Solution to the Hydrocoin steady state saltwater intrusion problem with a sloping boundary (Guo and Langevin 2002) Comparison of numerical solutions to the accepted benchmark solutions similar to those presented above is referred to as numerical code verification (Simpson and Clement 2003). Numerical code verification is very different from numerical code validation which involves a comparison of the numerical solutions with lab or field data (Holzbecher 1998). 18 19 2.2.2.1 Review of Density Coupled Flow and Transport Codes To represent a problem within a computational framework the physical representation to be modeled must be expressed in terms of the appropriate mathematical equations. Presently, there are many numerical tools available to predict the movement of contaminant plumes in the subsurface. These numerical tools solve the partial differential equations using various schemes such as finite element (SUTRA) or finite difference (HST3D, SEAWAT 2000) method. These tools are perhaps our best means for understanding how a contaminant plume will behave under various flow conditions. HST3D Model HST3D, Heat and Solute Transport Program, is a modified version of the parent code SWIP, Survey Waste Injection Program, which was developed to calculate the effects of liquid-waste disposal into deep saline aquifers. The code has the capability to simulate density and viscosity effects on groundwater flow due to solute and heat transport processes. As listed in the HST3D manual (Kipp 1997), this code may be used to simulate waste injection into saline aquifers, landfill-contaminant movement, seawater intrusion, brine disposal, heat storage in aquifers, and similar transport situations. HST3D uses a finite difference scheme which can be computationally expensive when attempting to capture the regions of low longitudinal and transverse dispersivity (Kipp 1997). For this scenario a fine grid is required in order to capture those effects which increase the computer-storage and computation-time requirements significantly. Another limitation of HST3D is that solutions of a scenario with a high viscosity or density contrast will vary from experimental results. Typically, the high viscosity contrast simulation will result in a fingering effect which is less pronounced than what is 20 actually observed in the lab. Similarly, for a high density contrast scenario in which a dense fluid overlies a less dense liquid the magnitude of the perturbations is less than what actually occurs. For a detailed development of the governing equations and a list of the assumptions incorporated into HST3D, refer to the HST3D user?s manual (Kipp 1997). SUTRA Model Solute transport simulation with SUTRA may be used for modeling variable density leachate movement, and for cross-sectional modeling of salt-water intrusion in aquifers at near-well or regional scales, with either dispersed or relatively sharp transition zones between fresh and saline waters (Voss 1984). SUTRA employs a finite element scheme to approximate the governing equations of a saturated-unsaturated transport model. SUTRA includes unsaturated flow and transport processes but is primarily intended for saturated systems as the numerical algorithm is not specialized for the non-linearities of unsaturated flow. SUTRA has the capability to model the two independent processes that include groundwater flow coupled with solute transport or flow coupled with thermal energy transport. Within the solute transport module, SUTRA can be used to model the effects of sorption and contaminant decay. SUTRA has the capability to model either in the cross-sectional or in the areal domain. SEAWAT 2000 Model The numerical model used in this study was SEAWAT 2000. SEAWAT 2000 (Langevin et al. 2003) is a widely used program that was developed to simulate three- dimensional, variable density, transient ground-water flow in porous media. The source code for SEAWAT was developed by combining a density modified version of 21 MODFLOW and MT3DMS codes into a single tool. The SEAWAT code was tested by simulating five benchmark problems involving density-driven groundwater flow. SEAWAT can be used in the fully-implicit finite difference mode. The governing equations used by SEAWAT are presented in Chapter 3. 2.2.3 Physical Models Since their development, numerical models have been used as the primary tool for predicting saltwater transport in density-driven flow scenarios. However, much has been gained from physical model studies with regards to conceptualization of a groundwater system. Physical studies have increased the importance of visualizing density-driven groundwater flow studies. Furthermore, it is a costly and difficult task to attain quality field scale data. For this reason physical laboratory scale data presents an attractive alternative to field scale studies and it can aid numerical modelers with code validation. Observation of a groundwater problem can be performed either in the field or in the laboratory. However, the laboratory setting provides an ideal setup for measurement of flow and transport variables. Also, the level of control in small scale laboratory experiments allows the simplifying assumptions to be made with greater confidence than at the field scale. Therefore, laboratory data is preferred in many instances to avoid the uncertainty of data collected in the field. Furthermore, laboratory experiments tend to be less time and resources demanding than field experiments. Physical modeling is of particular importance when analyzing coupled flow and transport groundwater problems that are difficult to model numerically. For example, the problem of plume stability in a density-driven flow (DDF) problem may be extremely 22 sensitive to the small perturbations which are introduced by numerical errors. For such complex problems, physical modeling is a necessary tool. 2.2.3.1 Review of Physical Model Studies Density-driven groundwater flow has been studied by groundwater researchers since Ghyben (1889) and Herzberg (1901) independently discovered that saltwater exists in the subsurface. Using a simple hydrostatic model they estimated that for every one foot of freshwater above sea level in the subsurface there exists 40 feet of freshwater below sea level (Fetter 1988). Since this discovery many researchers have attempted to simulate this variable density flow phenomena in the lab. Some of these works are summarized in the following paragraphs. One of the first studies investigating the mixing of two stratified fluids with different densities in horizontal motion was performed by List (1965). The purpose of List?s study was to observe the mixing and the interface behavior between the two fluids when the heavier fluid was situated above the lighter fluid. A ?sandbox? flow model of dimensions 250 cm long, 15 cm wide, and 35 cm deep was used for the laboratory experiments. An Ottawa Flint Sand with a mean diameter of 0.530 mm was used as the porous medium to simulate an isotropic and homogeneous aquifer. List discovered that if there is a dense fluid flowing horizontally above a less dense fluid then the flow would always be unstable. However, there exists a scenario in which the growth rate of the unstable waves was low enough such that the flow could be considered quasi-stable. List describes the mixing of fluids in a porous medium to be a combination of molecular diffusion processes inside the pore spaces and a macro-scale convective 23 mixing caused by the arbitrary flow through the porous medium. Since the convection caused by the flow was relatively much higher than the pore scale diffusion, the problem is referred to as a convection dominated problem. A hydrodynamic dispersion term is used to represent the cumulative effects of molecular diffusion and convection processes. List concluded that as the longitudinal velocity increases, the growth rate of unstable waves decreased. This is correct because an increase in velocity increases the rate of both longitudinal and lateral dispersion which will overcome the ability of the light fluid to enter the heavy fluid and vice versa. The instability growth rates will be high if the plume descends faster than it disperses. However, if the dispersion velocity is higher than the velocity of descent, the plume should be more stable. List observed density difference (??) up to 1% to have no affect on the lateral dispersion coefficient. Oostrom (1991) conducted another experimental study of dense aqueous phase leachate plume behavior which was conducted to observe the effects of altering the relative density difference, Darcy flux, and loading or leakage rate in a homogeneous unconfined aquifer model (Oostrom et al. 1992). Two ?sandbox? models referred to as flow containers A and B were used for flow visualization. The dimensions of the flow containers A and B were 205 ? 100 ? 7.8 (cm) and 167.5 ? 100 ? 5 (cm), respectively. The experiments conducted in flow container A were preliminary experiments for flow container B. Flow container B was more advanced because in addition to flow visualization there was a dual-energy gamma radiation technique that could be used to study changes in density. The study simulated leakage from a coastal lagoon situated above an unconfined aquifer. The porous media consisted of flintshot Ottawa sand with a 90% grain size diameter distribution between 0.25 mm and 1 mm. The hydraulic conductivity in flow container A and flow container B was found to be 62 and 36 m/day, respectively. The difference in the hydraulic conductivities was due to packing. A NaI and NaBr- solution was used to simulate the dense lagoon liquid. Fluorescent dye was mixed with the NaI and NaBr- to serve as an optical tracer in the experiments. NaBr- was the preferred tracer because it did not dull the fluorescent tracer. Each set of experiments consisted of the development of a neutral plume, the displacement of a neutral plume by a dense plume, and the displacement of a dense plume by a neutral plume. The neutral plume was released at a constant loading rate using a volume displacement peristaltic pump. The Darcy flux was altered to observe the effects on plume stability. It was observed that a large portion of the neutral plume was transported in the unsaturated zone. One experiment showed that the entire plume was transported in the capillary fringe. There was interestingly some transport in the uppermost layer which was dry due to evaporation. The neutral plume penetration depth (H dp ) could be predicted based on the following equation: sat x L dp H Q Q H ? ? ? ? ? ? ? ? ? = (2.1) H sat = the height of the saturated zone at x = 1.025 m (half-way point) * This equation considers advection only. The neutral plumes were then displaced by the dense plumes in which case the plume began to develop lobe shaped instabilities. The plume penetration increased as expected. This shows that an increase in density leads to an increase in the level of instabilities which increases the potential to contaminate larger regions of the aquifer. A noticeable 24 difference in plume behavior was observed on the front and rear sides of the tank. This suggests that perhaps the flow container was not level. Furthermore, the ambient groundwater was able to move around the instabilities which suggests that this was a quasi ? three dimensional study even though the tank width was relatively small. A non- dimensional ? 1 expression was developed in order to characterize the stability of a dense plume in a homogeneous porous medium with a constant loading rate. The ? 1 number represents a ratio of the density-driven flux to the horizontal Darcy flux as shown below: xsat qK / 0 1 ? ? ? ? ? ? ? ?? = ? ? ? (2.2) It was determined that all plumes with a ? 1 greater than 0.5 were unstable. This is a powerful ratio for conceptualizing the stability of dense plumes. An example presented in Oostrom (1991) states that a 10 g/L contaminant solution in a porous medium with a horizontal darcy flux, q x, of 1 m/day will produce an unstable plume. It can then be predicted that a 2 g/L contaminant solution in a porous medium with a flux of 0.2 m/day will behave similarly since the ratio is the same. The phenomenon of plume transport in the unsaturated zone was further investigated. Concentration was plotted as a function of time to obtain the breakthrough curves. Using a curve fitting technique, the coefficients of longitudinal hydrodynamic dispersion, D hl , were determined. The longitudinal dispersivity could then be determined from the following relationship: () V DD dhl L ? =? (2.3) Longitudinal dispersivity values, ? L , were found to be relatively small which indicates that advection transport was dominant. Molecular diffusion was responsible for less than 25 26 3 % of the longitudinal hydrodynamic mixing. Due to the small ? T values it was assumed that there was little mixing in the vertical direction. Molecular diffusion was found to account for 7.8% of the vertical mixing. It was observed that when the loading rate, Q L , was increased there was an observed decrease in the darcy flux, q x . This can be explained by the mounding that resulted as a result of the increased loading rate. Concentrations were observed using a gamma radiation technique which averaged the concentration over the whole thickness. This means that the actual mixing was probably overestimated when there was a lot of shear flow. Flow container B had the similar problem of container A in that the plume was concentrated against the wall of the container. A C/C 0 value was only observed near the very top of the container. The deeper the plume penetrated the lower the concentration became due to increased mixing. A similar plume was observed for plumes with similar ? 1 but only 50 % of the q x and Q L values. It was observed that the neutral plume tended to travel slightly faster than that of the stable dense plume. Instabilities in the dense leachate plume are more likely to develop in the lower part of the plume. When higher flow rates were used, the occurrence of instabilities was prevented. For all plumes the concentration gradient was steepest near the bottom. For stable plumes the concentration went from 1 to 0 in a matter of a few centimeters. Unstable plumes tended to have an initial deep penetration that was then stabilized by the freshwater flux. Over time the magnitude of instabilities decreased. Longitudinal dispersivities increased with distance from the source. In the capillary fringe the contaminants nearly moved with the same velocity as below the water table. Contaminant transport in the unsaturated zone was dominated by molecular diffusion vertically upwards from the saturated zone. Instabilities were generally large at first and then became smaller with time. Figure 2.4: Diagram of flow container (Oostrom 1991) Hayworth (1993) conducted both a physical and numerical study investigating the three-dimensional behavior of dense aqueous plumes. Dense plumes tended to become unstable under certain flow conditions. Unstable plume behavior can affect the trajectory and shape of the plume. An unstable plume also has more surface area which means that there is a large mixing zone with the ambient groundwater. Two flow containers were constructed for simulating homogeneous and isotropic phreatic aquifers. They are referred to as flow container A and flow container D with dimensions of 80 ? 40 ? 5 (cm) and 200 ? 100 ? 30 (cm), respectively. Flow visualization and direct fluid sampling techniques were used. Deionized (Millipore) water was used for the groundwater flow. Deionization produces water low in ionic, organic, and particulate contamination. Two different porous media were used in these 27 28 experiments. An Ottawa sand and industrial grade class beads (Potters, A-120). In flow container A, K sat was determined in situ to be 798 m/day and 62 m/day for the glass beads and the sand, respectively. In tank B, the K sat value was only determined for flow container D (K sat = 30.7 m/day). The difference in K sat from flow container A to flow container D was due to packing. Packing flow container D took an estimated 30 hours. A 3 cm and 12.5 cm capillary zone was observed for the glass beads and the sand, respectively. Darcy velocities ranged from about 11 to 20 m/day. Density differences ranged from 0.01 to 0.037. The dense solution consisted of fluorescin, NaBr, and NaN 3 . The sodium azide was used to kill bacteria. Unstable and stable dense plumes were observed. Hayworth also observed a plume referred to as quasi-stable according to List (1965). A quasi stable plume is characterized by the development of waves in the near source region of the plume along the plume- groundwater interface. These waves move through out the medium without significant growth or decay. The classification of quasi-stable suggests that classification of plumes with a ? 1 value near the transition range is subjective. There is no reason to assume that a distinct ? 1 number can be determined because there are many undefined variables. For example, small changes in fluid viscosity or small changes in fluid viscosity or heterogeneities may act to initiate unstable plume behavior. An unstable plume was observed for ? 1 values greater than 0.5 for the point source and for values greater than 0.3 for the line source. The point source plumes tended to remain more stable at higher values of ? 1 because of a destabilizing effect of forced penetration of ambient groundwater into the plume. For dense plumes emanating from a point source the ambient groundwater flow was allowed to move around the sinking body of the plume and was not forced to move through the plume. Thus pathways through the plume did not become desirable pathways until the ? 1 value was increased to 0.5. Rayleigh numbers were determined to not be the best way to determine plume stability as the number requires a value for transverse dispersivity which is difficult to determine. T satdp D KH Ra ? ?? )/( 0* ? = (2.4) The problem described in (Shackelford 1991) work involved wastes, consisting of coal ash and flue gas desulfurization sludge, which are typically disposed of in ponds or landfills. The experimentation involved a density-driven sinking plume experiment without a horizontal flow in a laboratory model. The model was ?? thick, filled with spherical glass beads (potters A-120 bead). Shackelford observed that different volumes of injected contaminant solution did not affect the rates of descent. However the smaller volume did not become unstable at the interface. Density difference however is significant with respect to the rate of descent. Plume size increased with increasing density due to larger settling velocity. Fingering was more prevalent at smaller densities. Fingers and their paths seemed random in nature, the plumes exhibited no symmetry. For a given density difference there is some small volume of plume that will be stable (Shackelford 1991). Shackelford did not observe a relationship between the volume injected and the number of fingers that developed. 29 Figure 2.5: Shackelford fingering experiment Paschke and Hoopes (1984) conducted a study in which a laboratory sand tank was used for validation purposes of an analytical model which predicts the trajectory, concentration, and boundary of the influenced contaminant plume. The model involved hydrodynamic dispersion, buoyancy induced vertical motion, and the ambient groundwater flow. Steady flow was applied in an isotropic, homogeneous, fully saturated aquifer of uniform, horizontal groundwater flow. The tank dimensions were 184 long ? 60 high ? 30 wide (cm). The porous medium consisted of a sand with a d 10 of 1.1 mm (d 10 = Hazen's effective grain size in mm, relative to which 10% of the sample is finer). The selected contaminant was a dissolved conservative tracer which is miscible with the ambient groundwater. The tracer NaNO 3 was preferred because it is less corrosive than NaCl; part of their tank was built with aluminum. The contaminant had the same viscosity but different density than the ambient groundwater. Lateral spreading due to density differences was assumed to be minimal compared to the lateral spreading due to 30 31 dispersion. Hydrodynamic dispersion in the x and y directions, D x and D y , were determined by fitting appropriate solutions of the mass conservation equation for a dilute, conservative substance. It was discovered that the plume would continue to move vertically downward and mix with the ambient groundwater until the plume becomes so dilute that downward movement becomes very small. Once this state is attained the flow is no longer density-driven and can be adequately modeled using the advection dispersion equations for non-buoyant fluids. Paschke and Hoopes (1984) were able to numerically model the movement of the groundwater plume with respect to concentration, trajectory, and boundary. Schincariol and Schwartz conducted a physical model study concerning plume stability in a homogeneous, layered, and a lenticular porous medium (Schincariol, Schwartz 1990). The study investigated the role of heterogeneity in hydraulic conductivity and instability development in promoting fluid mixing in density driven flow scenarios. The experiments were performed in a flow container of dimensions 106.8 ? 5 ? 71 (cm). Rhodamine was used as the optical tracer as its density was similar to that of the ambient water. For the homogeneous medium, it was observed that the highest concentrations of solute occurred in the lower portions of the plume, while the upper sections were more dilute. Increased plume penetration was observed in addition to longer and narrower lobes as plume concentrations were increased. Vertical mixing of the solute plume with the ambient groundwater was discovered to be dependent on the density of the source and the occurrence of instabilities. 32 Unstable convective mixing is an example of a nonlinear dynamic (chaotic) process (Scincariol and Schwartz 1990). Nonlinear systems are sensitive to small changes in the initial conditions. Schincariol and Schwartz (1990) found it difficult to replicate experiments due to the sensitivity of the system. However, a general observation of the experiments suggested that the plume profile in single density systems tended to have a steady state behavior near the source while the density driven plume problem tended to reach steady-state behavior further from the source, once the forces of free and forced convection dissipated each other. It was concluded that these systems become much more complex when layered and lenticular systems are introduced. Also, instabilities create a complex distribution of concentrations within the plume, with mass accumulation forming in the lower part of the plume. CHAPTER 3 DENSITY DRIVEN FLOW GOVERNING EQUATIONS The numerical code SEAWAT 2000 was used in this study to model the physical observations. A brief overview of the groundwater flow and solute transport equations as employed by the SEAWAT code are presented in this chapter. 3.1 Governing Equations For groundwater problems the potential energy can be quantified on a volume, mass, or weight basis. However, energy on a weight basis is the most commonly used for conceptual reasons because the units simplify to length. The quantification of energy on a weight basis is referred to as the hydraulic head. The hydraulic head is comprised of three components referred to as the pressure head, velocity head, and the elevation head. The hydraulic head is expressed mathematically as: Hydraulic Head = Pressure Head + Velocity Head + Elevation Head z g VP H ++= 2 2 ? (3.1) Typically, in groundwater problems the velocity head is relatively small and can be neglected. Hence, the hydraulic head in groundwater systems consists of a pressure head term and an elevation head term. The pressure head is defined as the amount of work required to bring a unit weight of water to a particular pressure. Likewise, the elevation 33 head is the amount of work required to bring a unit weight of water from a reference position to its actual elevation. 3.1.1 Equivalent Freshwater Head The first step in solving a density driven flow problem is to calculate the hydraulic head gradient in the system. However, the concept of hydraulic head is not straight forward for a scenario involving layered fluids with different densities. For example, in coastal aquifers the low density freshwater (density = ? f ) is underlain by a layer of high density saltwater (density = ? s ), as shown in Figure 3.1. The total hydraulic head (H f ) measured just above the interface would yield a different value than the total hydraulic head (H) measured just below the interface due to the density differences. To correct for this fact, SEAWAT uses an equivalent freshwater head approach. This approach helps to preserve continuity of the total hydraulic head across the interface. H f H Figure 3.1: A conceptual diagram of a coastal aquifer illustrating the layering of two fluids with different densities (Guo and Langevin 2002) The total hydraulic head, H f , measured at a location just above the freshwater/saltwater interface, as shown in Figure 3.1, can be expressed as the sum of the pressure head, h f , and the elevation, Z, above a reference datum: 34 35 H f = h f + Z (3.2) where, the freshwater pressure head, h f , can be expressed as: h f = P / (? f g) (3.3) Likewise the total hydraulic head, H, measured just below the interface can be expressed as: H = h+ Z (3.4) where, the pressure head, h, can be expressed as: h = P / (? s g) (3.5) Due to the different values of density (? f ? ? s ) there is a discrepancy in the value of total hydraulic head in the freshwater region, H f , and the value of total hydraulic head in the saltwater region, H. To correct for this inconsistency, the total hydraulic head measured in the freshwater and saltwater regions must be represented in terms of pressure, P: h f = P / (? f g) = H f ? Z (3.6) P = ? f g (H f ? Z) (3.7) Likewise, H = P / (? s g) = H ? Z (3.8) P = ? s g (H ? Z) (3.9) Setting (3.9) and (3.7) equal to each other yields: ? f g (H f ? Z) = ? s g (H ? Z) (3.10) Total hydraulic head, H f , in terms of total hydraulic head, H, can be expressed as: H f = (H? s ) / ? f - Z ( (? s - ? f ) / ? f ) (3.11) SEAWAT uses this equivalent freshwater head expression which is incorporated into the general form of the pressure-based groundwater flow equation. 3.1.2 Derivation of General Groundwater Flow Equation Groundwater flows through small pore size openings in the soil. However, the cross- sectional area through which the groundwater moves is often on a much larger scale. Therefore, groundwater flow is a microscopic process that must be scaled up to the macroscopic level. To accomplish this, a representative elementary volume (REV) is assumed. Predictions of groundwater flow through the subsurface can be attained by assuming conservation of mass, also referred to as the continuity principle, through a series of REVs. Essentially, an influx of mass into a REV must be equivalent to the sum of the outflux mass, changes in storage, and any source/sinks. Figure 3.2: Representative elementary volume (REV) In SEAWAT conservation of mass is mathematically expressed as: t qq s ? ? =+?? )( )( ?? ?? v =? source density (3.12) 36 The term on the right side of Equation 3.12 is equal to rate of mass accumulation stored in the REV. Using the chain rule to expand the rate of mass accumulation term, ttt ? ? + ? ? = ? ? ? ? ? ? ??)( (3.13) Assuming that the changes in porosity, ?, are restricted to those associated with the change of fluid pressure, P, it can be stated that: t P Pt ? ? ? ? = ? ? ?? (3.14) Assuming isothermal conditions the density is a function of fluid pore pressure and solute concentration, C: ),( CPf=? (3.15) Which can be mathematically represented as: t C Ct P Pt ? ? ? ? + ? ? ? ? = ? ? ??? (3.16) Substituting 3.14 and 3.16 into 3.13: t C Ct P Pt P Pt ? ? ? ? + ? ? ? ? + ? ? ? ? = ? ? ? ? ? ? ? ? ??)( (3.17) SEAWAT defines the coefficient of compressibility of a porous medium as: P? ? ? = ? ? ? )1( 1 (M -1 LT 2 ) (3.18) And the coefficient of compressibility of water is defined as: P? ? = ? ? ? 1 (M -1 LT 2 ) (3.19) Substituting equations 3.18 and 3.19 into 3.17: 37 [] t C Ct P t ? ? ? ? + ? ? +?= ? ? ? ?????? ?? )1( )( (3.20) Specific storage is defined as: [])1( ???? +?= p S (3.21) Substituting 3.21 into 3.20: t C Ct P S t p ? ? ? ? + ? ? = ? ? ? ?? ??)( (3.22) Substituting 3.22 into 3.12 yields the general form of the conservation of mass equation for the pressure-based variable-density ground-water flow in porous media: t C Ct P Sqq ps ? ? ? ? + ? ? =+??? ? ???? )( v (3.23) 3.1.3 Equivalent Freshwater Head Form of Groundwater Flow Equation SEAWAT modifies from a pressure based form to an equivalent head based form of the variable density ground-water flow equation. For flow in the horizontal plane the equation simplifies to: ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ?= g z Pk y P k x Pk q xz xy xx x ? ??? (3.24) Likewise, ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ?= g z P k y P k x P k q yzyyyx y ? ??? (3.25) In the vertical plane Darcy?s Law is expressed as, ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ?= g z Pk y P k x Pk q zz zy zx z ? ??? (3.26) 38 Substituting Equations 3.24, 3.25, and 3.26 into Equation 3.23: t C Ct P Sq kg z Pk y P k x Pk jg z P k y P k x P k ig z Pk y P k x Pk ps zz zy zx yzyyyx xz xy xx ? ? ? ? + ? ? =+ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ??? ? ??? ? ??? ? ? ??? ? ? ??? ? ? ? ? (3.27) The freshwater head can also be expressed as: g P zH f f ? += (3.28) Solving for P: ( )zHgP ff ?=? (3.29) Substituting (3.29) into equation 3.28, yields: () t C Ct H gSq k z H k y Hk x H k g j z Hk y Hk x Hk g i z H k y Hk x H k g f fps f f zz fzyf zx f f fyzfyyfyx f f f xz fxyf xx f ? ? ? ? + ? ? =+ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? ??? ? ???? ? ? ??? ?? ? ? ??? ?? ? ? ??? ?? ? 1 ? 1 ? 1 (3.30) The term S p g? = S f and f f gk K ? ? ?? ? = , substituting yields 39 () s f f f ffzzffzyffzx f ffyzffyyffyx f ffxzffxyffxx q t C Ct H S k z HK y HK x HK j z HK y HK x HK i z HK y HK x HK ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? ??? ? 1 ? 1 ? 1 (3.31) The term ? ? f can be considered equal to one. () s f f f f zz f zy f zx f f yz f yy f yx f f xz f xy f xx q t C Ct H S k z H K y H K x H K j z H K y H K x H K i z H K y H K x H K ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? ??? ? 1 ? 1 ? 1 (3.32) If we assume that the x, y and z axes are aligned with the principal permeability directions, then the above equation reduces to: () s f f f f zz f yy f xx q t C Ct H Sk z H Kj y H Ki x H K ? ? ? ? ? ? ? ? ? ? ? + ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ?+ ? ? ? ? ? ? ? ? ?+ ? ? ? ? ? ? ? ? ???? ? 1 ?? (3.33) Assuming homogeneous isotropic conditions for a two-dimensional groundwater problem the final equation can be expressed as: ( ) s f f ff x q t C Ct H Sk z H Kzi x H K ? ? ?? ? ? ? ? ? + ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ???? ?? (3.34) 40 3.1.4 Solute Transport Equation SEAWAT derives the solute transport equation by using the law of conservation of mass. SEAWAT begins the derivation by assuming a representative elementary volume (REV), shown in Figure 3.2. Concentration is a measure of the mass of a solute per unit volume. Therefore, the mass per unit volume of porous medium is equivalent to the product of porosity and concentration. Assuming porosity is a constant: x C x C ? ? = ? ? ? ? )( (3.35) For the transport of salt in an inert porous medium with negligible sorption capacity there are two dominant transport processes: advection and dispersion. Advection considers solute that is transported with the flow. The solute transport equation may be expressed in any flow direction. Rate of contaminant transport due to advection can be expressed as: Advective Mass Transport Rate CdAv x ?= (3.36) Hydrodynamic dispersion considers the processes of mechanical dispersion and molecular diffusion. These processes are lumped into one term referred to as the dispersion coefficient, D x , where D* represents molecular diffusion: *DvD xxx +=? (3.37) Dispersive Mass Transport Rate dA x C D x ? ? =? (3.38) The total solute transport rate per cross sectional area, F x , may then be expressed as: x C DCvF x xxxx ? ? ?= ?? (3.39) 41 For a three-dimensional problem, there will similarly be a F y and F z , (total solute transport in the y and z directions). To apply conservation of mass, the total mass entering the REV must equal the mass leaving the REV and the mass accumulation. This is expressed as: Mass entering the REV is: dxdyFdzdxFdzdyFREV zyxin ++= (3.40) Mass leaving the REV is: dxdydz z F Fdzdxdy y F Fyzx x F FREV z z y y x xout ? ? ? ? ? ? ? ? ++ ? ? ? ? ? ? ? ? ? ? ++?? ? ? ? ? ? ? ? ? ? += (3.41) The difference between entering and leaving mass can be expressed as: dxdydz z F y F x F z y x ? ? ? ? ? ? ? ? ? ? + ? ? + ? ? (3.42) The difference is equal to the amount of dissolved substance accumulated in the element and the effect of source/sink. Therefore, the rate of mass change in the rev may be expressed as: dxdydzCqdxdydz t C ss ? ? ? ?? (3.43) The complete conservation of mass equation may then be expressed: ss z y x Cq t C z F y F x F ? ? ? ?= ? ? + ? ? + ? ? ? (3.44) Substituting values for Fx, Fy, and Fz into (3.44) and dividing by ? yields: () () () ? ss zyxzyx Cq t C Cv z Cv y Cv xz C D zy C D yx C D x + ? ? = ? ? ? ? ? ? ? ? + ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? +? ? ? ? ? ? ? ? ? ? This is the form of the solute transport equation solved in this study. (3.45) 42 3.2 SEAWAT Coupling Procedure In SEAWAT, the flow and transport equations derived above are coupled using an iterative procedure (See Figure 3.3). Essentially, the groundwater flow equation is solved to attain the velocity field. The velocity is then passed to the transport equation which calculates solute transport. The density field is then updated and the heads are recalculated. This iterative procedure is continued until the density difference is less than some user-specified tolerance. This iterative computational process is summarized below. Figure 3.3: Flow chart of SEAWAT coupling procedure 43 44 CHAPTER 4 MATERIALS AND METHODS USED IN LABORATORY STUDY In this chapter the materials and methods used for each of the individual physical model studies are presented. The objectives and construction details of the flow container are presented in Section 4.1. Section 4.2 includes discussion on the experimental setup and peripheral equipments used during the study. The three experiments presented in this work include: 1) Combined beach-face and open-well infiltration experiment, 2) isolated beach face infiltration experiment, and 3) sinking plume experiment. The experimental setups for each of these experiments are presented in Sections 4.3, 4.4, and 4.5, respectively. Details concerning the dye, porous medium, dispersivity values, hydraulic conductivity measurements are also presented in this chapter. 4.1 Flow Container Design and Construction A flow container was constructed to physically model density driven flow in an unconfined, homogeneous, isotropic, and saturated porous media system. The first design objective was to develop a modular flow container which can allow a flow balance to be measured non-intrusively. A non-intrusive photographic measurement technique was used and hence the concentrations of the plumes in the physical model analyzed by 45 qualitative methods. The second goal of the flow container design was to reduce the time and resource demand of the experimentation process. This was achieved by creating an experimental setup that was easy to set up and required little maintenance during experimentation and preparation between trials. The dimensions of the flow container were selected as 63 ? 2.7 ? 30.5 (cm). At this size the flow container was large enough to capture the processes of concern, yet still manageable from a maintenance perspective. The third design goal was to enhance the visualization of contaminant plume behavior in the subsurface. This was facilitated by selecting a transparent material for the flow container and a translucent porous media. Backlighting was used to enhance the quality of the image. Details of flow container design, porous media, and peripheral equipment are described in the following figures. Typical construction details of a flow container are shown in Figure 4.1. For this study, all of the walls of the flow container were constructed of 6.35 mm thick Plexiglass?. However, for increased support while transporting the model, a thicker piece of Plexiglass? was used for the base of the flow container (1.3 cm). All Plexiglass? pieces were cut with a table saw to the desired size. A series of head control orifices were created by drilling holes into the left and right wall pieces of the flow container. Orifices were staggered along the left and right walls to allow greater head control in the chamber (See Figure 4.3). A small strip of Plexiglass? was glued to the front and rear walls of the inside of the tank approximately 5 cm from the left and right walls (not shown in flow container schematic). These Plexiglass? obstructions were necessary to hold the stainless steel mesh in place to prevent the porous medium from moving into the head chambers on either end of the flow container. A water-proofing silicon adhesive was then applied to the base of the flow container to connect the front and rear wall pieces to the base. The left and right boundary walls were also lined with a silicon adhesive which were connected to the base and the front and rear walls. A series of stainless steel screws were drilled to permanently seal and reinforce the flow container. The flow container was allowed to dry for 3 days before checking for leaks. Any leaks were sealed with a water proof silicon adhesive. Plexiglass? C-clamp reinforcements were applied to the top of the flow container for additional support and to reduce the affects of container bulging due to the weight of the porous medium. 3 4 2 3 2 4 1 Figure 4.1: Flow container Dimensions of Flow Container Pieces: 1) Base: 63?10 cm 2) Walls: 63?30.5 cm 3) Ends: 4?30.5 cm 4) US 16 Mesh Stainless Steel Screen: Approximately 2.7?30 (cm), the screen size was cut to fit 46 Stainless steel mesh Figure 4.2: Flow container Head Control Orifices Orifices used to drain porous medium at end of experiment Figure 4.3: Flow container 47 4.2 Experimental Setup The six main components of the experimental setup were: 1) Freshwater constant head supply reservoir (FWCHSR), 2) Saltwater constant head supply reservoir (SWCHSR), 3) Right constant head chamber (RCHC), 4) Left constant head chamber (SWCHC), 5) Porous medium chamber (PMC), 6) Dense contaminant source supply. A generic schematic of the main components Left Constant Head Cha Porous M From Pum Dense Sour 2 Valve 4 3 mber Figure 4.4: Experiment setup 48 of the experiment is shown in Figure 21. Valve 6 3 1 edium 5 4 Right Constant Head Chamber Peristaltic p, Supplies ce 49 The FWCHSR and SWCHSR each consisted of a five-gallon plastic cubic-shaped container with two outflows. The FWCHSR supplied the right constant head chamber and SWCHSR supplied the left constant head chamber as shown in Figure 4.4. A valve was used to control the flow rate exiting the FWCHSR and SWCHR through a ?? diameter tube. The flow from the FWCHSR supplied freshwater to the right constant head chamber (RCHC). Likewise, the SWCHSR supplied the saltwater to the left saltwater constant head chamber (LCHC). The SWCHSR was only used for experiments requiring a saltwater boundary condition (i.e. experiments simulating near-ocean conditions). The RCHC was used for controlling the boundary head level, H2, on the right side of the porous media. The head, H2, was controlled by a series of orifices lining the right wall of the RCHC. The orifices could be opened or closed with impermeable plugs to alter the head, H2. The saltwater head, H1, was controlled in a similar manner. For all flow experiments, H2 was higher than H1 in order to create a head gradient causing flow from the right boundary to the left boundary. The constant head chambers located on either side of the porous media were separated from the porous media flow chamber by a #16 stainless steel mesh screen provided by GeoTest Instrument Corp. Depending on the specific experiment (beach face infiltration or sinking plume) different components of the experimental setup were used. Therefore, the next two sections describe the peripheral equipment and the experimental setup for the combined beach face and open well, beach face infiltration, and sinking plume experiments, respectively. 4.3 Beach Face Infiltration with Well Experiment For the beach face infiltration with well experiment, the tsunami water entered the freshwater flow region by injection through an open well, in addition to the beach face infiltration. A 2-cm diameter impermeable glass tube, open at both ends, was inserted into the top of the porous media model to represent an open dug well as shown in Figure 4.5. Due to density differences, the saltwater from the left boundary intruded the porous medium, as described in the previous section, creating a wedge. Beach face infiltration was simulated by applying a volume of the denser saltwater to the surface (top) of the porous medium. Input to the well was achieved by a simultaneous and instantaneous filling of saltwater in the empty head space above the groundwater table in the well. BEACH FACE INJECTION OPEN WELL Figure 4.5: Generalized schematic of experimental setup illustrating the combined Source-2 and Source-3 saltwater discharges into an aquifer 4.4 Beach Face Infiltration Experiment This experiment required the presence of a salt wedge to study the interaction of the wedge with the tsunami water. The same flow container used for the sinking plume experiment was used for the beach face infiltration experiment. However, for this experiment, the left freshwater boundary was replaced with a saltwater boundary, to 50 represent the coastal interface, by injecting saltwater at a constant rate from a saltwater supply reservoir as shown in Figure 4.6. The saltwater from the boundary was allowed to intrude into the aquifer from the left to form a stable saltwater wedge. The steady initial condition was attained when the saltwater wedge reached a steady position and the lateral movement of the wedge ceased to penetrate further into the porous medium. Furthermore, this initial steady condition, representing the pre-tsunami scenario was attained by setting up a regional hydraulic gradient across the length of the model that would force the uncolored freshwater to flow from right to left. The inundation caused by the tsunami wave was simulated by instantaneously discharging approximately 500 mL of saltwater at the top surface of the model. The freshwater flow was continuously supplied to the model from the constant head boundary at the right side. To avoid any variation in density or chemical composition, the same batch of saltwater was used in a particular set of experiments to supply both the saltwater boundary and the tsunami water from above. Experiments were performed for different ambient groundwater flow rates. BEACH FACE INJECTION Figure 4.6 ? Generalized schematic of experimental setup illustrating saltwater contamination from Source-2 type region infiltrating through the beach face 51 52 4.5 Sinking Plume Experiment Sinking plume experiment was conducted in a thin transparent rectangular flow container as described in Section 4.1. For the sinking plume experiment, the only inflow into the system was supplied through the freshwater constant head chamber from the freshwater supply reservoir and into the porous media flow chamber (PMFC) due to a head gradient. The PMFC was filled with glass beads uniformly packed at 5 cm increments until the container was filled. In all experiments the ambient fresh groundwater flowed from right to left. Seawater was prepared by dissolving NaCl in de- ionized water. The average density of saltwater used in the experiments was 1025 kg/m 3 , which is representative of sea water. A small amount of red food coloring was added to the saltwater to differentiate the injected tsunami water from the freshwater. The saltwater was supplied as a continuous point source to the surface of the porous media through a continuous volume displacement peristaltic pump at a constant rate (Q L =8.9 mL/min). The flow rate was sufficiently low to prevent any ponding effects. The system was used to simulate a leaking saltwater pond. Three experiments were performed to study the effects of different ambient groundwater flow rates. The rates were controlled by varying the hydraulic gradient in the flow direction by adjusting the water levels in the end reservoirs. A drain tube with relatively low conveyance was connected near the bottom of the outflow chamber to allow the dense water which accumulates in the left boundary to exit. This was important to avoid changes in the head along the left boundary due to NaCl accumulation in the constant head chamber. Digital photographs were collected with a 7.1 mega-pixel camera. Figure 4.7: Generalized schematic of experimental setup illustrating saltwater contaminant from a saltwater plume infiltrating into a coastal aquifer 4.6 Porous Medium A laboratory simulation of realistic groundwater flow and transport scenarios is a difficult task due to the complexities of the porous media itself. Porosities and permeabilities will typically vary with position and direction. Ideal porous media, such as glass beads, can be used to minimize the physical complexities existing in the field. In addition, the use of glass beads as a porous media in physical groundwater flow studies enhances flow visualization due to its translucent property when lighting is introduced. The selection of a larger diameter porous media will typically lead large conductivity values (> 500 m/day). The porous medium selected for this study was an A-110 glass bead obtained from Potters Industries. The mean bead size diameter was 1.1 mm with a range of 1.0 - 1.2 mm. The bead size yielded a high value for hydraulic conductivity, but was acceptable for this experimental work since the goal of the study was to develop conceptual models 53 54 using easy to manage flow visualizations experiments. Other similar studies required a significantly larger time demand to setup and perform experiments. 4.7 Dye A red food dye (FDC 40&3) was used as an optical tracer to observe the infiltrated saltwater in both the sinking plume and tsunami experiment. A green food dye was used as an optical tracer for the saltwater wedge in the tsunami experiment. The food dye was selected because it had a density similar to that of freshwater and it was found to have no measurable sorption to the glass beads. Furthermore, the selected dyes were very soluble in water even at high concentrations. 4.8 Hydraulic Conductivity Measurements The hydraulic conductivity was measured using three separate techniques: 1) Ex-situ measurement using a constant head permeameter setup, 2) in-situ measurement of flow through the flow container using a known head difference (Oostrom et al. 1992) 3) and in-situ measurement of the velocity of a tracer injected into the system. Once the necessary datasets were attained, Darcy?s law was manipulated to determine values for hydraulic conductivity as described in the following paragraphs. 4.8.1 Technique 1 A glass column of inner diameter 4.8 cm and length 30 cm was uniformly packed to fill approximately half of the column length with the saturated porous medium. The column was submerged in a bucket of water to preserve saturation. A stainless steel mesh screen was applied to the bottom of the column to hold the porous media in place. The column was situated in a vertically upright position by means of a column stand (See Figure 4.8). A flow was then introduced into the top of the column and allowed to exit via the stainless steel screen at the bottom of the column and into the plastic bucket which was allowed to overflow into the drain. The vertical position of the column could be altered by adjusting the column stand. The column?s vertical position was modified until the water level in the column reached a steady state. After obtaining the head difference and the corresponding flow rate (Q) for a steady state condition the hydraulic conductivity could be calculated. The mean value was calculated by this approach to be approximately 1130 m/day. Figure 4.8: Hydraulic conductivity measurement with column setup 55 Figure 4.9: Photograph of column setup 4.8.2 Technique 2 This test was the first of two hydraulic conductivity tests performed in-situ in the flow container. The glass beads were uniformly wet-packed to preserve saturation in the system. The packing was conducted in layers, approximately 5-cm thick followed by tamping, to achieve homogeneity and isotropy in the porous medium. A constant head difference was established in the system by introducing a flow into the right chamber of the flow container. Overflow orifices allowed excess fluid to leave the system through the right boundary wall of the flow container. Likewise on the leftmost boundary wall, fluid was allowed to exit the system through orifices where the fluid could be collected for measurement. The flow was allowed to reach steady state for a period of 5 to 10 minutes. Once steady state was attained the flow rate through the system, Q, was measured with a graduated cylinder and stop watch. The excess fluid exiting through the right wall could also be measured to capture the flow balance. The hydraulic conductivity was then calculated by manipulating Darcy?s Law. Multiple experiments were completed and estimated values ranged from 1050 to 1150 with an average value of 1100 m/day. 56 4.8.3 Technique 3 The third technique for measuring hydraulic conductivity involved injecting a tracer of the same density as the ambient groundwater and observing its movement through the porous medium. The system was packed and the horizontal flow was established as prescribed in technique 2. However, instead of measuring the flow through the system, a conservative tracer was injected at the right boundary to track the velocity of the tracer through the system as shown in Figure 4.10. A stopwatch was used to measure the velocity of the tracer. The product of the groundwater velocity and the porosity is equal to the specific discharge or Darcy flux (q). Given the head gradient and the specific discharge the hydraulic conductivity could be calculated. The average value estimate using this approach was 1030 m/day. Note, this method suffers from qualitative interpretation of the advective front location, which can influence the estimated value. Figure 4.10: Tracer test for hydraulic conductivity 4.9 Porosity and Dispersivity Values The porosity was measured by filling a graduated cylinder of a known volume with glass beads. After packing the beads in the graduated cylinder in a similar manner to what would be performed in the actual flow container, the displaced water could be subtracted from the initial total volume of water to determine the volume of water present in the pores, commonly referred to as the volume of voids. By definition, porosity is the 57 ratio of the volume of voids to the total volume (? = V v /V T ). Porosity was measured multiple times and found to be 0.385. See Figure 4.11 for a description of this porosity measurement technique. V w = Volume of Water V D = Volume of Displaced Water V V = Volume of Voids = V w - V D V T = Total Volume of Water and Porous Media ? = Porosity = V v /V T Figure 4.11 : Porosity experiment Table 1: Porosity Measurement Trial 1 Trial 2 Volume of Voids, Vv 42 54 Total Volume of Water and Porous Media, VT 109 141 Porosity, n 0.385 0.383 Porosity Measurement Longitudinal dispersivity values were taken from (Goswami and Clement 2007), who performed a tracer test on the same porous medium and fitted the observed front. The values were attained for the same glass beads used in this study. The estimated value of longitudinal dispersivity was 1.0 mm, which is of the order of the grain size. Others have found similar small longitudinal dispersivity values of the order 1 to 2 mm, for uniform glass beads (Oswald 2004). 58 59 CHAPTER 5 RESULTS OF PHYSICAL MODEL EXPERIMENTS 5.1 Problem Definition Preliminary site investigation studies completed in Sri Lanka indicated that the tsunami waters entered and contaminated these unconfined aquifers via 1) direct infiltration from the land surface, 2) discharge through the open dug wells, 3) slow infiltration/leakage from shallow depressions and lagoons located along the coastal zone (Illangasekare et al. 2006). These contamination sources are identified in the conceptual model shown in Figure 5.1. Sources 1 and 2 represent the infiltration of saltwater over the beach face and the deposition of saltwater inside the open well, respectively. Source 3 represents the saltwater that traveled with the tsunami run-ups and was deposited in ponds and depressions at considerable distances inland. It is important to distinguish Source 1 from Source 3 due to the residence time of the saltwater in the aquifer. Source 1 is deposited on the beach face near the ocean hence has a relatively short residence time, while Source 3 contaminants are located behind sand dune formations where the saltwater remains pooled and represents a continuous source to the aquifer. 60 Figure 5.1: Cross-sectional view of the sources of contamination to a coastal aquifer (Villholth et al. 2005). Source-1 Source-3 Source-4 Source-2 Essentially, energy from wind traveling over the sea will decrease once it encounters the elevated land surface with obstructions such as vegetation (mangrove forests). When the wind energy is reduced, dust and sand particles are deposited which creates a topographical high point or mound, commonly known as a sand dune. Over time vegetation will grow on the sand dune which gives it more strength due to the roots. Sand dunes serve as a natural energy dissipater against storm waves, surges, and wind. However, if run-up levels exceed the height of the sand dune then the seawater becomes trapped on the other side of the sand dune where it develops into a pooled contaminant source (source 3) for coastal aquifers. Typically, sand dune formations run parallel to the coastline which means that the aquifer of the region can be conceptualized as a series of two-dimensional strips situated perpendicular to the sea. The profile of the plume can be assumed to be the same at any cross-section along the coast. Therefore, the plume movement can be conceptualized as a two-dimensional problem, as shown in Figure 5.2. Beach (Source 1) Open Wells (Source 2) Lagoon (Source 3) Ocean Saltwater Wedge (Source 4) Figure 5.2: Conceptual picture illustrating 2-D assumption and contaminant sources (Zheng, Bennett 2002). The low topography of the Sri Lankan coastline and the small head gradients made the region susceptible to saltwater intrusion from beneath the fresh groundwater. The level of intrusion was also impacted by the waves which increased the hydraulic head of the ocean, causing the saltwater wedge to move inland temporarily until the ocean receded to its normal level. The transient nature of the waves followed by the recession may have increased the width of the zone of dispersion between the fresh groundwater and saltwater in the subsurface. This source represents a fourth potential source of saltwater contamination to the Sri Lankan coastal groundwater resources. The combination of the contamination from open wells and infiltra surface coupled with the effects of saltwater intrusion from below created a unique variable-density flow problem. In the past, the effects of dense contaminant plumes emanating from surface sources and saltwater intrusion have been studied independently. However, in Sri Lanka, these two processes models of the behavior of the des have occurred simultaneously. Physical cribed sources were investigated in three separate 61 tion at the land?s 62 studies. First, a beach face infiltration and open well experiment is presented in Section 5.2. In the subsequent experiments, the beach face infiltration (Source 1) and the sinking plume emanating from a salt pond (Source 3) are isolated to understand how the sources behave independent of each other. Beach face infiltration (Source 1) and sinking plume (Source 3) experiments are presented in Sections 5.2 and 5.3, respectively. Further investigation of a Source 3 plume was performed via numerical modeling and a sensitivity analysis in the final sections of this chapter. 5.2 Beach Face Infiltration and Direct Injection Through an Open Dug Well Tsunami-scenario simulations were conducted to investigate the migration of infiltrated seawater along a beach face (Source 1) in an unconfined aquifer and also to gain a qualitative understanding of the fate of the saltwater discharged into the open dug wells (Source 2). Experiments were conducted in the rectangular flow container described in Chapter 4. The system was fully saturated with a small capillary fringe zone. A 2-cm diameter impermeable glass tube, open at both ends, was inserted into the model to represent an open dug well. Seawater was prepared by dissolving NaCl in de- ionized water. The average density of saltwater used in the experiment was similar to that of seawater (1025 kg/m 3 ). To avoid any variation in density or chemical composition, the same batch of saltwater was used in a particular set of experiments to supply both the seawater boundary in the left reservoir and the tsunami water from above. A small amount of food coloring was used to differentiate sea, fresh and tsunami waters for visualization of the saltwater migration. A schematic of the experimental setup is shown in Figure 4.5, in the previous section. The initial condition, representing the pre-tsunami scenario was simulated by setting up a regional hydraulic gradient across the length of the model that would force uncolored fresh water to flow from right to left (Refer to Figure 4.5). The solid-colored saltwater was allowed to intrude into the aquifer from the left to form a stable saltwater wedge. The inundation caused by the tsunami wave was simulated by evenly discharging approximately 500 mL of red-colored saltwater at the top surface of the model and also directly into the open well. Tsunami waters invaded the unconfined aquifer model by infiltration from the top of the water table boundary. Tsunami water also entered the fresh water flow region by injection through the open well. The fresh water flow was continuously supplied to the model from the constant head boundary at the right side. Photographs of tsunami experiment are shown in Figure 5.3. Tsunami Problem T = 0 min (Initial Condition) T = 2 min T = 4 min T = 6 min 63 T = 8 min T = 16 min T = 10 min T = 18 min After Well-Cleaning T = 12 min T = 19 min T = 14 min T = 21 min Figure 5.3: Photographs of the beach face infiltration (Source 1) and open well injection (Source 2) experiment 64 65 The initial steady state condition at time = 0 min is shown in the first photograph of Figure 5.3. The experimental results demonstrate that within a short period of time the tsunami water that flooded the well, descended as a large slug and contaminated the deeper aquifer. The infiltrated tsunami water (Source 1) eventually merged with seawater from the flooded well (Source 2). Another major observation was that the open well remained contaminated for a relatively long period. As observed in the figure, the bulk of the tsunami water injected into the well moved downward as a blob within about five minutes. However, a small amount of seawater remained trapped inside the well, and this remnant saltwater acted as a continuous source of saline water that discharged from the well. After about 19 minutes, a syringe was used to extract three storage volumes directly from the open well and this helped to remove the remnant saltwater from the well. 5.3 Beach Face Infiltration Experiment (Source 1) The beach face infiltration experiment represented a Source 1 type contamination. As discussed in the materials and methods section the tsunami source (Source-1) was created by instantaneously distributing a fixed volume of saltwater directly on the top boundary of the model (see Figure 4.6). This experiment was designed to isolate a Source 1 contaminant. Experiments were performed in a high and low flow system. The physical dataset shown in Figure 5.4 shows the transport of this dense saltwater slug over a time period of 11 minutes. The saltwater wedge intruded considerably less in the high flow experiment than in the low flow experiment. Under high flow conditions the saltwater slug was transported quickly through the model and formed a rather shallow contamination zone. On the other hand, when the ambient flow was low, the slug of saltwater migrated downward and contaminated larger volume of aquifer. There was also more intense fingering in the low flow case. An important observation made in these experiments was that tsunami water (red- colored salt water) always remained above the regional seawater interface. This should be expected because the density of the tsunami water sources decreased due to dilution as they migrated through the aquifer and mixed with freshwater. Therefore, the tsunami water would never be able to penetrate the regional seawater wedge. Instead, it will remain floating above the wedge and will eventually be advected by the regional freshwater flow along the wedge. 3 minutes 7 minutes 11 minutes High Flow Low Flow Figure 5.4: Beach face infiltration in a high (left photographs) and low (right photographs) flow groundwater system 66 67 5.4 Sinking Plume Experiment Sinking plume studies were conducted to investigate the effects of contaminant seawater deposited at some distance inland due to a tsunami or storm surge event. This study describes a Source 3 contaminant as shown in Figure 5.1. It is assumed that the contaminant source, i.e. flooded lagoon or inundated region, is located at a distance inland such that the effects of the seawater boundary can be ignored (Zhang et al. 2000). The experiments were conducted in the same flow container using a setup shown in Figure 4.7 in the previous section. As shown in Figure 4.7, An injection tube connected to a peristaltic volume displacement pump was situated at the top of the porous medium chamber. This pump supplied the dense saltwater which simulated the trapped tsunami water deposited above the shallow unconfined aquifer. Once a steady fresh groundwater flow from the freshwater constant head tank had been established, the peristaltic volume displacement pump was turned on and maintained at a constant loading rate, Q L . The dense saltwater was created by mixing NaCl with dyed deionized water. The density was monitored throughout the experiment with a hydrometer (ASTM-111H) to ensure that the saltwater was maintained at the accepted density of seawater (1025 kg/m 3 ). The migration of the dense plume was observed and recorded using digital photographs taken at specific time intervals. Photographs were taken throughout the experiment using a 7.1 megapixel camera and are presented in this section. Typically, photographs of the early time data were taken every minute until a quasi-stable state had been reached upon which the time interval between photographs was increased to five minutes. All experiments were run for at least 30 minutes. Fluxes were observed for each experiment to verify a volume flow balance in the system. Experiments were repeated for various head differences in the system. A summary of the hydraulic conditions used in these experiments are given in Table 2. Table 2: Summary of Sinking Plume Experiment Experiment ID Loading Rate Hf Hs Head Gradient Rin Rout Lout Flux Error mL/min cm cm % mL/sec mL/sec mL/sec % A-1 8.9 28.9 27.0 3.58 4.22 0.56 3.52 3.32 A-2 8.9 28.9 27.2 3.21 4.11 0.97 3.15 0.24 A-3 8.9 29.0 28.0 1.89 4.11 2.11 1.95 1.22 A-4 8.9 29.0 28.4 1.13 4.11 2.62 1.43 1.46 A-5 8.9 29.0 28.3 1.32 4.11 2.75 1.25 2.68 A-6 8.9 29.0 28.8 0.38 4.01 3.50 0.52 0.25 Sinking Plume Experiment The experiments are presented in the order of decreasing plume stability. The A-1 experiment had the largest head gradient and therefore was the most stable. The only experiment in which the head gradient did not correlate with plume stability was experiment A-5 in which the experiment was observed to have had a higher head gradient than A-4 but was found to be less stable. This could be attributed to human error of reading the left and right head levels. Figure 5.5 shows the result of experiment A-1 which had a very high head gradient. As the figure shows, the plume behavior appears to be stable since the body of the plume remains concentrated in the upper most region of the model. This means that the forced convection due to the flow is high compared to the free downward convection due to density differences. The contaminant moves relatively quickly through the model compared to the experiments with a lower head gradient. 68 Sinking Plume Experiment Stable A-1 T = 1 minute T = 3 min T = 5 min T = 7 min T = 9 min T = 11 min T = 13 min T = 15 min 69 T = 17 min T = 19 min T = 26 min T = 31 min Figure 5.5: A-1 Stable plume experiment Experiment A-2 had only a slightly smaller head gradient than experiment A-1 therefore the plume profiles are similar. As expected, the A-2 plume penetrates deeper into the model than the A-1 plume. The plume still exhibits a stable behavior because fingering and instabilities were not observed along the interface. Stable A-2 T = 1 min T = 3 min 70 T = 5 min T = 7 min T = 9 min T = 12 min T = 14 min T = 16 min T = 18 min T = 20 min 71 T = 25 min T = 30 min Figure 5.6: A-2 Stable Plume Experiment The gradient was reduced considerably in experiment A-3. Instabilities and fingering are on the verge of forming which is why this problem can be considered an intermediate case. In fact, this experiment was useful for determining the critical ? 1 number for our system. This was because the experiment represented a border case in which any decrease in head gradient would have caused noticeable fingering to develop. Another observation, was the length of time it took the plume to reach the left exit boundary was considerably less than the two prior experiments A-1 and A-2. Stable A-3 T = 1 min T = 3 min 72 T = 5 min T = 7 min T = 9 min T = 11 min T = 13 min T = 15 min T = 17 min T = 19 min 73 T = 24 min T = 29 min Figure 5.7: A-3 Stable Plume Experiment Experiment A-4 was the first of three unstable plume problems. Clearly, instability and fingering are generated along the interface. The fingers appear to have a similar shape but are randomly spaced. It appears that the size of the fingers tend to grow as the distance from the source is increased. Unstable A-4 T = 1 min T = 3 min T = 5 min T = 7 min 74 T = 9 min T = 12 min T = 14 min T = 16 min T = 18 min T = 20 min T = 25 min T = 30 min Figure 5.8: A-4 Unstable plume experiment 75 The head gradient was reduced further for experiment A-5. There was more fingering evident in Experiment A-5 due to the fact that the fingers were more closely spaced in Experiment A-5 than in A-4. The free convective force due to density now appears to be a major factor in the plume?s behavior. Unstable A-5 T = 1 min T = 3 min T = 5 min T = 7 min T = 9 min T = 12 min 76 T = 14 min T = 16 min T = 18 min T = 20 min T = 25 min T = 30 min Figure 5.9: A-5 Unstable plume experiment The head gradient was greatly reduced for the final experiment, A-6. This system is a free convection dominated problem. The behavior of the plume is highly unstable which is evident in the photographs shown in Figure 5.8. The plume sinks vertically downward in the shape of one big finger and accumulates on the bottom of the tank. This degree of instability represents the least desirable to encounter in the field environment as this 77 plume has the potential to contaminate larger regions of an aquifer than the previous experiments, A-1 through A-5. A-6 (Highly Unstable) T = 1 min T = 3 min T = 5 min T = 7 min T = 9 min T = 12 min 78 T = 14 min T = 16 min T = 18 min T = 20 min T = 24 min T = 28 min T = 35 min T = 40 min Figure 5.10: A-6 Unstable plume experiment 79 5.4.1 Discussions of Plume Stability Plume stability was largely dependent on the head gradient in the system. Experiments A-1, A-2, and A-3 all had high head gradients which corresponded to high groundwater fluxes. In these experiments the forced convection of the flow was much greater than the downward free convection. Therefore, the stability of the plume was high and the plume penetration was low. However, for the small head gradient systems, A-4, A-5, and A-6, the plume tended to become unstable which was manifest by the occurrence of fingering and deeper penetration of the plume into the aquifer. This can be attributed to the fact that the forced convection decreased while the free convection remained the same. The ratio of free to forced convection is referred to as a non-dimensional stability number, ? 1 : x sat q K ? ? ? ? ? ? ? ?? = 0 1 ? ? ? (5.1) This number represents a system dependent empirical parameter that can characterize the behavior of a plume in a contaminated aquifer. For our study, the ? 1 number was observed to be approximately 1.3 which represents the value where a plume transitions from stable to unstable. The plume?s behavior will always obey the trend of the ? 1 number in the same porous medium provided that the loading rate remains constant. A drawback of the ? 1 number is that it does not take the loading rate into consideration (Schincariol et al. 1994). Another stability number called the Modified Rayleigh Number 80 (G?ven et al. 1992), Ra*, should be used if the effects of loading rate are to be taken into account. The Rayleigh number is expressed as: T satdp D KH Ra ? ?? )/( 0* ? = (5.2) Table 3: Plume Stability Experiment ID Ksat Density Difference Flow Rate Area Darcy Flux Pi Number Stability cm/sec g/cc mL/sec cm^2 cm^3/(cm^2*sec) A-1 1.31 0.025 3.52 75.465 0.0466 0.70 Stable A-2 1.31 0.025 3.15 75.735 0.0416 0.79 Stable A-3 1.31 0.025 1.95 76.95 0.0253 1.29 Stable A-4 1.31 0.025 1.25 77.355 0.0162 2.03 Unstable A-5 1.31 0.025 1.43 77.49 0.0185 1.77 Unstable A-6 1.31 0.025 0.52 78.03 0.0067 4.91 Unstable Stability of Sinking Plume Experiment 5.4.2 SEAWAT 2000 Simulations The design for the rectangular mesh used in the simulations consisted of 6264 active cells, 108 cells in the X-direction, and 58 cells in Y-direction. The mesh was discretized at 0.5 cm in the X and Y directions. Dirichlet boundary conditions were implemented at the left and right boundaries as the head values did not fluctuate considerably during the experiments. A well was used to simulate the contaminant source at the top layer of the model. The single cell well was positioned 10 cells to the left of the right boundary. The source density was 1.025 g/mL with a loading rate of 0.055 mL/sec as measured in the experiments. All numerical simulations were run for a simulation time of 3000 seconds (50 minutes). 81 SEAWAT was used to simulate Experiment A-3 and Experiment A-6. The total head difference of the system was the only variable that was manipulated to match the numerical profiles with the physical data. A 1 mm human error was allowed for each head reading. Head values were not adjusted beyond 1 mm of the recorded head reading. The numerical solutions and physical datasets are presented in the figures below. The numerical solutions are presented graphically (left) next to the photographs of the physical experiments (right). The physical data was presented in the previous section. The four contours represent the 20, 40, 50, and 60% concentration countours. The grid for the numerical simulations was divided into 0.5 cm ? 0.5 cm cells. 1 minute 3 minutes 6 minutes 9 minutes 82 12 minutes Figure 5.11: Stable plume: comparison of numerical solutions with physical data (20% outer contour, 40%, 50%, 60% inner contour) 1 minute 3 minutes 6 minutes 83 9 minutes 12 minutes Figure 5.12: Unstable plume: comparison of numerical solutions with physical data (20% outer contour, 40%, 50%, 60% inner contour) The numerical solutions matched the physical model data reasonably well for the stable and unstable plume scenarios. The numerical model was extremely accurate in predicting the movement of the front for the stable plume. Interestingly, for the unstable plume the penetration depth was modeled reasonably well while the front of the plume tended to curl in the direction of the fresh water boundary for the numerical simulation which was not observed in the physical model. The exact plume profile as observed in the physical model was not expected for several reasons. The system is highly sensitive to the initial conditions (chaotic system) and therefore the plume profiles can not be replicated exactly in a physical model setting. Furthermore, numerical errors caused by the necessary approximations of mesh discretization will be introduced, which will affect the profile of the plume. It appeared from the qualitative data that the transition zone was relatively small for both cases. However, in order to accurately compare the transition zone generated numerically and physically, a quantitative technique is needed to measure the concentration within the plume. 84 5.4.3 Sensitivity Analysis A sensitivity analysis was performed on the parameters of dispersivity, density, head difference, and loading rate which were the main flow parameters affecting the concentration profiles. The sensitivity analysis was performed by changing the hydraulic parameters for the unstable plume, A-6, which represented the base case. The parameters for the base case simulation are shown in Table 4. A Modified Choleski solver was used initially but was found to produce more numerical dispersion and took a considerably longer time to run than the SSOR solver. Therefore, we switched to the SSOR to more efficiently solve the problem. Table 4: Parameter and Problem Description for Base Case Study Number of Cells in X-direction 108 Number of Cells in Y-direction 1 Number of Cells in Z-direction 58 Transport Step Size 0.1 sec Total Head Gradient 0.38% Hydraulic Conductivity 1.31 cm/sec Longitudinal and Horizontal Dispersivity 0.1 cm Source Loading Rate 8.9 mL/min Source Density 1.025 g/cc Density Concentration Slope 0.714 Convergence Parameter for Coupling Iterations (Density Change) 0.01 SEAWAT Base Case Numerical Simulation for Unstable Plume The sensitivity analysis was performed by altering the parameters of a selected base case scenario. The base case parameters are presented in Table 4. The base case used for the sensitivity analysis was the same as the numerical solution presented in the previous section, refer to Figure 5.12. 85 The first parameter investigated in the sensitivity analysis was the head gradient. The head gradient is a significant parameter affecting plume stability because it is directly related to the groundwater flux in the system. The first simulation (leftmost) of Figure 86 5.13 illustrates the base case scenario (head gradient = 0.38%). The gradient was increased to 0.94, 1.5, and 2.1% for subsequent simulations, respectively. These head gradient values were selected because it was observed in the physical model that the plume transitioned from an unstable to a stable case over a head difference range of about 1.1 cm. This head difference range corresponds to a head gradient range of 0.0038 to 0.021. The second parameter that was investigated was the source density. For the physical model study the source density was taken to be that of seawater (1025 kg/m 3 ). The base simulation is presented as the middle simulation in Figure 5.14. The numerical simulation was run for source density values of 1012.5 and 1050 kg/m 3 , representing half the density and twice the source density of the base case, respectively. It is understood that 1050 kg/m 3 is an unrealistically high value for source density but this was acceptable as the goal was to understand the degree of sensitivity of the parameter to the system. Next, the loading rate was altered to values corresponding to half (0.0275 mL/sec) and twice (0.110 mL/sec) that of the base case (0.055 mL/sec). These simulations are presented in Figure 5.15. The base case is presented as the middle simulation in the figure. The final simulation investigated the effects of the longitudinal and transverse (downward) dispersivity. The values of longitudinal and transverse dispersivity were taken to be equal for this analysis. The base case scenario is presented as the left most simulation in Figure 5.16. The dispersivity values were increased in the subsequent simulations from the base case value (? L = ? T = 0.1 cm) to values of 1.0, 2.0, and 10.0 cm presented as the second, third, and fourth datasets, respectively. Sensitivity Analysis of Plume Stability Total Head Gradient 0.0038 0.0094 0.015 0.021 1 minute 3 minute 6 minutes 9 minutes 12 minutes Figure 5.13: Sensitivity analysis of head gradient (20% outer contour, 40%, 50%, 60% inner contour) 87 Sensitivity Analysis of Plume Stability Source Density 1012.5 kg/m 3 1025 kg/m 3 1050 kg/m 3 1 minute 3 minutes 6 minutes 9 minutes 12 minutes Figure 5.14: Sensitivity analysis of source density (20% outer contour, 40%, 50%, 60% inner contour) 88 Sensitivity Analysis of Plume Stability Source Loading Rate Q L 0.0275 mL/sec 0.055mL/sec 0.11 mL/sec 1 minute 3 minutes 6 minutes 9 minutes 12 minutes Figure 5.15: Sensitivity analysis of source loading rate 89 (20% outer contour, 40%, 50%, 60% inner contour) Sensitivity Analysis of Plume Stability Longitudinal and Downward Dispersivity ? L = 0.1 cm ? L = 1.0 cm ? L = 2 cm ? L = 10 cm ? D = 0.1 cm ? D = 1.0 cm ? D = 2 cm ? D = 10 cm 1 minute 3 minutes 6 minutes 9 minutes 12 minutes Figure 5.16: Sensitivity analysis of dispersivityz (20% outer contour, 40%, 50%, 60% inner contour) 90 91 The general trends of the sensitivity analysis suggest that the parameters had a largely varying effect on the plume?s profile. In the first simulation it is observed that an increase of the head gradient from 0.0038 (base-case) to 0.021 was enough to transition the plume from unstable to stable. This is a similar gradient to what was required to transition the plume from stable to unstable in the physical datasets. The following simulations were not physically modeled so the discussion is limited to the observed trends of the numerical solutions. As the source density increases the plume tended to penetrate faster while the width of the plume is reduced (See Figure5.18). The faster penetration is expected as the free convection would increase with increasing density. Interestingly, if the source density is held constant and the loading rate is increased the width of the plume will increase while the penetration of the plume increases to a lesser degree. Shenoy (1994) concluded that the loading rate did play a role in plume stability but to a much less effect than other critical parameters as is observed in this sensitivity analysis. The dispersivity represents to what degree the plume will mix with the ambient groundwater. Mixing with the ambient groundwater is a major factor affecting plume stability. As mixing is increased the likelihood of generating a stable plume also increases. For the base-case scenario, (? L = ? T = 0.1 cm), the level of plume mixing is small. Therefore, the plume stability is also small and tended to penetrate deeper into the aquifer than the cases with higher dispersivity. The simulations with higher dispersivity penetrated much less and tended to spread more in the lateral direction. 92 CHAPTER 6 SUMMARY AND CONCLUSIONS On December 26, 2004, the country of Sri Lanka was devastated by an ocean- wide tsunami. The tsunami water greatly impacted the freshwater resources of the region which contaminated an estimated 40,000 open dug wells located near the coast. Just after the monsoon season, a decrease in the salinity level of the wells near the coast was observed (Villholth et al. 2005). In the months after the monsoon, the rate of salinity decline in the wells decreased and actually began to increase in some of the wells. In this study, the mechanisms of saltwater contamination due to a tsunami run-up into the Sri Lankan coastal aquifer were classified into three main categories: 1) Beach face infiltration, 2) Open well injection, and 3) Slow infiltration through inundated coastal lagoons and ponds (See Figure 5.1). These sources were studied simultaneously in a laboratory setting to develop some general hypotheses about transport and describe the coupled processes that transpires in the subsurface after a large-scale inundation event, such as a tsunami. The physical model data presented in this study attempted to describe the mixing and transport of the saltwater in a coastal aquifer discharged from multiple sources. The sources were simulated in the physical model as a single isolated source as well as multiple sources that simultaneously discharged dense water into the system. The 93 first experimental study simultaneously simulated beach-face infiltration and open- well sources. The physical experiments indicated that the tsunami water penetrated the porous medium relatively quickly from the beach face infiltration and from injection through an open well. The beach injection source penetrated the aquifer almost immediately with no observable run off. The well source penetrated the aquifer rapidly but migrated laterally due to the high groundwater flow. The porous media had a high permeability resulting in a flushing time of less than half an hour. There was very little mixing of the contaminant with the saltwater wedge. This was due to the mixing of the contaminant solution with the ambient groundwater flow. It took a much longer time for the contaminant that reached the wedge to leave the system. In the second laboratory experiment we isolated and studied the beach face infiltration source in a coastal environment. The experimental setup was similar to the first study but the open well was not used in this case. The experiment was performed in a high and low flow scenario which was attained by altering the head difference at the boundaries of the system. The first observation was that the size of the saltwater wedge changed depending on whether the system had a high or low flow initial condition. A noticeable increase in the amount of fingering was observed in the low flow experiment. There also appeared to be some saltwater mounding near the saltwater/freshwater interface for the low flow scenario. Some variability in the experiments was expected as the tsunami water appeared to be applied non-uniformly to the top of the model which may have affected the saltwater transport behavior. 94 The final experimental study employed a setup similar to ones described above to investigate the injection of saltwater into coastal aquifers via infiltration from the inundated ponds and lagoons. Therefore, the principal direction of the ambient groundwater flow can be represented as a uniform flow field toward the sea. Six trials were conducted to investigate the stability patterns of the evolving plume. The plume stability was determined to be highly sensitive to the head gradient. This can be partially attributed to the high conductance of the porous medium, since the flux through the system was found to be greatly impacted by slight changes (fractions of a cm) in the head gradient. The major observation of the sinking plume study was that there are two convective processes which are occurring, a downward convection due to density differences between the tsunami water and a horizontal convection due to the ambient groundwater flow. For a high flow regime, the plume appears to be carried quickly through the model and experienced minimal mixing with the ambient groundwater. For a low flow regime, the plume sunk vertically downward in a fairly uniform shape. Fingering was observed for the intermediate flow case which means that the system was in a quasi-stable state. The USGS numerical code SEAWAT was used to compare the numerical solutions with the actual plume profile as observed in the physical model under stable and unstable cases. The numerical solutions appeared to match reasonably well. Since we did not collect any quantitative concentration data from the physical model, a detailed analysis of the level of mixing simulated by the numerical model was not performed. The datasets could only be compared quantitatively. 95 A sensitivity analysis was performed by altering four key parameters that are expected to greatly influence the plume stability: 1) Horizontal groundwater flow 2) Relative density difference 3) Loading rate of contaminant solution, and 4) Dispersion The sensitivity analysis of plume stability was compared to a base case numerical solution. The base case corresponded to a numerical solution of a sinking plume that was found to have a reasonable match with the physical model data. Each of the sensitivity parameters investigated was found to have an effect on plume stability. However, our results show that loading rate is the most insensitive parameter. In our sensitivity analysis all four variables were perturbed by similar factors. While such a study provides relative sensitivity, it may not be indicative of what may transpire in the field setting since certain parameters may have a greater range of variability than others. However, these simulation results provided some general hypotheses about the flow and transport processes which occur following an inundation event. In this study, a conceptual basis was established for understanding the saline water transport processes from both surface infiltration and from direct injection into the wells. These experiments were performed in a highly idealized system. The data presented here is based on a preliminary conceptual study and should be interpreted as such. Additional studies are needed to create physical and/or numerical models that are more representative of the field case. Under field conditions, heterogeneities would play a larger role in the mixing of the variable density fluids. However, these 96 conceptual experiments are a prototype for future studies with the goal being to provide insight into the density driven transport processes that would have occurred as the result of a large-scale tsunami event. 97 REFERENCES Barlow, P. M. (2003). "Groundwater in Freshwater - Saltwater Environments of the Atlantic Coast." U. S. Geological Survey. Chapra, S.C., R. P. Canale (2002). Numerical Methods for Engineers with Software and Programming Applications, Mc Graw Hill. Custodio, E. (2002). Coastal Aquifers as Important Natural Hydrogeological Structures. Groundwater and Human Development, Bremerhaven, Germany Diersch, H.-J.G. and O. Kolditz (2002). ?Variable-density flow and transport in porous media: approaches and challenges.? Advances in Water Resources (25): 899-944. Elder, J. W. (1967). "Transient convection in a porous medium." Journal of Fluid Mechanics 27(3): 609-623. Fetter, C. W. (1988). Applied Hydrogeology. Columbus, Merrill Publishing Company. Freeze, R.A. and J.A. Cherry (1979). Groundwater. Prentice Hall, Inc. Guo, W. and C. D. Langevin (2002). User's guide to SEAWAT: A computer program for simulation of three dimensional variable-density groundwater flow. Goswami, R.R. and T.P. Clement (2007). Laboratory-scale Investigation of Saltwater Intrusion Dynamics (Draft Document). Civil Engineering. Auburn University. Auburn, Alabama. G?ven, O., J. H. Dane, W. E. Hill, J. G. Melville (1991). A review of the near-source physical behavior of subsurface contaminant plumes at electric utility sites. Auburn, Alabama, Electric Power Research Institute. Hanf, L. F. M., T.M. Poston, R.L. Dirkes (2005). Summary of the Hanford Site: Environmental Report for Calendar Year 2004. Richland, WA, Pacific Northwest Laboratory / U.S. Department of Energy. Hayworth, J.S. (1993). A physical and numerical model study of three-dimensional behavior of dense aqueous phase contaminant plumes in porous media. Civil Engineering. Auburn, Auburn University: 326. 98 Henry, H. R. (1964). "Interface between saltwater and freshwater in coastal aquifers." USGS Water Supply Paper(#1613C). Holzbecher, E. (1998). Modeling Density-Driven Flow in Porous Media. Berlin, Springer. Illangasekare , T., S. W. Tyler, T. P. Clement, K. G. Villholth, A. P. G. R. L. Perera, J. Obeysekera, A. Gunatilaka, C. R. Panabokke, D. W. Hyndman, K. J. Cunningham, J. J. Kaluarachchi, W. Yeh, M. T. van Genuchten, K. Jensen (2006). "Impacts of the 2004 tsunami on groundwater resources in Sri Lanka." Water Resources Research 42. Kipp Jr., K.L. (1997). Guide to the revised heat and solute transport simulator: HST3D - Version 2. Denver, CO, USGS Konikow, L. F. (1996). Numerical models of groundwater flow and transport. Vienna, Austria, International Atomic Energy Agency: 59-112. Konikow, W. E. Sanford, and P.J. Campbell (1997). "Constant Concentration boundary condition: Lessons from the Hydrocoin variable-density groundwater benchmark problem." Water Resources Research 33(10): 2253-2261. Langevin, C.D., Shoemaker, W.B., and Guo, Weixing, (2003). MODFLOW-2000, the U.S. 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): U.S. Geological Survey Open-File Report 03-426, 43 p. List, E. J. (1965). The stability and mixing of a density-stratified horizontal flow in a saturated porous medium. Civil Engineering. Pasadena, California, California Institute of Technology: 164. Liu, P. L., H. Fernando, B. E. Jaffe (2005). "Observations by the international tsunami survey team in Sri Lanka." Science 308: 1595. Oostrom, M. (1991). Behavior of dense leachate plumes in unconfined aquifer models. Agronomy and Soils. Auburn, Auburn University: 355. Oostrom, J. S. Hayworth, J. H. Dane, and O. Guven (1992). "Behavior of dense aqueous phase leachate plumes in homogeneous porous media." Water Resources Research 28(8): 2123-2134. 99 Oswald, S. and W. Kinzelbach (2004). ?Three-dimensional physical benchmark experiments to test variable-density flow models.? Journal of Hydrology 290 (1,2) 22- 42. Paschke, N. W., J. A. Hoopes (1984). "Buoyant contaminant plumes in groundwater." Water Resources Research 20(9): 1183-1192. Penny, M.-K. L., and C. Morton (2003). "Groundwater and microbial processes of Alabama coastal plain aquifers." Water Resources Research 39(11): 1320. Schincariol, R. A., F. W. Schwartz. (1990). "An experimental investigation of variable density flow and mixing in homogeneous and heterogeneous media." Water Resources Research 26(10): 2317-2329. Schincariol, R.A., E.E. Herderick, F.W. Schwartz (1993). "On the application of image analysis to determine concentration distributions in laboratory experiments." Journal of Contaminant Hydrology 12: 197-215. Schincariol , R.A., F. W. Schwartz, C.A. Mendoza (1994). "On the generation of instabilities in variable-density flow." Water Resources Research 30(4): 913-927. Shackelford, R. S. (1991). Physical Modeling of Liquids in Groundwater. Auburn, Auburn University. Shenoy, S. (1994). An experimental study of the three-dimensional behavior of dense aqueous phase contaminant plumes in an unconfined aquifer model. Civil Engineering. Auburn, Auburn University: 198. 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.? Water Resources Research (26): 17-31. UNESCO-IOC (2006). Tsunami Glossary. Paris, UNESCO. Villholth, K.G., P. H. A., P. Jeyakumar (2005). Tsunami impacts on shallow groundwater and associated water supply on the east coast of Sri Lanka. Colombo, Sri Lanka, International Water Management Institute: 1-68. Voss, C. (1984). A finite-element simulation model for saturated-unsaturated, fluid- density-dependent groundwater-flow with energy transport or chemically-reactive single- species solute transport, US Geological Survey: 409. 100 Wang, M. P. Anderson. (1982). Introduction to groundwater modeling: finite difference and finite element methods, Freeman and Co. Zhang, Q., R. E. Volker, D. A. Lockington (2000). "Influence of seaward boundary condition on contaminant transport in unconfined coastal aquifers." Journal of Contaminant Hydrology 49: 201-215. Zheng C., G.D. Bennett (2002). Applied Contaminant Transport Modeling ? 2nd ed. New York, NY, John Wiley & Sons, Inc.