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