EXPERIMENTAL AND NUMERICAL ANALYSIS OF VARIABLE-DENSITY FLOW AND TRANSPORT SCENARIOS Except where reference is made to the work of others, the work described in this dissertation is my own or was done in the collaboration with my advisory committee. This dissertation does not include proprietary or classified information. ____________________________ Rohit Raj Goswami Certificate of approval: ___________________________ Joel G. Melville Professor Civil Engineering ___________________________ Mark O. Barnett Associate Professor Civil Engineering ___________________________ T. Prabhakar Clement, Chair Professor Civil Engineering ___________________________ Jacob H. Dane Professor Agronomy and Soils ___________________________ George T. Flowers Dean Graduate School EXPERIMENTAL AND NUMERICAL ANALYSIS OF VARIABLE-DENSITY FLOW AND TRANSPORT SCENARIOS Rohit Raj Goswami A Dissertation Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy Auburn, Alabama December 19, 2008 iii EXPERIMENTAL AND NUMERICAL ANALYSIS OF VARIABLE-DENSITY FLOW AND TRANSPORT SCENARIOS Rohit Raj Goswami Permission is granted to Auburn University to make copies of this dissertation 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 iv VITA Rohit Raj Goswami was born to Shuchindra and Chanchal Goswami in the picturesque city of Chandigarh, India. In 2002, he received a Bachelor of Engineering degree in Civil Engineering from Punjab Engineering College affliated with Panjab University, Chandigarh. He joined the Ph.D. program at Auburn University in Fall 2003 under the guidance of Dr. T. Prabhakar Clement. During the course of his studies, he met Linzy Brakefield who joined the groundwater research group as a Masters candidate at Auburn. Linzy and Rohit are now married and have a beautiful daughter, Anamaia Kaye Goswami. v DISSERTATION ABSTRACT EXPERIMENTAL AND NUMERICAL ANALYSIS OF VARIABLE-DENSITY FLOW AND TRANSPORT SCENARIOS Rohit Raj Goswami Doctor of Philosophy, December 19, 2008 (M.C.E., Auburn University, Auburn, AL, 2007) (B. E., Panjab University, Chandigarh, India, 2002) 200 Typed Pages Directed by T. Prabhakar Clement Saltwater contamination has received a lot of attention in the past few years with natural disasters such as tsunamis and hurricanes causing widespread groundwater pollution. Groundwater modeling of saline contaminant plumes is complex due to the underlying density effects. The research presented here aims to improve the understanding of variable-density processes under both laboratory and field conditions. New benchmarking datasets were developed for saltwater intrusion (saltwedge) scenarios under controlled laboratory conditions. Three types of datasets- steady-state, transient and flux- were presented for the purpose of verifying the numerical codes. Results from SEAWAT, a widely used variable-density model, were tested against the acquired datasets and shown to match the experimental results. Other variable-density models (SUTRA, SWI) were also tested against the proposed datasets and produced close matches. vi Image analysis procedure was used to measure the spatial concentration profiles in the experimental tanks during variable-density experiments. A detailed error analysis method was developed to estimate accuracy in the application of image analysis procedure. Also, the proposed error estimation method was shown to be better than the often used methods available in literature. The image analysis procedure was improved based on the results from the analysis and a more accurate concentration dataset was obtained using the improved procedure. Variable-density experiments were conducted to obtain datasets with and without instabilities. The datasets were then modeled using SEAWAT with different numerical parameters. It was shown that MOC (method of characteristics), TVD (total variation diminishing) and regular finite difference methods are not suitable to be used for modeling variable-density scenarios when instabilities can develop as they lead to wrong predictions. Finite difference in conjunction with heterogeneity was found to be the best method. Variable-density flow and transport theory was applied to study two different phenomenon- a) migration of zero-valent iron nanoparticles (ZVI); and b) fate and transport of saltwater trapped in open wells after a saltwater contamination event such as a tsunami or a hurricane. vii ACKNOWLEDGEMENTS I would like to acknowledge my advisor, Dr. T. Prabhakar Clement, for giving me the opportunity to pursue graduate studies at Auburn. Dr. Clement taught me what I know about groundwater modeling and numerical methods, the study of which sparked my interest in pursuing a PhD. He provided support whenever I hit a brick wall and guided me through the maze of doctoral research. Thanks are also due to my committee members, Dr. Joel Melville, Dr. Jacob Dane and Dr. Mark Barnett for providing valuable suggestions at various times. I would also like to thank Dr. Puneet Srivastava for being the external reviewer of my dissertation. I feel indebted to the entire academic community at Auburn University for providing me with an enriching learning experience. It was an honor to study Soil Physics from Dr. Jacob Dane who made learning a joyful experience. I learned much about class organization and conductance from Dr. Mark Barnett whose chemistry courses were the most streamlined classes I have taken. Assisting Dr. Joel Melville in teaching numerical analysis and hydraulics were some of the better memories of my stay at Auburn and I miss that job. I was fortunate to take classes with Dr. Dongye Zhao who constantly challenged the students to do better. I would specially like to mention Dr. Sushil Raj Kanel, formerly a postgraduate scholar in our research group, for giving me the opportunity to conduct interdisciplinary research. Mr. Bharath Ambale, PhD candidate in the Department of Electrical Engineering, deserves thanks for his assistance viii during our collaborative research effort on image analysis. I would also like to thank my colleagues? Mr. Vijay Loganathan, Mr. Anjani Kumar, Mr. Matthew Hogan, Mr. Justin McDonald, Mr. Gautham Jeppu, Mr. Venkatraman Srinivasan, Mr. Che-An Kuo, Mrs. Linzy Brakefield-Goswami, Ms. Noemi Sanchez-Mendez, Dr. Feng He, and Dr. Xiong Zhong - for some stimulating intellectual discussions and for creating a cordial atmosphere. Thanks are also due to Mr. Dumindu Jayasekara and Mr. Sanjeewa Manamperi for their help during the fieldwork conducted in Sri Lanka. Dr. Christian Langevin, Dr. Karen Villholth, and Dr. Tissa Illangasekare gave me their time and support and contributed to my technical and professional development for which I am extremely grateful. Acknowledgement is also due to Dr. Elena Abarca who shared her expertise on various aspects of variable-density flow and transport modeling. I would also like to thank my family for providing much needed support and understanding. Words cannot describe my gratitude for my mother Mrs. Chanchal Goswami who has always believed in me and has been a constant source of inspiration. The two leading ladies of my life- my wife Linzy Brakefield-Goswami and my daughter Anamaia Kaye Goswami have given me much joy to cherish. Last but most importantly I would like to thank God for giving me the opportunity and the strength to complete this work. ix Style manual or journal used Auburn University manuals and guides for the preparation of theses and dissertations: Organizing the manuscript ? Water Resources Research format Computer software used Microsoft? Word? & Excel? 2007; EndNote?Web; Compaq Visual Fortran 6; MATLAB? 7.0; ModelViewer 1.1; SEAWAT 2000; GW Vistas 4 x TABLE OF CONTENTS LIST OF TABLES ............................................................................................................ xv LIST OF FIGURES ........................................................................................................ xvii CHAPTER 1. INTRODUCTION AND OBJECTIVES ................................................ 1 1.1 Background ................................................................................................................. 1 1.2 Stable Density Configuration ...................................................................................... 2 1.3 Unstable Density Configuration ................................................................................. 3 1.3.1 Purely density-driven systems ............................................................................ 4 1.3.2 Advection-Coupled density-driven systems ....................................................... 5 1.4 Mixed Density Configuration ..................................................................................... 5 1.5 Objectives and Organization ....................................................................................... 7 1.5.1 Objective-1 .......................................................................................................... 7 1.5.2 Objective-2 .......................................................................................................... 7 1.5.3 Objective-3 .......................................................................................................... 8 1.5.4 Objective-4 .......................................................................................................... 8 1.5.5 Objective-5 .......................................................................................................... 8 CHAPTER 2. LABORATORY-SCALE INVESTIGATION OF SALTWATER INTRUSION DYNAMICS ............................................................................................... 10 2.1 Introduction ............................................................................................................... 10 2.2 Laboratory Methods .................................................................................................. 14 2.2.1 Experimental Set-up.......................................................................................... 14 xi 2.2.2 Data Collection Methods .................................................................................. 18 2.3 Experimental Design ................................................................................................. 21 2.4 Experimental Results ................................................................................................ 22 2.5 Modeling Methods .................................................................................................... 28 2.6 Discussion of Modeling Results ............................................................................... 30 2.6.1 Steady-State Test .............................................................................................. 30 2.6.2 Flux Test ........................................................................................................... 32 2.6.3 Transient Test .................................................................................................... 32 2.6.4 Worthiness Analysis of the New Benchmark Datasets ..................................... 35 2.7 Conclusions ............................................................................................................... 39 CHAPTER 3. ESTIMATING ERRORS IN CONCENTRATION MEASUREMENTS OBTAINED FROM IMAGE ANALYSIS ....................................................................... 41 3.1 Introduction ............................................................................................................... 41 3.2 Error Analysis Using a Theoretical Test Problem .................................................... 47 3.3 Introducing Errors in the Theoretical Test Problem ................................................. 50 3.3.1 Introducing Calibration-Relationship Error ...................................................... 50 3.3.2 Introducing Calibration-Relationship Error ...................................................... 51 3.3.3 Introducing Both Errors .................................................................................... 51 3.4 Assessing the Limitations of Using Mass Balance Dispersion Coefficient Methods for Error Estimation .......................................................................................................... 52 3.4.1 Calculation of Mass Balance and Diffusion Coefficient .................................. 52 3.4.2 Comparison of Estimated Concentration Contours .......................................... 55 3.5 Statistics Based Error Estimation Method ................................................................ 60 3.5.1 Using Statistics to Quantify Calibration-Relationship Error ............................ 60 3.5.2 Using Statistics to Quantify Experimental Errors ............................................. 61 xii 3.5.3 Estimating Total Error and Fractional Error ..................................................... 62 3.5.4 Verification of Statistics-based Error Estimation Method ................................ 63 3.6 Experimental Analysis .............................................................................................. 66 3.6.1 Laboratory Set-up ............................................................................................. 66 3.6.2 Acquiring Calibration Dataset and Image Processing ...................................... 68 3.6.3 Designing Experiments Based on Error Analysis ............................................. 71 3.7 Conducting Tracer Experiment ................................................................................. 72 3.7.1 Computation of Mass Balance Error for the Experimental Data ...................... 74 3.8 Summary and Conclusions ....................................................................................... 75 CHAPTER 4. ANALYSIS OF NUMERICAL TECHNIQUES USED FOR SIMULATING INSTABILITIES OBSERVED IN VARIABLE-DENSITY FLOW EXPERIMENTS.. ............................................................................................................. 77 4.1 Introduction ............................................................................................................... 77 4.2 Laboratory Methods .................................................................................................. 79 4.2.1 Design of Variable-density experiments ........................................................... 80 4.3 Laboratory data ......................................................................................................... 81 4.3.1 Results of Rising Plume Experiment ................................................................ 82 4.3.2 Results of sinking plume experiment ................................................................ 83 4.4 Concentration data derived from Image Analysis .................................................... 84 4.5 Numerical Analysis ................................................................................................... 86 4.5.1 Numerical Model Details .................................................................................. 87 4.5.2 Homogeneous Simulations with MOC ............................................................. 89 4.5.3 Homogenous Simulations with finite-difference .............................................. 90 4.5.4 Simulations with TVD ...................................................................................... 94 4.5.5 Simulations with Heterogeneity ........................................................................ 96 4.5.6 Finite-difference with heterogeneity ................................................................. 96 xiii 4.6 Summary ................................................................................................................... 99 CHAPTER 5. TWO-DIMENSIONAL TRANSPORT CHARACTERISTICS OF SURFACE-STABILIZED ZERO-VALENT IRON NANOPARTICLES IN POROUS MEDIA??????. .................................................................................................. 102 5.1 Introduction ............................................................................................................. 102 5.2 Materials and Methods ............................................................................................ 104 5.3 Tracer Test .............................................................................................................. 106 5.4 Preparation of INP and S-INP................................................................................. 106 5.5 Modeling Tracer and S-INP Transport Processes. .................................................. 107 5.6 Results and Discussion ........................................................................................... 108 5.6.1 Stability of INP and S-INP in Vials ................................................................ 108 5.6.2 Transport of INP and S-INP in Porous Media ................................................ 109 CHAPTER 6. ANALYSIS OF SALTWATER CONTAMINATION OF OPEN DUG WELLS IN SRI LANKA AFTER THE 2004 TSUNAMI EVENT- PHYSICAL AND NUMERICAL MODELING EFFORTS ........................................................................ 114 6.1 Introduction ............................................................................................................. 114 6.2 Small-scale laboratory experiments ........................................................................ 118 6.3 Field Experiments ................................................................................................... 118 6.4 Collection of Laboratory Data and Numerical Model Validation .......................... 124 6.4.1 Tracer experiment ........................................................................................... 126 6.4.2 Numerical Modeling of the tracer data ........................................................... 127 6.4.3 Results from tracer experiment ....................................................................... 129 6.4.4 Saltwater injection experiment ....................................................................... 130 6.4.5 Numerical modeling of saltwater experiment ................................................. 133 6.4.6 Results from the saltwater experiment ............................................................ 134 6.5 Scenario analysis using the validated numerical model ......................................... 138 6.5.1 Details of the numerical model ....................................................................... 140 xiv 6.5.2 Results from the scenario analysis .................................................................. 142 6.6 Conclusions and Recommendations ....................................................................... 148 CHAPTER 7. SUMMARY, CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK ............................................................................................................ 150 7.1 Summary ................................................................................................................. 150 7.2 Suggestions for future work .................................................................................... 151 REFERENCES? ........................................................................................................... 153 APPENDIX-1 IMAGE ANALYSIS PROCEDURE ...................................................... 171 xv LIST OF TABLES Table 2-1. Steady-state Salt Wedge Dataset a ................................................................... 26? Table 2-2. Transient Advancing Salt-wedge Dataset a ..................................................... 27? Table 2-3. Transient Receding Salt-wedge Dataset a ........................................................ 28? Table 2-4. Comparison of Measured and Model-Predicted Freshwater Flows under Steady-State Conditions a ................................................................................................... 32? Table 3-1. Parameter values used in the theoretical test problem .................................... 48? Table 3-2. Percentage errors in mass balance and diffusion coefficient values obtained with different relationships. ?? represent overestimation and underestimation respectively with respect to the true values ...................................................................... 53? Table 3-3. Errors computed by the proposed statistical-method for the theoretical test problem ............................................................................................................................. 63? xvi Table 3-4. Errors computed by the proposed statistical-method for the experimental problem ............................................................................................................................. 71? Table 3-5. Errors in mass balance calculation for all experimental datasets ................... 75? Table 6-1. Physical parameters of the injected tracer fluid ........................................... 126? Table 6-2. Location of the center of the injection well and monitoring point for the tracer experiment....................................................................................................................... 127? Table 6-3. Parameters and grid discretization used for SEAWAT ................................ 128? Table 6-4. Locations of the center of the injection well and monitoring point for saltwater experiment ....................................................................................................... 131? Table 6-5. Physical parameters of the injected saltwater ............................................... 131? Table 6-6. Grid discretization used for the grid convergence analysis .......................... 133? Table 6-7. Parameters and grid discretization used for SEAWAT ................................ 134 Table 6-8. Scenarios simulated with the SEAWAT model ........................................... 139? xvii LIST OF FIGURES Figure 1-1. Natural groundwater flow and transport in a coastal unconfined aquifer ....... 2? Figure 1-2. Saltwater upconing due to over-pumping ....................................................... 3? Figure 1-3. Purely density-driven convective system, taken from Simmons [2002] ......... 5? Figure 1-4. Mixed density configuration after a saltwater contamination event. Taken from Hogan [2006] ............................................................................................................. 6? Figure 2-1. Conceptual diagram of the experimental set-up ............................................ 14? Figure 2-2. Details of the wedge delineation procedure .................................................. 17? Figure 2-3. Observed salt wedge profiles for the steady-state and transient experiments ....................................................................................................................... 23? Figure 2-4. Computation domain and boundary conditions used in the numerical model................................................................................................................................. 25? xviii Figure 2-5. Comparison of steady-state model results against experimental data ........... 31? Figure 2-6. Comparison of model generated (0.10, 0.50 and 0.90) isochlors with the experimentally generated saltwater wedge data ............................................................... 31? Figure 2-7. Comparison of model predicted transient intruding salt wedge locations against experimental data .................................................................................................. 34? Figure 2-8. Comparison of model predicted transient receding salt wedge locations against experimental data .................................................................................................. 34? Figure 2-9. Comparison of steady-state numerical solutions obtained from coupled and uncoupled SEAWAT simulations ..................................................................................... 38? Figure 2-10. Comparison of transient saltwater intrusion solutions obtained from coupled and uncoupled SEAWAT simulations ................................................................ 39? Figure 3-1. Concentration contours or the true solution for the theoretical problem ...... 48? Figure 3-2. Calibration data and the fitted polynomial relationships .............................. 49? Figure 3-3. 5 mg/L and 20 mg/L contours for C true and C estimate, EE . Part a shows the C estimate, EE and C true with corresponding mass balance error bound. Part b shows an xix enlarged view of a section of part a. Part c presents C estimate, EE with ?trial and error? bounds of 10% and 15% for 20 mg/L and 5 mg/L contours, respectively. ...................... 57? Figure 3-4. 5 mg/L and 20 mg/L contours, C true and C estimate, CRE+EE , for linear, quadratic and cubic relationships. Only the mass balance error bounds are shown. ........................ 59? Figure 3-5. Illustration of concentration-intensity function and its relationship to standard deviation in intensity and the corresponding concentration value ..................... 62? Figure 3-6. Concentration error bounds based on statistical analysis. a presents the error bounds for linear relationship with the only CRE. b presents the error bounds for only EE. c and d show error bounds for linear and cubic relationships with effects of both CRE and EE. ............................................................................................................. 66? Figure 3-7. Layout of the experimental set-up................................................................. 67? Figure 3-8. Calibration data and the fitted intensity-concentration relationship for the experimental system .......................................................................................................... 70? Figure 3-9. Plots of different components of fractional error with respect to concentration ..................................................................................................................... 72? Figure 3-10. a) digital images from the experiment, b) filled concentration distribution xx obtained from the IA technique, c) 1 ml/l (red) and 3.5 ml/l (blue) contours. Numbers- 1,2 and 3 are for different times after injection ............................................... 74? Figure 4-1. Experimental results from the rising plume experiment at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection ........................................................ 83? Figure 4-2. Experimental results from the sinking plume experiment at- a) 0 min, b) 2 min, and c) 5 min - after injection ............................................................................. 84? Figure 4-3. Concentration profiles obtained from the rising plume experiment at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection ......................................... 86? Figure 4-4. Concentration profiles obtained from the sinking plume experiment at- a) 0 min, b) 2 min, and c) 5 min - after injection .............................................................. 86? Figure 4-5. Conceptual description of- a) rising plume, and b) sinking plume experiments ....................................................................................................................... 87? Figure 4-6. Results from MOC simulations of the rising plume experiment with homogeneous packing at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection............................................................................................................................. 90 Figure 4-7. Results from MOC simulations of the sinking plume experiment with xxi homogeneous packing at- a) 0 min, b) 2 min, and c) 5 min - after injection .................... 90? Figure 4-8. Results from upstream finite-difference simulations of the rising plume experiment with homogeneous packing at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection .............................................................................................. 92? Figure 4-9. Results from central-in-space finite-difference simulations of the rising plume experiment with homogeneous packing at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection .............................................................................................. 92? Figure 4-10. Results from upstream finite-difference simulations of the sinking plume experiment with homogeneous packing at- a) 0 min, b) 2 min, and c) 5 min - after injection............................................................................................................................. 93? Figure 4-11. Results from central-in-space finite-difference simulations of the sinking plume experiment with homogeneous packing at- a) 0 min, b) 2 min, and c) 5 min - after injection .................................................................................................................... 93? Figure 4-12. Numerical modeling results obtained with TVD for rising plume experiment with homogeneous packing at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection .............................................................................................. 95? Figure 4-13. Numerical simulation results obtained with TVD for sinking plume xxii experiment with homogeneous packing at- a) 0 min, b) 2 min, and c) 5 min .................. 95? Figure 4-14. Numerical modeling results for rising plume experiment with heterogeneous packing and upstream finite-difference method at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection ........................................................................ 98? Figure 4-15. Numerical modeling results for rising plume experiment with heterogeneous packing and central finite-difference method at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection ........................................................................ 98? Figure 4-16. Numerical modeling results for sinking plume experiment with heterogeneous packing and upstream finite-difference at- a) 0 min, b) 2 min, and c) 5 min- after injection ..................................................................................................... 99? Figure 4-17. Numerical modeling results for sinking plume experiment with heterogeneous packing and central finite-difference at- a) 0 min, b) 2 min, and c) 5 min- after injection ..................................................................................................... 99? Figure 5-1. Conceptual diagram of the flow container .................................................. 105? Figure 5-2. Vials containing INP and S-INP at various times: (a) and (b) INP and S-INP after 1 min, (c) and (d) after 10 min, (e) and (f) after 2 h, (g) INP after 2 days, and (h) S-INP after 60 days ............................................................................................ 109? xxiii Figure 5-3. Transport of tracer, pristine INP, and S-INP in the flow container ............ 110? Figure 5-4. Comparison of experimental and numerical modeling results for tracer-dye transport. Numerical results show 1 and 10% shaded contour levels. .......... 111? Figure 5-5. Comparison of experimental and numerical modeling results for S-INP transport. Numerical results show 1 and 10% shaded contour levels. ........................... 112? Figure 6-1. Conceptual diagram showing the three possible sources of saltwater infiltration into the unconfined aquifer (modified from Villholth et al. [2005]) ............ 116? Figure 6-2. Tsunami Impacted Field Site in Sri Lanka .................................................. 116? Figure 6-3. Selected field site shown on- a) the complete map of Sri Lanka and b) a map of the Kalpitiya peninsular region. Taken from Google Earth ? ........................... 120? Figure 6-4. a) Setting up field equipment and b) the final experimental set up ............ 120? Figure 6-5. Wells installed at the Field Site ................................................................... 121? Figure 6-6. Cross-sectional profile (Wells 1a and 1b) ................................................... 121? xxiv Figure 6-7. Change in salinity over time for the 14 day inject-and-pump experiment, for well 1a (injection well). Pumping events are shown by dotted lines ........................ 123? Figure 6-8. Change in salinity over time for inject-and-monitor experiment. Breakthrough profile from well 1b (monitoring well) .................................................... 124? Figure 6-9. Front view of the three-dimensional flow container ................................... 125? Figure 6-10. Experimental set-up for the tracer experiment .......................................... 127? Figure 6-11. Normalized breakthrough curve obtained from the physical tracer experiment....................................................................................................................... 129? Figure 6-12. Numerical modeling of the tracer experiment- breakthrough curves obtained at different depths ............................................................................................. 130? Figure 6-13. Conductivity probe shown here with the screen attached around the element ............................................................................................................................ 131? Figure 6-14. Setup for the saltwater injection experiment ............................................. 132? Figure 6-15. Breakthrough curves obtained from the physical experiments. ................ 132? xxv Figure 6-16. Breakthrough curves obtained from Grid # A- Very Fine. Layer 22 represents a depth of 11 cm. Similarly layers 24, 26 and 28 represent depths of 12 cm, 13 cm and 14 cm respectively. ............................................................................ 136? Figure 6-17. Breakthrough curves obtained from Grid # B- Fine. Layers 11-14 represent depths of 11- 14 cm respectively..................................................................... 136? Figure 6-18. Breakthrough curves obtained from Grid # C- Coarse. Layers 11-14 represent depths of 11- 14 cm respectively..................................................................... 137? Figure 6-19. Breakthrough curves obtained from Grid # D- Very Coarse. Layers 5-7 represent depths of 10 cm, 12 cm and 14 cm respectively. ............................................ 137? Figure 6-20. Breakthrough curves obtained from all grid discretizations at corresponding layers and comparison with experimental dataset ................................... 138? Figure 6-21. Conceptual model of a saltwater plume moving through an aquifer after a saltwater flooding event. The saltwater originates from the ground surface (orange part) and the well (red part). a. The initial situation. b. After the saltwater plume has sunk into the aquifer ...................................................................................... 140? Figure 6-22. The movement of saltwater in the aquifer (Scenario 1a): a) immediately after the tsunami, and b) after 5 days. 6 well volumes injected. .................................... 143? xxvi Figure 6-23. The movement of saltwater into the aquifer (Scenario 1b): a) immediately after the tsunami, and b) after 5 days. 12 well volumes injected. .............. 143? Figure 6-24. Effect of early pumping and subsequent re-entry of saltwater from the surface source under medium vulnerability conditions (Scenario 2) .............................. 145? Figure 6-25. Effect of early pumping under high vulnerability condition (high hydraulic conductivity and shallow well) (Scenario 3). Ingress of the surface source is very rapid..................................................................................................................... 146? Figure 6-26. Effect of late pumping (Scenario 4). The re-entry of surface-derived saltwater into the well after pumping is almost immediate ............................................ 147 1 CHAPTER 1. INTRODUCTION AND OBJECTIVES 1.1 Background Recent catastrophic events such as the 2004 tsunami and hurricane Katrina have brought the saltwater invasion problem to the forefront. Furthermore, the forecasted rise in the sea level due to global warming is expected to cause additional large-scale saltwater intrusion problems throughout the world. Saltwater can enter natural groundwater systems through various environmental transport processes. The most significant processes include saltwater intrusion from open oceans [Herbert et al., 1988; Reilly and Goodman, 1985], saltwater up-coning from natural halite deposits [Younes et al., 1999], contamination from highly saline water from landfills, and saltwater inundation during catastrophic events such as tsunamis [Illangasekare et al., 2006]. Density differences play a significant role in controlling the fate and transport of constituents in these processes. Based on the density configuration, saltwater contamination problems can be broadly classified into three distinct types: a) stable configuration, b) unstable configuration, and c) mixed systems involving both stable and unstable density configurations [Clement et al., 2004; Hogan, 2006]. 2 1.2 Stable Density Configuration The classic saltwater intrusion problem is an excellent example of a stable density configuration [Reilly and Goodman, 1985]. Conceptual model of natural groundwater flow in a coastal unconfined aquifer undergoing saltwater intrusion is shown in Figure 1- 1. During the saltwater intrusion process, the saltwater (heavier water) will always be transported beneath the freshwater (lighter water) and hence the interface between the two water bodies will be stable. Mixing across the interface is expected to be small since the heavier saltwater will remain at the bottom and the lighter freshwater will remain floating on the top of the saltwater. Therefore, in coastal aquifers, the presence of a saltwater interface will force the floating freshwater to exit via a narrow discharge zone along the beach boundary, as shown in Figure 1-1. Understanding freshwater transport processes occurring in the vicinity of a stable interface has an enormous implication for sustainable management of coastal groundwater resources worldwide. Figure 1-1 Natural groundwater flow and transport in a coastal unconfined aquifer 3 Salt upconing near pumping wells is another example of stable density configuration. As shown in Figure 1-2, in this scenario the lighter freshwater remains over heavier saltwater but the saltwater tends to rise below the pumping well and mix with the freshwater. Such uplifting of saltwater can also be caused by freshwater flowing over saltwater in a geologically pressurized aquifer units [Younes et al., 1999] and by excessive pumping in freshwater wells near the saltwater interface [Todd and Mays, 2005]. Figure 1-2 Saltwater upconing due to over-pumping 1.3 Unstable Density Configuration In unstable density configurations heavier saltwater exists on top of lighter freshwater. However, the heavier saltwater will have the natural tendency to sink into freshwater under the influence of gravity. Such tendency to penetrate and mix from 4 above can give rise to instabilities in the freshwater-saltwater interface. The unstable density configuration problems can be further sub-divided into two categories: a) purely density-driven systems and b) advection coupled density-driven systems. 1.3.1 Purely density-driven systems Purely density-driven convective systems occur in aquifer regions where the ambient groundwater flow is extremely low. The low flow allows the denser water to sink into the aquifer almost entirely under the influence of density. This type of system is likely to be encountered in very low permeability aquifers with a saline boundary on the top. Since such scenarios are rare to find in field cases, the use of this system is limited to theoretical analysis and for developing benchmark cases for testing variable-density numerical codes [Elder, 1967b; Simmons et al., 1999]. Figure 1-3 shows an example of a purely density-driven system, taken from a laboratory experiment conducted by Simmons [2002]. A saltwater source placed at the top boundary discharged high density solution into the fresher ambient water. In the absence of any ambient groundwater flow, only density differences govern the movement of denser saltwater into the freshwater. This leads to the development of instabilities or fingers in flow systems where density differences are high. Figure 1-3 presents such a system where density differences are sufficiently high to produce instabilities. 5 Figure 1-3 Purely density-driven convective system, taken from Simmons [2002] 1.3.2 Advection-Coupled density-driven systems When heavier saltwater sinks through freshwater in the presence of considerable ambient groundwater flow, the resulting system can be categorized as a mixed convective system. List [1965] studied these types of systems in detail both experimentally and theoretically, and provided a stability analysis. Coupled convective systems may or may not develop instabilities depending upon system parameters such as density difference, hydraulic gradient, permeability, etc. More recently, Schincariol [1990; 1994; 1997] ,Oostrom [1992b] and Hogan [2006] reinvestigated such systems by varying parameters such as flow conditions. Coupled convective systems are closer to real life scenarios and have direct application to many problems, such as the tsunami problem. 1.4 Mixed Density Configuration In some field scenarios, one could encounter both stable and unstable density configurations within the same aquifer. Catastrophic events such as tsunamis and 6 hurricanes are examples of such scenarios. The wave run-up from tsunamis and hurricanes can discharge large quantities of saltwater over low-lying coastal regions giving rise to an unstable density configuration [Illangasekare et al., 2006]. The coastal aquifers where such catastrophes occur have classic saltwater intrusion problem (salt- wedge) as discussed in Section 1.1 above. Therefore, in such aquifers both the unstable and stable density configurations can exist simultaneously. An experimental dataset taken from Hogan [2006] presenting a laboratory-scale simulation of the tsunami problem is shown in Figure 1-4. The figure presents the unstable density configuration with the red colored ?tsunami? saltwater discharged over colorless ambient freshwater. Green colored saltwater represents the classic saltwater intrusion scenario in the aquifer model. Figure 1-4 Mixed density configuration after a saltwater contamination event. Taken from Hogan [2006] 7 1.5 Objectives and Organization Review of literature information indicated that there are several gaps in our understanding of the invasion of stable and unstable saltwater fronts into saturated groundwater systems. This study focuses on five specific objectives and each of the objectives, which are covered in separate chapters, attempts to focus on a specific research need or hypothesis. Since the chapters are organized in journal format, the literature review pertaining to the objectives are covered at the beginning of the respective chapter. 1.5.1 Objective-1 Our review indicated that currently physical experimental datasets are not available to benchmark transient saltwater intrusion models. Therefore, our first objective was to conduct a well- controlled laboratory experiment to develop a benchmarking dataset for testing saltwater intrusion models. Details of this experimental study are discussed in Chapter 2. A summary of the benchmark problems discussed in the chapter has been published in Water Resources Research [Goswami and Clement, 2007]. 1.5.2 Objective-2 Our second objective is to develop laboratory methods to extract quantitative concentration data from digital photographs using image analysis techniques. Details of the work are summarized in Chapter 3. A summary of this analysis has been accepted for publication in Vadose Zone Journal [Goswami et al., 2008]. 8 1.5.3 Objective-3 Our third objective is to apply the image analysis techniques discussed in the previous chapter to develop a benchmarking dataset involving unstable-configuration flow scenarios with and without instabilities. The goal is to establish a novel dataset for rigorously testing the numerical capabilities of variable-density flow codes. Details of this study are summarized in Chapter 4. A summary of these benchmark problems has been submitted for publication to the journal Ground Water. 1.5.4 Objective-4 Our fourth objective is to apply variable-density flow theory to analyze transport issues associated with a remediation methodology involving stabilized iron nano-particles (S-INP). The remediation methodology was developed by Dr. Sushil Raj Kanel, former post-doctoral fellow with our research group. Since major focus of this effort was to study the effects of density on S-INP transport the author of this dissertation contributed significantly to the work especially in the experimental and numerical modeling part. Details of this study are provided in Chapter 5. A summary of this work was published in Environmental Science & Technology [Kanel et al., 2008]. 1.5.5 Objective-5 Our fifth objective is to apply variable-density flow theory to study the contamination scenarios occurred during the 2004-tsunami event. This work involved field, laboratory, and numerical studies. The field work was completed at site on the eastern coast of Sri Lanka. Details of these efforts are presented in Chapter 6. Part of this work has been used to help develop well cleanup guidelines endorsed by the World 9 Health Organization [Villholth, 2008]. In addition, information provided in this chapter is being organized to prepare a paper for the journal Ground Water. The final chapter provides summary and conclusions with some suggestions for future work. 10 CHAPTER 2. LABORATORY-SCALE INVESTIGATION OF SALTWATER INTRUSION DYNAMICS 2.1 Introduction The management of saltwater intrusion into coastal aquifers is one of the most challenging environmental management problems faced by water resources planners worldwide. The intrusion of salt water into groundwater aquifers is normally prevented by the ambient freshwater flux discharging towards the ocean. However, overexploitation of coastal aquifers has lowered groundwater levels and reduced freshwater flux, and this has led to severe saltwater intrusion problems in several metropolitan areas [USGS, 2000]. Furthermore, catastrophic events such as tsunamis and hurricanes can inject saltwater into local aquifers and contaminate large volumes of freshwater reserves [Illangasekare et al., 2006]. Therefore, understanding the mixing dynamics of saltwater within freshwater aquifer systems is an important research problem. 11 Modeling saltwater flow in freshwater systems requires numerical codes that can solve density-coupled flow and transport equations. The performance of these codes is often validated by solving a set of benchmark problems. The most commonly used benchmark problem for testing saltwater intrusion codes is the Henry problem. This problem is based on an analytical solution developed by Harold Henry in his PhD dissertation [Henry, 1960]. Henry considered a salt-wedge transport problem in a rectangular, saturated, two-dimensional, confined porous media domain. He adapted a mathematical solution developed by Poots [1958], which was originally used for modeling heat transfer processes, to solve his steady-state, saltwater intrusion problem. Over the years, the original Henry problem has undergone several revisions since numerical modelers were unable to recreate the exact analytical solution for a variety of reasons [Bues and Oltean, 2000; Croucher and OSullivan, 1995]. Segol [1994] reinvestigated Henry?s analytical expression and pointed out some errors originally made by Henry and proposed appropriate corrections. These corrections finally allowed investigators to exactly match their numerical results against the revised Henry solution [Segol, 1994; Simpson and Clement, 2004]. However, one of the limitations of using the Henry problem to benchmark density-coupled codes is that the original profiles presented by Henry are relatively insensitive to density-coupling effects [Simpson and Clement, 2003]. Simpson and Clement [2004] formulated and solved a new Henry problem to address this issue. The Henry problem is a steady-state problem hence it cannot be used to test the performance of transient solutions. Therefore analysts often use the Elder problem and/or the Hydrologic Code Intercomparison (HYDROCOIN) problem to further test 12 their numerical codes. The standard form of the Elder problem used in the groundwater literature is a theoretical problem. It was formulated by Voss and Souza [1987] by modifying the heat transport parameters used by Elder [1967b] to reflect variable density flow in a porous media system. Over the years, researchers have found that the Elder problem is an inherently unstable problem due to its density configuration. The Elder problem has also been criticized for the lack of reproducibility due to its high level of sensitivity to numerical parameters such as grid discretization level and time step size [Kolditz et al., 1998]. Other benchmark problems such as the salt dome problem or the HYDROCOIN problem, Level 1, Case5 [Inspectorate, 1986], are essentially inter- comparison of numerical solutions . Simmons et al. [1999] presented a salt lake experimental problem for testing density-dependent groundwater models. This problem involves evaporation of dense brine overlying less dense fluid; this is a hydrodynamically unstable system and hence would naturally lead to downward convection of salt fingers. Numerical solution to the problem was compared against a Hele-Shaw experiment. The authors were able to qualitatively match the experimental plume patterns by introducing some artificial noise in the concentration data that corresponded to 1% of the total salinity difference between boundary and background fluid concentrations. In the past, researchers have attempted to study saltwater intrusion processes in Hele-Shaw cells [Bear and Dagan, 1964; Mualem and Bear, 1974; Vappicha, 1975]. However, these viscous analog experiments fail to capture the true transport patterns in the presence of a porous medium. In the case of Bear and Dagan?s [1964] work, the flow measurements were relatively crude and there were some difficulties in controlling 13 the freshwater inflow. Furthermore, as pointed out by Segol [1994], it is difficult to maintain a constant interspace width between the walls of a Hele-Shaw model, and this could lead to errors in obtaining a consistent set of permeability and specific discharge values. Therefore these analog model results have not been matched precisely against numerical solutions. Zhang et al. [2001] conducted experimental and numerical investigation to study the fate and transport of a contaminant plume in the vicinity of a steady-state, saltwater boundary. Their goal was to investigate the consequence of simplifying the seaward boundary condition by neglecting the seawater density and tidal variations. They concluded that neglecting the density effects would underestimate solute discharge rates and predict incorrect solute migration pathways. Oswald and Kinzelbach [2004] presented the results of a salt pool experiment as a three-dimensional benchmark problem. This work described a benchmarking case study where one knows the exact experimental results and hence can objectively evaluate the quality of the numerical solution. However, in this study, the authors reported that the numerical simulation of this problem required a very fine grid. They also encountered some difficulties in obtaining grid converged solution. Furthermore, these experiments simulated a hypothetical salt upconing scenario that has very little similarity to the classic saltwater intrusion problem. Therefore experimental work with a more realistic seaward boundary condition is needed to evaluate the performance of saltwater intrusion models. The objective of this work was to complete a laboratory investigation to generate experimental data sets for various saltwater intrusion scenarios. Our goal was to generate quantitative data sets to describe the transport patterns of an intruding and receding salt 14 wedges under different hydraulic gradient conditions, and compare the data against numerical solutions. Figure 2-1 Conceptual diagram of the experimental set-up 2.2 Laboratory Methods 2.2.1 Experimental Set-up Experiments were conducted in a rectangular flow tank. The tank was constructed using 6 mm thick Plexiglass?. The internal dimensions of the porous media region are 53 cm (length) ? 2.7 cm (width) ? 30.5 cm (height). A conceptual diagram of the experimental setup is shown in Figure 2-1. A digital picture of the flow tank used in this study is shown in Figure 2-2a. The flow tank was divided into three distinct chambers: a central flow chamber containing the porous medium, and two constant-head chambers 15 containing saltwater and freshwater. The constant head chambers are 5 cm long and are separated from the porous media chamber by two fine screens of mesh size US # 16. Uniform silica beads of average diameter 1.1 mm, obtained from Potter Industries Inc., were used as the porous medium. Self adhesive measurement tapes were pasted on the sides and at the bottom of the tank to allow direct measurements (see Figure 2-2). The experiments were recorded using Nikon Coolpix digital camera in a high-resolution mode. The digital data allowed us to zoom-in and observe small-scale variations occurring at the millimeter scale. Constant amounts of freshwater and saltwater flow from two large constant-head controlled reservoirs were injected into the saltwater and freshwater chambers, respectively. The saltwater and freshwater heads in the model were controlled by opening the appropriate overflow outlet as shown in Figure 2-1. About 40 gallons of saltwater was prepared in a large tank by dissolving commercial salt in deionized water. In the published literature, others have used various dyes including food color [Kalejaiye and Cardoso, 2005], potassium permanganate / fluorescein [Oostrom et al., 1992b], Rhodamine WT liquid [Schincariol and Schwartz, 1990] to track the movement of dense saltwater plumes in laboratory-scale aquifer models. In this study, we added about 100 ml of nonsorbing, red food color to 40 gallons of salt solution. We performed batch sorption experiments using dye solutions of various concentrations to study the interaction of the dye with the porous medium. A HP spectrophotometer was used to monitor the concentration levels with time. The data indicated that the concentrations remained at the initial value in all the batch experiments. These sorption test results 16 indicated that the dye used in the experiment will not sorb or react with the porous medium (silica beads) used in this study. The porosity of the uniform porous medium was measured using both volumetric and gravimetric methods, and the average value of the porosity was 0.385. The density of the salt solution was measured using the ASTM 111H hydrometer, and its value was 1.026 g/ml. Multiple density measurements were made throughout the experiment to monitor variations in the saltwater density level. The saltwater solution used in all the experiments was taken from the initially prepared 40-gallon batch; this helped to avoid any minor density or color variations between the experiments. 17 Figure 2-2 Details of the wedge delineation procedure c a b de 18 2.2.2 Data Collection Methods Prior to the experiment, the tank was packed with wet porous medium under fully saturated conditions; this procedure helped us avoid air entrapment. The packing was done in layers of about 5 cm. After distributing the porous media in each layer, the tank and the porous medium were tamped to achieve homogeneous packing conditions. The hydraulic conductivity of the porous medium was estimated by setting up a uniform flow field through the system and measuring the hydraulic gradient and the corresponding volumetric discharge. The flow through the system under various gradient conditions was measured, and the average in-situ hydraulic conductivity value was subsequently calculated using the Darcy?s law. The overall conductivity of laboratory-scale aquifer models would depend on packing conditions; hence an in-situ method is the most appropriate method for estimating the average hydraulic conductivity value [Oostrom et al., 1992b]. Using the in-situ approach, the average hydraulic conductivity value of the flow tank was estimated to be 1050 (?25) m/day. A tracer test was conducted by instantaneously injecting a slug of tracer into the freshwater chamber; this allowed a long pulse of tracer to develop along the entire height of the box. The pulse was allowed to transport across the flow domain. The homogeneity of the packing was ensured by visually verifying the uniformity of the pulse as it migrated through the system. Laboratory observations indicated very little spreading, and the tracer front was relatively sharp. A tracer transport simulation with a longitudinal dispersivity value of 1 mm, which of the order of the average grain diameter, was able to predict the average spreading observed in the tracer experiment. Such small dispersivity values (of about a millimeter) have been previously reported for similar 19 uniform porous media systems [Oswald and Kinzelbach, 2004].The transverse dispersivity value was assumed to be 1/10th of the longitudinal dispersivity [Johannsen et al., 2002]. Before starting the saltwater intrusion experiment, the system was allowed to transmit freshwater from right to left at a fixed gradient condition. Excess amount of freshwater water was injected into the right-hand-side constant-head chamber, and the overflow outlet was adjusted to maintain a constant head h f = 26.70 cm. This allowed the freshwater to transmit from the right chamber to the left chamber. This transmitted water exited the system through the overflow outlet in the left constant-head chamber; and the outlet level was adjusted to maintain the head h s = 25.50 cm. All these head values were measured from the bottom of the tank. After establishing steady freshwater flow, the saltwater intrusion process was initiated by inserting the saltwater supply tube from the reservoir into the left chamber (as shown in Figure 2-1). The dense saltwater rapidly flushed the freshwater in the left chamber, and then the dense water started to invade the porous medium. Considerable mixing of saltwater and freshwater flows was observed during this initial invasion period. The system was allowed to evolve for about an hour until it reached the steady-state condition. The steady-state condition was confirmed by ensuring that there was no observable change in the location of interface and also no measurable change in the freshwater discharge transmitted through the system. Figure 2-2a shows the steady-state distribution of saltwater wedge observed under the initial steady-state condition, referred as steady state 1 (also as SS-1). Figure 2-2b shows the close-up details of the freshwater discharge zone, which is a sharp transition boundary between the fresh uncolored water in the porous media chamber and the dense 20 colored water in the saltwater chamber. Laboratory observations indicated that there was very little mixing between these two regions. As the uncolored freshwater flow approached the saltwater boundary, the flow turned upward and moved almost vertically along the mesh screen. This freshwater flux eventually floated over the saltwater in the left chamber and was flushed out through the overflow outlet. The saltwater injection rate into the left chamber (which was supplied from the external saltwater supply reservoir) was adjusted to rapidly flush the floating freshwater (or any mixed, less dense water) from the system. In our experiment, the rate of new saltwater injected into the chamber was maintained at 6.6 ml/s, and the volume of saltwater chamber was 334 ml; this yielded an average residence time of 52 seconds. This short residence time allowed the system to maintain a constant density fluid by refreshing the saltwater in the left chamber every 52 seconds. Also the net rate of freshwater discharged into the saltwater tank was relatively small, and it ranged from 0.5 to 1.4 ml/s; these values were much less than the fresh saltwater supply rate of 6.6 ml/s. This operational strategy along with our innovative outlet design mitigated the effects of mixing and dilution at the boundary and helped us establish a constant-head saltwater boundary condition. Details shown in Figure 2-2b clearly demonstrate that our experimental design was able to prevent mixing between saltwater and freshwater flows along the entire freshwater discharge zone. Figure 2-2c shows the saltwater wedge and illustrates how the scale attached to the tank was utilized to establish a rectangular grid pattern. The grid was used to make quantitative estimates of the wedge locations. Note the bottom of the porous media tank is located approximately 0.5 cm below the visible region (which is behind the scale); therefore the toe of the salt wedge extended slightly beyond the visible region that could 21 be observed in these digital photographs. Figure 2-2d and Figure 2-2e show the finer details of the mixing zone (or the dispersion zone) at the inlet and at the toe, respectively, of the salt wedge. The color variations in these figures indicate that the dispersion zone is relatively narrow and is estimated to be about 1 cm wide. Therefore, the wedge delineation line shown in these figures (which is assumed to be the 0.5 isochlor) has an error in the range of ?0.5 cm. 2.3 Experimental Design The saltwater intrusion experiments completed in this study included three phases involving distinct steady-state conditions. In the first phase, a sharp steady-state wedge (designated as SS-1) was established. Details of this initial, steady-state phase were discussed above (also shown in Figure 2-2). This steady-state condition was considered as the initial condition of our experimental design. After completing all necessary measurements, the freshwater head in the right chamber was instantaneously lowered to force a milder gradient which resulted in reduced freshwater flow. This allowed the salt wedge to advance into the freshwater system. This transient saltwater intrusion phase is referred in this manuscript as the ?advancing-front condition?. The wedge was allowed to migrate until it reached another steady-state condition, which is referred to as the ?steady-state 2 condition (SS-2)?. Flux measurements were made under this second steady-state condition. The head in the freshwater tanks was then instantaneously raised to push the wedge back toward the saltwater boundary. This transient, receding phase is referred in the manuscript as the ?receding-front condition?. The system finally reached the ?steady state 3 condition (SS-3)?. 22 Under each steady-state condition, measurements were made at several points to quantify the freshwater supply rate from the external freshwater reservoir, saltwater supply from the external saltwater reservoir, freshwater overflow from the right chamber, and the combined saltwater and freshwater overflow from the left chamber. The net amount of freshwater flux transmitted through the system was quantified by subtracting the freshwater overflow rate from the freshwater supply rate and also by subtracting the saltwater overflow rate from the saltwater supply rate. The steady-state freshwater flux values estimated from these two distinct flow balances were always within 5%. 2.4 Experimental Results The data sets collected from the saltwater intrusion experiments completed in this study are summarized in Figure 2-3. The data include: (1) a steady-state dataset (Figure 2-3 a), (2) a transient dataset under the advancing-front condition (Figure 2-3 b), and (3) a transient dataset under the receding-front condition (Figure 2-3 c). Note that the time levels shown in Figure 2-3b are measured from SS-1 until the system reached SS-2) (took 80 min to reach SS-2 from SS-1), and the time levels shown in Figure 2-3c are measured from SS-2 until the system reached SS-3 (took 60 min to reach SS-3 from SS-2) The first steady-state condition SS-1, presented in Figure 2-3a (also discussed in detail in Figure 2-2), was the initial condition used in our study. The freshwater and saltwater heads used to simulate this condition were, 26.70 cm and 25.50, respectively. Analysis of the digital data indicated that the toe of wedge was approximately located at X T = 15 cm (see Figure 1 for the definition). Note that this value was estimated at a height of 0.5 cm above the aquifer bottom since the region below this level is not visible in the photograph. At the steady-state-1 condition, the elevation of the saltwater interface 23 (Z W ) at the saltwater boundary (at X = 0) was 13 cm. The flow measurements were made at both sides of the tank to estimate the freshwater flux through the system. Using these data the net freshwater flow transmitted through the system was estimated to be 1.42 cm 3 /sec. 3a. Steady State Dataset SS-1 (0 min) SS-2 (80 min/ 0 min) SS-3 (60 min) 3b. Transient Advancing Front Dataset 5 min (from SS-1) 15 min 55 min 3c. Transient Receding Front Dataset 10 min (from SS-2) 15 min 25 min Figure 2-3 Observed salt wedge profiles for the steady-state and transient experiments 24 To initiate transient salt wedge intrusion conditions, the freshwater head was instantaneously lowered from 26.70 cm to 26.20 cm. The transient data were recorded at regular time intervals as the wedge advanced from SS-1 to SS-2. Figure 2-3b shows the transient data after 5, 15, and 55 min of salt-wedge transport. Laboratory observations indicated that when the head was lowered, the wedge rapidly invaded the system; however, the rate of invasion decreased with time. The system attained close to steady- state conditions after about 60 min, and the picture taken after 80 min of transport is reported as the second steady state data set SS-2. Detailed analysis of the SS-2 data indicated that the toe of the second steady-state salt wedge (X T ) was approximately located at 38 cm (which was estimated at the height of 0.5 cm above the bottom) away from the saltwater boundary. The elevation of the saltwater interface (Z W ) at the left boundary (at X = 0) was 22 cm. Flow measurements were made under this steady-state condition and using these data the net freshwater flow transmitted through the system was estimated to be 0.59 cm 3 /sec. To initiate the salt-wedge receding condition, the freshwater head was instantaneously raised from 26.2 cm to 26.55 cm. The transient data were recorded at regular intervals as the wedge receded from SS-2 (which was the initial condition for this receding front experiment) to the final steady-state condition SS-3. Figure 2-3c shows the observed data after 10, 15 and 25 min of transport. The system attained close to steady-state conditions in about 35 min, and the picture taken after 60 min is presented as the third steady-state dataset SS-3. Analysis of the steady-state data indicated that that the toe of the salt wedge (X T ) was approximately located at 18 cm ( measured at a height of 0.5 cm above the bottom) away from the saltwater boundary. The elevation of the 25 saltwater interface (Z W ) at the left boundary (at X = 0) was 15.5 cm. Using the measured flow data, the net freshwater flux transmitted through the system was estimated to be 1.19 cm 3 /sec. The digital data and the flow measurements collected in this study helped us assemble multiple datasets which can be used for benchmarking the performance of saltwater intrusion models. This include (1) the locations of three steady-state saltwater wedges, (2) locations of three transient saltwedges under intruding conditions, and (3) locations of three transient salt wedges under receding conditions. These datasets are summarized in Table 2-1through Table 2-3. In addition, freshwater fluxed transmitted through the system were also measured under the three steady-state conditions. A detailed summary of the physical problem along with the initial and boundary conditions used to develop these experimental datasets are presented in Figure 2-4. Figure 2-4 Computation domain and boundary conditions used in the numerical model q z = 0 q z = 0 p = ? f gh f p = ? s gh s c = 0 ?c/ ?x = 0, if v < 0 ?c/ ?z = 0 ?c/ ?z = 0 c = c s , if v > 0 x z h s =25.5 cm h f =26.2 cm (SS-2) 53 cm 26 cm Fresh water Salt water h f =26.55 cm (SS-3) h f =26.7 cm (SS-1) 26 Table 2-1 Steady-state Salt Wedge Dataset a a X is distance along the x-direction (cm); Z w is the elevation of the wedge (cm) SS-1 SS-2 SS-3 X Z w X Z w XZ w 0 13.2 0 22 0 15.5 0.2 12.5 0.4 21.5 0.6 14.5 1 11.5 1.1 20.5 1.5 13.5 1.9 10.5 1.9 19.5 2.3 12.5 2.8 9.5 3 18.5 3.3 11.5 4 8.5 4.3 17.5 4.4 10.5 5 7.5 5.7 16.5 5.5 9.5 6.3 6.5 7.4 15.5 6.6 8.5 7.5 5.5 9 14.5 7.7 7.5 8.7 4.5 10.7 13.5 8.8 6.5 10.1 3.5 12.6 12.5 10.1 5.5 11.6 2.5 14.4 11.5 11.6 4.5 13.1 1.5 16.2 10.5 13.1 3.5 15.1 0.5 18.1 9.5 14.7 2.5 19.9 8.5 16.2 1.5 21.9 7.5 18.2 0.5 23.9 6.5 26 5.5 28.3 4.5 30.5 3.5 32.9 2.5 35.2 1.5 37.6 0.5 27 Table 2-2 Transient Advancing Salt-wedge Dataset a 5 min 15 min 55 min X Z w X Z w XZ w 0 19.3 0 20.6 0 21.6 0.3 18.5 0.6 19.5 0.6 20.5 0.9 17.5 1.4 18.5 1.4 19.5 1.5 16.5 2.2 17.5 2.4 18.5 2.2 15.5 3.2 16.5 3.7 17.5 3 14.5 4.4 15.5 5.3 16.5 4.1 13.5 5.6 14.5 7 15.5 5.1 12.5 6.9 13.5 8.5 14.5 6.2 11.5 8.2 12.5 10.1 13.5 7.4 10.5 9.4 11.5 11.7 12.5 8.6 9.5 10.7 10.5 13.4 11.5 9.7 8.5 12.2 9.5 15.2 10.5 11 7.5 13.7 8.5 17.2 9.5 12.3 6.5 15.3 7.5 19.1 8.5 13.7 5.5 17 6.5 21.2 7.5 15.2 4.5 18.7 5.5 23.3 6.5 16.8 3.5 20.6 4.5 25.4 5.5 18.4 2.5 22.6 3.5 27.5 4.5 20.1 1.5 24.8 2.5 29.7 3.5 21.5 0.5 27 1.5 32 2.5 29 0.5 34.5 1.5 35.9 1 a X is distance along the x-direction (cm); Z w is the elevation of the wedge (cm) 28 Table 2-3 Transient Receding Salt-wedge Dataset a 10 min 15 min 25 min X Z w X Z w XZ w 0 18.8 0 18.3 0 16.7 0.7 17.5 0.2 17.5 0.6 15.5 1.6 16.5 0.7 16.5 1.5 14.5 2.5 15.5 1.6 15.5 2.5 13.5 3.6 14.5 2.6 14.5 3.6 12.5 4.8 13.5 3.7 13.5 4.7 11.5 6 12.5 4.8 12.5 5.8 10.5 7.5 11.5 6 11.5 6.9 9.5 9 10.5 7.35 10.5 8 8.5 10.5 9.5 8.65 9.5 9.1 7.5 12.1 8.5 10 8.5 10.2 6.5 13.8 7.5 11.4 7.5 11.5 5.5 15.4 6.5 12.9 6.5 13 4.5 17.1 5.5 14.4 5.5 14.5 3.5 18.9 4.5 15.9 4.5 16 2.5 20.7 3.5 17.45 3.5 17.7 1.5 22.1 2.5 19.15 2.5 19.4 0.5 24.7 1.5 21 1.5 26.6 0.5 23.1 0.5 a X is distance along the x-direction (cm); Z w is the elevation of the wedge (cm) 2.5 Modeling Methods The finite-difference model SEAWAT [Langevin et al., 2003] was used to simulate the steady-state and transient experiments completed in this study. In addition, we also used the finite element model SUTRA [Voss and Souza, 1987] and a modified version of the MODFLOW model with the sharp interface package Sea Water Intrusion (SWI) [Bakker and Schaars, 2003] to simulate the experiment results. The goal of the modeling exercise was to test whether the experimental data are consistent with the predictions made by these widely used numerical models. In section below, however, 29 we present the results obtained from SEAWAT. Simulation results from other numerical models results were almost identical to the SEAWAT results. The numerical description of experimental setup involves a rectangular domain of 53 cm ? 26 cm, as shown in Figure 2-4. Constant-head boundary conditions were forced at both fresh and saltwater sides. No flow boundary conditions were forced at the top and bottom of the computational domain. The location of the salt wedge was delineated by interpolating the 50% saltwater concentration points predicted by the numerical model. All the simulations reported in this paper were completed using the SEAWAT code. A uniform grid of dimensions ?x = ?z = 0.5 cm was used in these simulations. Convergence tests were performed using 1 cm and 0.25 cm grids. The results of our convergence study indicated that the 0.5 isochlor predicted on various grids were almost identical. However, the width of the dispersion zone was sensitive to grid resolution. The average width of the dispersion zone predicted using the 0.5 cm mesh was about a centimeter, which was close to the width observed in our experiment. The initial condition used was a quiescent aquifer. The freshwater and saltwater boundary conditions were forced upon this system, and the system was allowed to reach the first steady-state condition. The SEAWAT simulations employed the fully-implicit finite difference package with the highly accurate TVD advection package. The hydraulic head tolerance level was set at 10 -4 cm, concentration tolerance level was 10 -8 g/L, the density-concentration slope was set at 700, salt concentration was set at 0.0371 g/mL, and the value of time step ?t was set at 1 sec. All other numerical parameters were at the default values set within the GW Vistas model of the domain. 30 2.6 Discussion of Modeling Results In this section we compare our experimental data against model predictions and use the results to establish three distinct tests for evaluating the performance of saltwater intrusion models. In the first test, identified as the steady-state test, we compare the saltwater interface locations predicted by the numerical model against three steady-state datasets. In the second test, identified as the flux test, we compare the steady state fluxes estimated by the numerical model against the measured steady-state fluxes. In the third test, identified as the transient test, we compare the transient numerical model results against two sets of experimental data collected under transient conditions. 2.6.1 Steady-State Test To test the code performance under steady-state conditions, we compare the salt wedge locations (0.5 isochlors) predicted by SEAWAT against the experimental observations in Figure 2-5. The results show an excellent match between the model results and the experimental data. The SEAWAT simulations exhibited some minor differences at the exit boundary, where there was a sharp transition between the freshwater and saltwater regions. Figure 2-6 compares the steady-state 2 data along with the 0.1, 0.5, and 0.9 isochlors predicted by the numerical model. This figure shows that the width of the dispersion zone predicted by the numerical model was about 1 cm, which was similar to the level of dispersion observed in our experiments (see Figure 2-2). All the experimental data points for 0.5 isochlor are distributed close to this narrow 1 cm dispersion zone. 31 0 5 10 15 20 25 30 0 5 10 15 20 25 30 35 40 45 50 X Value (cm) Z Va lue (cm) Steady State 1 Steady State 2 Steady State 3 SEAWAT Figure 2-5 Comparison of steady-state model results against experimental data 0 5 10 15 20 25 30 0 5 10 15 20 25 30 35 40 45 50 X Value (cm) Z Value (cm) Steady State 2 0.10 Isochlor 0.50 Isochlor 0.90 Isochlor Figure 2-6 Comparison of model generated (0.10, 0.50 and 0.90) isochlors with the experimentally generated saltwater wedge data 32 Table 2-4 Comparison of Measured and Model-Predicted Freshwater Flows under Steady-State Conditions a SS-1 SS-2 SS-3 Experimental Data 1.42 0.59 1.19 SEAWAT 4.46 0.59 1.13 Q f (No Salt Wedge) 1.94 1.12 1.69 a Flow rates are reported in cubic centimeters per second 2.6.2 Flux Test In Table 2-4, we compare the model-predicted freshwater flows against the measured flux data. For reference purposes, we have also provided the flow rate transmitted in an equivalent freshwater system (Q f ) at different gradient conditions. These values were computed using Darcy?s law after ignoring the density effects. The results shown in Table 2-4 demonstrate that the numerical models closely predicted the measured flux values. As expected, the net freshwater flow transmitted in the presence of a salt wedge is less than the flow values computed for an equivalent freshwater flow system. 2.6.3 Transient Test As part of the transient test, we first compare the experimental data from our advancing salt wedge experiment against numerical predictions. These results are presented in Figure 2-7. The advancing salt wedge simulation employed the output of steady-state-1 (SS-1) as the initial condition and the wedge was allowed to intrude until 33 the second steady condition was reached after 80 minutes of simulation. The figure shows that SEAWAT was able to closely predict the location of the invading front. In Figure 2-8 we compare the model results against the laboratory dataset collected for the receding salt wedge. The numerical study employed the steady-state 2 (SS-2) profile as initial condition, and the wedge was forced to move back until the steady-state 3 (SS-3) condition was reached after 60 minutes of simulation. The results shown in the figure indicate a reasonable match between model predications and the experimental data. However, SEAWAT predictions for early time periods tend to lag behind the data. While conducting the receding saltwedge experiment, we did encounter some difficulties in maintaining a constant freshwater head level during the initial phase. This resulted in minor head variations during the first 5 min of this particular experiment. Some of the observed initial lags could perhaps be related to this experimental problem. The model predictions, however, did reproduce the observed trends, and furthermore, the predictions matched the data very well at later times. To the best of our knowledge, this is the first time experimental data obtained from a laboratory-scale aquifer model are presented for a Henry-type problem, and the data are successfully compared against model-predicted saltwedge profiles and flux values. 34 0 5 10 15 20 25 30 0 5 10 15 20 25 30 35 40 45 50 X Value (cm) Z Va l u e ( c m ) 5 min 15 min 55 min SEAWAT Figure 2-7 Comparison of model predicted transient intruding salt wedge locations against experimental data 0 5 10 15 20 25 30 0 5 10 15 20 25 30 35 40 45 50 X Value (cm) Z Va lu e ( c m ) 10 min 15 min 25 min SEAWAT Figure 2-8 Comparison of model predicted transient receding salt wedge locations against experimental data 35 2.6.4 Worthiness Analysis of the New Benchmark Datasets Simpson and Clement [2003] proposed a coupled versus uncoupled strategy to evaluate the worthiness of benchmark problems used for testing density-coupled flow models. They pointed out that in a typical density-coupled flow problem, a dense front will be transported into a less dense region due to the forced convection generated from the boundary forcing and due to the free convection generated from internal density variations. A standard numerical solution to the density-coupled flow problem (identified as the ?coupled solution?) will account for both of these transport processes. Simpson and Clement [2003] devised an ?uncoupled solution? by disabling the density-coupling step between the flow and transport solutions. This uncoupled solution can be conceptually viewed as a method to isolate the transport due to internal free convection caused by density-dependent flow processes from the transport due to forced convection caused by external forcing (boundary conditions). By directly comparing the coupled and uncoupled solutions, one could quantify the sensitivity of a benchmark problem to density-coupling effects. Simpson and Clement [2003] used this analysis to evaluate the worthiness of Henry and Elder problems and concluded that the original Henry problem is a relatively less worthy benchmark for testing the performance of density-coupled flow simulators. Simpson and Clement (2004) later used this approach to test and develop a more worthy Henry problem. More recently, Abarca et al. (2006) compared coupled and uncoupled solutions to analyze diffusive and dispersive versions of Henry problems and concluded that solutions to the dispersive problem are better benchmarks. Dentz et al. (2006) completed a detailed analytical analysis of Simpson and Clement?s [2003] uncoupled approach (identified in their study as the pseudocoupled model). They 36 examined its accuracy and validity for modeling a range of variable density flow and transport conditions. Here we employed Simpson and Clement?s [2003] coupled/uncoupled strategy to analyze the worthiness of the benchmark problems developed in this study. The numerical code SEAWAT was used in this analysis. The coupled solution results presented in this section are identical to the standard SEAWAT simulation results presented in the previous section. Generation of the uncoupled solution results, however, required some modifications to SEAWAT input files. The modifications were as follows: (1) we set the density slope parameter in the model to zero to nullify the density-coupling effect; and (2) we computed the equivalent freshwater head distribution in the saltwater chamber and forced these values at the saltwater boundary nodes. The method used for computing the equivalent freshwater head distribution of a column of salt water is described by Langevin et al. [2003]. The above two modifications allowed SEAWAT to solve the problem in the uncoupled mode where the density variations are ignored within the flow domain but not at the boundaries. We first completed the worthiness analysis for the steady-state experimental system. Figure 2-9 compares the three steady-state profiles generated using the coupled and uncoupled SEAWAT simulations. The difference between the predicted extent of the salt-wedge toe penetration lengths (X T ) estimated from the coupled and uncoupled solutions are approximately 9, 29 and 12 cm for SS-1, SS-2 and SS-3, respectively. These results indicate that when the freshwater flow is increased, the profiles become relatively insensitive to the density coupling effects. This is because when the flow (or hydraulic gradient) is increased, the profiles become relatively less sensitive to density- 37 coupling effects; these results are consistent with the theoretical scaling arguments made by Simpson and Clement (2004). When the freshwater flow (or hydraulic gradient) is high, the transport in the system is dominated by forced convection driven by the boundary conditions. Hence the difference between coupled and uncoupled profiles under high flow conditions are small, as indicated in SS-1 and SS-3 results. Therefore, the low freshwater flow dataset (SS-2) is the worthiest dataset for testing steady-state saltwater intrusion models. Finally, we completed the coupled/uncoupled analysis to test the worthiness of the transient saltwater intrusion experimental problems. Figure 2-10 compares the transient results for the intruding salt wedge generated using coupled and uncoupled SEAWAT simulations. The differences between the predicted extent of the toe penetration lengths (X T ) simulated by the two solutions are approximately 11, 17, and 28 cm at 5, 15, and 55 min, respectively. The figures also show that the overall differences between the coupled and uncoupled profiles increase with increase in the transport time. These results clearly indicate that the transient simulations are highly sensitive to density-coupling effects. It is also interesting to note that the uncoupled simulations, which essentially describe the pressure propagation (or forced convection) process, reached close to steady-state conditions within 15 minutes. However, the true coupled solution to this problem, which accounted for both forced and free convection, required over an hour to reach steady- state. These time estimates indicate that in this system the timescale required for the boundary-controlled forced convection processes to reach steady-state is about 15 min, and the timescale required for the internal free convection processes to reach stead-state is about an hour. This mismatch between the timescales associated with the forced and 38 free convection processes amplifies the overall sensitivity of the problem to density variations. Hence this transient problem is a highly worthy problem for testing density- coupled flow codes. 0 5 10 15 20 25 30 0 5 10 15 20 25 30 35 40 45 50 X Value (cm) Z Valu e ( c m) Coupled SS-1 Coupled SS-2 Coupled SS-3 Uncoupled SS-1 Uncoupled SS-2 Uncoupled SS-3 Figure 2-9 Comparison of steady-state numerical solutions obtained from coupled and uncoupled SEAWAT simulations 39 0 5 10 15 20 25 30 0 5 10 15 20 25 30 35 40 45 50 X Value (cm) Z V a l u e (cm ) Coupled 5 min Coupled 15 min Coupled 55 min Uncoupled 5 min Uncoupled 15 min Uncoupled 55 min Figure 2-10 Comparison of transient saltwater intrusion solutions obtained from coupled and uncoupled SEAWAT simulations 2.7 Conclusions The steady-state Henry solution, which was developed in 1960, has almost exclusively served as the benchmark for testing saltwater intrusion models for over four decades. In this work, we provide a set of new benchmark datasets for testing saltwater intrusion models. This new benchmark is a more robust alternative to the original Henry problem since it considers saltwater transport under both steady-state and transient conditions. We present multiple experimental datasets that include a steady-state wedge dataset, a steady-state flux dataset, and two types of transient salt-wedge datasets. We used the saltwater intrusion model SEAWAT to successfully simulate these datasets. We completed a worthiness analysis to quantify the sensitivity of the system to density- 40 coupling effects. The results of this analysis indicated that the transient datasets are more worthy benchmarks for testing the performance of density-coupled, numerical formulations. The experimental setup described in this study provides a novel approach for simulating the saltwater boundary condition in a laboratory-scale saltwater intrusion model. The results demonstrate that the mixing of saltwater and freshwater within the saltwater reservoir was minimal, and hence the system was able to simulate the constant- head saltwater boundary condition. The experimental results also show that the transition zone between the saltwater and freshwater regions is sharp under both steady-state and transient transport conditions. In the published literature, researchers have provided numerical solutions to Henry-type problems under transient conditions [e.g. Frind. 1982; Simpson and Clement, 2004]. However, none of the transient numerical solutions have been compared against experimental observations or against analytical solutions. Therefore the transient saltwater intrusion dataset is an important contribution. Furthermore, other commonly used transient, density-coupled flow problems, such as the Elder problem or the salt lake problem, involve an unstable interface that evolves to produce complex fingering patterns. Numerical solutions to these problems yield nonunique results, and hence it is difficult to solve these problems in a consistent manner. On the contrary, the transient, density-coupled flow problem considered in this study involves transport of a stable interface which can yield a unique, stable solution. Therefore the new transient dataset provides a better alternative for testing the transient behavior of density-dependent flow codes. 41 CHAPTER 3. ESTIMATING ERRORS IN CONCENTRATION MEASUREMENTS OBTAINED FROM IMAGE ANALYSIS 3.1 Introduction In the published literature, both qualitative and quantitative methods have been used to study the behavior of porous media systems. Qualitative methods have largely consisted of visual observations of the plume movement [Goswami and Clement, 2007]. Quantitative methods have varied from sophisticated non-invasive techniques such as gamma radiation [Oostrom et al., 1992a] and magnetic resonance imaging (MRI) [Oswald and Kinzelbach, 2004], and invasive techniques that use fluid extraction from ports and pressure transducers. Equipment needed to set up and operate techniques such as gamma radiation and MRI are expensive. Also, these techniques have several limitations. For example, gamma radiation is limited to capturing single-point measurements at a given time and therefore cannot measure transient phenomenon at a large spatial scale. Similarly, MRI cannot be used for large physical models as the available equipment can accommodate only small experimental set-ups. In the past, experimentalists used sampling ports extensively for measuring in situ system properties [Silliman et al., 1998]. However, these flow systems are often sensitive to perturbations that can be potentially triggered by the ports. Therefore, scientists studying flow experiments have started utilizing non-invasive techniques for measuring the concentration field. 42 One of the rapidly evolving non-invasive techniques is image analysis (IA). IA records an optical image property (such as pixel intensity, hue, color, etc.) and relates it to a system property such as concentration or water content. Recent advances in the field of electronics and computing have provided more readily available equipments such as cameras, scanners and computers to acquire and process digital images in any laboratory setting. In the past decades, various researchers have demonstrated the application of different types of IA methods for subsurface experiments. These IA methods generally employ the following steps: a. A calibration dataset is obtained to relate the recorded property of an image (e.g., intensity) with a known system property (e.g., concentration). The primary properties obtained from a digital image are the pixel intensity values for red, green and blue colors, commonly known as the RGB data. The RGB data can be used to determine various other image properties such as hue, saturation etc. b. The calibration dataset is processed to reduce errors. Error reduction is accomplished using various types of image processing techniques and software. c. Regression analysis is conducted on the calibration dataset to obtain a relationship between the system property (e.g., concentration) and the acquired image property (e.g., intensity). This relationship will be referred to as ?calibration-relationship? in this work. Since the focus of this study is to measure concentrations, we only consider problems associated with concentration-intensity relationships. 43 d. The intended physical experiment is conducted and images (also known as experimental datasets) are obtained. These images are processed using the same error reduction techniques that were used for processing the calibration data. e. System properties (such as concentration contours) are estimated from the experimental datasets using the calibration-relationship. The general steps outlined above are followed in all IA applications. However, the specific methods used for conducting IA may slightly vary due to the differences in experimental setups and data acquisition procedures. Schincariol and Schwartz [1990] and Schincariol et al. [1993] developed their IA technique to measure concentrations based on obtained pixel intensity and optical density of the images reflected from the experimental tank. This process is referred to as the light reflection method (LRM) [Oostrom et al., 2007]. Swartz and Scwartz [1998] developed a LRM procedure, to obtain concentration data for their experiments, based on the technique developed by Schincariol et al. [1993]. A major assumption in the LRM method is that concentration is distributed uniformly across the width of the experimental tank. Oostrom et al. [2007] note that such an assumption may not be correct for flow experiments. Therefore, light transmission method (LTM) may be a better method than LRM. In LTM, the properties of the image transmitted through the experimental tank are obtained thereby capturing an averaged profile of the system property. Tidwell and Glass [1994] measured saturation profiles in unsaturated flow experiments by applying an LTM procedure. Darnault et al. [1998] used hue values obtained from their LTM process to estimate water and oil contents in two-phase flow experiments. LTM process was also used by Niemet and Selker [2001] to obtain different relationships between fluid 44 saturation and transmitted light intensity. Fluorescent methods have also been employed in IA studies. Mainhagu et al. [2007] used laser induced fluorescence (LIF) as the optical technique to measure fluid concentrations in a Hele-Shaw cell, but they only report selected vertical profiles of measured concentrations instead of complete spatial two- dimensional measurements. The authors state that LIF is a better technique than LTM and attribute the improvements to higher pixel resolution and color visualization, both of which are properties of the camera and not of the LIF technique. Huang et al. [2002] used UV light and fluorescent dye to obtain concentration distributions in a saturated flow experiment. The two major components of an IA technique are: (a) obtaining and processing the digital images, and (b) determining calibration-relationship between the measured system property and the image property (e.g., concentration-intensity relationship). Many sources of errors may be introduced while obtaining and processing digital images. The most common source of error is lighting non-uniformity, which can be both spatial and temporal [Detwiler et al., 1999; Tidwell and Glass, 1994]. Porous media experiments are also subject to uneven light transmittance, as pointed out by Bridge et al. [2006]. There may also be errors in image acquisition and processing due to hardware and software problems [Hansen et al., 2006]. Evaluation of the best-fit function for the calibration data would also involve errors; this is because, the functional relationship is unknown a priori and is simply selected based on regression statistics. In this work, we categorize all these errors into two categories: 45 1. Calibration-Relationship Error (CRE), which is the error introduced by the approximate relationship used to correlate the image property and measured system property 2. Experimental Error (EE), which is the error in the image property (e.g., intensity) due to various experimental factors such as lighting non-uniformity, light scattering, hardware problems, etc. Both CRE and EE must be quantified for the selected IA technique to determine the reliability of the technique and to determine the error associated with the measurements. In the published literature, researchers have primarily used three methods to verify their IA techniques and to quantify the errors. The first method involves comparison of the measurements obtained from IA with those obtained from analytical solutions and/or from other experimental methods (X-Ray attenuation, extraction ports etc). The second method involves comparison of global mass balance errors in the measurements made by IA. The third method is based on calculating dispersion coefficients from the measurements made by IA and comparing them with generally accepted values reported in literature. Tidwell and Glass [1994] compared results from IA with those obtained from X- Ray absorption, gravimetric method and gamma densiometry for saturation experiments. Darnault et al. [2001] compared fluid content measurements made by IA with those obtained from X-Ray synchrotron to validate the IA method. Swartz and Schwartz [1998] compared the concentrations obtained from IA to those obtained from extracted samples and estimated a generic underestimation error of 7%. However, one can clearly observe from their plotted results that both underestimation and overestimation exist in 46 their IA data. Oltean et al. [2004] used LTM in a Hele-Shaw cell and presented a thorough analytical discussion of the optical behavior of their experimental system. However, their results lack detailed explanation of the reported differences between theoretically estimated concentrations and those calculated based on the output signal from the digital images. Although the authors stated that higher concentration measurements had more errors, which were attributed to lighting non-uniformity, the errors were not quantified. Niemet and Selker [2001] measured saturation levels in an unsaturated-flow experiment and quantified errors using the standard deviation of transmitted light intensity. The reported errors in saturation were based on the root mean square error (RMSE) in mass balance. Huang et al. [2002] estimated errors based on deviations of their experimental results from an analytical solution. Several studies have used mass balance error as a general test for verifying the consistency between the laboratory datasets and the computed concentrations [Bridge et al., 2006; Darnault et al., 1998; Tidwell and Glass, 1994]. McNeil et al. [2006] compared longitudinal and transverse dispersion coefficients obtained from application of their IA process with those obtained from analytical methods and also with those reported in the published literature. All available error quantification methods have certain limitations. For example, comparison of measurements obtained by IA against analytical solutions is limited since analytical solutions are not available for most experimental conditions. Calculation of mass balance and dispersion coefficients requires use of the entire range of measured values. Consequently, these comparisons might indicate the presence global errors, but they cannot be used to quantify local error associated with a specific concentration value. 47 However, no one has clearly demonstrated these limitations. This is because all previous analyses were completed using experimental data for which the exact error structure and the ?true? solution are not available. Hence, the first objective of this study is to develop a theoretical test dataset that have known an error structure and a true solution. The dataset will be used to quantify the limitations of existing error analysis methods. Our second objective is to develop a general error quantification method to estimate the local error in a selected system property (e.g., a specific concentration contour value) and demonstrate its use by using the proposed theoretical dataset and by using a new experimental dataset. This chapter is organized into four sections following the present section. In the first section we develop a theoretical test problem and use it to demonstrate the shortcomings of the commonly used mass-balance and dispersion-coefficient methods to measure errors. In the second section, a statistics-based error estimation method is proposed and the theoretical problem was used to demonstrate the robustness of the method. In the third section, we provide details of our laboratory setup and IA technique along with an experimental dataset. In the fourth section, we describe the application of the error estimation process to the experimental dataset. 3.2 Error Analysis Using a Theoretical Test Problem We first develop a theoretical test problem, which will be used to evaluate the robustness of mass balance and dispersion coefficient comparisons for error estimation. The test problem considers diffusion of a point source in a two-dimensional domain. The following analytical solution to the problem was used to simulate the concentration profiles [Fischer et al., 1979]: 48 ? ? ? ? ? ? ? ? ??= tD z tD x DDt M zxC zx zx 44 exp 4 ),( 22 ? (1) where C(x,z) is the concentration at a given point (x,z) in the domain, t is the time from the start of point-injection, M is the mass injected in the domain, and D x and D z are the diffusion coefficients in x and z directions respectively. Values of parameters used for calculating the concentration distribution are summarized in Table 3-1. The solution domain was set to 250 mm ? 250 mm in the x and z directions, respectively. The predicted plume concentration distribution is given in Figure 3-1. As shown in the figure, the concentrations of the plume vary within the range of 0-30 mg/l. Table 3-1 Parameter values used in the theoretical test problem Parameter M (gm) t (sec) D x (m 2 /sec) D z (m 2 /sec) Value 100 1800 4.00?10-7 4.00?10-8 Figure 3-1 Concentration contours or the true solution for the theoretical problem 49 The analytical concentrations obtained were converted into pixel intensity values using the following concentration-intensity relationship: ()() ii BCAI ??= exp1 (2) where I i is the pixel intensity corresponding to a given concentration value C i , The value of A and B were set to 255 and 0.04 respectively. Parameters A and B were chosen to obtain a smooth relationship between concentrations of 0 mg/L to 30 mg/L and the corresponding pixel intensity range of 0 to 255 (8 bit image intensity range). The concentration-intensity relationship given above will be referred to as the ?true- relationship? in this analysis. A hypothetical ?calibration dataset? was generated by computing a set of intensity values using Equation (2) for a selected set of concentration values. These calibration data points are shown in Figure 3-2 (solid squares). Figure 3-2 Calibration data and the fitted polynomial relationships 50 3.3 Introducing Errors in the Theoretical Test Problem 3.3.1 Introducing Calibration-Relationship Error We assume the structure of the concentration-intensity relationship is unknown and use a set of polynomial functions to fit the hypothetical calibration dataset; this introduced calibration relationship error (CRE) into our test problem. In the literature, polynomials have been used to obtain the calibration relationships [Huang et al., 2002; Persson, 2005a; b]. These studies used R 2 values to select the best-fit calibration- relationships. Huang et al. [2002] and Swartz and Schwartz [1998] obtained a linear calibration-relationship with R 2 values of 0.978 and 0.98, respectively. Yu and Schwartz [1999] selected a third-order polynomial calibration-relationship and report the R 2 values as 0.98. Niemet and Selker [2001] used a fourth-order calibration-relationship with the R 2 value of 0.999. In this study, we fitted linear, quadratic and cubic relationships to our hypothetical calibration dataset. The values of best-fit model parameters for selected polynomials and the corresponding R 2 values are summarized below: 973.0;1555.0: 2 == RICLinear (3 a) 998.0;0722.000053333.0: 22 =+= RIICQuadratic (3 b) 999.0;1069.010468.7105359.2: 22536 =+???= ?? RIIICCubic (3 c) The fitted relationships (along with the calibration data points) are shown in Figure 3-2. The estimated R 2 values for all the fitted curves were within the range of R 2 values reported in the literature. The fitted polynomial relationships (Equation 3 a-c) were then used to estimate concentration distributions. Since none of these fitted relationships is the exact function, 51 the estimated concentrations will have some error, which is referred in this text as CRE. The concentrations obtained from these calibration-error affected datasets will be identified as C estimate, CRE . 3.3.2 Introducing Calibration-Relationship Error Experimental error (EE) is generally associated with the randomly distributed errors in the acquired image property [Tidwell and Glass, 1994]. Such random errors introduce speckle-like noise in pixel intensities [Jain, 1988]. Hansen et al. [2006] suggest that the MATLAB? function randn can be used for simulating noisy images. We introduced speckle-type noise by using the MATLAB? function randn with a standard deviation of 2% of the maximum intensity level (approximately 5 pixel intensities). The simulated noise was added to the intensity profile corresponding to the theoretical concentration data, thereby generating a ?noisy intensity? profile. We then used the true relationship (given by Equation 2) to estimate the noise-influence on concentration distributions. Because of the introduced noise, the estimated concentrations will have some error, which is referred in this text as EE. The concentrations obtained from this noisy dataset will be identified as C estimate, EE . 3.3.3 Introducing Both Errors In the above sections, we generated two datasets that either have CRE or EE. In real experimental datasets, both EE and CRE will exist together. To reflect such real case scenarios, we generated datasets that introduced both errors simultaneously by combining the use of the fitted polynomial relationships and the noisy intensity profile. Since the 52 resulting concentration profiles will be affected by both CRE and EE, they will be identified as C estimate, CRE+EE . In summary, we simulated seven scenarios by introducing various combinations of known errors. Three scenarios involving CRE were simulated using the three fitted polynomial relationships; one scenario involving EE was simulated using the noisy intensity profile; and finally, three scenarios involving both CRE and EE were simulated using the three fitted relationships and the noisy intensity profile. 3.4 Assessing the Limitations of Using Mass Balance Dispersion Coefficient Methods for Error Estimation 3.4.1 Calculation of Mass Balance and Diffusion Coefficient The estimated concentration datasets obtained after introducing various types of errors were processed to calculate the values of mass balance error and diffusion coefficient. These values were calculated using the moment analysis method [Fischer et al., 1979] : t D i i 2 2 ? = , wher, (4 a) i i i i M M ?? ?= ,0 ,22 (4 b) i i i M M ,0 ,1 =? (4 c) ??? = dVCM tzxi ,,,0 (4 d) ??? = dVCiM tzxi ,,,1 (4 e) 53 ??? = dVCiM tzxi ,, 2 ,1 (4 f) The transcript i represents both x and z directions. D i is the diffusion coefficient; ? i 2 is the variance of the concentration distribution; ? i is the location of the center of mass in the direction i; M 0,i ,M 1,i and M 2,i are the zeroth, first and second moments of the concentration distribution, respectively; and M 0,i provides an estimate of the total mass in the system. Table 3-2 Percentage errors in mass balance and diffusion coefficient values obtained with different relationships. ?? represent overestimation and underestimation respectively with respect to the true values Scenario / Relationship Used Mass Balance Error (%) Estimated Diffusion Coefficient (m 2 /s) Error in Diffusion Coefficient (%) D x D z D x D z Linear 16.48 ? 4.57?10-7 4.61?10-8 14.47 ? 15.37 ? Quadratic 02.72 ? 3.70?10-7 3.72?10-8 07.33 ? 06.83 ? Cubic 00.14 ? 3.99?10-7 4.02?10-8 00.14 ? 00.57? Noise- Exp. Error 00.26? 3.99?10-7 4.79?10-8 00.16 ? 19.77? Linear + Noise 16.56 ? 4.59?10-7 5.12?10-8 14.89 ? 28.11 ? Quadratic + Noise 01.85 ? 3.81?10-7 5.26?10-8 04.68 ? 31.69 ? Cubic + Noise 00.17 ? 3.99?10-7 4.26?10-8 00.03 ? 06.54? 54 The values of the estimated mass balance error and diffusion coefficient are summarized in Table 3-2. The symbols ?? are used to indicate overestimation or underestimation, with respect to the true values. The influences of CRE are reflected in the values reported in the first three rows, EE in the fourth row, and the effects of both errors in the last three rows. From these values, it appears that the quadratic and cubic relationships provide lower error than the linear relationship when only CRE is present. On the other hand, introduction of EE has no effect on mass balance and we get near perfect mass balance for the noisy intensity profile. This indicates that random noises in the dataset will not be quantified by the mass balance error. In terms of overestimation and underestimation of values, no definite correlation can be observed between the estimated value of mass balance error and diffusion coefficient values. There appears to be a general trend that indicates that when the mass balance errors decreased the error in the diffusion coefficients also decreased. However, there are some exceptions to this general trend; for example, for the case of cubic relationship overestimation in mass result in underestimation of diffusion coefficient value. Overall, a consistent quantitative estimate of error cannot be made by applying mass-balance and diffusion-coefficient analyses. Percentage errors in mass balance can be calculated for laboratory datasets since the true total mass in the experimental system is generally a known quantity. However, diffusion/dispersion coefficient value for a system is an unknown quantify and hence true percentage errors cannot be calculated. Therefore, in the published literature, the mass balance error is commonly used as the quantitative indicator of the system error, whereas diffusion and dispersion coefficient value is only used as a qualitative indicator. 55 3.4.2 Comparison of Estimated Concentration Contours We assumed the percentage error in mass balance as a quantitative estimate of error in the computed concentration values. We then compared two concentration contours [5 mg/L (low) and 20 mg/L (high)] obtained from C estimate,EE and C estimate, CRE + EE against the C true values. The mass balance error was used to compute the error bounds using the relationship: () truetruebounderror CErrorBalanceMassCC ??= %1 (5) Figure 3-3a presents the 5 mg/L and 20 mg/L concentration contours obtained from the dataset influenced by EE. It can be seen that the C estimate,EE values have speckle-type noise that are distributed around C true values. The figure also has error bound contours estimated using Equation (5). Since the overall mass balance error is very low (<1%) the C estimate,EE and C error bound cannot be clearly differentiated in the figure. Therefore, we have presented an enlarged section of Figure 3-3 a in Figure 3-3 b, which shows the blue contours (C error bound ) being very close to the red contours (C true values) but not bounding the green contours (C estimate,EE ). It is clear is from Figure 3-3 a and Figure 3-3 b that C error bound contours do not bound the true concentration values. This implies that the mass balance error is not a good indicator of error in the estimated concentration contour values. We employed a trial-and-error approach to generate two contours that will bound the observed spread around the predicted 5 mg/L and 20 mg/L contours. These bounding contours are expected to give a better estimate of error around each of these predicted contours. We estimated that a value of ?10% error will bound the 20 mg/L contour (corresponding to contours of 22 mg/L and 18 mg/L) and ?15% error will bound the 5 56 mg/L contour (corresponding to contours of 4.25 mg/L and 5.75 mg/L). These ?trial and error? bounds are shown in Figure 3-3c along with C estimate and C true . Note that these trial- and-error method based contours fully bound the true contours and therefore it can be concluded that the error associated with 5 mg/L and 20 mg/L contours are approximately ?15% and ?10%, respectively. These values are considerable higher than the mass balance error, which is less than 1%. These results indicate that the mass-balance method fails to correctly predict the local errors associated with a contour level. 57 Figure 3-3 5 mg/L and 20 mg/L contours for C true and C estimate, EE . Part a shows the C estimate, EE and C true with corresponding mass balance error bound. Part b shows an enlarged view of a section of part a. Part c presents C estimate, EE with ?trial and error? bounds of 10% and 15% for 20 mg/L and 5 mg/L contours, respectively. a b c 58 The 5 mg/L and 20 mg/L concentration contours that are influenced by both CRE and EE (C estimate, CRE + EE ) are compared against C true in Figure 3-4. Figure 3-4 a through c, show the contours generated using linear, quadratic, and cubic relationships, respectively. These figures also show the corresponding C error bound contours, which were computed using the mass balance method. Figure 3-4a, which uses the linear relationship, indicates maximum differences. It can be seen in Figure 3-4a that C estimate, CRE + EE contour plumes (for both 5 mg/L and 20 mg/L) are consistently larger than C true values. Also, the mass balance based error bounds cannot correctly envelop C estimate, CRE + EE . Both error sources cause discrepancies between C estimate, CRE + EE and C true values. CRE may cause a significant over or under estimation of concentration, as can be observed in Figure 3-2; note, the linear relationship is shown to overestimate at both 5 mg/L and 20 mg/L. This causes a bias in the estimation process and this effect can be clearly seen in Figure 3-4a, where C estimate, CRE + EE contours are skewed towards one side of C true . Since C error bound is drawn around C true values, this bias also causes the C estimate, CRE + EE contours to be skewed towards either the lower or the higher values of C error bound . For quadratic (Figure 3-4b) and cubic (Figure 3-4c) relationships, the value of CRE is not high as can be observed from Figure 3-2. The speckle-like noise (EE) causes a spread in the contour values similar to those observed in Figure 3-3. From these observations, it can be concluded that the discrepancy between C estimate, CRE + EE and C true values is primarily caused by EE. The above analysis demonstrates that the errors associated with predicated concentration contours are not fully reflected in the mass balance analysis. For example, as shown in Figure 3-3, EE cannot be captured by the mass balance error. Therefore, it can be concluded that mass balance is a necessary but not a sufficient condition. This 59 implies that large mass balance inconsistencies simply indicate presence of errors in the IA technique, however good mass balance does not guarantee its correctness. Furthermore, mass balance and diffusion coefficient errors cannot predict the value of the local errors associated with concentration contours. Therefore, there is a need to develop a robust method that can quantify these local errors. Figure 3-4 5 mg/L and 20 mg/L contours, C true and C estimate, CRE+EE , for linear, quadratic and cubic relationships. Only the mass balance error bounds are shown. b c a 60 3.5 Statistics Based Error Estimation Method Statistical methods can be used to estimate the margin of error for a given set of measurements if large number of measurements are available [Taylor, 1997]. Typical IA datasets involves large number of measurements hence one can use standard statistical method to quantify error. In the sections below we use simple statistical methods to estimate CRE and EE errors and combine them to evaluate the total error. 3.5.1 Using Statistics to Quantify Calibration-Relationship Error As described before, the appropriate calibration-relationship for an observed intensity-concentration dataset is obtained through regression analysis. Typically, R 2 values are used to assess the quality of the regression and the best possible fit is then used for further processing. However, during this regression step one could also calculate the root mean square residual error for the fit, which can provide an indirect measure for CRE. If the selected relationship fits the data well then the residual error () ,,i measured i estimated CC? obtained from regression analysis will be randomly distributed around zero. Under such conditions, the overall error associated with the calibration step can be estimated using the relation [Taylor, 1997]: () df CC N i estimatedimeasuredi iprelationshc 2 1 ,, , ? = ? =? (6) Where ? c,relationship is an estimate for CRE computed directly from regression statistics; N is the number of measurements in the calibration dataset; p is the number of parameters used in the relationship; df is the degrees of freedom of the relationship given by N-p. It 61 should be noted that ? c,relationship represents the average residual error along the entire relationship and it is also known as root mean square error (RMSE). 3.5.2 Using Statistics to Quantify Experimental Errors To quantify EE, one has to devise an approach to quantify the noise in the estimated property (e.g., concentration). This noise is directly related to the standard deviation of the measured image intensity [Hansen et al., 2006]. If the standard deviation in the measured intensity data is known, then the standard deviation in the concentration can be estimated using the relationship [Taylor, 1997]: Ierimentc dI dC ?? = exp, (7) where, ? c,experiment is the error in concentration due to noise in pixel intensity values (image property); dC dI is the slope of the concentration-intensity relationship; and ? i is the randomly distributed error in the recorded pixel intensity field. The physical implication of this relationship is illustrated in Figure 3-5. It can be seen from the figure, that if the observed value of ?i, which represents the standard deviation in pixel intensity corresponding to a known concentration, then the derivative of the calibration relationship at a given concentration value can be used to compute ? c,experiment (using Equation 7). It is important to note that the slope of the relationship changes with concentration and also the slope can be negative or positive depending on the type of IA technique used. Therefore, the absolute (modulus) value of the product of slope and ? i should be used. 62 Figure 3-5 Illustration of concentration-intensity function and its relationship to standard deviation in intensity and the corresponding concentration value 3.5.3 Estimating Total Error and Fractional Error To compute total error due to both CRE and EE, the individual errors can be added in quadrature [Taylor, 1997]: erimentchprelationshc totalc exp, 2 , 2 , ??? += (8) The total error is also a function of concentration values. Also, a low total error value at a low concentration or a high total error at a high concentration may not provide insight into the relative magnitude of the error. Therefore, a quantity known as the fractional error, C totalc, ? , is a better value to quantify errors [Taylor, 1997]. Catania et al. [2008] also reported a fractional error analysis in their study. The fractional error normalizes Intensity ? C,Experiment Standard Deviation: ? I Slope dC dI Concentration 63 the total error with respect to the selected concentration level. Changes in the value of fractional error will directly indicate the magnitude of error in the IA technique. Here, we use the concept of fractional error can be used for two purposes: 1) to estimate errors in a given IA experiment, and 2) to design a better IA experiment. In the sections below we use the theoretical dataset and later an experimental dataset to illustrate the first use. We also design and conduct an IA experiment to demonstrate the second use. Table 3-3 Errors computed by the proposed statistical-method for the theoretical test problem Calibration C 0.5 5 10 20 30 Estimated C True 0.5 5 10 20 30 ? relationship, true 0.00 0.00 0.00 0.00 0.00 ? experiment, true 0.51 0.61 0.74 1.11 1.66 ? total, true 0.51 0.61 0.74 1.11 1.66 Estimated C linear 0.78 7.19 13.07 21.84 27.71 ? relationship, linear 2.05 2.05 2.05 2.05 2.05 ? experiment, linear 0.79 0.79 0.79 0.79 0.79 ? total, linear 2.20 2.20 2.20 2.20 2.20 Estimated C cubic 0.54 5.01 9.86 20.07 30.02 ? relationship, cubic 0.11 0.11 0.11 0.11 0.11 ? experiment, cubic 0.54 0.62 0.79 1.25 1.69 ? total, cubic 0.55 0.63 0.81 1.26 1.69 3.5.4 Verification of Statistics-based Error Estimation Method Our previous analysis indicated that both the mass balance method and the dispersion coefficient method fail to give proper error bounds on the estimated concentration contour levels. The proposed statistics-based method can be considered to be a superior method if the errors estimated by the method can bound the C true values. 64 The statistics-based method discussed above was used to quantify the error in our theoretical problem. Table 3-3 presents the values of CRE (? c,relationship ), EE (? c,experiment ) and total error (? c,total ) for the true, linear and the cubic relationships. In case of the true relationship CRE is not present therefore the true error is the same as EE. It can be observed from Table 3-3 that, as expected, the CRE value was low for the cubic relationship as compared to the linear relationship. Also CRE yields a constant value (RMSE) for all fitted relationships, whereas EE is a constant only when the linear relationship was used. This is because when the concentration-intensity relation is non-linear, EE is expected to have a different value at different concentrations. The most important observation that can be made from Table 3-3 is that ? c,total when added to C true appears to predict the C estimate values quite closely. To verify whether ? c,total calculated in Table 3-3 can successfully predict C estimate values, we calculated and plotted the error bounds (C error bound ) along with C estimate and C true . The C error bound can be calculated using the relationship: totalctruebounderror CC , ??= (9) Figure 3-6 presents C error bound , C estimate and C true . We selected four simulation scenarios to compare the contours. The first scenario (Figure 3-6 a) uses the linear relationship and it represents maximum deviation due to CRE. The second scenario, shown in Figure 3-6 b, uses the noise-only dataset and it represents EE. The third and fourth scenarios show the combined effect of the CRE and EE for the linear and cubic relationships and are presented in Figure 3-6 c and d, respectively. It can be observed from the figures that 65 C error bound enclose C true with good accuracy at both the higher and lower concentration values. This further confirms the results shown in Table 3-3. a b c 66 Figure 3-6 Concentration error bounds based on statistical analysis. a presents the error bounds for linear relationship with the only CRE. b presents the error bounds for only EE. c and d show error bounds for linear and cubic relationships with effects of both CRE and EE. 3.6 Experimental Analysis The purpose of this section is to illustrate the application of the error analysis techniques proposed in this paper using a laboratory dataset. Furthermore, the experiments themselves were designed using one of the error estimation approaches, with a goal to illustrate the power of the statistics-based approach to design better IA experiments. 3.6.1 Laboratory Set-up Figure 3-7 shows a conceptual diagram of our experimental set up. The flow tank was constructed using 6 mm thick PlexiglasTM. The dimensions of the flow regime were: 225 mm (length) ? 200 mm (height) ? 12 mm (width). Constant head chambers were made on either side of the flow region. The chambers were separated from the porous media region by placing US #16 mesh screens on either side. Overflow outlets d 67 were provided in the constant head chambers to control the head gradient across the tank. Clear silica glass beads, with mean diameter of 1.1 mm, were used as the porous media. The tank was packed under saturated conditions in layers of 5 cm with frequent tamping to avoid air entrapment and layering effects. Porosity of the homogeneously packed porous media was measured by both volumetric and gravimetric methods and found to be 0.385. The tank was placed on a leveled platform. Flow Tank CCD Camera Translucent Sheet Leveled Tripod Lighting Figure 3-7 Layout of the experimental set-up As shown in Figure 3-7 we employed backlighting and used light transmittance method (LTM). We have chosen LTM since it provides a better estimate of averaged concentration profile across the width of the tank [Oostrom et al., 2007]. Generic compact fluorescent bulbs (CFL) bulbs were used to light the flow tank. Non-uniform lighting has been cited as a major cause of problems in IA and researchers have suggested use of diffused lighting to reduce such problems [Theodoropoulou et al., 2003; Zinn et al., 2004]. We diffused the light source by placing translucent sheets between the light source and the tank. The experiments were conducted in a darkened environment during late nights to avoid interference from external daylight sources. 68 A 6 MP CCD (charge coupled device) canon digital camera was used to capture the images. General specifications of the camera can be obtained from the online manual [Canon, 2008]. The camera was mounted on a leveled tripod and placed at a distance of approximately 500 mm from the front of the flow tank. Images were captured manually at regular intervals. Care was taken not to disturb the camera. The spatial resolution of obtained images was 0.0267 mm2 per pixel. The digital resolution of the camera was 8 bits per channel (red, blue and green- RGB) providing 256 (0-255) intensity levels at each pixel. A sample image of the porous media tank was acquired and processed before the experiment. The spectral distribution of RGB components indicated that the intensity of red and green signals were more than blue. Dyeing water with a color component will lead to an increase in the intensity of that color and reduction in the intensity of the other components. We selected Red FD&C #4 (liquid) as our dye. Food dyes (FD&C) have been used successfully in other IA applications [Darnault et al., 1998; Detwiler et al., 1999]. 3.6.2 Acquiring Calibration Dataset and Image Processing Different dye solutions with concentrations ranging from 5 ml of dye per liter (5 mL/L) to 0.1 mL/L were prepared. Each of these solutions was flushed through the tank to record the constant image intensity field corresponding to a known concentration. Approximately 10 pore volumes were flushed for each concentration to ensure that the entire tank contained the intended solution. This method for obtaining calibration dataset was used by Schincariol et al. [1993] and later adopted by other researchers including 69 Dong and Selvadurai [2006], Huang et al. [2002], McNeil et al. [2006] and Swartz and Schwartz [1998]. We observed large variations in the intensity field at the edges of our calibration datasets. This was due to edge effects. Previous researchers have made such observation and have removed the edges in their analysis and recommended to conduct experiments away from the edges [Huang et al., 2002]. Parker et al. [2006] cropped their image edges during image processing step. We adopted this approach and cropped the edges of our images. To reduce noise in our image intensities, we employed median filtering using the MATLAB? function medfilt2. Median filtering has been used in the past to reduce errors [Schincariol et al., 1993]. Median filtering is a well accepted image processing technique for removing speckle noises [Jain, 1988]. Also, errors introduced by pixel shifting (due to camera movement) can be removed by median filtering [Vanderborght et al., 2002]. However, it important to note that median-filtering removes noise at the expense of image resolution. We employed median filtering with pixel widths of 3, 5, 7 and 9 pixels. The spatial resolution of the original digital images obtained by our camera was 0.0267 mm 2 per pixel. Therefore, with the largest pixel width of 9, the resulting spatial resolution would not be worse than 1.4 mm 2 , which is of the order of our porous media grain size. Employing median filtering decreased the standard deviations in the intensity data; however, it also decreased the mean pixel intensity level. This is an expected outcome because of the median filtering algorithm. We selected the pixel width of 9 to employ median filtering on our dataset. The standard deviation in image 70 intensities was decreased to a constant value of about 5 pixel intensities (about 2% of our maximum intensity value of 255) for all the datasets. Obtaining an appropriate calibration-relationship is an important next step in IA. We chose a power-law function to fit our observed concentration-intensity data points. Using regression analysis we determined the best-fit parameter values and the resulting relation is: 6988.022190 98.1 ??= ? IC (10) Also, 98.2 2.43936 ? ??= I dI dC (1) Figure 3-8 presents a comparison between the calibration data and the fitted relationship. Dye concentration values above 4 mL/L provided an indiscernible change in intensity and hence these values were not used. Figure 3-8 Calibration data and the fitted intensity-concentration relationship for the experimental system 71 Table 3-4 Errors computed by the proposed statistical-method for the experimental problem Measured C 0.10 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 Estimated C 0.04 0.59 1.06 1.56 1.97 2.40 2.94 3.47 4.08 ? relationship 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 ? experiment 0.04 0.09 0.15 0.22 0.28 0.35 0.44 0.54 0.66 ? total 0.09 0.13 0.17 0.23 0.29 0.36 0.45 0.55 0.67 3.6.3 Designing Experiments Based on Error Analysis Various possible errors were estimated using the statistical approach for the entire range of the calibration-relationship (0 mL/L to 4 mL/L). Table 3-4 provides the values of errors due to calibration relationship, EE and the total error at different values of concentration. It can be observed that the system is dominated by EE. The data indicates that the absolute value of total error increases with increase in concentration. We calculated fractional error ,C total C ? to normalize these variations. Figure 3-9 shows the changes in fractional error with respect to concentration. It is interesting to note that fractional error is high at very low concentrations and then it decreases rapidly and reaches fairly constant value. Figure 3-9 also shows the individual effects of fractional CRE and fractional EE on total fractional error. It can be observed that at lower concentrations, the high values of total fractional errors are the result of high values of fractional CRE and fractional EE. However, at higher concentrations, fractional EE controls the total fractional error. This is an important observation that will help select appropriate design concentration range for conducting the physical experiment. 72 Figure 3-9 Plots of different components of fractional error with respect to concentration 3.7 Conducting Tracer Experiment If we assume the concentration range in the physical experiment to be between 0 mL/L and 4 mL/L, then from the fractional error curves we can estimate the error associated with specific contours a priori. For example, if we want to contour a normalized concentration value of 10% (0.4 mL/L) and 15% (0.6 mL/L), the corresponding fractional error will be about 0.20 and 0.27, respectively (from Figure 3-9). This implies that we will be making an error of 20% and 27% in estimating the 10% and 15% contours, respectively. The error can be reduced considerably if we select a concentration range that has lower total fractional error. It can be observed from Figure 73 3-9 that concentrations between 0.5 mL/L and 4.0 mL/L have a near-constant fractional error of 0.15. We selected this concentration range for the present analysis. We conducted all our transport experiments in an ambient flow field that had a background concentration of 0.5 mL/L. The tracer transport experiment involved simulation of a dense tracer plume sinking in a uniform flow field. We injected 10 ml of saltwater tracer containing a dye concentration of 4.0 mL/L of dye (and having a density of 1.025 mg/L) into the tank. The digital images of the migrating saltwater plume recorded at regular intervals. All the experimental conditions remained identical to those used for obtaining the calibration dataset. The acquired digital images of the physical experiments were cropped and filtered, similar to the processes used for the calibration images. The calibration relationship (Equation 10) was used to compute the dye concentrations in the tank. The digital images of the tracer plume captured at different times are presented in Figure 3-10 a (a1-a3). The plume concentrations (color filled contours) estimated by IA (scale shown in the figure) are given in Figure 3-10 b (b1-b3). It is evident from these two figures that the computed concentration contours qualitatively match the digital image of the tracer plume. Two line contours corresponding to the concentrations1 mg/L and 3.5 mg/L are also presented in Figure 3-10 c (c1-c3). We can use the fractional error data shown in Figure 3-9 to directly quantify the error associated with 1 mg/L and 3.5 mg/L contours levels, which are approximately 15%. Use of the fractional error plot not only helped us to calculate the error level, but also helped us to design a working concentration range for the experimental system. 74 Figure 3-10 a) digital images from the experiment, b) filled concentration distribution obtained from the IA technique, c) 1 ml/l (red) and 3.5 ml/l (blue) contours. Numbers- 1,2 and 3 are for different times after injection 3.7.1 Computation of Mass Balance Error for the Experimental Data The computed mass balance errors for the three experimental datasets (plumes recorded at 3 different times) are provided in Table 3-5. The percentage mass balance error is below 5%, except for the last dataset which was affected by boundaries. The high value of mass balance error in the last dataset (~10%) is attributed to the edge effects a1 a2 a3 b1 b2 b3 c1 c2 c3 75 from the exit boundary. It should be noted that based on fractional error calculations the error associated with this IA system is expected to be about 15%. Clearly, the mass balance calculations underestimated this error, which further indicates that mass balance are not appropriate for estimating error in individual concentration measurements. Table 3-5 Errors in mass balance calculation for all experimental datasets Experimental Dataset 1 2 3 Mass Balance Error (%) 1.2 4.8 10.2 3.8 Summary and Conclusions A method for using IA techniques to measure concentration values was described along with a detailed discussion of errors associated with the application of this technique. The errors were categorized into two broad categories that include calibration-relationship error (CRE) and experimental error (EE). The CRE is introduced into an IA method due to the lack of a priori knowledge of the structure of the relationship that describes the concentration-intensity data (also known as the calibration data). EE is introduced due to various random experimental errors associated with lighting intensity and other image capturing/processing techniques. We developed a theoretical test problem where CRE, EE, and both CRE and EE were systematically introduced. The two error quantification methods that are widely used in the literature-- the mass balance error method and the dispersion coefficient method, were used to quantify the effects of the errors in theoretical test problem. Our analyses show that the dispersion approach can only be used for making certain qualitative statements about errors. Furthermore, the mass balance approach fails to correctly quantify the local error associated with specific contour levels. This is primarily due the fact that the mass 76 balance approach cannot capture EE. Therefore, it is concluded that while mass balance closure is a necessary condition it is not a sufficient condition to guarantee the correctness of IA measurements. This implies that large mass balance inconsistencies simply indicate presence of errors in the IA technique, however good mass balance closure does not imply the results are correct. We propose an alternate, statistics-based approach to evaluate errors associated with IA techniques. The approach was first used to successfully quantify the error in our theoretical test problem. A fraction error analysis, which was derived from the proposed approach, allows one to not only quantify errors but to also help design an optimal concentration range for the experiment. A tracer transport experiment involving the movement of a dense tracer plume in a uniform flow field was designed using this approach. The fractional error analysis indicated that the use of 0.5 mg/L of background concentration would considerably reduce error in the experimental system. Three sets of digital images were captured for the dense-tracer transport experiment. The images were processed and the concentration contours for the system were generated using the proposed approach. The results show that the proposed approach can successfully generate concentration contours from a captured digital image. This study is unique and provides a comprehensive assessment of errors associate with IA techniques. A theoretical problem with a known error structure was used for the first time for error analysis. Use of theoretical problem allowed us to unequivocally identify the limitations of existing error quantification methods. The proposed statistics- based error quantification approach is shown to be a better alternative to quantify IA errors. 77 CHAPTER 4. ANALYSIS OF NUMERICAL TECHNIQUES USED FOR SIMULATING INSTABILITIES OBSERVED IN VARIABLE- DENSITY FLOW EXPERIMENTS 4.1 Introduction Several types of variable-density experimental datasets, collected under controlled laboratory conditions, are available for benchmarking numerical codes. These datasets can be classified into two broad categories- 1) systems involving stable density- configuration, and 2) systems involving unstable density-configuration. Stable density- configurations involve a lighter fluid overlying a denser fluid. Examples of stable density-configurations include saltwedge experiments that have considered the transport patterns of saltwater wedge into a saturated aquifer [Goswami and Clement, 2007], and the saltpool experiment that considered the upward movement of saltwater under pressurized conditions [Oswald and Kinzelbach, 2004]. Unstable density-configuration involves denser fluid overlying a lighter fluid. Examples of unstable density- configuration systems include dense sinking plumes [Hogan, 2006; Oostrom et al., 1992b; Simmons et al., 1999; Simmons et al., 2002], rising plumes [Brakefield, 2008], and large-scale saltwater contamination events such as tsunami [Illangasekare et al., 2006]. 78 Numerical modeling of unstable density-configuration problems is considered to be more difficult than stable density-configuration problems. This is because, unstable density-configurations can lead to generation of instabilities or fingering patterns during migration of the plume [List, 1965; Wooding et al., 1997a; Wooding et al., 1997b]. The fingering patterns observed in such experiments are difficult to be replicated numerically and require introduction of artificial perturbations [Schincariol et al., 1994; Simmons et al., 1999]. In general, two methods have been used to simulate fingering patterns in variable-density numerical models without introducing artificial perturbations in the flow or concentration field [Liu and Dane, 1997]. The first method requires the use of particle-tracking algorithm, such as the method of characteristics, to solve the transport equation. The second method is to introduce heterogeneity in the porous media characterization [Schincariol et al., 1997; Schincariol, 1998; Simmons et al., 2001]. Onset of instabilities in density coupled flow systems are often analyzed using dimensionless parameters [List, 1965; Oostrom et al., 1992b]. However, recent experimental and theoretical studies have shown that dimensionless numbers have limited ability to predict the occurrence of instabilities in variable-density flow scenarios [Menand and Woods, 2005; Simmons et al., 2001; Simmons, 2005]. Therefore, groundwater modelers and practitioners are primarily dependant on numerical simulations to estimate if instabilities would develop in a given variable-density flow scenario. However, there are no guidelines available in the literature that aid in evaluating the applicability of different numerical approaches to correctly predict the generation of instabilities. Therefore, the primary objective of this work is to provide an analysis of the numerical approaches available for correctly estimating development of 79 instabilities observed in variable-density experiments. Dispersion has a substantial effect on the generation of instabilities [Menand and Woods, 2005]. Therefore, our secondary objective is to test the ability of the various numerical approaches to correctly evaluate the extent of the - mixing or dispersion zone in unstable density-configuration systems. This study includes both experimental and numerical tasks. For the experimental part, we completed a set of unstable density-configuration laboratory experiments to generate two distinctly different datasets- one with instabilities and another without instabilities. We obtained both datasets under similar porous media conditions and density differences but with different boundary and initial conditions. Image analysis techniques were used to obtain concentration data. In the numerical part of the study, we tested both the particle-tracking and the heterogeneity approach to model the experimental datasets. Method of characteristics (MOC), total variation diminishing (TVD), and finite-difference solvers available within the USGS code SEAWAT [Langevin and Guo, 2006] were used to model the two datasets using identical numerical parameters. We also evaluated the performance of the same solvers to predict instabilities under both homogeneous and heterogeneous hydraulic conductivity fields. The results of the study are used to develop an understanding of the advantages and disadvantages of using these numerical techniques for simulating unstable density- configuration scenarios. 4.2 Laboratory Methods The dimensions of the flow tank used in this study are: 225 mm (length) ? 200 mm (height) ? 12 mm (width). Overflow outlets provided in the constant head chambers controlled the head gradient across the tank. Clear Silica glass beads (mean diameter 1.1 80 mm) were used as the porous media. The laboratory setup and experimental procedures are similar to those described in Goswami et al [2008]. The porous medium was packed under saturated conditions in layers of 5 cm and disturbed frequently during packing to avoid air entrapment and to reduce layering effects. The porosity was measured by both volumetric and gravimetric methods and found to be approximately 0.385. Saltwater used for injection was prepared by dissolving NaCl in water. Density of the prepared saltwater was measured to be 1.025 g/mL (salt concentration of 0.035 gm/mL) by ASTM # 111H hydrometer and by gravimetric method. A 6 MP Canon digital camera providing 8 bit data (0-255 intensity levels) was used to capture images manually. An image analysis procedure was used to extract concentration contours from the digital images. Details of the image analysis procedure are provided in Goswami et al [2008] and Appendix-I. The average saturated hydraulic conductivity of the porous medium was calculated from flux measurements during the calibration experiment and found to be approximately 11.57 mm/sec (1000 m/day). A similar flux measurement method was used by Goswami and Clement [2007] to calculate in-situ hydraulic conductivity value. 4.2.1 Design of Variable-density experiments Two unstable density-configuration experiments, related to different practical problems, were completed. The first experiment was intended to simulate the movement of treated wastewater effluents injected into a saltwater aquifer. Such a scenario is encountered in Florida during the operation of deep-well injection wells for disposing treated wastewater effluents and aquifer storage and recovery systems in saline aquifers. For the physical experiment, the flow tank was filled with saltwater and then freshwater was introduced in the lower part of the flow tank through an injection port. This 81 experiment was conducted in the absence of ambient flow to create optimal conditions for generating instabilities. We will refer to this experiment as the ?rising plume? experiment. The second experiment was designed to simulate the injection of a saltwater slug into a freshwater aquifer. Such a scenario represents the contamination of a well during a saltwater inundation event such as a tsunami or hurricane [Carlson et al., 2008; Illangasekare et al., 2006]. Ambient freshwater flow conditions were established by introducing a hydraulic gradient in the experimental tank from left to right. Saltwater was then injected in the upper part of the flow tank. The flow field was designed to facilitate migration (sinking) of the saltwater slug without any instabilities. This experiment will be referred as the ?sinking plume? experiment in this work. 4.3 Laboratory data For the rising plume experiment, the flow tank was flushed with more than 10 pore volumes of saltwater (solution density of 1.025 g/mL, salt concentration of 0.035 gm/mL) resulting in a completely saline system. Ambient flow conditions were eliminated by setting zero head gradient across the tank. About 5 mL of red-colored freshwater (containing dye concentration of 4.0 mL/L and solution density of 1.000 g/mL) was injected into the tank with a syringe. The injection port was located approximately 114 mm from the left end and 17 mm from the bottom of the tank. The freshwater was injected slowly to avoid disturbing the packing in the vicinity of the injection port. It took almost 10 seconds to inject the 5 mL of freshwater into the tank. Digital images of the migrating freshwater plume were acquired at different times during the experiment. The acquired images were cropped to minimize edge effects that can 82 introduce errors in the datasets [Parker et al., 2006]. Cropped images of the freshwater plume captured at four different times are presented in Figure 4-1 (a-d). 4.3.1 Results of Rising Plume Experiment Figure 4-1a was taken just after injection was completed which was approximately 10 seconds after the start of injection. Note the shape of the injected freshwater bubble is almost circular at the end of injection. Immediately after injection, the density contrast between the injected freshwater and ambient saline water started the transport of the freshwater bubble towards the top of the tank. At the same time, the unstable density-configuration of denser water on top of the fresher water also introduced instabilities in the migrating front of the freshwater plume. These instabilities represented by the fingering effect can be seen emerging at the top of the main body of the plume in the second image (Figure 4-1b) captured 2 minutes after injection. The fingers develop further during the upward migration of the plume in the tank. As the freshwater migrates it also mixes with the ambient saline water causing the fingers to disperse. Such dispersion smoothens the profile of the plume as the fingers coalesce [Schincariol et al., 1994]. Figure 4-1c shows the plume 4 minutes after injection; in this figure, although the fingers appear to be longer than those observed in Figure 4-1 b they are connected to each other due to the dispersive mechanism. The plume approaches a conical shape as it migrates upwards through the tank which can be seen in Figure 4-1d (taken after 6 minutes) where the top of edge of the plume has almost reached the top of the tank. 83 Figure 4-1 Experimental results from the rising plume experiment at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection 4.3.2 Results of sinking plume experiment The flow tank was repacked for conducting the sinking plume experiment. Freshwater with a density of approximately 1.000 g/mL and red dye concentration of 0.5 mL/L was used to pack the tank. In this study, a light red background was used to enhance the quality of the image processing procedure (for image processing details see Chapter 3, also [Goswami et al., 2008]). Ambient flow conditions were created from left to right by forcing a head gradient of 4 mm across the tank. About 5 mL of a saltwater tracer containing a dye concentration of 4.0 mL/L (combined solution density of 1.025 g/mL, salt concentration of 0.035 gm/mL) was injected in the upper left portion of the tank. The injection port was located approximately 54 mm from the left end and 145 mm from the bottom of the tank. Similar to the first experiment, it took almost 10 seconds to inject the 5 mL of saltwater into the tank. Images of the saltwater plume were also acquired and cropped at the edges. The final processed digital images of the saltwater plume captured at three different times are presented in Figure 4-2 (a-c). a b c d 84 Figure 4-2 Experimental results from the sinking plume experiment at- a) 0 min, b) 2 min, and c) 5 min - after injection Figure 4-2a was taken just after injection was completed. The shape of the injected saltwater is initially circular, similar to that obtained in the rising plume experiment. With time, the heavier saltwater moves downwards and the ambient freshwater flow simultaneously transports the saltwater towards the right boundary. As seen in the rising plume experiment, one would expect to see instability in the movement of the migrating saltwater plume since the density configuration of this system is unstable. However, due to the presence of ambient flow, the generation of instabilities is suppressed. The ambient flow enhances mixing and the increased dispersion decreases the density contrast between the injected and the ambient fluid density which is the driving force for creating instabilities. The dispersed plume attains an elliptical shape, as shown in Figure 4-2b and Figure 4-2c. 4.4 Concentration data derived from Image Analysis The image analysis procedure discussed in Appendix-I was used to derive salt concentration values from the physical experiments discussed above. The salt concentration datasets obtained from the rising and the sinking plume experiments are plotted in Figure 4-3 and Figure 4-4 and will be referred to as datasets 1 and 2 in this a b c 85 work. Results from the rising plume experiment (Figure 4-3) show distinct fingers/instabilities confirming the results presented in Figure 4-1. It should be noted that the plume has a narrow mixing zone just after injection (Figure 4-3 a). As the plume migrates, the instabilities are generated at the leading front of the plume. Also the instabilities decrease as the plume migrates upwards through the porous media which can be explained by the decreasing concentration gradients and a more dispersed plume with the progression of time. It can be seen from Figure 4-3 d that not much freshwater (blue) remains in the system since the plume appears to be mostly consisting of salt concentrations in the range of 0.02 to 0.03 gm/mL. Some preferential flow appears to have slightly skewed the transport towards the right in Figure 4-3 b) and towards the left in Figure 4-3 c. However, the overall transport pattern of the rising plume experiment displays a symmetrical pattern. The concentration profiles obtained from the sinking plumes experimental dataset (Figure 4-4) do not show instabilities. Furthermore, it should be noted that the saltwater plume has a sharp mixing zone just after injection (Figure 4-4a). High salt concentrations (~0.03 gm/mL) appear to be present in large sections of the plume at later times (Figure 4-4 b and c). We will refer back to these mixing zone observations while comparing these experimental results against numerical predictions. 86 Figure 4-3 Concentration profiles obtained from the rising plume experiment at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection Figure 4-4 Concentration profiles obtained from the sinking plume experiment at- a) 0 min, b) 2 min, and c) 5 min - after injection 4.5 Numerical Analysis Numerical simulation of variable-density flow scenarios requires the use of coupled flow and transport models. However even with the use of coupled models, typical fingering patterns observed in the physical experiments are difficult to reproduce numerically. As discussed earlier, here we will focus on the two types of approaches to simulate instabilities in variable-density experiments. As the first approach we will a b c d a b c 87 explore the use of numerical perturbations generated by different MT3DMS solver options including MOC, finite-difference, and TVD. As the second approach we will explore the use of perturbation generated by physical heterogeneities in the hydraulic conductivity field. SEAWAT [Langevin and Guo, 2006] was used to model the physical experiments. Our objective is to analyze the performance of numerical codes based on two criteria- a) generating appropriate level of instabilities, and b) reproducing the appropriate level of dispersion/mixing (as obtained from image analysis). 4.5.1 Numerical Model Details As the first step in developing the numerical model we set up conceptual descriptions of the physical experiments. Figure 4-5 a and b present the conceptual diagrams representing both experiments. Figure 4-5 Conceptual description of- a) rising plume, and b) sinking plume experiments The problem domain was considered to be two-dimensional since the width (12 mm) of the tank was much less than the other two dimensions (225 mm and 180 mm). Separate models were developed for both experiments with horizontal and vertical dimensions of a b 88 225 mm and 180 mm respectively. Since our goal is not to quantitatively match the physical experiments but only to analyze the ability of different numerical approaches to correctly predict the generation of instabilities in the experiments, we conducted all numerical simulations with a uniform grid of 1 mm, which is about the order of the mean grain size. The grid discretization was the same in both horizontal and vertical directions. The boundary conditions shown in Figure 4-5 were used in the numerical model. Initial conditions with fully-saline and fully-fresh water conditions were used for the rising plume and the sinking plume experiments, respectively. Two cells at the upper left and right corners of the rising plume model were set to constant zero pressure throughout the simulation. The zero pressure cells allowed a small amount of fluid to leave the system during and after injection and similar conditions were used in the modified version of the Elder problem solved by Voss and Souza [1987]. Porous media packing was considered to be isotropic. Hydraulic conductivity distribution was considered to be homogeneous in the first set of simulations, and was assumed to be heterogeneous with a lognormal distribution in the second set of simulations. Porosity and average hydraulic conductivity values were set at the measured value of 0.385 and 11.57 mm/sec, respectively (refer Section 4.2). Longitudinal dispersivity values was set at 0.1 mm and the value of transverse dispersivity was set to be 1/10 th of the longitudinal dispersivity [Goswami and Clement, 2007; Johannsen et al., 2002]. Three stress periods were used to simulate- a) steady-state conditions before the experiment, b) injection of freshwater/saltwater, and c) migration of injectant through the tank. Time stepping was set to an initial minimum of 0.1 seconds with a maximum value of 1 second; the transient values were also limited by setting a courant number of 0.5. All other parameters for 89 SEAWAT were set at default values set within GW Vistas [Rumbaugh and Rumbaugh, 2007]. 4.5.2 Homogeneous Simulations with MOC MOC, a particle-tracking algorithm, was the first solver used to simulate the experiments. The simulation results for the rising plume experiment are presented in Figure 4-6. As shown in the figure, MOC was able to reproduce the complex instability generation process (or fingering patterns) exhibited in the experiment. It can also be observed (Figure 4-6a) that the simulated mixing zone between freshwater and saltwater concentrations is narrow at the beginning of rising plume experiment, as shown in dataset 1 (Figure 4-3a). Therefore MOC satisfied both of our model performance criteria- generation of instabilities and reproducing the appropriate level of dispersion. MOC simulation results from the sinking plume model are presented in Figure 4-7. Interestingly, the simulation results show fingering patterns whereas no such fingering or instabilities were observed in the physical experimental data shown in Figure 4-4. Thus, MOC failed to satisfy the first criterion. Next, we tested the ability of the finite-difference solver to simulate the experiments. 90 Figure 4-6 Results from MOC simulations of the rising plume experiment with homogeneous packing at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection Figure 4-7 Results from MOC simulations of the sinking plume experiment with homogeneous packing at- a) 0 min, b) 2 min, and c) 5 min - after injection 4.5.3 Homogenous Simulations with finite-difference Both the upstream and central options were used in finite-difference to simulate the experiments with a homogeneous hydraulic conductivity distribution. Results from both these solvers exhibited significant differences; therefore, we have summarized both of these simulation results in Figure 4-8 through Figure 4-11. Unlike the MOC a b c d a b c 91 simulations, both finite-difference methods completely failed to simulate instabilities in the rising plume. Both the upstream and central schemes also fail to reproduce the general mixing zone in the rising plume experiment. The upstream scheme incorrectly simulates a large mixing zone at the beginning of the experiments (Figure 4-8 a and Figure 4-10 a). However, the upstream method closely estimated the volume of freshwater (blue) left in the system at end of the experiment (Figure 4-3d and Figure 4-8d). On the other hand, it failed to estimate the amount of saltwater (red) left in the plume at the end of the sinking plume experiment (Figure 4-4d and Figure 4-10d). Results from the central scheme also show similar anomalous behavior. The mixing zone after injection, however, was correctly simulated -for both the rising and sinking plume problems (Figure 4-9 a and Figure 4-11 a). Also, the amount of saltwater (red) remaining at the end of the sinking plume experiment was estimated to be close to that in the physical experiment (Figure 4-10 d and Figure 4-4 d). However, the central scheme falsely estimated a large amount of freshwater (blue) to be present in the system at the end of the rising plume experiment (Figure 4-9 d and Figure 4-3 d). The observed anomalies are primarily caused by numerical dispersion which is different for the upstream and central weighting schemes [Woods, 2004]. The numerical dispersion also suppressed the instabilities [Schincariol et al., 1994]. Therefore, we simulate the experiments using the total-variation-diminishing (TVD) scheme, a higher order finite- difference scheme that reduces numerical dispersion effects and hence is capable of tracking sharp-fronts [Zheng and Bennett, 2002]. 92 Figure 4-8 Results from upstream finite-difference simulations of the rising plume experiment with homogeneous packing at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection Figure 4-9 Results from central-in-space finite-difference simulations of the rising plume experiment with homogeneous packing at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection a b c d a b c d 93 Figure 4-10 Results from upstream finite-difference simulations of the sinking plume experiment with homogeneous packing at- a) 0 min, b) 2 min, and c) 5 min - after injection Figure 4-11 Results from central-in-space finite-difference simulations of the sinking plume experiment with homogeneous packing at- a) 0 min, b) 2 min, and c) 5 min - after injection a b c a b c 94 4.5.4 Simulations with TVD TVD simulations have been shown to reproduce sharp transition (or mixing) zones in stable density-configuration experiments [Goswami and Clement, 2007]. Simulation results obtained using the TVD solver for our unstable density configurations are presented in Figure 4-12 and Figure 4-13. For initial time steps, the rising plume results (Figure 4-12 a and b) are similar to those obtained from the central finite- difference scheme (Figure 4-9 a and b); the model simulated a small mixing zone and a large finger. Thereafter, for the next two time steps, the TVD results are able to simulate complex instability patterns (Figure 4-12 c and d). However, the instabilities simulated by the TVD scheme do not correspond to those observed in the physical experiments (Figure 4-3). The major differences lie in the migration patterns of the central finger, which is faster and larger in the TVD simulation when compared to the physical experimental data. However, if we ignore these differences, TVD simulations were able to reproduce the instabilities observed in the rising plume experiment. The TVD simulation results for the sinking plume problem closely estimated the small mixing zone observed at the beginning (Figure 4-4 a and Figure 4-13 a). In addition, the amount of saltwater (red) at the end of the simulation (Figure 4-4 c and Figure 4-13 c) is similar to those observed in the experimental data. However, the sinking plume simulations show instabilities which were not observed in the physical experiments. In summary, with a homogeneous hydraulic conductivity distribution, both MOC and TVD solvers were able to generate instabilities in the rising plume experiment. However, they also generated spurious instabilities in the sinking plume experiment. While the finite-difference methods (upstream and central) correctly predicted the sinking 95 plume experiment without any instabilities they failed to generate the instabilities observed in the rising plume experiment. Therefore, all these methods- MOC, TVD and finite-difference- failed to satisfy the first criterion of our analysis. To overcome this limitation, we employed a heterogeneous hydraulic conductivity field resolved the two problems using all three solvers. Figure 4-12 Numerical modeling results obtained with TVD for rising plume experiment with homogeneous packing at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection Figure 4-13 Numerical simulation results obtained with TVD for sinking plume experiment with homogeneous packing at- a) 0 min, b) 2 min, and c) 5 min-after injection b c d a b ca 96 4.5.5 Simulations with Heterogeneity As described earlier, in the literature, heterogeneous hydraulic conductivity fields have been used to simulate fingering patterns observed in variable-density plumes. Although our tank was packed carefully to obtain a homogeneous flow field, some level of small scale heterogeneities (of the order of 1% variation) are difficult to eliminate. We used the turning bands method (by employing the computer code TUBA) [Zimmerman and Wilson, 2005] to generate the desired lognormally distributed heterogeneous conductivity distribution. Different hydraulic conductivity realizations were generated using a correlation length of 10 mm and a small scale variability of 1%. All other parameters within TUBA were set at default values. Except for conductivity values, all other numerical model parameters used in the heterogeneous simulations were identical to those used in the homogeneous simulations. 4.5.6 Finite-difference with heterogeneity As both MOC and TVD solvers had failed the first criterion in the homogeneous field itself, we started our analysis with the heterogeneous field by using the finite- difference method. Results of numerical simulations obtained using both upstream and central finite-difference methods are presented in Figure 4-14 through Figure 4-17. Modeling results from both methods show complex fingering patterns for the rising plume experiment, as can be seen in Figure 4-14 and Figure 4-15. The complexity in the generation of fingers appears to be more in the case of central method when compared to the upstream method. Simulation results from the sinking plume experiment are presented in Figure 4-16 and Figure 4-17. Both upstream and central methods provide a good qualitative match with the sinking plume data. Unlike the MOC and TVD results 97 for the homogeneous case, no spurious instabilities were generated in the numerical simulation of the sinking plume problem. However, similar to the homogeneous simulations, the upstream method produced more dispersion and hence a larger mixing zone than those observed in the experimental datasets. Although the numerical dispersion appears to be less for the central method, it cannot be totally avoided. The numerical dispersion problems can be reduced by using MOC or TVD schemes. Therefore, we applied both these methods with a heterogeneous field with the goal to reduce numerical dispersion effects. However, as expected, both MOC and TVD schemes introduced spurious oscillations in the sinking plume experiment. The general trend of for the sinking plume simulation results was similar to those obtained using the homogenous hydraulic conductivity field. We are not presenting the simulation results for MOC and TVD with heterogeneity since the results were similar to those obtained before (refer Figure 4-6, Figure 4-7 (MOC); Figure 4-12, Figure 4-13 (TVD)). In summary, the application of MOC and TVD schemes resulted in the failure to satisfy both criteria. The difference in the case of heterogeneity simulations was an increased level of complexity in the generation of instabilities as compared to the homogenous case. This is expected since heterogeneities would introduce more perturbations to the flow field thereby increasing the chances of introducing more instabilities. Also, as expected, the mixing zones were smaller in both MOC and TVD simulations. 98 Figure 4-14 Numerical modeling results for rising plume experiment with heterogeneous packing and upstream finite-difference method at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection Figure 4-15 Numerical modeling results for rising plume experiment with heterogeneous packing and central finite-difference method at- a) 0 min, b) 2 min, c) 4 min, and d) 6 minutes- after injection b c d a b c d a 99 Figure 4-16 Numerical modeling results for sinking plume experiment with heterogeneous packing and upstream finite-difference at- a) 0 min, b) 2 min, and c) 5 min- after injection Figure 4-17 Numerical modeling results for sinking plume experiment with heterogeneous packing and central finite-difference at- a) 0 min, b) 2 min, and c) 5 min- after injection 4.6 Summary Two unstable density-configuration physical experiments were conducted. In the first experiment the freshwater was injected into a saline aquifer resulting in the upward migration of the freshwater bubble. This experiment was identified as the rising b ca b ca 100 plume experiment. In the second experiment, a saltwater bubble was injected into an ambient freshwater flow field. This experiment was identified as the sinking plume experiment. Experimental conditions set to allow generation of instabilities in the rising plume experiment and none in the sinking plume experiment. Image analysis (IA) was used to extract concentration distributions from both experimental datasets. Two types of numerical modeling approaches were used to simulate the experiments. In the first approach, MOC (method of characteristics), a particle-tracking method was used to simulate both experiments. A lognormal heterogeneity distribution for the hydraulic conductivity field was used in the second approach in conjunction with the same solvers. MOC results indicated instabilities in either both experiments or in neither experiment based on the parameters used. These results suggest that MOC is not appropriate for simulating variable-density flow problems since it could not correctly predict the generation of instabilities in both physical experiments. Finite-difference and TVD (total-variation-diminishing) methods were also used to simulate the experiments with homogenous conditions. While TVD results were similar to those from MOC in that they produced instabilities in both experiments, the finite-difference methods failed to generate instabilities. However, instabilities could be predicted correctly in the physical experiments using the finite-difference methods by introducing heterogeneity in the hydraulic conductivity field. MOC and TVD solvers failed to correctly predict instabilities again with a heterogeneous field. It should be noted that the uniqueness of this work lies in the combined use of the two physical experiments for conducting the numerical analysis. Taken individually, we can successfully model the experiments using different numerical approaches. 101 Instabilities in the rising plume experiment are well predicted by using any of the MOC or the TVD solvers or the finite-difference solver with a heterogeneous hydraulic conductivity field. Similarly, the sinking plume experiment can be modeled successfully with the finite-difference method. However the options seem to decrease when we combine the two experiments together and try to solve them using the same numerical approach. From our results, it appears that the finite-difference method used in conjunction with a heterogeneous hydraulic conductivity distribution is the optimal way to model an unstable density-configuration variable-density scenario. However, overall the finite-difference methods failed to reproduce the narrow mixing zone due to numerical dispersion. This analysis presents a need to improve the current state-of- science in numerical modeling of unstable density-configuration flow problems in particular and variable-density flow problems in general. 102 CHAPTER 5. TWO-DIMENSIONAL TRANSPORT CHARACTERISTICS OF SURFACE-STABILIZED ZERO-VALENT IRON NANOPARTICLES IN POROUS MEDIA 5.1 Introduction The motivation for this work was derived from the possibility of field application of a remediation technology developed by Dr. Sushil Raj Kanel, formerly a postdoctoral researcher at Auburn University. The basic premise is to develop a fundamental understanding of iron nanoparticle delivery methods by studying its transport behavior in a two-dimensional aquifer model, which was described in previous sections. Synthesis, preparation and analysis of nanoparticles and the related chemical analysis were performed by Dr. Kanel. As a joint lead investigator, the author of this dissertation conducted the physical nano-particle transport experiments and the numerical modeling studies. Iron nanoparticles are used in a variety of areas for magnetic/electronic, catalytic, and biomedical applications [Huber, 2005]. In the environmental area, nanoscale iron materials have been widely researched to explore their potential for treating contaminated soil and groundwater [Li et al., 2006]. Among available iron nanoparticles, zero-valent iron nanoparticles (INP) have attracted significant interest due to their ability to reduce a variety of environmental contaminants. For example, INP have been found to degrade 103 chlorinated hydrocarbons such as trichloroethene (TCE), tetrachloroethene (PCE), and carbon tetrachloride [Nurmi et al., 2005; Zhang, 2003]. In addition, environmental contaminants such as perchlorate [Cao et al., 2005], nitrate [Yang and Lee, 2005], and metals such as Cr(VI) [Manning et al., 2007; Ponder et al., 2000], lead, nickel, mercury [Li et al., 2006], arsenic [Kanel et al., 2005; Kanel et al., 2006],and U(VI) [Burghardt and Kassahun, 2005] can be transformed using INP. The INP can also produce hydroxyl radicals in the presence of oxygen to oxidize a variety of organic contaminants such as carbothioate herbicide/molinate [Joo et al., 2004] and benzoic acid [Joo et al., 2005]. Despite its high reactivity, the natural tendency of INP to aggregate, due to its magnetic properties [Phenrat et al., 2007], may severely limit our ability to deliver INP into deep porous media formations [Nurmi et al., 2005]. To overcome this limitation, various surface modification and particle stabilization strategies have been developed by using different types of additives such as surfactant (Tween-20) [Kanel et al., 2007], poly acrylic acid (PAA) [Schrick et al., 2004], carboxymethyl cellulose (CMC) [He et al., 2007], cellulose acetate [Wu et al., 2005] starch [He and Zhao, 2005], noble metals [Elliott and Zhang, 2001], and oil emulsions [Quinn et al., 2005]. A majority of these studies used batch experiments to demonstrate the additive?s potential to stabilize the INP. However, the transport dynamics of stabilized INP can only be tested under dynamic flow conditions. Only a few studies have explored the transport behavior of S- INP in soil columns. Schrick et al. [2004] studied PAA-stabilized INP and its reactivity in a glass burette ; Kanel et. al [2007] studied surfactant (Tween 20) stabilized INP and PAA-stabilized INP in a sand column and in a glass-bead packed column. They also studied the reactivity of various forms of stabilized INP for removing arsenic species. All 104 of the above INP transport studies were limited to one?dimensional analysis. To the best of our knowledge, transport of S-INP under two-dimensional flow conditions has not been reported in the literature. Furthermore, there have been no studies on numerical modeling of the observed transport characteristics of INP in groundwater systems. In this study, we hypothesize that two-dimensional physical models can be used to unravel the multidimensional transport dynamics of S-INP, which may be influenced by small density gradients. We use a novel experimental setup to demonstrate the importance of density effects while injecting nanoparticles into saturated aquifer formations. We compare the two-dimensional transport data of S-INP and INP plumes against a tracer plume to demonstrate the efficiency of the stabilization process. Finally, we use the numerical model SEAWAT to test whether the observed S-INP plume can be conceptually modeled as a density-driven conservative plume. 5.2 Materials and Methods All the chemicals used in the experiments were reagent grade. Chemicals such as NaBH4 and PAA were obtained from Sigma-Aldrich Chemical Co. (Sigma-Aldrich, St. Louis, MO). Ferrous iron (FeSO4.7H2O), was obtained from Fisher Chemical Company (Fisher Scientific, Fairlawn, NJ). The porous media selected for this study was A-110 silica beads obtained from Potters Industries (Malverne, PA). The mean bead diameter was 1.1 mm with a variation of 0.1 mm. The porous medium properties were estimated using methods reported in our previous work [Goswami and Clement, 2007]. The average porosity of the packed system was estimated to be 0.385. The average hydraulic conductivity was estimated to be 1050 m/day from in situ flow and head measurements. The value of longitudinal dispersivity was estimated to be 1 mm from tracer experiments. 105 Nonreactive dye (FD&C Red # 40) was used as an optical tracer in all the experiments. Transport characteristics of this dye have been verified in previous experiments where it has been used to visualize the movements of nonreactive solutes [Goswami and Clement, 2007]. Ultrapure (18?cm) deionized water purified by a Barnstead purification system was used to prepare all nanoparticles suspensions. A two-dimensional flow container, shown in Figure 5-1, was used as the physical model to conduct experiments. The dimensions of the flow container are: 50 cm (length) ? 2 cm (width) ? 28.5 cm (height). Two chambers (5 cm wide) were built at the two ends (left and right) of the container to set constant-head boundary conditions. A series of overflow orifices were drilled in the left and right chambers of the flow container. These overflow orifices allowed excess fluid to drain from the system and also controlled the water level (head) in the chamber. In all the experiments the ambient freshwater flowed from right to left by establishing a head difference of 0.7 cm (gradient of 1.4%) between the head chambers. The flow was allowed to reach steady state for a period of 10 min before starting the injection experiments. Figure 5-1 Conceptual diagram of the flow container 106 5.3 Tracer Test After establishing steady-state flow conditions, 20 mL of freshwater colored with an optical tracer (red dye) was injected into the porous media to characterize the movement of freshwater in the physical system. The location of the injection point was approximately 15 cm from the right end and 16cmfrom the bottom of the inner dimensions of the flow container. It took approximately 16 s to manually inject 20 mL of tracer (at the rate 1.25 mL/s). The transport of the tracer was recorded for about 15 min by taking high resolution digital pictures at regular intervals. Similar data collection techniques were previously employed to study the migration patterns of dense plumes in porous media systems (24). 5.4 Preparation of INP and S-INP INP and S-INP were synthesized using previously reported methods with minor modifications [Kanel et al., 2005; Kanel et al., 2007]. In this study, 3.25 g of FeSO4?7H2O and 3.05 g of polyacrylic acid (PAA), (C3H4O2) mol wt.: 1800 g/mol, were dissolved in 100 mL of deionized water. The metal salts were reduced by adding 2 g of NaBH4 dissolved in 50 mL deionized water under nitrogen environment. The total concentration of nano-Fe obtained was 4 g/L, which was measured using an atomic absorption spectrophotometer (Spectra AA 220 FS). The density of S-INP suspension was 1.036 g/cm 3 and the pristine INP suspension was 1.030 g/cm3. The surface area of INP and S-INP preparations synthesized using this approach is expected to be in the range of 15-30 m 2 /g [Kanel et al., 2005]. INP and S-INP were freshly prepared prior to each transport experiment. Similar to the tracer test, 20 mL of 4 g/L INP and S-INP were injected separately into the system using the procedure described earlier. Digital images 107 were taken every minute for a total time of 14 min to record the movement of the nanoparticle plumes. 5.5 Modeling Tracer and S-INP Transport Processes. Computer modeling of the data sets obtained from the experiments described above requires the capability to solve both the groundwater flow and transport equations. Our characterization studies indicated that the INP and S-INP solutions are denser than water and hence flow and transport had to be solved simultaneously in a coupled mode. In the groundwater literature, it is a common practice to solve the flow and transport problems independently by assuming that changes in the concentration field have no effect on the flow field [Zheng and Bennett, 2002]. This assumption may not hold for our experiments since the small-scale density gradients introduced by the dense INP plume may influence the groundwater flow patterns. Accurate numerical simulation of density- dependent flow systems requires the ability to solve variable-density flows. Several types of variable-density modeling codes have been reported in the literature. Among available codes, the most widely used numerical codes are the U.S. Geological Survey public- domain codes SUTRA [Voss and Souza, 1987] and SEAWAT [Langevin and Guo, 2006]. We selected SEAWAT for simulating our experimental results. The variable-density flow and transport code SEAWAT uses a modified form of the MODFLOW code to solve the groundwater flow equations and the MT3DMS code [Zheng and Bennett, 2002] to solve the transport equations and to simulate the associated density-coupling effects. We chose a uniform grid size of 0.5 cm in our numerical model. Constant head boundaries were set up on the right and left end of the numerical model to describe the head gradient used in the experiments. The model discretization and other general 108 transport parameters are summarized in Table 4-1. Three stress periods were employed to simulate (1) the initial steady-state before injection, (2) injection of tracer or nanoparticles in the domain, and (3) transport of tracer and nanoparticles through the domain. To obtain the steadystate flow conditions, the first stress period was run with constant head boundaries for 200 s. The second stress period was run for 16 s; during this period, 20 mL of tracer or nanoparticles solution was injected into the system through a well assuming a constant flow rate of 1.25 mL/s. The well was turned off in the third stress period to simulate the transport of tracer or nanoparticles. Since the transport was highly advection dominated in our experiments, we used the total variation diminishing (TVD) technique to solve the advection part of the transport equation. This technique helps to minimize numerical dispersion effects [Zheng and Bennett, 2002]. 5.6 Results and Discussion 5.6.1 Stability of INP and S-INP in Vials Figure 5-2 shows the digital pictures of the vials containing pristine INP and S- INP dispersed in the aqueous phase. The pictures show that INP started to aggregate and settled in the vials within 1?10 min (Figure 5-2 a and c). The INP suspension was well segregated and the particles settled at the bottom of vial after 2 days (Figure 5-2 g). However, S-INP was well dispersed and remained as a homogeneous gray-colored solution for up to 60 days (Figure 5-2 h). This qualitative data demonstrates that the stabilizer (PAA) kept the particles suspended. 109 Figure 5-2 Vials containing INP and S-INP at various times: (a) and (b) INP and S-INP after 1 min, (c) and (d) after 10 min, (e) and (f) after 2 h, (g) INP after 2 days, and (h) S- INP after 60 days 5.6.2 Transport of INP and S-INP in Porous Media The transport patterns of the tracer, pristine INP and S-INP observed at different time intervals are shown in Figure 3. The pictures in the first column (time- 0 min) of Figure 5-3 present the location of the plumes just after injection into the porous medium. The other pictures in the figure were taken at the following times: 4, 7, and 10 min. The data indicates that the colored freshwater tracer dispersed and moved horizontally and reached the left boundary in approximately 10 min (see Figure 5-3a). However, S-INP moved vertically downward toward the bottom of the flow container as it migrated horizontally toward the left boundary (see Figure 5-3b). The downward movement of S- INP particles is due to the higher density of S-INP (1.036 g/cm 3 ) compared with water (1 g/cm 3 ). Pristine INP, on the other hand, showed no transport even after 10 min, as indicated by Figure 5-3c. The pristine INP is expected to have positive charge at neutral pH; whereas, the porous medium will be negatively charged. Therefore, INP can attach to the porous medium and become immobilized because of charge interactions. 110 Furthermore, there are other physical processes that can immobilize both INP and S-INP in porous media systems. It is reported in the literature that INP transport can be affected by three distinct transport mechanisms including diffusion, interception, and sedimentation [Schrick et al., 2004; Tufenkji and Elimelech, 2004]. Lecoanet et al. [2004; Lecoanet and Wiesner, 2004] reported that the transport of different nanomaterials (fullerene, single wall nanotube, alumoxane, C60, and ferroxane) in a glass bead column was influenced by both hydrophobic and hydrophilic interactions, blocking of deposition sites, and changes in the attachment efficiency due to ionic strength and steric stabilizations. Figure 5-3 Transport of tracer, pristine INP, and S-INP in the flow container Interestingly, our S-INP transport data did not show much retardation when compared to the tracer data. The higher mobility of S-INP particles can be explained by the association of the hydrophobic part of the PAA with INP and the orientation of the polar head group toward the aqueous phase [Kanel et al., 2007]. Furthermore, unlike INP (which will be positively charged), the S-INP suspension will be negatively charged at neutral pH since it was stabilized using an anionic polymer [Schrick et al., 2004]. Hence, one can expect very little interaction between S-INP and the porous medium (glass 111 beads). In this study, we hypothesized that the observed S-INP transport was primarily controlled by advection and dispersion processes that are coupled to small-scale density gradients. In the section below we provide numerical calculations to test the validity of this hypothesis. Figure 5-4 Comparison of experimental and numerical modeling results for tracer-dye transport. Numerical results show 1 and 10% shaded contour levels. 112 Figure 5-5 Comparison of experimental and numerical modeling results for S-INP transport. Numerical results show 1 and 10% shaded contour levels. The experimental observations were simulated using the variable-density flow code SEAWAT. The results obtained from the numerical model (1 and 10% shaded contour levels) are compared with the experimental results in Figure 5-4 and Figure 5-5. Figure 5-4 compares SEAWAT results against the tracer data and the results show that the model was able to capture the transport patterns very well. Figure 5-5 compares SEAWAT results against S-INP transport data. The model simulations employed a retardation factor value of unity. Therefore, other than the standard advection dispersion processes, the density coupled transport effects is the only additional process included in the simulation. The figure shows that the numerical model was able to accurately predict the sinking effects of the plume and also the overall shape of the plume. These results indicate that our assumption of modeling S-INP transport as a conservative dense plume 113 was indeed correct, and small-scale density gradients can play a significant role in transporting S-INP. In this study we have provided two-dimensional transport data which show that the standard INP is virtually immobile, whereas the PAA-stabilized INP can be transported without any significant retardation. To the best of our knowledge, there is no data set available in the published literature that demonstrates the transport behavior of stabilized nanoparticles in a two-dimensional setting. The two-dimensional data set also indicates that density-driven flow can play an important role in transporting certain classes of nanoparticles into deeper aquifer regions. Since these density-driven sinking effects cannot be fully discerned from conventional one-dimensional column experiments, it is important to characterize nanoparticle transport in multidimensional systems. Our modeling results show that the density-coupled groundwater flow model SEAWAT can be used to predict the movement of S-INP, and hence, the model can potentially be used as a tool for designing nanoparticle injection experiments. However, it is important to note that the physical and chemical heterogeneities inherently present in field soils can interact with S-INP resulting in filtration and clogging of the particles. Therefore, carefully designed feasibility studies should be completed using site soils prior to any field-scale INP injection experiments. 114 CHAPTER 6. ANALYSIS OF SALTWATER CONTAMINATION OF OPEN DUG WELLS IN SRI LANKA AFTER THE 2004 TSUNAMI EVENT- PHYSICAL AND NUMERICAL MODELING EFFORTS 6.1 Introduction The 2004 Asian tsunami destroyed more than 85000 homes and displaced several more thousand residents [Illangasekare et al., 2006]. In addition to the loss of life and property, the tsunami also caused a massive influx of saltwater into the unconfined aquifers underlying the impacted region. Seawater contains high percentage of chloride concentrations; therefore, even small amount of seawater is sufficient (about 2%) to raise the chloride levels of freshwater beyond the recommended standard values [EPA, 2007]. The saltwater influx caused by the tsunami has contaminated several regional aquifers located along the east coast of Sri Lanka. Villholth et. al [2005] conducted a field study after the tsunami and identified three possible sources from which saltwater could have entered the underlying aquifers. These sources, as shown in Figure 6-1, are: (1) infiltration of seawater through the vadose zone during the residence time of tsunami waters on the land surface, (2) entrapment and subsequent leaching of saltwater from low-lying areas or surface depressions, and (3) direct injection through the open dug wells. 115 The open dug wells, found in the areas affected by the tsunami, have a diameter of 1 to 2 meters and a depth of 3-8 meters. These wells serve as the major source of freshwater in most parts of rural Sri Lanka. A typical tsunami-impacted field site with multiple open dug wells is shown in Figure 6-2. In Sri Lanka, the tsunami contaminated an estimated 60000 open wells with saltwater resulting in a crisis for freshwater access [UNEP, 2005]. In an effort to desalinize the wells, various relief agencies started to pump from these wells with an aim to decrease the salinity levels [Illangasekare et al., 2006]. At the time of the tsunami event no guidelines were available for cleaning open wells contaminated with saltwater, therefore, the wells were repeatedly pumped. However, these efforts did not lead to any significant improvement as the salinity levels and the wells remained saline even months after the event [Villholth et al., 2005]. Such field observations have brought to the forefront the need for understanding the well salination and restoration processes which can help develop better methods for managing open well contamination from large scale inundation events like hurricanes and tsunamis. The objective of this study is to address this need. 116 Figure 6-1 Conceptual diagram showing the three possible sources of saltwater infiltration into the unconfined aquifer (modified from Villholth et al., [2005]) Figure 6-2 Tsunami Impacted Field Site in Sri Lanka 117 As shown in Figure 6-1, large scale inundation events can create three types of saltwater contamination sources in a coastal aquifer. Individually, each of the three sources described in the conceptual diagram (Figure 6-1) are potentially unstable plume scenarios since in all these cases the denser saltwater overlies the freshwater in the aquifer [Illangasekare et al., 2006]. List [1965] presented an experimental and theoretical analysis of unstable plume scenarios caused by different variations of Source Type-2. According to his study, there are various parameters that govern the stability of dense aqueous plumes in porous media such as the ambient groundwater flow, the solute concentration, dispersivity etc. Later, various laboratory and theoretical studies have attempted to explain the relationship between stability and these transport parameters [Elder, 1967a; b; Oostrom et al., 1992b; Schincariol et al., 1994; 1997; Wooding et al., 1997a; Wooding et al., 1997b]. However, many of these parameters associated with an unstable dense plume are subjective in definition and hence are very hard, if not impossible, to estimate in the field. During a tsunami event, all the sources presented above were introduced into the groundwater system concurrently, leading to a complex scenario with plumes from differing sources interacting with each other. This interaction makes it difficult to isolate various plumes from their sources. Presence of all these complexities has created a need to study all three source types individually and together. In this work, we focused on Source Type-3 (saltwater inundation through open wells) and conducted careful field and laboratory experiments together with numerical modeling to better understand the transport dynamics of saltwater from this source configuration. 118 6.2 Small-scale laboratory experiments Small-scale laboratory experiments were conducted to first qualitatively analyze the mechanisms of groundwater contamination in a saltwater flooding scenario. The experiments were performed in a two-dimensional sand box model with silica glass beads as the porous media. The complete details of the experimental set-up can be obtained from Hogan [2006]. A summary of the results from these experiments, reported in Illangasekare et. al[2006], show that the saltwater plume emanating from the well moved through the system quite rapidly. An important observation was that although a majority of the saltwater migrated vertically through the well quite rapidly, the well remained contaminated by a small amount of saltwater which is trapped inside the well. Illangasekare et al. [2006] argued based on simple scaling calculations that saltwater should have traveled through the freshwater in Sri Lanka in a time frame of the order of days. However, field observations have later indicated that the residual salinity levels remained high for several months after the tsunami event [Villholth et al., 2005]. In order to better understand the system, we completed a saltwater transport experiment during the Summer of 2006 at an uncontaminated coastal site in Sri Lanka. 6.3 Field Experiments The selected field site was located on the west coast of Sri Lanka in the Kalpitiya peninsular region as presented in Figure 6-3. This field site was unaffected by the tsunami and its geology was very similar to that of the east coast areas which were severely affected by the tsunami. The experiments were conducted near the sea boundary on a pristine beach. Figure 6-4 presents the experimental setup. A platform was constructed to put the outlet in the storage reservoir at a higher elevation than the well 119 opening. The storage tank was manually filled with seawater carried in buckets. The valve on the storage tank was then opened to allow seawater to discharge into the well under gravity simulating a tsunami-type inundation. Two wells were constructed at the field site. The upstream well was used as the injection well where seawater was discharged. The second- downstream well was used as a monitoring well where the migration profile of the seawater injected into the first well was recorded. Figure 6-5 shows the well field used for this study. In this work, we specifically focus on a saltwater injection experimental study conducted using wells 1a and 1b shown in the figure. Figure 6-6 presents a cross-sectional view providing a detailed description of the depth and location of screening of the wells related to this work (Wells 1a and 1b). The experiment was started by injecting three well volumes of saltwater (salinity ~ 50.7 mS/cm), obtained from the nearby ocean, into the injection well (Well 1a). In-situ salinity values were recorded in the wells using a hand-held salinity meter (YSI 30) throughout the experiment. Pumping was conducted by withdrawing one well volume from the injection well using a hand-pump attached to the well. The results from the field experiments are presented and discussed below. 120 Figure 6-3 Selected field site shown on- a) the complete map of Sri Lanka and b) a map of the Kalpitiya peninsular region. Taken from Google Earth ? Figure 6-4 a) Setting up field equipment and b) the final experimental set up Valve Storage tank Hose Well Platform b Platform Well Storage Tank a a b 121 Figure 6-5 Wells installed at the Field Site F.W. Level Screening S.W. Level 1.5 m 1.8 m 2.7 m Well 1 a Well 1 b Sandy aquifer Silt lenses .15 m 20 m Saltwater Interface Figure 6-6 Cross-sectional profile (Wells 1a and 1b) Figure 6-7 shows the salinity profile in Well 1 a; note the pumping events are represented by the dotted lines. The different lines shown in the figure represent the salinity values recorded at the different depths within the well over the period of time (note time is given on log scale). Salinity data for the four shallower depths (7 feet-15 feet) show similar trends and they break away from the salinity profiles observed at the deeper depths (18-19 feet). This observation indicates that the injected saltwater sunk Well 1 a Well 1 b Well 2 a Well 3 a Well 2 b Well 3 b 122 rapidly within the well and reached the deeper parts of the well within a short time. Also, under the high saltwater head condition (during injection), the- top level screens would have allowed an easy pathway for the saltwater to flow out of the well and into the ambient water. After allowing the initial flow out of the wells, these screens would also have acted as conduits allowing freshwater and saltwater to mix. The fate and transport of saltwater injected into Well 1a was observed in the monitoring well (Well 1 b). Figure 6-8 shows the salinity measurements obtained from Well 1b. It can be clearly seen from Figure 6-8 that within few hours (approximately 4 hrs) after saltwater injection, measurable amounts of salinity were observed in the monitoring well. The salinity increased almost exponentially over the next 4 days and peaked at a maximum salinity level of 25 mS/cm. After attaining this peak, the salinity decreased consistently to a relatively constant value of 8 mS/cm during the 30 day monitoring period. An important observation to be made here is that the salinity levels in the monitoring well (well 1 b) remained almost constant at depths shallower than 15 feet. We can observe salinity values only at deeper depths (15 feet-19 feet). Even at these deeper depths, the salinity values are almost the same, implying that as the plume moved through the monitoring well it mixed with the ambient water. Therefore, we can conclude that screening the entire depth of the monitoring well caused the averaging of the salinity profile. The most important observation from this experimental dataset is the long length of time taken by the saltwater (breakthrough curve) to move through the monitoring well. From the experimental results, a plausible conclusion can be drawn that the presence of upstream saltwater source (Well 1 a) could have resulted in the long term salinity recordings in the monitoring well (Well 1 b). 123 Our experimental efforts, in the small-scale laboratory tank as well as in the field site suggested that the saltwater would move rapidly through the coastal groundwater system. However, as stated before, monitoring data from tsunami-impacted sites showed that the wells had high salinity levels even an year after the tsunami event [Villholth et al., 2005]. This implies that other saltwater sources from Type-1 and 2 (Figure 6-1) may be present concurrently with open well source Type-3. Also, there is possibility that the downstream wells got contaminated by the saltwater slug discharged from upstream wells. We wanted to investigate these complex field-scale mutually interactive scenarios using a numerical modeling study. In order to benchmark the numerical simulations, we generated a laboratory dataset which considered a problem similar to the field-scale experimental problem. 0 10000 20000 30000 40000 50000 60000 0.1 1 10 100 1000 Time (Hrs) S p e c i f i c C onduc t i v i t y ( u S / c m ) 7 feet 10 feet 12 feet 15 feet 18 feet 19 feet Figure 6-7 Change in salinity over time for the 14 day inject-and-pump experiment, for well 1a (injection well). Pumping events are shown by dotted lines 124 0 5000 10000 15000 20000 25000 30000 0 100 200 300 400 500 600 700 800 Time (hr) Spe c i f i c Conduc t i v i t y ( u S/ c m ) 7 feet 10 feet 12 feet 15 feet 18 feet 19 feet Figure 6-8 Change in salinity over time for inject-and-monitor experiment. Breakthrough profile from well 1b (monitoring well) 6.4 Collection of Laboratory Data and Numerical Model Validation A three-dimensional flow container was constructed for conducting the physical experiments. The flow container was packed with silica glass beads of average diameter 1.1 mm using the procedure given in Goswami and Clement [2007]. The dimensions of the porous media in the flow container after packing were: 42 cm (length) ? 24.5 cm (width) ? 27 cm (height). Figure 6-9 shows the front view of the flow container after packing. A head gradient of approximately 0.6 cm was created between the two head chambers from left to right to initiate flow across the porous media. The head gradient and the flow were measured regularly to confirm steady-state has been reached at which time the flow measurements became constant. The hydraulic conductivity of the porous 125 media in the flow container was calculated from the head and flow measurements using Darcy?s law. Two sets of experiments were conducted in the physical model: 1) tracer experiment, and 2) saltwater experiment. A slug of water with known physical parameters (conductivity, density, temperature, T.D.S and salinity) was injected at a constant flow rate, using a peristaltic pump, into an upstream well near the left head chamber. The resulting plume was monitored using a conductivity probe at a downstream point near the right head chamber. The locations of the upstream wells and the downstream monitoring points varied for the two sets of experiments. Figure 6-9 Front view of the three-dimensional flow container Conductivity probe (Hach SensION 5?) was used to obtain the quantitative measurements of the movement of the tracer and saltwater at the monitoring points. The probe can measure four physical parameters of the fluid: temperature, conductivity, total dissolved solids (TDS) and salinity. The last three values vary dependent upon the temperature and are therefore corrected for a change in temperature from a known Constant head tank Left head chamber Right head chamber Porous Media 126 standard temperature. In our case a non-linear relationship was chosen to correct for the temperature-varying parameter values. The probe was calibrated using a standard NaCl solution. 6.4.1 Tracer experiment Tracer experiment was conducted to obtain an experimental dataset that could help in calibrating the numerical model. A small quantity of salt was added to freshwater to prepare the tracer solution. Note that the salt was used here as a tracer and the quantity was sufficiently small that it did not change the density significantly. The physical parameters of the tracer fluid were measured using the conductivity probe and are provided in Table 6-1. A plastic tube, open at both ends, with a diameter of 1.5 cm was modified to be used as the injection well by attaching a screen to the bottom for preventing the beads from entering the well. The conductivity probe was also modified slightly by attaching a screen around the probe element to prevent the beads from contacting the element. The location co-ordinates of the injection well and monitoring point for the tracer experiment are provided in Table 6-10. Figure 6-10 shows the experimental set-up for this set of experiments. Table 6-1 Physical parameters of the injected tracer fluid Conductivity Temperature TDS Salinity 12620 ?S/cm 20.4 ?C 6.31 g/l 7.2 ppt 127 Table 6-2 Location of the center of the injection well and monitoring point for the tracer experiment Point X (from left Y (from Z (from water Injection well 4.5 cm 12 cm 6 cm Monitoring point 34.5 cm 11.5 cm 7 cm Figure 6-10 Experimental set-up for the tracer experiment Tracer fluid was pumped into the injection well at a constant flow rate of 0.625 mL/s for duration of 12 minutes. Time of the experiment was taken to start at the start of injection. Fluid parameters were measured at the monitoring point at an interval of at least 1 minute from the start of the injection to 35 minutes into the experiment. The experiment was repeated twice to ensure repeatability. 6.4.2 Numerical Modeling of the tracer data SEAWAT [Langevin and Guo, 2006] was used to simulate the experiments. The parameters and grid discretization used in the model are provided in Table 6-3. Three stress periods were employed to simulate: 1) steady-state condition before injection; 2) injection of tracer through the injection well; and 3) monitoring the migration of tracer a b Conductivity probe Injection well Monitoring point 128 through the flow container. The injection well was simulated using a well boundary condition within SEAWAT with a constant injection rate of 0.625 mL/s for the duration of the second stress period. The monitoring point was simulated using an observation well package provided within SEAWAT/MODFLOW. Table 6-3 Parameters and grid discretization used for SEAWAT MODEL Parameter Value GENERAL Grid discretization (x) 2 cm Grid discretization (y) 2 cm Grid discretization (z) 1 cm Number of columns 21 Number of rows 12 Number of layers 24 Number of stress periods 3 Length or stress periods 600 s, 720 s, 1680 s Porosity 0.39 Hydraulic conductivity 0.99 cm/s FLOW Constant head in left chamber 24.9 cm Constant head in left chamber 24.3 cm Head convergence criterion 0.001 TRANSPORT Advection package TVD Concentration convergence 1e-8 Transport time step 1 s Longitudinal dispersivity 0.1 cm Transverse dispersivity 0.01 cm Vertical dispersivity 0.01 cm 129 6.4.3 Results from tracer experiment The datasets obtained from the physical and numerical experiments for the tracer were normalized to compare the results. Note, when normalized, the conductivity and the salt concentration breakthrough datasets should provide a comparable signal. We can plot the breakthrough profiles of the normalized concentration (from the numerical data) and conductivity or T.D.S (from the physical data) and compare both the results to examine if the there is a good enough match. Figure 6-11and Figure 6-12 show the breakthrough curves from the physical and numerical experiments respectively for the tracer study. The numerical graph has been shown for the breakthrough at five different layers. Since the probe covers multiple layers, we have to select the dataset which provides the best match. From the normalized values, it seems that the most acceptable values for the probe dataset are from the 4 th layer. However, the time at which breakthrough occurred, and both the rising and the falling limbs of the curves matched well for all the layers. 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 5 101520253035 C ondu c t i v i t y : C / C o Time (min) Figure 6-11 Normalized breakthrough curve obtained from the physical tracer experiment. 130 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 5 101520253035 C/ C o Time (min) C/Co 3 C/Co 4 C/Co 5 C/Co 6 C/Co 7 Figure 6-12 Numerical modeling of the tracer experiment- breakthrough curves obtained at different depths 6.4.4 Saltwater injection experiment A set of experiments was conducted by instantaneously injecting a slug of saltwater in to the physical model. The locations of both the injection and monitoring points were modified from the tracer experiment and are provided in Table 6-4. The physical parameters of the injected saltwater are given in Table 6-5. As in the tracer experiment, Hach Sension 5 conductivity probe was used for monitoring the physical parameters. Note, in these experiments, the conductivity probe was directly installed into the porous medium to record the measurements. However, the conductivity probe was slightly modified as shown in Figure 6-13 by putting a screen around it to prevent beads 131 from coming into contact with the probe element. Henceforth, the words ?monitoring point? or ?monitoring probe? shall be used to refer to the conductivity probe. Table 6-4 Locations of the center of the injection well and monitoring point for saltwater experiment Point X (from left Y (from Z (from water Injection well 9 cm 13 cm 6 cm Monitoring point 29 cm 13 cm 10 cm Table 6-5 Physical parameters of the injected saltwater Conductivity Temperature TDS Salinity 56400 ?S/cm 25.5 ?C 28.2 g/l 37.5 ppt Figure 6-13 Conductivity probe shown here with the screen attached around the element 132 Figure 6-14 Setup for the saltwater injection experiment The experiments were conducted in the physical model by injecting 50 ml of saltwater of known physical parameters (conductivity, density, temperature, T.D.S and salinity) at the location specified in Table 1. It took nearly 40 seconds to inject 50 ml of saltwater into the porous medium. The resulting plume was monitored using the conductivity probe at the downstream location. The experiment was repeated to ensure reproducibility and both datasets (with normalized concentration) are shown in Figure 6-15. 0.00 0.05 0.10 0.15 0.20 0.25 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 Time (min) N o r m a l i z e d C o nd uc t i v i t y Exp # 1 Exp # 2 Figure 6-15 Breakthrough curves obtained from the physical experiments. 133 6.4.5 Numerical modeling of saltwater experiment The variable density code SEAWAT was used for numerically modeling the experimental results. Four grid levels were chosen to perform the grid convergence analysis. The details of the grid levels are discussed in Table 6-6. The parameters employed within the numerical model are provided in Table 6-7. Three stress periods were employed to simulate: 1) steady-state condition before injection; 2) injection of saltwater through the injection well; and 3) monitoring the breakthrough of the saltwater tracer through the flow container. The injection well was simulated using a well boundary condition within SEAWAT with a constant injection rate of 1.66 ml/s (50 ml in 30 s) for the duration of the second stress period. The monitoring points were simulated using the observation well package provided in SEAWAT/MODFLOW. Table 6-6 Grid discretization used for the grid convergence analysis Grid Level ? x ? y ? z A- Very Fine 0.5 cm 0.5 cm 0.5 cm B- Fine 1 cm 1 cm 1 cm C- Coarse 2 cm 2 cm 1 cm D- Very Coarse 2 cm 2 cm 2 cm 134 Table 6-7 Parameters and grid discretization used for SEAWAT MODEL Parameter Value GENERAL Grid discretization (x) Dependent on Grid # Grid discretization (y) Grid discretization (z) Number of columns Dependent on Grid # Number of rows Number of layers Number of stress periods 3 Length or stress periods 600 s, 30 s, 1260 s Porosity 0.39 Hydraulic conductivity 0.99 cm/s FLOW Constant head in left chamber 24.9 cm Constant head in left chamber 24.3 cm Head convergence criterion 0.001 TRANSPORT Advection package TVD Concentration convergence 1e-8 Transport time step 1 s Longitudinal dispersivity 0.1 cm Transverse dispersivity 0.01 cm Vertical dispersivity 0.01 cm 6.4.6 Results from the saltwater experiment The breakthrough curves obtained from the numerical simulations of the saltwater experiments are compared here. Figure 6-16 through Figure 6-19 present the normalized concentration calculated at the monitoring point (14 cm depth). It is important to note here that the monitoring probe is approximately 3 cm deep and requires saltwater throughout that depth to measure the conductivity of the fluid. From the figures, it is 135 apparent that the most appropriate layer for comparison with the physical dataset is the 11-12 th layer with a peak normalized conductivity of approximately 0.2. In the case of very fine mesh (Level A) layer 12 has a better match than layer 11 whereas in the case of the fine (Level B) and coarse mesh (Level C) layer 11 has a better match than layer 12. In the case of mesh level D, very coarse mesh, layer 10 has a better match than layer 12. Overall, our results suggest that the mesh discretization can provide results to accuracy of the grid size on the vertical axis. For example, mesh level D has a grid discretization of 2 cm on the z-axis, and it provides a match within 2 cm of the most accurate layer (Layer 10 -12). Similarly, the coarse mesh level C and fine mesh level B provide accuracy to within 1 cm accuracy (Layer 11-12) and the very fine mesh also within 1 cm (Layer 11- 12). Figure 6-20 presents the breakthrough curves obtained at the 11-12 th layer for all four numerical experiments and compares them with the experimental data. A qualitative analysis of Figure 6-20 indicates that the best match is obtained for the 11 th layer of the fine mesh level B or the 12 th layer of the very fine mesh level A. However, the very fine mesh size took almost seven times as long to compute one simulation when compared to the fine mesh (135 minutes v/s 18 minutes). 136 0.00 0.05 0.10 0.15 0.20 0.25 0.30 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 N o r m a l i z e d C onc ( C / C o) Time (s) Layer 22 Layer 24 Layer 26 Figure 6-16 Breakthrough curves obtained from Grid # A- Very Fine. Layer 22 represents a depth of 11 cm. Similarly layers 24, 26 and 28 represent depths of 12 cm, 13 cm and 14 cm respectively. 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 N o r m a l i z e d C o nc ( C / C o) Time (s) Layer 11 Layer 12 Layer 13 Layer 14 Figure 6-17 Breakthrough curves obtained from Grid # B- Fine. Layers 11-14 represent depths of 11- 14 cm respectively. 137 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 Time (s) N o r m aliz e d C o n c ( C /C o ) Layer 12 Layer 13 Layer 14 Layer 11 Figure 6-18 Breakthrough curves obtained from Grid # C- Coarse. Layers 11-14 represent depths of 11- 14 cm respectively. 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 N o r m a liz e d C onc ( C /C o) Time (s) Layer 5 Layer 6 Layer 7 Figure 6-19 Breakthrough curves obtained from Grid # D- Very Coarse. Layers 5-7 represent depths of 10 cm, 12 cm and 14 cm respectively. 138 Experimental vs Numerical 0.00 0.05 0.10 0.15 0.20 0.25 4.0 6.0 8.0 10.0 12.0 14.0 Time (min) N o r m al iz ed V a lu es Very Fine- Layer 12 Fine- Layer 11 Coarse- Layer 11 Very Coarse- Layer 5 Exp # 1 Exp # 2 Figure 6-20 Breakthrough curves obtained from all grid discretizations at corresponding layers and comparison with experimental dataset 6.5 Scenario analysis using the validated numerical model We used validated numerical model to simulate various saltwater migration scenarios related to the well contamination problem. The goal is develop a fundamental understanding to the transport processes involving multiple sources and use this knowledge to develop guidance for well cleaning activities. Our project sponsor International Water Management Institute (IWMI) was specifically interested in developing a concise guidance document that was intended to be endorsed and distributed by WHO. Based on sponsor?s suggestions, the modeling effort was completed with the following specific objectives: 1. Corroborate the field and laboratory findings reported in earlier sections 139 2. Support the development of well clean up guidelines 3. Visualize the saltwater flooding problem within groundwater aquifer 4. Simulate the effect of pumping of wells in an area flooded with saltwater Under natural condition, coastal aquifers will have stable salt-wedge beneath the freshwater. The contamination event, caused by the tsunami caused an ?unnatural? and unstable situation with saltwater lying on top of fresh water. Figure 6-21 shows a conceptual model for a saltwater contamination event such as a tsunami event and its fate after a certain amount of time t. We assume that the tsunami waters entered the aquifer as bulk recharge/infiltration evenly from the top of the aquifer and also as point discharges via open wells, as shown in Figure 6-21. To support the development of the guidelines and to visualize the saltwater flooding processes, four different scenarios were simulated (Table 6-8). Table 6-8 Scenarios simulated with the SEAWAT model Scenario # Purpose Results 1 To visualize the movement of saltwater from a flooding event into the aquifer and wells Figure 6-22 and Figure 6-23 2 To simulate the effect of early pumping of a well under medium vulnerability conditions Figure 6-24 3 To simulate the effect of early pumping of a well under high vulnerability conditions Figure 6-25 4 To simulate the effect of late pumping of a well under medium vulnerability conditions Figure 6-26 140 Figure 6-21 Conceptual model of a saltwater plume moving through an aquifer after a saltwater flooding event. The saltwater originates from the ground surface (orange part) and the well (red part). a. The initial situation. b. After the saltwater plume has sunk into the aquifer 6.5.1 Details of the numerical model To simulate the conceptual processes depicted in Figure 6-21, SEAWAT was set up in three dimensions with parameter values representative of sandy coastal aquifer conditions in Sri Lanka. A base case simulation was used to visualize the processes and Initial Condition After some time 141 simulate a medium vulnerable condition (Scenario 1 and 2, Table 6-8). This setup was slightly modified to describe the other scenarios (Scenario 3 and 4, Table 6-8). The dimensions of the three-dimensional region simulated in the base case were: 100 m (length perpendicular to the shoreline) ? 50 m (deep) ? 5 m (width along the shore line). The grid discretization was at a uniform size of 1 m in all three directions. A constant hydraulic gradient of 1% was simulated between the left (fresh) and right (sea) boundaries. The tsunami wave was simulated using the Generalized Head Boundary (GHB) condition provided within MODFLOW. The head for the GHB was kept at 5 m above the average head in the aquifer to simulate the wave height, and the density of saltwater injected was fixed at 1.025 gm/ml, representing seawater density. The hydraulic conductivity was assumed to be constant and isotropic at 100 m/day. The open well in the model was located at 60 m from the sea-boundary and had a depth of 3 m. The well was hydraulically in contact with the surrounding aquifer through the bottom only. If simulated, the saltwater wedge at the site would be small and located very deep near the right boundary and hence would have negligible effect on the flow field. Therefore, the model was simulated without incorporating the effects of the saltwater wedge. The bottom boundary (at 50 m depth) was a no-flow boundary. Three ?stress periods? were employed to simulate: 1) steady-state flow conditions prior to the tsunami event (2 days of simulation), 2) tsunami inundation event (0.5 hours) and 3) subsequent saltwater movement (20 days). Two simulations were run in the base case, with well injection rates of 6 and 12 well volumes, to explore effects of variation in the well injection on the overall plume movement. 142 A summary of the numerical model parameters for the base case and the modifications for the other scenarios are provided in the technical report submitted to IWMI. 6.5.2 Results from the scenario analysis 6.5.2.1 Scenario 1 The objective of the first numerical experiment was to simulate the transport described by the model and to visualize the sinking of the dense plumes discharged from the top of the aquifer and from the well along a cross-section of the aquifer normal to the ocean boundary, as depicted in Figure 6-21. The results from the simulations were post-processed using GW Vistas. Figure 6-22 and Figure 6-23 present the results from the first and second simulation, respectively, with different volume of saltwater injection into the well. It is seen how the saltwater moves as an irregular blanket from the surface into the aquifer and that the saltwater from the well moves like a separate plume into the aquifer. It is evident from the figures that a higher injection rate (Scenario 1b) gives rise to a bigger plume under the well and also forces the plume to sink deeper into the aquifer. However, the well injection rate has no effect on the plumes from the surface source which remains the same in both simulations. It is interesting to note that the saltwater that discharged from the well migrated rather rapidly into the aquifer (within 5 days) and hence it will be difficult to pump out this contamination source. However, the exact transport pathways and times simulated should be viewed as indicative only, as the model has not been calibrated against field observations 143 Figure 6-22 The movement of saltwater in the aquifer (Scenario 1a): a) immediately after the tsunami, and b) after 5 days. 6 well volumes injected. Figure 6-23 The movement of saltwater into the aquifer (Scenario 1b): a) immediately after the tsunami, and b) after 5 days. 12 well volumes injected. a b a b 144 6.5.2.2 Scenarios 2 and 3 The objective of the second and third numerical experiment was to understand the effects of pumping under medium and high vulnerability conditions. Vulnerability is here governed by the hydraulic conductivity and the depth of the well. In Scenario 2, the hydraulic conductivity is the same as in the base case (Scenario 1), while in Scenario 3 it is doubled (100 m/day to 200 m/day). The well depth was decreased from 1.5 m (Scenario 2) to 0.5 m (Scenario 3). Five stress periods were employed to simulate: 1) steady-state conditions prior to the tsunami event (2 days of simulation), 2) tsunami inundation event (0.5 hours), 3) saltwater movement after the tsunami event (5 days), 4) a pumping event (1 hour), and 5) subsequent saltwater migration through the system (100 days). Two well volumes were pumped out during stress period 4. Simulated saltwater concentration levels against time at the well for scenario 2 are shown in Figure 6-24. The relative concentrations, i.e. the concentrations relative to the initial tsunami water concentration are shown. It can be seen that the concentration drops sharply immediately after the tsunami. This should be expected based on results already shown in Figure 6-22 and Figure 6-23, which indicated that the saltwater injected into the well sunk rather rapidly from the well. However, as time progresses, the saltwater infiltration from the surface source reaches and contaminates the well. The peak salt concentration predicted at the well due to the contamination from the surface is about an order of magnitude less than the initial saltwater concentration. Figure 6-24 describes four distinct stages after the well is contaminated with the tsunami water. Time zero represents a very high concentration level that will be seen in the well immediately after the tsunami event. This will be followed by a period where the 145 concentration will fall rapidly as the well water sinks into the aquifer. Pumping during this period will further remove some salinity. Finally, the figure shows the arrival of the surface-infiltrated saltwater source into the well, which can last for a substantial period of time, depending on the physical characteristics of the aquifer and the infiltration source. 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 2040608010120 Time (days) Relative Con c Pumping Figure 6-24 Effect of early pumping and subsequent re-entry of saltwater from the surface source under medium vulnerability conditions (Scenario 2) For Scenario 3, where the vulnerability of the system and the well is increased, it is expected that the overall transport in the system will be faster and hence the infiltrated water will reach the well more rapidly. Furthermore, decreasing the depth of the well would make it more vulnerable to direct contamination from the surface-source; hence the concentrations at the well are expected to be higher than those predicted in the previous experiment. 146 The simulated breakthrough obtained in the well is show in Figure 6-25. In a high conductivity system, the saltwater concentration is not decreased to zero as it occurred for the low conductivity system (Scenario 2- Figure 6-24). The relative concentration level had a much higher peak concentration of about 0.3 for Scenario 3 (Figure 6-26). This can be primarily attributed to the shallow depth of the depth. Also, overall the breakthrough curve has a shorter residence time in the well for this experiment, which is due to the faster transport velocities that resulted from the higher conductivity field. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 2040608010 Time (days) Relati ve Co n c Pumping Figure 6-25 Effect of early pumping under high vulnerability condition (high hydraulic conductivity and shallow well) (Scenario 3). Ingress of the surface source is very rapid 147 6.5.2.3 Scenario 4 The objective of Scenario 4 was to investigate the effect of late pumping on the well water salinity. The conditions were similar to those in Scenario 2, except the well is pumped a second time - after the arrival of the infiltration source, at about 5 weeks after the tsunami event. Also, the depth of the well was decreased from 1.5 meter to 1.0 meter to observe breakthrough at a shallower region. The relative predicted saltwater concentration levels against time are shown in Figure 6-26. The generated profiles and the associated transport processes were very similar to the ones shown for Scenario 2 (Figure 6-24) until the second pumping event. It can be seen from the figure that the new pumping event simulated after 5 weeks decreases the concentration in the well for a short period. However, after an immediate decrease, the concentration comes right back as the ambient saltwater reenters the well. The modeling results clearly demonstrate the futility of delayed pumping events. 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 2040608010 Time (days) Re l a t i v e C o n c 2nd Pumping1st Pumping Time between pumpings Figure 6-26 Effect of late pumping (Scenario 4). The re-entry of surface-derived saltwater into the well after pumping is almost immediate 148 6.6 Conclusions and Recommendations The results of our physical and numerical experiments were used to develop the well cleaning guidelines document [Villholth, 2008]. Some general conclusions and recommendations are provided here. a. Laboratory experiments have shown that an isolated open well, affected by saltwater inundation, can be cleaned effectively by pumping. However, multiple pumping events were required in the field experiment to bring down the salinity levels to drinkable limits (~1000 ?S/cm) b. The movement of saltwater from the upstream source to a hydraulically connected well was captured in the small-scale and field experiments. Breakthrough curves obtained from the monitoring wells in these experiments shows that salinity levels increase to a peak value and then decrease. c. By continuously monitoring the open wells, we can estimate the time period when the salinity levels start to decline. This would represent the trailing (decreasing) part of the breakthrough curve indicating that the net salinity delivered from the upstream sources is beginning to decrease. d. It is important to realize that over-pumping can cause the well to be hydraulically linked to other sources of saline water (e.g., salt-wedge or the ocean). Therefore excessive pumping should not be conducted to avoid such problems. e. A quick indicator of over-pumping is the increase in salinity level of the extracted water. Therefore, salinity levels of the water pumped out from the well should be continuously monitored. A sample of water (enough for getting a good measurement from the salinity meter) should be checked for increase in salinity 149 after a fraction (0.3 or 0.4) of well volume has been pumped out. If the salinity increases then pumping should be stopped immediately. f. It is recommended that only one well volume be removed after the saltwater contamination event. The extracted water should be discharged away from the well so as to avoid recirculation of contaminated water in the well. g. Any further pumping (more than one well volume) may act to hydraulically link the well to upstream salinity sources. h. In most cases, the best technique for cleaning the well is to leave it for natural restoration. Depending on climatic factors (such as rainfall) and hydraulic properties of the soil, it may take more than a year to obtain background freshwater conditions in the wells. 150 CHAPTER 7. SUMMARY, CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK 7.1 Summary The goal of this research effort is to investigate the basics of various types of variable-density flow and transport problems that are relevant to environmental engineering and management. Chapters 2-4 involve basic research into benchmarking variable-density flow problems. Benchmarking datasets in the form of steady-state and transient stable density-configuration experiments were provided for testing variable- density flow codes in Chapter 2. It was shown that the presented benchmarking datasets were better than the often-used Henry problem. Quantitative techniques were not needed to estimate concentrations in experiments presented in Chapter 2 since the transition zone between freshwater and saltwater was narrow. However, variable-density experiments can have larger transition zone which may make it hard to estimate concentration distributions in the physical experiments. Therefore, an image analysis (IA) technique was developed for measuring concentrations in porous media experiments. A detailed analysis of error-estimation in the application of IA was conducted and is presented in Chapter 3. It was shown that the generally accepted methods of estimating errors, namely the mass-balance and dispersion coefficient comparison, have limitations. A statistics-based technique was also presented to quantify the errors in application of IA to porous media systems. In Chapter 4 we apply the image analysis methods to conduct a 151 set of unstable density-configuration experiments. These physical experiments offer a unique combination of datasets involving transport with and without instabilities. Such a combination offers a unique opportunity to rigorously test the applicability of numerical approaches. Results from numerical modeling showed that none of the option available in the USGS code SEAWAT can accurately simulate the two experimental datasets. Therefore, Chapter 4 highlights the need to develop and/or identify the appropriate numerical approaches that should be used to simulate variable-density flow problems. Chapter 5 and 6 discuss two applied problems involving variable-density flow scenarios. Chapter 5 describes the experimental results of stabilized iron nano-particle (S-INP) transport in a porous media tank. Variable-density flow was found to play a dominant role in the transport mechanism of the S-INP. Numerical modeling using variable-density approaches provided a good match with experimental results demonstrating a novel use of variable-density flow codes. Chapter 6 outlines the work conducted through a joint project effort with the International Water Management Institute (IWMI) to study the effect of saltwater contamination on well salination. As a part of this effort, field experiments were conducted in Sri Lanka at a pristine coastal aquifer. Physical experiments were also conducted in the laboratory in a small-scale porous media tank. The results of the field, laboratory, and modeling studies were used together to develop a set of guidelines of managing the contamination of open wells after large scale inundation events such as hurricanes and tsunamis. 7.2 Suggestions for future work Based on the ideas developed in this dissertation the following research efforts can be pursued for future work: 152 1. There is a lack of combined heat and solute transport variable-density flow experiments in the literature. Such experiments can be conducted using the techniques presented here. Similarly, other novel flow and transport experiments can be conducted with the experimental setup described here, for example, nano- particle transport experiments discussed in Chapter 5. 2. Filtering algorithms used in image analysis (IA) techniques are based on averaging techniques. Such averaging can spuriously disperse the concentration profile observed in the physical experiment. IA is a powerful tool and its use in porous media experiments may have widely accepted in the future. Therefore, it is important to study the effect of filtering algorithms on resulting concentrations and their subsequent use and interpretation. 3. It was suggested in Chapter 4 that heterogeneity can accurately predict the generation of instabilities in variable-density experiments. However, the extent of variability to be used in generating heterogeneous fields was not explored. There is a need to analyze the effect of heterogeneity on generation of instability using results from the combined problem of rising and sinking plume experiments. 4. Although heterogeneity can correctly estimate the generation of instabilities it fails to provide a narrow transition zone due to artificial numerical dispersion. Similarly, MOC (method of characteristics) and TVD (total variation diminishing) methods also have limitations as discussed in Chapter 4. Thus, there is a general need to develop better numerical methods or approaches to simulate variable- density flow problems. 153 REFERENCES Bakker, M., and F. Schaars (2003), The Sea Water Intrusion (SWI) Package Manual, version 0.2. Athens, Georgia: University of Georgia. http://www.engr.uga.edu/~mbakker/swi.html. Bear, J., and G. Dagan (1964), Moving interface in coastal aquifers, ASCE Journal of the Hydraulics Division, 90(HY4), 193-216. Brakefield, L. K. (2008), Physical and Numerical Modeling of Buoyant Groundwater Plumes, MS Thesis, 107 pp, Auburn University, Auburn. Bridge, J. W., et al. (2006), Noninvasive quantitative measurement of colloid transport in mesoscale porous media using time lapse fluorescence imaging, Environmental Science & Technology, 40(19), 5930-5936. Bues, M. A., and C. Oltean (2000), Numerical simulations for saltwater intrusion by the mixed hybrid finite element method and discontinuous finite element method, Transport in Porous Media, 40(2), 171-200. 154 Burghardt, D., and A. Kassahun (2005), Development of a reactive zone technology for simultaneous in situ immobilisation of radium and uranium, Environmental Geology, 49(2), 314-320. Canon, U. (2008), http://www.usa.canon.com/consumer/controller?act=ModelInfoAct& fcategoryid=145&modelid=12913. Cao, J. S., et al. (2005), Perchlorate reduction by nanoscale iron particles, Journal of Nanoparticle Research, 7(4-5), 499-506. Carlson, D. A., et al. (2008), Storm-damaged saline-contaminated boreholes as a means of aquifer contamination, Ground Water, 46(1), 69-79. Catania, F., et al. (2008), Assessment of quantitative imaging of contaminant distributions in porous media, Experiments in Fluids, 44(1), 167-177. Clement, T. P., et al. (2004), Understanding the dynamics of groundwater flow and reactive transport patterns near a saltwater boundary, NSF proposal, edited. Croucher, A. E., and OSullivan (1995), The Henry problem for saltwater intrusion, Water Resources Research, 31(7), 1809-1814. 155 Darnault, C. J. G., et al. (1998), Visualization by light transmission of oil and water contents in transient two-phase flow fields, Journal of Contaminant Hydrology, 31(3-4), 337-348. Darnault, C. J. G., et al. (2001), Measurement of fluid contents by light transmission in transient three-phase oil-water-air systems in sand, Water Resources Research, 37(7), 1859-1868. Detwiler, R. L., et al. (1999), Measurement of fracture aperture fields using transmitted light: An evaluation of measurement errors and their influence on simulations of flow and transport through a single fracture, Water Resources Research, 35(9), 2605-2617. Dong, W. J., and A. R. S. Selvadurai (2006), Image processing technique for determining the concentration of a chemical in a fluid-saturated porous medium, Geotechnical Testing Journal, 29(5), 385-391. Elder, J. W. (1967a), Steady free convection in a porous medium heated from below, J. Fluid Mech., 27(1), 29-48. Elder, J. W. (1967b), Transient convection in a porous medium, J. Fluid Mech., 27(3), 609-623. 156 Elliott, D. W., and W. X. Zhang (2001), Field assessment of nanoscale biometallic particles for groundwater treatment, Environmental Science & Technology, 35(24), 4922- 4926. EPA, U. (2007), Current National Recommended Water Quality Criteria, edited. Fischer, H. B., et al. (1979), Mixing in inland and coastal waters, Academic Press, Inc. Goswami, R. R., and T. P. Clement (2007), Laboratory-scale investigation of saltwater intrusion dynamics, Water Resources Research, 43(4). Goswami, R. R., et al. (2008), Estimating error in Concentrations Measurements obtained from Image Analysis, Vadose Zone Journal, in press. Hansen, P. C., et al. (2006), Deblurring Images: Matrices, Spectra and Filtering, Society for Industrial and Applied Mathematics. He, F., and D. Y. Zhao (2005), Preparation and characterization of a new class of starch- stabilized bimetallic nanoparticles for degradation of chlorinated hydrocarbons in water, Environmental Science & Technology, 39(9), 3314-3320. He, F., et al. (2007), Stabilization of Fe-Pd nanoparticles with sodium carboxymethyl cellulose for enhanced transport and dechlorination of trichloroethylene in soil and groundwater, Industrial & Engineering Chemistry Research, 46(1), 29-34. 157 Henry, H. R. (1960), Salt intrusion into coastal aquifers, Ph.D. Thesis, Columbia University. Herbert, A. W., et al. (1988), Coupled groundwater flow and solute transport with fluid density strongly dependent on concentration, Water Resources Research, 24(10), 1781- 1795. Hogan, M. B. (2006), Understanding the flow and mixing dynamics of saline water discharged into coastal freshwater aquifers, MS Thesis, 114 pp, Auburn University, Auburn. Huang, W. E., et al. (2002), Physical modelling of solute transport in porous media: evaluation of an imaging technique using UV excited fluorescent dye, Water Research, 36(7), 1843-1853. Huber, D. L. (2005), Synthesis, properties, and applications of iron nanoparticles, Small, 1(5), 482-501. Illangasekare, T., et al. (2006), Impacts of the 2004 tsunami on groundwater resources in Sri Lanka, Water Resources Research, 42(5). 158 Inspectorate, S. N. P. (1986), HYDROCOIN- an international project for studying groundwater hydrology modeling srategies, Level 1 Final Report, Stockholm. Jain, A. K. (1988), Fundamentals of digital image processing, Prentice Hall. Johannsen, K., et al. (2002), The saltpool benchmark problem - numerical simulation of saltwater upconing in a porous medium, Advances in Water Resources, 25(3), 335-348. Joo, S. H., et al. (2004), Oxidative degradation of the carbothioate herbicide, molinate, using nanoscale zero-valent iron, Environmental Science & Technology, 38(7), 2242- 2247. Joo, S. H., et al. (2005), Quantification of the oxidizing capacity of nanoparticulate zero- valent iron, Environmental Science & Technology, 39(5), 1263-1268. Kalejaiye, B. O., and S. S. S. Cardoso (2005), Specification of the dispersion coefficient in the modeling of gravity-driven flow in porous media, Water Resources Research, 41. Kanel, S. R., et al. (2005), Removal of arsenic(III) from groundwater by nanoscale zero- valent iron, Environmental Science & Technology, 39(5), 1291-1298. 159 Kanel, S. R., et al. (2006), Arsenic(V) removal kom groundwater using nano scale zero- valent iron as a colloidal reactive barrier material, Environmental Science & Technology, 40(6), 2045-2050. Kanel, S. R., et al. (2007), Transport of surface-modified iron nanoparticle in porous media and application to arsenic(III) remediation, Journal of Nanoparticle Research, 9(5), 725-735. Kanel, S. R., et al. (2008), Two dimensional transport characteristics of surface stabilized zero-valent iron nanoparticles in porous media, Environmental Science & Technology, 42(3), 896-900. Kolditz, O., et al. (1998), Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models, Advances in Water Resources, 21(1), 27-46. Lajeunesse, E., et al. (1999), Miscible displacement in a Hele-Shaw cell at high rates, Journal of Fluid Mechanics, 398, 299-319. Langevin, C. D., et al. (2003), MODFLOW-2000, the U.S. Geological Survey modular ground-water model: Documentation of the SEAWT-2000 version with the variable- density flow processes (VDF) and the integrated MT3DMS Transport Processes (IMT). U.S. Geological Survey Open-File Report 03-426. 160 Langevin, C. D., and W. X. Guo (2006), MODFLOW/MT3DMS-based simulation of variable-density ground water flow and transport, Ground Water, 44(3), 339-351. Lecoanet, H. F., et al. (2004), Laboratory assessment of the mobility of nanomaterials in porous media, Environmental Science & Technology, 38(19), 5164-5169. Lecoanet, H. F., and M. R. Wiesner (2004), Velocity effects on fullerene and oxide nanoparticle deposition in porous media, Environmental Science & Technology, 38(16), 4377-4382. Li, L., et al. (2006), Synthesis, properties, and environmental applications of nanoscale iron-based materials: A review, Critical Reviews in Environmental Science and Technology, 36(5), 405-431. List, E. J. (1965), The stability and mixing of a densty-stratified horizontal flow in a saturated porous medium, 164 pp, California Institute of Technology. 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(1-4), 126-142. 161 Mainhagu, J., et al. (2007), Measurement by laser induced fluorescence on miscible density driven flows in a Hele-Shaw cell: settings and preliminary results, Comptes Rendus Mecanique, 335(2), 105-112. Manning, B. A., et al. (2007), Spectroscopic investigation of Cr(III)- and Cr(VI)-treated nanoscale zerovalent iron, Environmental Science & Technology, 41(2), 586-592. McNeil, J. D., et al. (2006), Quantitative imaging of contaminant distributions in heterogeneous porous media laboratory experiments, Journal of Contaminant Hydrology, 84(1-2), 36-54. Menand, T., and A. W. Woods (2005), Dispersion, scale, and time dependence of mixing zones under gravitationally stable and unstable displacements in porous media, Water Resources Research, 41(5). Mualem, Y., and J. Bear (1974), The shape of the interface in steady flow in a stratified aquifer, Water Resources Research, 10(6), 1207-1215. Niemet, M. R., and J. S. Selker (2001), A new method for quantification of liquid saturation in 2D translucent porous media systems using light transmission, Advances in Water Resources, 24(6), 651-666. 162 Nurmi, J. T., et al. (2005), Characterization and properties of metallic iron nanoparticles: Spectroscopy, electrochemistry, and kinetics, Environmental Science & Technology, 39(5), 1221-1230. Oltean, C., et al. (2004), Transport with a very low density contrast in Hele-Shaw cell and porous medium: Evolution of the mixing zone, Transport in Porous Media, 55(3), 339- 360. Oostrom, M., et al. (1992a), Dispersivity Values Determined from Effluent and Nonintrusive Resident Concentration Measurements, Soil Science Society of America Journal, 56(6), 1754-1758. Oostrom, M., et al. (1992b), Behaviour of dense aqueous phase leachate plumes in homogenous porous media, Water Resources Research, 28(8), 2123-2134. Oostrom, M., et al. (2007), A review of multidimensional, multifluid, intermediate-scale experiments: Flow behavior, saturation imaging, and tracer detection and quantification, Vadose Zone Journal, 6(3), 610-637. Oswald, S. E., and W. Kinzelbach (2004), Three-dimensional physical benchmark experiments to test variable-density flow models, Journal of Hydrology, 290, 22-42. 163 Parker, L., et al. (2006), Observations of gas flow in porous media using a light transmission technique, Water Resources Research, 42(5). Persson, M. (2005a), Accurate dye tracer concentration estimations using image analysis, Soil Science Society of America Journal, 69(4), 967-975. Persson, M. (2005b), Estimating surface soil moisture from soil color using image analysis, Vadose Zone Journal, 4(4), 1119-1122. Phenrat, T., et al. (2007), Aggregation and sedimentation of aqueous nanoscale zerovalent iron dispersions, Environmental Science & Technology, 41(1), 284-290. Ponder, S. M., et al. (2000), Remediation of Cr(VI) and Pb(II) aqueous solutions using supported, nanoscale zero-valent iron, Environmental Science & Technology, 34(12), 2564-2569. Poots, G. (1958), Heat transfer by laminar free convection in enclosed plane gas layers, Quarterly Journal of Mechanics and Applied Mathematics, 11(3), 257-273. Quinn, J., et al. (2005), Field demonstration of DNAPL dehalogenation using emulsified zero-valent iron, Environmental Science & Technology, 39(5), 1309-1318. 164 Reilly, T. E., and A. S. Goodman (1985), Quantitative-Analysis of Saltwater Fresh-Water Relationships in Groundwater Systems - a Historical-Perspective, Journal of Hydrology, 80(1-2), 125-160. Rumbaugh, J., and D. Rumbaugh (2007), Groundwater Vistas, edited. 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(10), 2317-2329. Schincariol, R. A., et al. (1993), On the Application of Image-Analysis to Determine Concentration Distributions in Laboratory Experiments, Journal of Contaminant Hydrology, 12(3), 197-215. Schincariol, R. A., et al. (1994), ON THE GENERATION OF INSTABILITIES IN VARIABLE-DENSITY FLOW, Water Resources Research, 30(4), 913-927. Schincariol, R. A., et al. (1997), Instabilities in variable density flows: Stability and sensitivity analyses for homogeneous and heterogeneous media, Water Resources Research, 33(1), 31-41. 165 Schincariol, R. A. (1998), Dispersive mixing dynamics of dense miscible plumes: natural perturbation initiation by local-scale heterogeneities, Journal of Contaminant Hydrology, 34(3), 247-271. Schrick, B., et al. (2004), Delivery vehicles for zerovalent metal nanoparticles in soil and groundwater, Chemistry of Materials, 16(11), 2187-2193. Segol, G. (1994), Classic groundwater simulations: proving and improving numerical models, Prentice Hall Inc., New Jersey. Silliman, S. E., et al. (1998), The use of laboratory experiments for the study of conservative solute transport in heterogeneous porous media, Hydrogeology Journal, 6(1), 166-177. Simmons, C. T., et al. (1999), On a test case for density-dependent groundwater flow and solute transport models: The salt lake problem, Water Resources Research, 35(12), 3607- 3620. Simmons, C. T., et al. (2001), Variable-density groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges, Journal of Contaminant Hydrology, 52(1-4), 245-275. 166 Simmons, C. T., et al. (2002), Laboratory investigation of variable-density flow and solute transport in unsaturated-saturated porous media, Transport in Porous Media, 47(2), 215-244. Simmons, C. T. (2005), Variable density groundwater flow: From current challenges to future possibilities, Hydrogeology Journal, 13(1), 116-119. Simpson, M. J., and T. P. Clement (2003), Theoritical analysis of the worthiness of the Henry and Elder problems as benchmarks of density-dependent groundwater flow models, Advances in Water Resources, 26, 17-31. Simpson, M. J., and T. P. Clement (2004), Improving the worthiness of the Henry problem as a benchmark for density-dependent groundwater flow models, Water Resources Research, 40(1), W01504, doi:01510.01029/02003WR002199. Swartz, C. H., and F. W. Schwartz (1998), An experimental study of mixing and instability development in variable-density systems, Journal of Contaminant Hydrology, 34(3), 169-189. Taylor, J. R. (1997), An introduction to error analysis, Second Edition ed., University Science Books. 167 Theodoropoulou, M. A., et al. (2003), A new visualization technique for the study of solute dispersion in model porous media, Journal of Hydrology, 274(1-4), 176-197. Tidwell, V. C., and R. J. Glass (1994), X-Ray and Visible-Light Transmission for Laboratory Measurement of 2-Dimensional Saturation Fields in Thin-Slab Systems, Water Resources Research, 30(11), 2873-2882. Todd, D. K., and L. W. Mays (2005), Groundwater Hydrology, John Wiley & Sons. Tufenkji, N., and M. Elimelech (2004), Correlation equation for predicting single- collector efficiency in physicochemical filtration in saturated porous media, Environmental Science & Technology, 38(2), 529-536. UNEP, U. N. E. P. (2005), After the Tsunami: Rapid Environmental Assessment (available at http://www.unep.org/tsunami/reports/TSUNAMI_report_complete.pdf), edited. USGS (2000), Groundwater resources for the future- Atlantic Coastal Zone, USGS Fact Sheet, 085-00. Vanderborght, J., et al. (2002), Imaging fluorescent dye concentrations on soil surfaces: Uncertainty of concentration estimates, Soil Science Society of America Journal, 66(3), 760-773. 168 Vappicha, V. N. (1975), Steady and Transient Interfaces in Sea Water Intrusion into Coastal Aquifers. Ph.D. Thesis, Indian Institute of Technology, Bombay, 202. Villholth, K. G., et al. (2005), Tsunami impacts on shallow groundwater and associated water supply on the east coast of Sri Lanka, IWMI internal report. Villholth, K. G. (2008), Cleaning wells after seawater flooding- Emergency guidelines. Voss, C. I., and W. R. Souza (1987), Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater-saltwater transition zone, Water Resources Research, 23(10), 1851-1866. Weisbrod, N., et al. (2004), Migration of saline solutions in variably saturated porous media, Journal of Contaminant Hydrology, 72(1-4), 109-133. Weissberger, and Rossiter Physical methods of chemistry. Wooding, R. A., et al. (1997a), Convection in groundwater below an evaporating salt lake .1. Onset of instability, Water Resources Research, 33(6), 1199-1217. Wooding, R. A., et al. (1997b), Convection in groundwater below an evaporating salt lake .2. Evolution of fingers or plumes, Water Resources Research, 33(6), 1219-1228. 169 Woods, J. A. (2004), Numerical accuracy of variable-density groundwater flow and solute transport simulations, PhD Dissertation, 213 pp, Univeristy of Adelaide, Adelaide. Wu, L., et al. (2005), Preparation of cellulose acetate supported zero-valent iron nanoparticles for the dechlorination of trichloroethylene in water, Journal of Nanoparticle Research, 7(4-5), 469-476. Yang, G. C. C., and H. L. Lee (2005), Chemical reduction of nitrate by nanosized iron: Kinetics and pathways, Water Research, 39(5), 884-894. Younes, A., et al. (1999), Modeling variable density flow and solute transport in porous medium: 2. Re-evaluation of the salt dome flow problem, Transport in Porous Media, 35(3), 375-394. Yu, Z. B., and F. W. Schwartz (1999), Determining concentration fields of tracer plumes for layered porous media in flow-tank experiments, Hydrogeology Journal, 7(2), 236- 240. Zhang, Q., et al. (2001), Influence of seaward boundary condition on contaminant transport in unconfined coastal aquifers, Journal of Contaminant Hydrology, 49(3-4), 201-215. 170 Zhang, W. X. (2003), Nanoscale iron particles for environmental remediation: An overview, Journal of Nanoparticle Research, 5(3-4), 323-332. Zheng, C., and G. D. Bennett (2002), Applied Contaminant Transport Modeling, 2nd ed., John Wiley and Sons, Inc, New York. Zimmerman, D. A., and J. L. Wilson (2005), Description of and user's manual for TUBA: A computer code for generating two-dimensional random fields via the turning bands method, edited, SETA, Inc., 27 Ponderosa Drove, Suite 100, Cedar Crest, New Mexico 87008, New Mexico. Zinn, B., et al. (2004), Experimental visualization of solute transport and mass transfer processes in two-dimensional conductivity fields with connected regions of high conductivity, Environmental Science & Technology, 38(14), 3916-3926. 171 APPENDIX-1 IMAGE ANALYSIS PROCEDURE Four major steps were involved in development and application of the IA technique to measure concentration. In the first step we obtained calibration datasets for relating image properties with concentration. A filtering procedure was then selected and applied to reduce errors in the calibration datasets in the second step. In the third step, a relationship was established that correlates image properties to concentration based on the calibration datasets. Finally, the developed IA technique was applied to variable- density flow experiments to obtain concentration measurements. Description of the physical experiments is also provided in this section. Calibration Datasets Solutions containing different concentrations of the red dye ranging from 4 mL of dye per liter of water (4mL/L) to 0.25 mL/L were prepared. Each of these dye solutions was flushed through the tank to obtain a constant image intensity field. More than 10 times the tank pore volume was used for each flush to ensure that the entire tank was filled with the intended dye concentration. This process of obtaining the calibration datasets was originally proposed by Schincariol et al. [1993] and later adopted by others [Huang et al., 2002; McNeil et al., 2006; Swartz and Schwartz, 1998]. Both the experiments were conducted on different days and as such the conditions could have varied slightly. 172 Therefore, calibration datasets were obtained for both the experiments. A filtering algorithm was applied to reduce the noise in the obtained images. Median filtering was selected and applied on the acquired digital images to reduce the standard deviations in image intensities. Median filtering has been used in the past in the IA process [Schincariol et al., 1993]. Also errors introduced by pixel shifting (due to camera movement) can be removed by median filtering [Vanderborght et al., 2002]. Median filtering replaces the intensity values in each pixel of a unit region of the image with the statistical median intensity value of all pixels in that region. By applying this process over the entire spatial domain of the image, standard deviations are reduced. To calculate the spatial extent of the unit region, it is important to know the spatial resolution of the original image data before applying filtering. The spatial resolution of the digital images obtained by our camera was 0.0118 mm 2 per pixel. Calibration Curve In the past, researchers have related the captured image property and targeted fluid property by using curve-fitting procedures. [Schincariol et al., 1993] used a polynomial approximation to find an equation relating optical density (image property) to the fluid concentration (of Rhodamine dye) and porous media (glass beads) size. [McNeil et al., 2006] utilize a IA technique similar to [Schincariol et al., 1993] and used the power law function to obtain a relationship between optical density and image intensity. Weisbrod et al. [2004] used a polynomial series fit to relate saturation (fluid property) with light transmission (image property). Yu and Schwartz [1999] chose a third order polynomial fit from linear, quadratic, fourth order and exponential series fit to equate 173 image properties (pixel intensity, optical density) with fluid concentration. Huang et al. [2002] used a spline relationship between pixel intensity and fluid concentration. A few researchers have also used optical laws as a basis to derive a relationship between the observed image property and the measured fluid property. Oltean et al. [2004] used the Bouguer-Beer-Lambert?s Law (B-B-L) [Weissberger and Rossiter] and calculated a theoretical relationship between concentration and intensity. However, their theoretical analysis was limited to Hele-Shaw cells and is not applicable to porous media. Furthermore, their theoretical relationship did not match well with the experimental values at higher normalized concentration values. They attribute the non-confirming behavior of the experimental values to possible non-uniform lighting and employed quadratic fitting to obtain a better relationship. Lajeunesse et al. [1999] also used the B- B-L for relating pixel intensity with dye concentrations in a Hele-Shaw experiment. Calibration-relationship given by the power-law was obtained using the curve- fitting toolbox available in MATLAB?. Cropped images of the variable-density experiment were also filtered using the same pixel widths yielding three distinct filtered intensity profiles. Conversion of the intensity profiles to concentration was then carried out using the calibration-relationship. The concentration obtained using the calibration relationship is the dye concentration. Concentration of salt was calculated from the concentration of dye by assuming that salt and dye move together. Salt concentrations range from 0 kg/m 3 to 0.0357 kg/m 3 ml/l corresponding to dye concentrations of 4 ml/l to 0 ml/l. Therefore, we get the following relationship between dye concentration and salt concentration: 174 We used the power-law function to relate our experimental dataset with concentrations: C dye = K 1 ? I n + K 2 (1) Regression analysis was conducted using the calibration dataset to determine the values of K 1 , K 2 and n. Since the flow was advection-dominated it was assumed that dye migrates at the same rate as the salt. Therefore a direct linear correlation was derived between the salt concentration (C salt ) and dye concentration (C dye ) to obtain the concentration of salt from the experiments.