Development and Application of One Dimensional Multi-component Reactive Transport Models by Jagadish Torlapati 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 May 4, 2013 Keywords: reactive transport, genetic algorithms, parallel computing, OpenMP Copyright 2013 by Jagadish Torlapati Approved by T. Prabhakar Clement, Committee chair, Groome Jr. Endowed Professor of Civil Engineering Ahjeong Son, Committee member, Assistant Professor of Civil Engineering Jose G Vasconcelos, Committee member, Assistant Professor of Civil Engineering Mark Dougherty, Committee member, Associate Professor of Biosystems Engineering ii Abstract The overarching goal of this dissertation is to develop reactive transport models and explore their applications to groundwater remediation problems. The primary focus of this dissertation is aimed at developing models that can support laboratory studies investigating remediation strategies, as this is an important intermediate step before the remediation methods can be scaled up to apply at field sites. As a part of this research effort, a comprehensive, one- dimensional, multi-component reactive transport model, RT1D, which can be used for simulating biochemical and geochemical reactive transport problems, has been developed. The code can be run within the standard Microsoft EXCEL Visual Basic platform and it does not require any additional software tools. The capabilities of the tool were illustrated by solving several benchmark problems taken from the literature that have varying levels of reaction complexity. These literature-derived benchmarks were used to highlight the versatility of the code for solving a variety of practical reactive transport problems. This model was subsequently applied to a published experimental dataset that described bioaugmentation processes to remediate PCE-DNAPL trapped in a fracture system. A mathematical framework was first formulated to model the bioremediation processes in a PCE contaminated single fracture system augmented with Dehalococcoides Sp. (DHC). The mathematical framework describes multi-species bioreactive transport processes that include bacterial growth and detachment dynamics, biodegradation of chlorinated species, competitive inhibition of various reactive species, and the loss of daughter products due to back-partitioning iii effects. Two sets of experimental data, available in Schaefer et al. (2010b), were used to calibrate and test the model. The simulation results indicate that the yield coefficient and the DHC maximum utilization rate coefficient were the two important process parameters. A detailed sensitivity study was completed to quantify the sensitivity of the model to variations in these two parameter values. The proposed model provides a rational mathematical framework for simulating remediation systems that employ DHC bioaugmentation for restoring chlorinated solvent contaminated groundwater aquifers. While calibrating the DHC bioaugumentation model, several inefficiencies related to the use of trial and error methods for parameter estimation were identified. In order to improve the efficiency of the parameter estimation process, a parallel genetic algorithm (PGA) was developed to automate the parameter estimation process. The performance of the PGA was tested by solving four benchmark problems that have published experimental data or analytical/numerical solutions. Benchmarking results indicate that the PGA estimated parameters are close to the true parameters. A shared memory parallel computing platform that utilized OpenMP FORTRAN was used to demonstrate the speedup of the code on a four processor desktop Pentium computer. The parallelized code showed linear speedup with increasing number of processors. The PGA routines used in this study are generic and can be easily adapted to solve parameter estimation problems in other environmental modeling applications. iv Acknowledgments I would like to acknowledge my advisor Dr. Prabhakar Clement for his continued support and guidance. I would also like to acknowledge my committee members Dr. Ahjeong Son, Dr. Jose Vasconcelos and Dr. Mark Dougherty for agreeing to be on my committee and for their valuable time and suggestions in shaping this dissertation into its final form. I would also like to acknowledge Dr. Jacob Dane, Dr. Joel Melville and Dr. Alice Smith for being excellent teachers and making learning exciting. I would also like to thank Dr. Ming-Kuo Lee for serving as an external reader of this dissertation. His comments greatly helped in ironing out the finer details of my dissertation. I would also like to thank the Civil Engineering Department of Auburn University for the Teaching Fellowship that funded my research over the past few years. I would like to thank Dr. Charles E Schaefer of Shaw Environmental Engineering for his support and input in development of bioaugmentation modeling framework. I would also like to thank my friends and colleagues at Auburn and away who helped me grow as a person. I would also like to thank my family for their support and patience. Finally, I would like to thank my wife for her unconditional support and unwavering affection over the years. v Table of Contents Abstract ......................................................................................................................................... ii Acknowledgments ....................................................................................................................... iv List of Tables ................................................................................................................................ x List of Illustrations ...................................................................................................................... xii Chapter 1 INTRODUCTION ....................................................................................................................... 1 1.1 Introduction to groundwater contamination .............................................................. 1 1.2 Numerical modeling of reactive transport problems ................................................. 3 1.3 Parameter estimation using parallel genetic algorithm (GA) .................................... 4 1.4 Research Objectives .................................................................................................. 6 Chapter 2 BENCHMARKING A VISUAL-BASIC BASED MULTI-COMPONENT ONE- DIMENSIONAL REACTIVE TRANSPORT MODELING TOOL ............................... 8 2.1 Review of existing reactive transport models ............................................................ 8 2.2 Model development and numerical solution ............................................................ 12 2.2.1 Transport module ............................................................................................ 16 2.2.1.1 Explicit advection scheme .................................................................. 16 2.2.1.2 Explicit TVD scheme .......................................................................... 17 2.2.1.3 Implicit finite difference method for dispersion scheme .................... 18 2.2.1.4 Fully implicit numerical solution for advection-dispersion ................ 20 vi 2.2.2 Reaction module ............................................................................................. 20 2.2.2.1 Kinetic type problem........................................................................... 20 2.2.2.2 Geochemistry equilibrium type problem ............................................ 21 2.3 Testing the model ..................................................................................................... 22 2.3.1 Pure advection ................................................................................................. 22 2.3.2 Advection dispersion modules ........................................................................ 24 2.3.2.1 Explicit advection and implicit dispersion .......................................... 24 2.3.2.2 Fully implicit advection dispersion scheme ........................................ 25 2.3.2.3 TVD advection and implicit dispersion scheme ................................. 27 2.3.3 Advection dispersion and reaction modules ................................................... 28 2.3.3.1 Explicit advection and implicit dispersion .......................................... 29 2.3.3.2 Fully implicit advection dispersion scheme ........................................ 30 2.3.3.3 TVD advection and implicit dispersion .............................................. 32 2.4 Benchmarking RT1D .............................................................................................. 34 2.4.1 First order sequential degradation ................................................................... 34 2.4.2 Four species coupled sequential first order degradation ................................. 35 2.4.3 Four component decay chain .......................................................................... 38 2.4.4 Modified Monod kinetics for TCE bioaugmentation ...................................... 40 2.4.5 Rate-limited sorption reaction in porous media .............................................. 43 2.4.6 Microbial transport and growth under denitrification conditions ................... 47 2.4.7 Carbon Tetrachloride Biodegradation ............................................................. 50 2.4.8 Geochemical transport involving a constant capacitance model .................... 56 2.4.9 Multiple Sequential Batch Reactor ................................................................. 60 vii 2.4.10 Ion exchange ................................................................................................. 64 2.5 Summary and conclusions ...................................................................................... 65 Chapter 3 MODELING Dehalococcoides Sp. AUGMENTED BIOREMEDIATION IN A SINGLE FRACTURE SYSTEM ................................................................................................... 67 3.1 Introduction ............................................................................................................. 67 3.2 Experimental method and governing equations ....................................................... 70 3.3 Results and discussion ............................................................................................. 75 3.3.1 Testing the batch kinetic model against Schaefer et al. (2009b) model predictions ...................................................................................................... 75 3.3.2 Calibration of the reactive transport model .................................................... 77 3.3.3 Testing the reactive transport model .............................................................. 80 3.3.4 Sensitivity Analysis ........................................................................................ 83 3.3.4.1 Model response to variations in the yield coefficient ......................... 84 3.3.4.2 Model response to variations in the DHC maximum utilization rate constant (q).......................................................................................... 86 3.4 Summary and Conclusions ...................................................................................... 88 Chapter 4 ASSESSMENT OF PARALLEL GENETIC ALGORITHMS FOR CALIBRATING ONE- DIMENSIONAL MULTI-COMPONENT REACTIVE TRANSPORT MODELS ...... 90 4.1 Background .............................................................................................................. 90 4.2 Methodology ? General Steps in Genetic Algorithm ............................................. 95 4.3 Genetic Algorithm Implemented in this Study ........................................................ 96 4.4 Parallelization of the Genetic Algorithm ................................................................. 97 4.5 Details of the Numerical Model used for Fitness Calculation ................................. 99 viii 4.6 Results and Discussion .......................................................................................... 101 4.6.1 Benchmark Problem-1: Parameter estimation in a rate-limited sorption problem ......................................................................................................... 101 4.6.2 Benchmark Problem-2: Parameter estimation in a sequential decay problem103 4.6.3 Benchmark Problem-3: Parameter estimation for a TCE Biodegradation model............................................................................................................. 105 4.6.4 Benchmark Problem-4: Parameter estimation for a Carbon tetrachloride bioremediation problem ................................................................................ 108 4.6.5 Scalability of Parallel GA ............................................................................. 112 4.6.6 Sensitivity to Initial Population and Generations ......................................... 114 4.6.7 Sensitivity to the Bounds of the Unknown Parameter .................................. 116 4.7 Conclusions ............................................................................................................ 118 Chapter 5 SUMMARY AND CONCLUSIONS ...................................................................................... 120 5.1 Summary and Conclusion ...................................................................................... 120 5.2 Recommendations for future work ....................................................................... 121 References ............................................................................................................................... 124 Appendix A1 ........................................................................................................................... 139 A1.1 Understanding the spreadsheet ............................................................................ 139 A1.2 Input spreadsheet label details ............................................................................ 145 A1.2.1 Transport module .................................................................................... 145 A1.2.2 Kinetic reaction module .......................................................................... 146 A1.2.3 Geochemistry equilibrium module.......................................................... 147 A1.2.4 Kinetic reaction parameters .................................................................... 148 A1.2.5 Geochemistry equilibrium parameters (without transport) ..................... 149 ix A1.2.6 Geochemistry equilibrium parameters (with transport) .......................... 151 Appendix A2 ........................................................................................................................... 152 x List of Tables Table 2.1: Parameters for low and high Peclet number simulations .......................................... 29 Table 2.2: Model parameters used in problem 2 obtained from Quezada et al. (2004) .............. 37 Table 2.3: Model parameters used in problem 3 obtained from Bauer et al. (2001) .................. 40 Table 2.4: Model parameters regressed from batch experiments in Schaefer et al. (2009b) ...... 42 Table 2.5: Model parameters used for problem?5 ...................................................................... 46 Table 2.6: Model parameters used for problem?6 ...................................................................... 49 Table 2.7: Model parameters used for problem? 7 ..................................................................... 54 Table 2.8: Stoichiometry of chemical reactions (tableau) for the problem?8 ............................ 59 Table 2.9: Model parameters used for problem?8 ...................................................................... 60 Table 2.10: Model parameters used for problem?9 .................................................................... 62 Table 2.11: Stoichiometry of chemical reactions (tableau) used for problem?9 ........................ 63 Table 2.12: Stoichiometric matrix for the simulation of Test Problem ? 10 .............................. 64 Table 3.1: Summary of numerical model parameters used for Experiment A ........................... 75 xi Table 3.2: Bio-augmentation parameters regressed from batch experiments in Schaefer et al (2009a). ........................................................................................................................... 76 Table 3.3: Calibrated Monod parameters for the Experiment A ................................................ 79 Table 3.4: Summary of numerical model parameters used for Experiment B ............................ 81 Table 4.1: Comparison of the PGA estimated parameters with true solutions along with their bounds and percentage error for benchmark problem -1 .............................................. 102 Table 4.2 Comparison of the PGA estimated parameters with true solutions along with their bounds for benchmark problem ? 2 .............................................................................. 104 Table 4.3: Comparison of the PGA estimated parameters with true solutions along with their bounds for benchmark problem ? 3 .............................................................................. 107 Table 4.4: Model parameters used for benchmark problem ? 4 ............................................... 111 Table 4.5: Comparison of the PGA estimated parameters with true solutions along with their bounds for benchmark problem ? 4 .............................................................................. 112 Table 4.6: Simulation times for all the benchmark problems for different number of processors in seconds .......................................................................................................................... 113 Table 4.7: High and low bounds for the four difference scenarios along with the GA estimated value .............................................................................................................................. 117 xii List of Illustrations Figure 1.1: Schematic diagram of different sources of groundwater (Canada, 2011) ................. 2 Figure 2.1: Illustration of the simulation options available in RT1D ......................................... 14 Figure 2.2: RT1D results for different Courant numbers for the explicit advection scheme with v=1 cm/day, dx=1 and duration of 20 days ................................................................... 23 Figure 2.3: RT1D results for different Courant numbers for the TVD advection scheme with v=1 cm/day, dx=1 and duration of 20 days ............................................................................ 23 Figure 2.4: RT1D results for low Peclet number simulations with varying Courant number using the explicit advection and implicit dispersion scheme (v=1cm/day, D=0.2 cm2/day, dx=0.1cm) ....................................................................................................................... 24 Figure 2.5: RT1D results for high Peclet number simulations with varying Courant number using the Explicit advection and Implicit dispersion scheme (v=1cm/day, D=0.05 cm2/day, dx=0.1cm) ....................................................................................................................... 25 Figure 2.6: RT1D results for low Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=1cm/day, D=0.2 cm2/day, dx=0.1cm)26 Figure 2.7: RT1D results for high Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=1cm/day, D=0.05 cm2/day, dx=0.1cm)??????????????. .............................................................. 26 xiii Figure 2.8: RT1D results for low Peclet number simulations with varying Courant number using the TVD advection and implict dispersion scheme (v=1cm/day, D=0.2 cm2/day, dx=0.1cm) ....................................................................................................................... 27 Figure 2.9: RT1D results for high Peclet number simulations with varying Courant number using the TVD advection and implict scheme (v=1cm/day, D=0.05 cm2/day, dx=0.1cm) ...... 28 Figure 2.10: RT1D results for low Peclet number simulations with varying Courant number using the explicit advection and implict dispersion scheme (v=1 cm/day, D=10 cm2/day, dx=5 cm, k=7.0E-4 day-1, T=3000 days) ....................................................................... 30 Figure 2.11: RT1D results for high Peclet number simulations with varying Courant number using the explicit advection and implict dispersion scheme (v=0.4 cm/day, D=0.08 cm2/day, dx=0.4 cm, k=7.5E-2 day-1, T=50 days) .......................................................... 30 Figure 2.12: RT1D results for low Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=1 cm/day, D=10 cm2/day, dx=5 cm, k=7.0E-4 day-1, T=3000 days) ................................................................................. 31 Figure 2.13: RT1D results for high Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=0.4 cm/day, D=0.08 cm2/day, dx=0.4 cm, k=7.5E-2 day-1, T=50 days) ......................................................................... 32 Figure 2.14: RT1D results for low Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=1 cm/day, D=10 cm2/day, dx=5 cm, k=7.0E-4 day-1, T=3000 days) ................................................................................. 33 Figure 2.15: RT1D results for high Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=0.4 cm/day, D=0.08 cm2/day, dx=0.4 cm, k=7.5E-2 day-1, T=50 days) ......................................................................... 33 Figure 2.16: RT1D simulation results for problem 1using two different k values ..................... 35 Figure 2.17: Comparison of RT1D simulation results with analytical solutions for problem 2 . 38 xiv Figure 2.18: Comparison of RT1D simulation results with analytical solutions for problem 3 . 39 Figure 2.19: Comparison of RT1D simulation results with the published model results for problem 4 ........................................................................................................................ 42 Figure 2.20: Comparison of RT1D simulations with analytical solutions for problem?5: (a) aqueous concentration for different mass transfer coefficients; (b) solid phase concentrations for ?=0.00015; (c) solid phase concentrations for ?=0.015; and (d) solid phase concentrations for ?=2 .......................................................................................... 45 Figure 2.21: Comparison of the RT1D results with the analytical solutions for problem?5 with a decay rate constant value of 0.03 day-1 ........................................................................... 46 Figure 2.22: Comparison of the RT1D model results (solid line) with the published model results (dotted line) and published data (dots) for benchmark problem?2: (a) effluent nitrate for low substrate conditions; (b) effluent biomass for low substrate conditions; (c) effluent nitrate for high substrate conditions; and (d) effluent biomass for high substrate conditions ........................................................................................................................ 50 Figure 2.23: Comparison between the published model data (dots) and the RT1D simulations (line) for benchmark problem?3: (a) carbon tetrachloride after 4 days; (b) acetate after 13 days; (c) nitrate after 4 days; (d) mobile KC after 7 days; (e) mobile KC after 11 days; (f) mobile KC after 14 days; (g) mobile KC after 4 days; and (h) mobile KC after 13 days. (Note: KC is the strain of the mobile bacteria) .......................................................... 55-56 Figure 2.24: Comparison of the RT1D results (solid lines) with the published model results (dots) for the benchmark problem?4: (a) simulation results for Case-I and (b) simulation results for Case-V ........................................................................................................... 58 Figure 2.25: Comparison of the RT1D results with the published model results for problem-9..63 Figure 2.26: Results for example problem 10 a) Breakthrough profile for Na+ b) Breakthrough profile for Mg2+ c) Breakthrough profile for Ca2+ .......................................................... 65 xv Figure 3.1: Batch model results: Comparison of proposed model results against the model results reported in Schaefer et al. (2009b) .................................................................................. 76 Figure 3.2: Model calibration results. Continuous lines represent the model results and the symbols represent Experimental-A data ......................................................................... 80 Figure 3.3: Model validation results. Continuous lines represent the model results and the symbols represent Experimental-B data ......................................................................... 82 Figure 3.4: Comparison of modeled and observed PCE concentration levels ............................ 83 Figure 3.5: Model response to an decrease in the yield value (by 0.5) for Experiment-A ......... 85 Figure 3.6: Model response to an increase in the yield value (by 5) for Experiment-A ............. 86 Figure 3.7: Model response a decrease in the q value (by 0.5) for Experiment-A ..................... 87 Figure 3.8: Model response to an increase in the q value (by 2) for Experiment-A ................... 88 Figure 4.1: Illustration showing the flow of a parallel genetic algorithm .................................. 99 Figure 4.2: Comparison of results from PGA estimated parameters and the true parameters for benchmark problem ? 1 .................................................................................................. 103 Figure 4.3: Comparison of results from PGA estimated parameters and the true parameters for benchmark problem ? 2 ................................................................................................ 105 Figure 4.4: Comparison of results from PGA estimated parameters and the true parameters for benchmark problem ? 3 ................................................................................................ 108 Figure 4.5: Comparison of results from PGA estimated parameters and the true parameters for benchmark problem ? 4 ................................................................................................ 112 xvi Figure 4.6: Speedup for all benchmark problems for different number of processors ............. 114 Figure 4.7: Fitness of the best solution for each generation by varying initial population for each benchmark problem ........................................................................................................ 115 Figure 4.8: Comparison of different scenarios of the high and low bounds with the true solution for the benchmark problem ? 1 ..................................................................................... 117 Figure A1.1: Spreadsheet layout for RT1D parameters to generate reaction-specific input template ......................................................................................................................... 140 Figure A1.2: Example input template for a kinetic-type reactive transport ............................. 142 Figure A1.3: Example input template for a geochemistry equilibrium coupled with transport.143 1 Chapter 1 INTRODUCTION 1.1 Introduction to groundwater contamination Groundwater is an important source of water supply used by populations all around the world (Todd, 1980). In the United States, groundwater accounts for approximately twenty one percent of the annual water supply budget and hence it is considered as an important natural resource (Perlman, 2011). Figure 1 shows the schematic diagram of different sources of groundwater. Inadvertent discharge of harmful contaminants including metals and organics into groundwater aquifers poses a significant threat to this resource. Groundwater systems could be contaminated by leachates emanating from several anthropogenic sources including landfills, mines, leaking underground storage tanks (LUSTs), and other industrial waste sites. Also, the extensive use of various forms of chlorinated solvents for dry cleaning and metal degreasing activities has resulted in widespread contamination of groundwater and soil systems (Coleman et al., 2002). In addition to heavy metal and chlorinated solvent issues, contamination of aquifers by petroleum products released from LUSTs has been reported at several field sites in different continents including North America, Europe, and Australia (Lu et al., 1999; Moreau, 1987; Prommer et al., 1998). In the US alone, about 10 to 20% of the estimated total of 2 million underground storage tanks are expected to be leaking (Atlas and Cerniglia, 1995). Therefore, contamination of groundwater by LUSTs and other sources poses a significant threat to groundwater quality. 2 Figure 1.1 Schematic diagram of different sources of groundwater (Canada, 2011) Once contaminated, groundwater aquifers require some type of remediation strategy to restore its water quality to safe drinking water levels. The type of remediation system used and the associated costs would depend on the type of contaminant, extent of contamination and the local geological conditions. Typically, most contaminated field sites are first remediated using some form of conventional pump-and-treat systems. However, due to solubility limitations, conventional pump-and-treat methods have been ineffective at several sites. Therefore, in recent years, engineers have attempted to use innovative in-situ technologies such as the bioremediation methods to transform the contaminants into non-toxic daughter products (Beeman and Bleckmann, 2002; Clement et al., 2004). Active and passive (or natural attenuation) bioremediation methods can be used to treat petroleum and chlorinated solvent plumes (Clement, 2011). Design and application of these bioremediation methods require tools that can model the fate and transport of the contaminants and the associated site-specific biogeochemical reactions. 3 The modeling step is particularly important when natural attenuation methods are employed for managing petroleum and/or chlorinated solvent plumes (Clement et al., 2000; Lu et al., 1999; Prommer et al., 1998; Rolle et al., 2008). 1.2 Numerical modeling of reactive transport problems Reactive transport problems can be broadly classified into either geochemical equilibrium or kinetic problems based on the nature of the chemical processes involved in the remediation methods employed. Chemical kinetics describes the rate of reaction in a fast chemical reaction whereas geochemical equilibrium describes the speciation at equilibrium for a slow chemical reaction. Currently, there are several models capable of simulating multi-species multi- dimensional reactive transport processes in groundwater for both kinetic and geochemical equilibrium problems. MT3DMS developed by (Zheng and Wang, 1999) is capable of simulating three-dimensional advective-dispersive multi-species transport processes. Multispecies bioreactive transport in one-dimensional soil columns has been numerically modeled by several researchers (Clement et al., 2004; Clement et al., 1996; Clement et al., 1997; Schaefer et al., 2009b; Yu and Semprini, 2004; Zysset et al., 1994). Bioplume III (Rifai et al., 1998) is a two- dimensional, finite difference model for simulating both aerobic and anaerobic biodegradation of hydrocarbons in groundwater in addition to advection, dispersion, sorption and ion exchange. RT3D (Clement et al., 1998) combines a multispecies sequential dechlorination, biodegradation and first-order decay processes in groundwater in three-dimensional domain. In addition to these models, there are several models capable of simulating geochemistry equilibrium problems. MINEQL (Westall, 1976) reduces the chemical equilibrium problem into a set of nonlinear equations that can be solved by Newton-Raphson iteration scheme. MICROQL-I (Westall, 1979a) computes chemical equilibrium in aqueous systems without sorbed or solid phase 4 species. MICROQL-II (Westall, 1979b) includes additional routines for solving adsorption equilibria in aqueous systems using constant capacitance, diffused layer, stern layer and triple layer adsorption models (Kumar, 2006). However, most of these models are three-dimensional models developed for field scale evaluations with added complexity for multi-dimensional parameters. This level of complexity is not required for simple laboratory experiments. The column experiments done in laboratories can be simplified as one-dimensional problems. Furthermore, there are very few models that are capable of solving both geochemistry and kinetic problems in the same model. Therefore, it is necessary to have a simple, user-friendly, accessible, comprehensive yet robust modeling tool that is capable of simulating a variety of chemical, biochemical and geochemical processes. Therefore, one of the goals of this study is to present a comprehensive one-dimensional modeling tool that using Visual Basic for Applications (VBA) in Excel that is capable for simulating a variety of bio-geochemical problems that can be used by laboratory researchers. 1.3 Parameter estimation using parallel genetic algorithm (GA) The multi-component advection-dispersion reaction equation explains the fate and transport of a contaminant in groundwater. The reaction part of the equation involves several kinetic parameters that are unique to the contaminant or remediation process. However, many of these reaction parameters used in the study of these experiments are theoretical and could be difficult to estimate using experimental studies. They may require several laboratory experiments isolating each chemical compound to estimate these parameters (Massoudieh et al., 2008; Prommer et al., 1998; Schaefer et al., 2009a; Schaefer et al., 2009b; Singh et al., 2008). The methods currently used for estimating the reaction parameters are trial and error-based methods 5 or use software like CXTFIT (Toride et al., 1995). However, the complexity of the problem increases as the number of parameters to be estimated increases. Modern optimization techniques such as genetic algorithm (GA) can be used to estimate these unknown parameters. GA follows the evolutionary concept of survival of the fittest for mathematical optimization using genetic recombination (Holland, 1975). GA searches through a solution space until it converge to a global minima or maxima. GA has been used in groundwater hydrology and hydrogeology for the parameter estimation with successful results (Babbar and Minsker, 2006; Singh et al., 2005; Sinha and Minsker, 2007; Wang, 1997). However, they are computationally intensive during the fitness calculation and they have to be run for long periods of times to find the global minima. This computational expense can be a significant drawback when compared to other parameter estimation techniques. In order to make the GAs computationally efficient, parallel computing techniques can be utilized. The concurrency in fitness calculations makes GA an ideal candidate for parallel computing (Cant?-Paz, 1998). Most loops in GA can be computed independently and this makes it an embarrassingly parallel problem. This means that there could be little to no communication between the processors. Based on the architecture, the type of parallel computers can be categorized into either distributed memory or shared memory computers. The new computers available in the market are equipped with multi-core processors and a shared memory parallel computing language could be used to optimize existing serial algorithms. Currently there are several studies involving parallel genetic algorithms for parameter estimation in groundwater problems in both shared memory and distributed memory architectures (Babbar and Minsker, 2006; He et al., 2007; Sinha and Minsker, 2007). However, there are no parallel implementations in OpenMP FORTRAN in a desktop environment. Therefore, in this effort, we explored the use 6 of a parallel genetic algorithm (PGA) for estimating the kinetic parameters in different types of reactive transport problems. 1.4 Research Objectives The aim of this dissertation is to develop a set of comprehensive tools for laboratory researchers for modeling multi-component reactive transport problems and also provide a tool to estimate the model parameters from laboratory experiments. The first objective of this dissertation study is to develop a comprehensive one- dimensional multi-component reactive transport tool that is capable of simulating both kinetic and geochemical equilibrium problems. The model is also able to run both kinetic and geochemical equilibrium problems in a batch mode as well as couple with transport problems. This tool is developed using Visual Basic for Application (VBA) in Microsoft Excel without any additional software to be installed. This model was validated using a variety of kinetic and geochemistry problems published in the literature. The second objective of this study is to apply this model to a laboratory experiment involving the bioaugmentation of chlorinated ethenes using Dehalococcoides Sp. in a single fracture system. We have developed a mathematical framework to simulate the bioaugmentation of PCE-DNAPL and estimated the parameters using a trial and error process for a low flow experimental data. These parameters were validated using the experimental data for the high flow experimental data. The third objective of this study is to develop a GA based parameter estimation tool. The GAs are computationally intensive search procedures. Therefore, the GA was optimized to run in parallel on a multicore desktop computer using the shared memory parallel computing language 7 OpenMP using a parallel genetic algorithm (PGA). This PGA was used to estimate the model parameters for four problems whose true parameters are already known. This dissertation consists of five chapters. The first chapter (the current chapter) provides a basic introduction to this research, lists the objectives, and provides a brief summary of each chapter. The second chapter focuses on the first objective, with the outcome of this effort already published in Computers and Geosciences (Torlapati and Clement, 2012b). The third chapter focuses on the second objective and the outcome of this effort has been published in Groundwater Monitoring and Remediation (Torlapati et al., 2012). The fourth chapter focuses on the third objective and the outcome of this effort were communicated to Environmental Modelling and Software (Torlapati and Clement, 2012a). The fifth chapter provides a brief summary of the key outcomes of this research and offers recommendations for future work. 8 Chapter 2 DEVELOPMENT OF A VISUAL-BASIC BASED MULTI-COMPONENT ONE- DIMENSIONAL REACTIVE TRANSPORT MODELING TOOL 2.1 Review of existing reactive transport models Laboratory-scale experiments are also routinely used to develop a better understanding of various biogeochemical transport processes expected to occur under field conditions. Both batch and column studies have been employed for establishing the feasibility of proposed remediation methods (Schaefer et al., 2009b; Schaefer et al., 2010b). Numerical models are also routinely utilized at this feasibility assessment stage to interpret the laboratory data and to develop a better understanding of underlying treatment processes (Clement et al., 1998; Phanikumar et al., 2002; Torlapati et al., 2012). The modeling exercises can greatly help the scaling and design steps that are required for deploying field-scale remediation technologies. Currently, there are several models available that are capable of simulating multi-component, multi-dimensional reactive transport processes. Zheng and Wang (1999) developed MT3DMS, which is capable of simulating three-dimensional advective-dispersive multi-component transport processes. Clement et al. (1998) developed a reactive transport code RT3D, which is based on MT3DMS, that can simulate bioreactive transport scenarios involving kinetic reactions. Prommer et al. (1998) combined MT3DMS with PHREEQC (Parkhurst and Appelo, 1999) to simulate both equilibrium and kinetic reactions. UTCHEM (de Blanc et al., 1996) and 9 BIOPLUME-II (Rafai et al., 1987) are two generic reactive transport models that are capable of simulating bioreactive transport processes. BIOPLUME-II uses a modified version of the USGS code MOC and can simulate aerobic biodegradation of petroleum plumes. A later version, known as BIOPLUME-III (Rafai et al., 1998), can simulate both aerobic and anaerobic reactions involved in petroleum biodegradation. BIOCHLOR (Aziz, 2000) is an EXCEL-based tool which implements a sequential decay analytical solution described in Sun and Clement (1999), Sun et al. (1999), and Clement (2001) to simulate natural attenuation processes occurring at chlorinated solvent contaminated sites. However, BIOCHLOR is an analytical model and is limited by the capabilities of the underlying solution procedure; some of these limitations are explained in Quezada et al. (2004), Srinivasan and Clement (2008) and Srinivasan et al. (2007). Most biological processes that degrade organic contaminants such as hydrocarbons and chlorinated solvents are kinetic-limited reactions, and they are described using a set of ordinary differential equations (ODEs). Within a numerical reactive transport formulation, these ODEs are normally implemented as a reaction package using the operator-split strategy and are independently solved by an ODE solver (Clement et al., 1998). The geochemical processes (which mediate the fate and transport of inorganic contaminants such as metals) on the other hand are mostly equilibrium-controlled reactions that require solution to a set of non-linear algebraic equations. These non-linear equilibrium equations can be solved by an independent geochemistry routine and they can also be integrated into a transport formulation using the operator-split strategy (Cederberg et al., 1985; Prommer et al., 1998). There are several computer codes available in the literature that can be used to solve chemical speciation problems. WATEQ was one of the first geochemical models that uses an iteration scheme to solve the system of non-linear equations; the code, however, cannot explicitly handle heterogeneous 10 reactions such as precipitation and dissolution processes (Truesdell and Jones, 1974). MINEQL is a commonly used model that employs the well-known tableau approach to define and solve the chemical equilibrium problem. MINEQL is capable of handling reactions such as mineral precipitation and dissolution (Westall, 1976). MICROQL-I is a simplified version of the MINEQL code which can solve chemical equilibrium without sorbed or solid phase species (Westall, 1979a). MICROQL-II is an updated version that includes routines for modeling adsorption equilibrium using constant capacitance, diffuse layer, and triple layer adsorption models (Westall, 1979b). The chemical speciation code MINTEQA2 is the most comprehensive software (also developed based on the ideas espoused in MINEQL and MICROQL codes) for modeling dilute aqueous solutions (Allison et al., 1990). PHREEQC is another widely used geochemical code that can perform a variety of geochemical speciation calculations (Parkhurst and Appelo, 1999). PHAST is a three-dimensional reactive-transport model derived from coupling the geochemical model PHREEQC with the solute-transport model HST3D (Kipp, 1987). The flow and transport model is restricted to constant density, constant temperature, and saturated ground-water flow conditions. Chemical reactions considered in HST3D include mineral and gas equilibrium, ion exchange, surface complexation, solid solutions, and kinetic reactions (Parkhurst, 2004). Most of the multi-dimensional, numerical transport codes discussed above require considerable experience and expertise to model coupled reactive transport problems. Therefore, in the published literature, laboratory researchers have developed several simpler one- dimensional reactive transport tools to model column-scale datasets. For example, Miller and Benson (1983) developed a numerical model CHEMTRN to simulate the transport of solutes in saturated porous media. The model simulated advection, dispersion, ion exchange, and formation 11 of aqueous-phase complexes. Engesgaard and Kipp (1992) developed a one-dimensional geochemical code to simulate precipitation-dissolution and oxidation-reduction reactions and used it to model pyrite oxidation processes at a field site in Denmark. Zysset et al. (1994) presented a numerical model for describing reactive transport processes occurring within a biofilm. (Clement et al., 1996) developed a one-dimensional model to simulate bioremediation patterns occurring near an injection well. Clement et al. (1997) developed a one-dimensional code to simulate bacterial transport and denitrification processes observed in a column experiment. Prommer et al. (1998) developed a one-dimensional numerical model for predicting biodegradation occurring at a petroleum hydrocarbons site. Islam and Singhal (2002) presented a one-dimensional multi-component reactive transport model coupled with geochemical equilibrium reactions to simulate the interactions between the microbial redox reactions and inorganic geochemical reactions. (Amos et al., 2009) developed a numerical model to study the enhanced dissolution of PCE-DNAPL in presence of dechlorination bacterial cultures. Clement et al. (2004) developed a code for modeling DNAPL-dissolution and rate-limited sorption occurring in a biologically reactive one-dimensional porous media system. Schaefer et al. (2009b) and Torlapati et al. (2012) developed one-dimensional models which were used to simulate laboratory studies that explored the effects of bioaugmentation on chlorinated solvent contaminants. Unfortunately, most of the one-dimensional tools discussed above, which are primarily developed for solving a specific research problem, have little or no documentation. Moreover, none of these codes are user friendly tools that can be used by other laboratory researchers. Also, these tools support either a kinetic formulation or an equilibrium formulation; none of these models provide a flexible framework, such as an EXCEL interface, which would allow 12 users to fine tune the code to their specific needs. The objective of this effort is to develop a comprehensive, one-dimensional reactive transport model within a user-friendly, EXCEL-based Visual Basic environment. Our goal is to provide a unified EXCEL tool that can be easily adapted by others to model laboratory-scale experiments involving different types of biological (kinetic) and geochemical (equilibrium) reactions. In this paper, we present the details of this tool and demonstrate its use by solving five benchmark problems. The benchmark problems illustrate the characteristics of a variety of bio-geochemical problems that would be of interest to a broad range of environmental scientists. The problems are described in detail to provide a comprehensive benchmarking database which can be used for testing other reactive transport codes. 2.2 Model development and numerical solution Reactive transport problems in porous media systems could be mediated by either kinetic or equilibrium processes. Kinetic models are used to describe relatively slow chemical reactions, whereas equilibrium reactions are used to describe fast chemical reactions. Several reaction transport codes are available in the literature, but they can only handle kinetic type reactions (e.g., RT3D) or can handle equilibrium type reactions (Cederberg et al., 1985; Parkhurst and Appelo, 1999). RT1D was designed to provide a unified platform for simulating transport problems involving both geochemical and kinetic reactions. However, it is important to note that the current version can either simulate a set of pure kinetic reactions or a set of pure geochemical reactions, but one cannot mix both types of reactions. The capabilities of RT1D model, which are designated as ?simulation options,? are summarized in Figure 2.1. The model currently supports four different ?simulation options? 13 involving the use of two types of reaction modules: a kinetic module and an equilibrium module which can simulate kinetic and equilibrium reactions, respectively. There are four types of simulation options available within RT1D. The first simulation option can be used to solve batch kinetic problems. The second option can be used to solve one-dimensional reactive transport problems involving kinetic reactions. The third option can solve batch equilibrium problems. The fourth option can solve one-dimensional reactive transport problems involving geochemical equilibrium reactions. If the kinetic module is selected, the user should also provide a problem- dependent reaction package. The kinetic module supports several standard reaction models that are already coded within a set of pre-programmed reaction packages. In addition to these preprogrammed packages, RT1D also supports a user-defined reaction package via which the user can input any type of kinetics. The geochemical module, on the other hand, does not require a problem-dependent reaction package, and instead uses a MICROQL-based chemistry package to solve the equilibrium problem. The information required for formulating a specific geochemical problem are input using the standard tableau. The multi-component one-dimensional reactive transport model developed in this study, designated as RT1D, solves a coupled set of advection-dispersion-reaction equations for a total of ?n? components. The model simulates the transport of ?m? mobile components that are either fully or partially coupled to a set of ?n-m? immobile components. The reactions between these components could be mediated by biological/geochemical kinetic reactions, or geochemical equilibrium reactions. The governing set of equations solved by the model can be written in a general form: 14 2i i i i 2i i iC C C ?VD=- + +t R x R x R where i = 1, 2,3? m (2.1) j jS=?t where j = (m+1), (m+2), (m+3),? (n) (2.2) Figure 2.1 Illustration of the simulation options available in RT1D where V is the velocity (m/day); D is the hydrodynamic dispersion coefficient (m2/day), Ci is the aqueous phase concentration (mg/L) of a mobile component i; Sj is the solid phase concentration (mg/mg) of an immobile component j; m is the total number of mobile 15 components; n is the total number of components (note that a total of (n-m) immobile components are numbered sequentially after numbering all the mobile components); Ri is the linear retardation factor of the ith mobile component [ d?KR1 , where ? is the bulk density (mg/L); ? is the porosity; Kd is the linear sorption constant (L/mg)]; and ?i and ?j are the reactions involving mobile and immobile components, respectively. The expressions used for ?i and ?j terms would vary depending on the type of reactions involved in the system. Note the immobile component equations do not have advection dispersion terms, but will have reaction terms that will be coupled to some of the mobile component reaction terms. Also, the mobile- component reaction terms themselves could be coupled to each other. The coupled set of reactive transport equations, represented by equations (1) and (2) are solved using the operator split approach (Clement et al., 1998; Torlapati and Clement, 2012b). Using this approach, the governing set of transport equations can be written as: ii i CCV=-t R x (2.3) 2ii 2iCCD=t R x (2.4) ii i dC ?=dt R (2.5) j jdS=?dt (2.6) In the numerical code, the advection terms in all the mobile components (equation 2.3) are first solved using an advection solver module. Next, the advected concentrations ( ) are C 16 used to solve for dispersion terms (equation 2.4) in the mobile component transport equations using a dispersion solver module. Finally, the dispersed concentrations ( ) are used to solve a set of coupled reaction terms involving both mobile and immobile components, represented by equations (2.5) and (2.6). For transport problems involving kinetic reactions, the reaction part of the transport equations (equations 2.5 & 2.6) would yield a set of coupled ODEs. These ODEs, referred as the kinetic reaction package, are solved using an ODE solver. For transport problems involving geochemical equilibrium reactions, the reaction terms would yield a set of coupled non-linear equations. These non-linear equilibrium equations are represented using the tableau approach (Westall, 1976) and are solved using the geochemical equilibrium solver, MICROQL, developed by (Westall, 1979a, b). 2.2.1 Transport module The advection module provides two explicit solver options: a total variation diminishing (TVD) solver and an explicit finite difference solver that uses backward difference approximation. The advected concentrations are then used to solve the dispersion equation, within the dispersion module, using the implicit finite difference method. In addition to these explicit-implicit solvers, there is also a fully-implicit option that solves the advection-dispersion terms together using a fully-implicit approach. 2.2.1.1 Explicit advection scheme The advection part of the transport equation can be solved using the explicit backward difference approximation as shown below: C 17 n + 1 n n ni i i i- 1 i C - C C - CV=-? t R ? x (2.7) where n+1iC is the concentration of the component at the current time step at the current node; niC is the concentration of the mobile component at the previous time step at the current node and ni-1C is the concentration of the mobile component at the previous time step at the preceding node. After further simplification, we can solve for the concentration of mobile component at the current node ( n+1iC ) as shown below: n + 1 n n ni i i- 1 i i CrC C C CR (2.8) where V?tCr ?x is known as the grid Courant number. 2.2.1.2 Explicit TVD scheme Numerical dispersion is a major concern while solving the advection dominated problems. RT1D includes a robust total variation diminishing (TVD) scheme that minimizes numerical dispersion errors. Details of this scheme are given below. Using the Taylor series expansion, the standard Lax-Wendroff (LW) scheme for the advection term can be written as (Leveque, 2002): 2n + 1 n n n n n n i i i + 1 i - 1 i + 1 i i - 1C r C rC = C - ( C - C ) + ( C - 2 C + C )22 (2.9) The above equation can be rearranged and written in a flux-balance format as: nnn + 1 n i+ 1 /2 i-1 /2ii F -FCC =- t ?x (2.10) 18 Where n n n ni+ 1 /2 i i+ 1 iVF = V C + ( 1 - C r ) ( C - C )2 and n n n ni- 1 /2 i- 1 i i- 1VF = V C + ( 1 - C r ) ( C - C )2 respectively. Note the flux terms defined above consist of a lower and higher order flux terms. The lower order flux term in ni+1/2F is iVpC and the higher order term in ni+1/2F is i+ 1 iV (1 -Cr )(p C -p C )2 . Furthermore in TVD schemes, a flux limiter will be used to minimize the potential numerical oscillations induced by the higher order term as shown below: n n n ni+ 1 /2 i i+ 1 iVF = V C + ( 1 - C r ) ( C - C ) * ?2 (2.11) Different types of flux limiters are available in the literature and in this study we have used the Van-leer flux limiter given as (Leveque, 2002): ?+??=1+? (2.12) where nni i-1 nni-1 i-2C -C?=C -C Farthing and Miller (2001) investigated the adaptive-stencil and finite volume schemes to capture sharp fronts and shocks in advective-dispersive transport. They observed that the Lax- Wendroff scheme performed better in the presence of a flux limiter. 2.2.1.3 Implicit finite difference method for dispersion scheme The dispersion part of the transport equation can be numerically discretized using a central difference approximation. 19 n + 1 n n + 1 n + 1 n + 1i i i- 1 i i+ 1 i 2C C C 2 C CR = Dtx (2.13) Where Ci-1 is the concentration at the previous node for the current time step, Ci+1 is the concentration at the next node for the current time step. The above equation can be further simplified as follows n n + 1 n + 1 n + 1i i- 1 i i+ 1 i i i ? 2 ? ?C C C 1 CR R R (2.14) Where 2Dt?=x . Assembling equation (14) on a node-by-node basis would yield a following tri- diagonal matrix of the form: n+1 01 n+1 22 n+1 33 n+1 44 n+1 55 n+1 nn C1 0 0 0 0 0 0 C da b c 0 0 0 0 C d0 a b c 0 0 0 C = d0 0 a b c 0 0 C d0 0 0 a b c 0 C .... .. .. .. .. .. . .. d0 0 0 0 0 0 1 C (2.15) where i ?a=R , i 2?b= 1R , i ?c=R and nid=-C ; Co is the concentration at the boundary node. The above matrix can be solved using a tridiagonal matrix solver to solve for all the unknown concentrations at the new time level. 20 2.2.1.4 Fully implicit numerical solution for advection-dispersion In this option, we solve the advection and dispersion together implicitly. We use a central difference approximation for the advection term and the numerically discretized form for the advection-dispersion equation is as follows: n + 1 n n + 1 n + 1 n + 1 n + 1 n + 1i i i+ 1 i- 1 i- 1 i i+ 1 2ii C - C C - C C - 2 C + CVD= - +? t R 2 ? x R ? x (2.16) The above equation can be further simplified as follows: n n + 1 n + 1 n + 1i i i - 1 i i i + 1C ( R . ) = ( ? + 1 )C ( 2 R . )C (1 ? )C (2.17) where 2?x= D?t and V?x?=2D . Expanding the equation (2.17) for all the nodes would yield a tridiagonal matrix similar to (2.15). For this problem, the values of a, b, c and d are given as?+1 , i2-R. , 1-? and nii-C (R .?) , respectively. 2.2.2 Reaction module 2.2.2.1 Kinetic type problem For a multi-species, the reaction part shown in equation (2.5) & (2.6) simplifies into a set of ordinary differential equations (ODE). These ODEs are referred to as a reaction package. The set of ODEs described within a kinetic reaction package can be solved using two different ODE solvers: a standard 4th order Runge-Kutta (RK) solver, or a more robust Runge-Kutta-Fehlberg 21 (RKF) solver (Chapra and Canale, 1998). The RK solver uses a constant reaction time step, whereas the RKF solver will automatically subdivide the reaction time step into sub steps to minimize the local error. 2.2.2.2 Geochemistry equilibrium type problem The geochemical equilibrium reaction problem is formulated in the form of a tableau that represents the interactions between all the components and species involved in the chemical system. As defined by Westall (1976), a species is a chemical entity of interest present in the system whereas a component is a basic building block used for forming various chemical species in the system. The stoichiometric relationship between the components and the species can be represented in the form of a matrix known as the tableau. The chemical speciation problem, defined by the tableau, is solved using an EXCEL-VBA version of MICROQL code. The details of the numerical solution schemes employed by MICROQL are discussed in Westall (1979b). Within RT1D, the transport equations that involve geochemical (or equilibrium) reactions are solved using an approach proposed by Cederberg et al. (1985), which is slightly different from the approach used for the solving the transport equations involving kinetic reactions. As discussed in Cederberg et al. (1985), first the aqueous concentration of a component of interest will be transported using the transport module. The aqueous concentration of a component at a particular node is calculated by subtracting the concentrations of all the sorbed species associated with that component from its total concentration. After the transport step, the updated (advection-dispersed) aqueous component concentration is added back to the sorbed concentrations of the respective component to compute the total component concentration at that node. This total component concentration is then transferred to MICROQL to solve the geochemical speciation problem. The equilibrated species concentrations are used to update the 22 values of aqueous component concentrations for the next time transport step. Further details of this transport algorithm are discussed in Cederberg et al. (1985). 2.3 Testing the model The explicit finite difference schemes used in the advection modules of RT1D program are constrained by certain stability condition criteria. To ensure that the results generated by the numerical models are stable and error-free, researchers have used Courant number (Vdx/(R.dt)) and Peclet number (Vdx/D) to determine the grid size and time step for the simulations. In the following sections, we have tested the stability of our numerical models for different Courant and Peclet numbers. For explicit schemes, the Courant number should be less than 1 to obtain oscillation free results. If the Courant number exceeds 1, numerical oscillations are observed near the advective front. In addition, to obtain good quality solution the Peclet number should be set below 2. 2.3.1 Pure advection In this section, we tested the Explicit and TVD schemes for different Courant numbers by setting the value of dispersion to 0. Simulations were performed for a one-dimensional column of 50 cm length. The pore velocity was about 1 cm/day and the simulations were performed for duration of 20 days. The grid size was set to 1 cm and the time step was varied to generate Courant numbers of 1, 0.5, 0.1 and 0.01. A constant boundary condition of 1 mg/L was supplied at the inlet for the complete duration of the experiment. The results of the simulations for the Explicit and TVD schemes are shown in Figure 2.2 and 2.3 respectively. It can be observed from these figures that at Courant number 1, the results from both the advection schemes produce sharp advective fronts. However, the numerical dispersion comes into effect as the Courant number is decreased. 23 The numerical dispersion is comparatively less in case of TVD schemes than the explicit schemes. Figure 2.2 RT1D results for different Courant numbers for the explicit advection scheme with v=1 cm/day, dx=1 and a duration of 20 days Figure 2.3 RT1D results for different Courant numbers for the TVD advection scheme with v=1 cm/day, dx=1 and a duration of 20 days 24 2.3.2 Advection dispersion modules The program was tested for a high and low Peclet numbers of 0.5 and 2 respectively for varying Courant numbers 0.01, 0.1, 0.5 and 1. A hydrodynamic dispersion coefficient of 0.05 cm2/day was used to perform high Peclet number simulations and a hydrodynamic dispersion coefficient of 0.2 cm2/day was used to perform low Peclet number simulations. A grid size of 0.1 cm and a pore velocity of 1 cm/day were used to perform the simulations. The results from these simulations were compared against the analytical solutions presented in van Genuchten and Alves (1982). 2.3.2.1 Explicit advection and implicit dispersion Simulations were performed with a time step of 0.1, 0.05, 0.01 and 0.001 for both high and low Peclet numbers to generate Courant numbers of 1, 0.5, 0.1 and 0.01 respectively. Figure 2.4 and 2.5 show the results for low and high Peclet number simulations respectively. Figure 2.4 RT1D results for low Peclet number simulations with varying Courant number using the explicit advection and implicit dispersion scheme (v=1cm/day, D=0.2 cm2/day, dx=0.1cm) 25 Figure 2.5 RT1D results for high Peclet number simulations with varying Courant number using the Explicit advection and Implicit dispersion scheme (v=1cm/day, D=0.05 cm2/day, dx=0.1cm) It was observed from the figures that the explicit scheme showed numerical dispersion in the presence of low Peclet number and this numerical dispersion decreased with increase in the Peclet number and the simulation results were closer to the analytical solutions. 2.3.2.2 Fully implicit advection dispersion scheme Simulations were performed with a time step of 0.1, 0.05, 0.01 and 0.001 for both high and low Peclet numbers using the fully implicit advection dispersion scheme to generate Courant numbers of 1, 0.5, 0.1 and 0.01 respectively. Figure 2.6 and 2.7 show the results for low and high Peclet number simulations respectively. 26 Figure 2.6 RT1D results for low Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=1cm/day, D=0.2 cm2/day, dx=0.1cm) It was observed from the simulations that the fully implicit scheme performed a lot better than the explicit scheme. However there was some numerical dispersion when the Peclet number was high. Figure 2.7 RT1D results for high Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=1cm/day, D=0.05 cm2/day, dx=0.1cm) 27 2.3.2.3 TVD advection and implicit dispersion scheme Simulations were performed with a time step of 0.1, 0.05, 0.01 and 0.001 for both high and low Peclet numbers using the TVD advection and implicit dispersion scheme to generate Courant numbers of 1, 0.5, 0.1 and 0.01 respectively. Figure 2.8 and 2.9 show the results for low and high Peclet number simulations respectively. It was observed from the results that the TVD scheme was consistent with the analytical results for all the Courant and Peclet numbers. It is highly recommended that the users of the RT1D use TVD scheme for accurate results. Figure 2.8 RT1D results for low Peclet number simulations with varying Courant number using the TVD advection and implict dispersion scheme (v=1cm/day, D=0.2 cm2/day, dx=0.1cm) 28 Figure 2.9 RT1D results for high Peclet number simulations with varying Courant number using the TVD advection and implict scheme (v=1cm/day, D=0.05 cm2/day, dx=0.1cm) 2.3.3 Advection dispersion and reaction modules In this section, we present the results to test the program?s stability under varying Peclet and Courant numbers in the presence of a first order decay constant for a single component. The simulations were compared with analytical solutions. For the high Peclet number simulations, the data from the first component of Bauer et al. (2001) was used and for low Peclet number simulations, the date from the first component of Quezada et al. (2004) was used. Simulations were performed for different Courant numbers of 0.01, 0.1, 0.5 and 1.0. The parameters for these simulations are presented in Table 2.1. 29 Table 2.1: Parameters for low and high Peclet number simulations Parameter High Low Length (cm) 40 3000 Time (days) 50 3000 dx (cm) 0.4 5 Velocity (cm/day) 0.4 1 Dispersion Coefficient (cm2/day) 0.08 10 Retardation 1 5.3 Decay constant (1/day) 0.075 7.00E-04 2.3.3.1 Explicit advection and implicit dispersion Simulations were performed with a time step (dt) of 1, 0.5, 0.1 and 0.01 for high Peclet number simulations and 26.5, 13.25, 2.65 and 0.265 days for low Peclet number simulations to generate Courant numbers of 1, 0.5, 0.1 and 0.01 respectively. Figure 2.10 and 2.11 show the results for low and high Peclet number simulations for the explicit advection and implicit dispersion respectively. It was observed from the figures that explicit scheme performed well with both high and low Peclet numbers. There was some numerical dispersion in the presence of high Peclet number. However, this is negligible. 30 Figure 2.10 RT1D results for low Peclet number simulations with varying Courant number using the explicit advection and implict dispersion scheme (v=1 cm/day, D=10 cm2/day, dx=5 cm, k=7.0E-4 day-1, T=3000 days) Figure 2.11 RT1D results for high Peclet number simulations with varying Courant number using the explicit advection and implict dispersion scheme (v=0.4 cm/day, D=0.08 cm2/day, dx=0.4 cm, k=7.5E-2 day-1, T=50 days) 31 2.3.3.2 Fully implicit advection dispersion scheme Simulations were performed with a time step of 1, 0.5, 0.1 and 0.01 for high Peclet number simulations and 26.5, 13.25, 2.65 and 0.265 days for low Peclet number simulations to generate Courant numbers of 1, 0.5, 0.1 and 0.01 respectively. Figure 2.12 and 2.13 show the results for low and high Peclet number simulations respectively. It was observed from the figures that explicit scheme performed well with both high and low Peclet numbers. There was some numerical dispersion in the presence of high Peclet number. However, this is negligible. Figure 2.12 RT1D results for low Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=1 cm/day, D=10 cm2/day, dx=5 cm, k=7.0E-4 day-1, T=3000 days) 32 Figure 2.13 RT1D results for high Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=0.4 cm/day, D=0.08 cm2/day, dx=0.4 cm, k=7.5E-2 day-1, T=50 days) 2.3.3.3 TVD advection and implicit dispersion Simulations were performed with a time step of 1, 0.5, 0.1 and 0.01 for high Peclet number simulations and 26.5, 13.25, 2.65 and 0.265 days for low Peclet number simulations to generate Courant numbers of 1, 0.5, 0.1 and 0.01 respectively. Figure 2.14 and 2.15 show the results for low and high Peclet number simulations respectively. 33 Figure 2.14 RT1D results for low Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=1 cm/day, D=10 cm2/day, dx=5 cm, k=7.0E-4 day-1, T=3000 days) Figure 2.15 RT1D results for high Peclet number simulations with varying Courant number using the fully implicit advection dispersion scheme (v=0.4 cm/day, D=0.08 cm2/day, dx=0.4 cm, k=7.5E-2 day-1, T=50 days) 34 There was some numerical dispersion for very low courant numbers in the high Peclet number simulations but this was negligible and the results from RT1D matched the analytical solutions well. 2.4 Benchmarking RT1D 2.4.1 First order sequential degradation A generalized reactive transport equation for a single component first order decay is present below. 2 2C C CR = -V D k Ct x x (2.18) where C is concentration of the mobile component, k is the first order decay constant (T-1). The above equation simulates a tracer when the decay constant is set to zero. This reaction package could be written in the code form as follows: dydt(1)=1/R(1)*RC(1)*Conc(1) Where dydt(1) is the ordinary differential equation for the mobile component 1, R(1) is the retardation factor for component 1 and RC(1) is the user-set reaction parameter (k), Conc(1) is the aqueous concentration of the mobile component. In order to test this reaction package, we considered a test column of 40 m length with contaminant transporting through the column for 50 days with a pore velocity of 0.4 m/day and 0.08 m2/day of hydrodynamic dispersion coefficient. We used a first order decay constant of 0.075 day-1 and a tracer simulation was also done with the decay constant set to zero and the results were compared against the analytical solutions presented in Quezada et al. (2004). The results from these simulations are presented in 35 Figure 2.16. It can be observed from the figure that the simulation results from the RT1D model are able to match the results from the analytical solutions of both the tracer as well as the decay constant. Figure 2.16 RT1D simulation results for problem-1 using two different k values 2.4.2 Four species coupled sequential first order degradation Quezada et al. (2004) presented the analytical solutions for a four species coupled sequential first order degradation reactions. The reaction equations simulate the transport and coupled decay of four mobile components. The governing equations are as follows: 21 1 1 1 1 12C C CR = - V + D - k Ct x x (20) 22 2 2 2 c 2 / c 1 c 2 / c 1 1 1 2 2 c 2 / c 3 c 2 / c 3 3 32C C CR = - V + D + F Y k C - k C + F Y k Ct x x (21) 36 23 3 3 3 c 3 / c 1 c 3 / c 1 1 1 c 3 / c 2 c 3 / c 2 2 2 3 32C C CR = - V + D + F Y k C + F Y k C - k Ct x x (22) 24 4 4 4 c 4 / c 2 c 4 / c 2 2 2 c 4 / c 3 c 4 / c 3 3 3 4 42C C CR = - V + D + F Y k C + F Y k C - k Ct x x (23) where Ri is the retardation factor; V the seepage velocity (LT-1), D is the hydrodynamic dispersion coefficient (LT-2), Y is the yield coefficient, F is the fraction, k (T-1) is the first order decay constant. The yield coefficient determines the number of moles of a component converted into its daughter product. For example, the term Yc4/c2 determines the number of moles of component 2 converting into component 4. The yield coefficient is set to 1 for all components and this means there is a complete conversion of the component at each time step. The fraction governs the amount of total degraded component converting from one component to another. For example, Fc4/c2 determines the fraction of total degraded component 2 converting into component 2. This dependence on fractions causes the coupling effects between different components. Also, the reaction equations themselves contain the terms for other components causing coupling effects. The above reaction equations could be written in code format as follows: dydt(1) = -RC(1) * Conc(1) / R(1) dydt(2) = (RC(11) * RC(5) * RC(1) * Conc(1) - RC(2) * Conc(2) + RC(12) * RC(6) * RC(3) * Conc(3)) / R(2) dydt(3) = (RC(13) * RC(7) * RC(1) * Conc(1) + RC(14) * RC(8) * RC(2) * Conc(2) - RC(3) * Conc(3)) / R(3) dydt(4) = (RC(14) * RC(9) * RC(2) * Conc(2) + RC(16) * RC(10) * RC(3) * Conc(3) - RC(4) * Conc(4)) / R(4) 37 where dydt(i) is the expression for ODE for component i; i could be either 1, 2, 3 or 4; R(i) is the retardation factor for the component i; RC(i) is the user-set reaction parameters, Conc(i) is the concentration of the component i. The model parameters for this problem are presented in Table 2.2 and the results from this simulation are presented in Figure 2.17. It can be observed from the figure that the results from the RT1D simulations were able to match the analytical solutions well. Table 2.2: Model parameters used in Test Problem 2 obtained from Quezada et al. (2004) Column Length, L (m) 40 Dispersivity (m) 0.2 Velocity (m d-1) 0.4 Singularity parameter (?) 0.1 R1 1 R2 2 R3 3 R4 4 k1 (days-1) 0.075 k2 (days-1) 0.05 k3 (days-1) 0.02 k4 (days-1) 0.045 Yield(all) 1 Fc2/c1 0.75 Fc3/c1 0.25 Fc3/c2 0.5 Fc4/c2 0.5 Fc2/c3 0.9 Fc4/c3 0.1 Boundary for Species 1 (mol l-1) 1.0 Boundary for Species 2-4 (mol l-1) 0 Total Time (d) 50 38 Figure 2.17 Comparison of RT1D simulation results with analytical solutions for problem-2 2.4.3 Four component decay chain Bauer et al. (2001) presented analytical solutions for the transport of a decay chain for in homogenous porous media. The reaction equations presented in the paper were used as a built-in reaction package for the RT1D. The reaction equations are as follows 21 1 1 1 1 12C C CR = - V + D - k Ct x x (2.23) 22 2 2 2 1 1 1 2 2 22C C CR = - V + D + k R C - k R Ct x x (2.24) 23 3 3 3 2 2 2 3 3 32C C CR = - V + D + k R C - k R Ct x x (2.25) 24 4 4 4 3 3 3 4 4 42C C CR = - V + D k R C - k R Ct x x (2.26) 39 where Ri is the retardation factor; V the seepage velocity (LT-1), D is the hydrodynamic dispersion coefficient (LT-2), ki (T-1) is the first order decay constant. Further details of the model are available in Bauer et al. (2001). The code format for this reaction package is as follows: dydt(1) = -RC(1) * Conc(1) dydt(2) = (RC(1) * Conc(1) * R(1) - RC(2) * Conc(2) * R(2)) / R(2) dydt(3) = (RC(2) * Conc(2) * R(2) - RC(3) * Conc(3) * R(3)) / R(3) dydt(4) = (RC(3) * Conc(3) * R(3) - RC(4) * Conc(4) * R(4)) / R(4) The model parameters for this problem are presented in Table 2.3 and the results from this simulation are available in Figure 2.18. It can be observed from the figure that the RT1D simulations were able to match the analytical solutions well. Figure 2.18 Comparison of RT1D simulation results with analytical solutions for problem 3 40 Table 2.3: Model parameters used in problem 3 obtained from Bauer et al (2002) Column Length, L (m) 3000 Dispersivity (m) 10 Velocity (m d-1) 1 R1 5.30 R2 1.90 R3 1.20 R4 1.30 K1 (days-1) 7.5E-4 K2 (days-1) 5.0E-4 K3 (days-1) 4.5E-4 K4 (days-1) 3.8E-4 Boundary for Species 1 (mg l-1) 100 Boundary for Species 2-4 (mg l-1) 0 Total Time (d) 3000 2.4.4 Modified Monod kinetics for TCE bioaugmentation Schaefer et al. (2009b) conducted batch experiments to study the bioaugmentation of TCE. They used Modified monod kinetics to model these reactions. The reaction package for this model is as follows: (note that there are no advection and dispersion terms as this is a batch system) TC E TC E TC E TC E TC E TC E d C q X C1=-d t R C + K (2.27) D CE D CE D CE T CE T CE D CE T CE T CE T CET CE D CE D CE T CE d C q X C q X C11=- d t R R C + KCC + K 1 + I (2.28) 41 V C V C V C D C E D C E V C D C ET C E D C E T C E V C V C D C E D C E T C E D C E T C E d C q X C q X C11= - + d t R RC C CC + K 1 + + C + K 1 + I I I (2.29) T CE T CE D CE D CE V C V C T CE T CE T CE D CE V CT CE T CE D CE D CE D CE V C V C T CE T CE D CE q C q C q Cd X 1 1 1= Y X + + d t R C + K R RC C CC + K 1 + C + K 1 + + I I I (2.30) where Ci (mM) and X (cells/L) are the concentration of ith compound and biomass respectively; i can be either TCE, DCE and VC; qi is the maximum biomass utilization rate, Ki is the half velocity coefficient of the compound, I is the competition coefficient, Ri is the retardation due to the presence of air-gap. For further details about the model and the experimental methods, refer to Schaefer et al. (2009b). The code form of this reaction package is given below: mTCE = ((RC(4) * Conc(1)) / (Conc(1) + RC(1))) mDCE = ((RC(5) * Conc(2)) / (Conc(2) + RC(2) * (1 + (Conc(1) / RC(7))))) mVC = ((RC(6) * Conc(3)) / (Conc(3) + RC(3) * (1 + (Conc(1) / RC(7)) + (Conc(2) / RC(8))))) dydt(1) = (-1 / R(1)) * Conc(5) * mTCE dydt(2) = -Conc(5) * (mDCE / R(2) - mTCE / R(1)) dydt(3) = -Conc(5) * (mVC / R(3) - mDCE / R(2)) dydt(4) = Conc(5) * mVC / R(3) dydt(5) = (RC(10) * Conc(5)) * (mTCE / R(1) + mDCE / R(2) + mVC / R(3)) To simplify the equations, we have defined three different variables (mTCE, mDCE and mVC) to define the Monod terms for each component. The monod parameters used in this 42 simulations are presented in Table 2.4 and the results from the RT1D simulations and it?s comparison with the model simulations from Schaefer et al. (2009b) is shown in Figure 2.19. It can be observed from the figure that the RT1D simulations were able to predict the concentration trends exactly. Table 2.4: Model parameters regressed from batch experiments in Schaefer et al (2009) Component Initial Condition (mM) Boundary Condition (mg/L) Kd (L/Kg) R K (mM) q (mmol L-1 (cell h)-1) I (mM) DCE 0.123788 0.103157 0.070 1.340 2.00E-03 7.00E-13 5.20E-03 VC 0 0 0.016 1.078 1.40E-02 1.40E-12 1.00E+06 Ethene 0 0 - 1.000 - - - DHC (Immobile) 0 0 - - - - - DHC (mobile) 0 0 - 1.000 - - - Figure 2.19 Comparison of RT1D simulation results with the published model results for problem-4 43 2.4.5 Rate-limited sorption reaction in porous media In case of non-equilibrium conditions, the sorption is controlled by a mass-transfer coefficient (?). When the ? value is really low, the plume acts like tracer because there is no sorption and when the ? value is really high, the plume acts like a retarded plume due to the linear sorption. This kind of rate-limited sorption kinetics can be modeled using the following reaction package: 2 2 d C C C SVD ? C - k Ct x x K (2.31) d dS ?? SC-dt ?K (2.32) where C is the concentration of the aqueous phase component (mg/L); S is the concentration of the solid phase component (mg/mg); ? is the bulk density (mg/L); ? is the porosity; Kd is the linear sorption constant (L/mg); k is the first order decay constant (day-1); and ? is the mass transfer coefficient (day-1). Clement et al. (1998) used a similar type of formulation to model rate-limited reactions, although their study ignored the first order decay term. This reactive transport problem involves two components: a mobile component that represents the aqueous phase concentration (C), and an immobile component that represents the solid phase concentration (S). Using the operator split strategy, the reaction kinetics for this problem can be formulated as: d d C S? C - kCd t K (2.33) 44 d dS ?? SC-dt ?K (2.34) Equations (2.33) and (2.34) are coded into a reaction package. The column was assumed to be initially clean and the left hand boundary condition was fixed at 1 mg/L. The porosity of the column was assumed to be 0.3, the bulk density of the porous media (?) was set to 1600 g/L, and the sorption constant (Kd) was set at 1.875 x 10-4 L/g. The model was run using three different mass transfer coefficients (?): 0.00015, 0.015, and 2 (day-1). Other model parameters used in this benchmark problem are summarized in Table 2.5. Note when the value of mass transfer coefficient is low, the solute is expected to behave like a tracer with R=1; on the other hand, when the mass transfer coefficient is high, the solute is expected to behave like a retarded plume with R=2. The scenarios in-between these two extreme conditions would result in rate limited, non-equilibrium transport conditions. Toride et al. (1993) presented a set of analytical solutions for transport equations involving non-equilibrium sorption and first-order decay terms. Valocchi and Werth (2004) developed a web-based Java applet to implement the analytical solutions developed by Toride et al. (1993). This Java applet was used to benchmark the results of the RT1D code. The definition of model parameters used in the analytical solution vary slightly from the model definitions described above; in order to compare the results, the mass transfer coefficient to be used in the analytical solutions must be calculated using the formula d ??=?K where ?` is the mass transfer coefficient used in the analytical solution. Figure 2.20a shows the aqueous phase concentrations simulated by RT1D for different values of mass transfer coefficients; the figure also shows the analytical solution results. Similar results for solid phase concentrations are shown in Figures 2.20b, 2.20c, and 2.20d. Simulations 45 were also completed using a constant decay co-efficient (k) of 0.03 day-1 and Figure 2.21 compares the numerical results with analytical results. The total computer time required for solving this benchmark problem was about 18 seconds. The figures show that the results from the RT1D simulations were able to match the analytical results. Furthermore, as expected, Figure 2.20a shows that the aqueous phase concentration profile was retarded by a factor of R=2, when the mass transfer coefficient was set to an arbitrarily high value. Figure 2.20 Comparison of RT1D simulations with analytical solutions for problem?5: (a) aqueous concentration for different mass transfer coefficients; (b) solid phase concentrations for ?=0.00015; (c) solid phase concentrations for ?=0.015; and (d) solid phase concentrations for ?=2 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 0 5 10 15 20 25 30 35 40 A q u e o u s c o n c e n tr a ti o n ( m g /L ) D i s t a n c e ( c m ) ? = 2 ( R T 1 D ) ? = 0 .0 1 5 ( R T 1 D ) ? = 0 .0 0 0 1 5 ( R T 1 D ) ? = 2 ( A n a l y t i c a l ) ? = 0 .0 1 5 ( A n a l y t i c a l ) ? = 0 .0 0 0 1 5 ( A n a l y t i c a l ) ( a ) 0 . 0 0 E + 0 0 2 . 0 0 E - 07 4 . 0 0 E - 07 6 . 0 0 E - 07 8 . 0 0 E - 07 1 . 0 0 E - 06 1 . 2 0 E - 06 1 . 4 0 E - 06 1 . 6 0 E - 06 0 10 20 30 40S o li d C o n c e n tr a ti o n ( m g /m g ) D i s t a n c e ( c m ) ? = 0 .0 0 0 1 5 (A n a l y t i c a l ) ? = 0 .0 0 0 1 5 (R T 1 D ) ( b ) 0 . 0 0 E + 0 0 2 . 0 0 E - 05 4 . 0 0 E - 05 6 . 0 0 E - 05 8 . 0 0 E - 05 1 . 0 0 E - 04 1 . 2 0 E - 04 0 10 20 30 40S o li d C o n c e n tr a ti o n ( m g /m g ) D i s t a n c e ( c m ) ? = 0 .0 1 5 ( A n a l y t i c a l ) ? = 0 .0 1 5 ( R T 1 D ) ( c ) 0 . 0 0 E + 0 0 4 . 0 0 E - 05 8 . 0 0 E - 05 1 . 2 0 E - 04 1 . 6 - 04 2 . 0 0 E - 04 0 10 20 30 40S o li d C o n c e n tr a ti o n ( m g /m g ) D i s t a n c e ( c m ) ? = 2 ( A n a l y t i c a l ) ? = 2 ( R T 1 D ) ( d ) 46 Table 2.5: Model parameters used for problem?5 Length (cm) 40 Total time (days) 50 ?x 0.4 ?t 0.01 Pore Velocity (cm/day) 0.53 Longitudinal Dispersion Coefficient (cm2/day) 0.08 # Mobile Species 1 # Immobile Species 1 Figure 2.21 Comparison of the RT1D results with the analytical solutions for problem?5 with a decay rate constant value of 0.03 day-1 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 0 10 20 30 40 C o n c e n t r a t i o n ( m g / L ) L e n g t h ( c m ) ? = 2 ( R T 1 D ) ? = 2 ( A n a l y t i c a l ) ? = 0 .0 0 0 1 5 ( A n a l y t i c a l ) ? = 0 .0 0 0 1 5 ( R T 1 D ) 47 2.4.6 Microbial transport and growth under denitrification conditions Clement et al. (1997) studied the effects of denitrifying conditions on the growth and transport of bacteria in a porous media column under two substrate loading conditions. A numerical model was developed to generate the breakthrough profiles of bacterial cells and substrates. A first order attachment and detachment model was used to describe the exchange processes between mobile and immobile-phase bacterial cells. This benchmark problem considered three mobile components namely nitrate, acetate and aqueous-phase bacteria, and one immobile component namely immobile bacteria. The reaction package used in the problem is given below: N N sNad C r X ?= -r X -d t n (2.35) AsA Aa rX ?dC = -r X -d t n (2.36) a d e sX a a t ad X K X ?r X K Xd t n (2.37) s a t aX s d e sd X n K Xr X K Xdt ? (2.38) where CN, CA, Xa and Xs are concentrations (mg/L) of nitrate, acetate, aqueous-phase bacteria and immobile-phase bacteria (mg/mg), respectively. The parameters Kat (day-1) and Kde (day-1) are the attachment and detachment coefficients of mobile and immobile phase bacteria, respectively; n is the porosity of soil; and ? is the bulk density of the soil (mg/L). The rate expression is the nitrate utilization rate described using Monod kinetics as: Nr 48 N AN m a x N N A A C Crq K C K C,where q max is the maximum nitrate utilization rate (mg nitrate/ mg biomass-day), KN is the half saturation coefficient for nitrate (mg/L); KA is the half saturation coefficient for acetate (mg/L). The specific utilization rate of acetate ( ) and biomass growth rate ( ) are given by the expressions: and X X/N N dr Y r K , where YA/N and YX/N are the yield coefficients for acetate and biomass, respectively, and Kd is the cell decay rate coefficient (day-1). A finite difference grid of size 1 cm and a time step of 0.001 day were used in this problem. Other model parameters used are summarized in Table 2.6. Further details of the experiments are available in Clement et al. (1997). The total amount of computer time required for solving this benchmark problem was about 590 seconds. Figure 2.22a-d compare RT1D simulation results with published model results and data available in the literature. Figures 2.22a and 2.22c compare effluent concentrations of nitrate at different times and Figures 2.22b and 2.22d show mobile phase bacteria concentrations in the effluent. The results show that the RT1D model simulations closely matched published data. Ar Xr A A/N Nr Y r 49 Table 2.6: Model parameters used for problem?6 Pore Velocity (cm/day) 1890.91 Length (cm) 50 Longitudinal dispersion coefficient (cm2/day) (D) 1890.91 Porosity (n) 0.44 Bulk density (?) (mg/l) 1.56E6 Time (days) 15 Microbial decay rate (day-1) (Kd) 0.06 Attachment coefficient (day-1) (Kat) 288 Detachment coefficient (day-1) (Kde) 0.32 Distribution coefficient (L/mg) 3.9E-7 Half saturation coefficient: (mg/L) Acetate (KA) Nitrate (KN) 1.20 0.66 Maximum specific nitrate utilization rate (mg nitrate/ mg biomass-day) (qmax) 7.21 Yield: Acetate (mg acetate/mg NO3) (YA/N) Biomass (mg biomass/mg NO3) (Yx/N) 0.84 0.13 Initial condition (mg/L): Acetate (CA) Nitrate (CN) Mobile bacteria (XM) Immobile bacteria (mg/mg) (XIM) 0 0 1.0E-15 3.0E-07 Boundary condition (Low Substrate) (mg/L): Acetate (CA) Nitrate (CN) Mobile bacteria (XM) 5.0 5.5 0 Boundary condition (High Substrate) (mg/L): Acetate (CA) Nitrate (CN) Mobile bacteria (XM) 48.0 58.0 0 50 Figure 2.22 Comparison of the RT1D model results (solid line) with the published model results (dotted line) and published data (dots) for benchmark problem?2: (a) effluent nitrate for low substrate conditions; (b) effluent biomass for low substrate conditions; (c) effluent nitrate for high substrate conditions; and (d) effluent biomass for high substrate conditions 2.4.7 Carbon Tetrachloride Biodegradation Phanikumar et al. (2002) developed a bioremediation model to predict carbon tetrachloride (CT) degradation processes observed in sequential column experiment. In this study, we have used one of their experiments, identified as once-fed (OF) column, as a benchmark problem. The laboratory experiment used a 200 cm long column fitted with an 11-cm long slug injection zone at a distance of 34 cm away from the column inlet. The injection zone was fitted with an inlet and an outlet to circulate flow within this zone. This injection-extraction setup was used to 0 1 2 3 4 5 6 0 5 10 15 20 25 E ff lu e n t N O 3 (m g /L ) T i m e ( d a y s ) N i t r a t e ( R T 1 D ) N i t r a t e ( M o d e l ) N i t r a t e ( D a t a ) (a) 0 . 0 0 0 1 0 . 0 0 1 0 . 0 1 0 . 1 1 0 5 10 15 20 25 E ff lu e n t B io m a s s C o n c e n tr a ti o n (m g /L ) T i m e ( d a y s ) A q u e o u s B i o m a s s ( R T 1 D ) A q u e o u s B i o m a s s ( M o d e l ) A q u e o u s B i o m a s s ( D a t a ) (b ) 0 10 20 30 40 50 60 0 5 10 15 20 25 E ff lu e n t N O 3 (m g /L ) T i me ( d a y s ) N i t r a t e ( R T 1 D ) N i t r a t e ( M o d e l ) N i t r a t e ( D a t a ) (c ) 0 . 0 0 1 0 . 0 1 0 . 1 1 10 0 5 10 15 20 25 E ff lu e n t B io m a s s C o n c e n tr a ti o n (m g /L ) T i m e ( d a y s ) A q u e o u s B i o m a s s ( R T 1 D ) A q u e o u s B i o m a s s ( M o d e l ) A q u e o u s B i o m a s s ( D a t a ) (d ) 51 inoculate the column with nutrients and mobile bacteria for about 16 minutes by circulating flow at rate of 20 ml/min. The inoculation step was completed just once at the beginning of the experiment. It was assumed that the inoculation step completely replaced the initial contents of the slug injection zone and hence the concentrations in the inoculant solution were used as the initial conditions for the 11 cm zone. Table 3 summarizes the details of the boundary and initial conditions used in this problem for the entire column. The transport problem considered four mobile components: carbon tetrachloride, acetate, nitrate, and mobile-phase bacteria; and two immobile components: sorbed carbon tetrachloride and immobile-phase bacteria. The reaction package used for modeling this bioremediation problem, as provided in Phanikumar et al. (2002) is : 'd C T C T M I M d C T C T?fK dC ??1 - k C ( X +X ) - 1- f K C S? dt ? (2.39) a m a x a na M I M a dC ? M MR - ( X X )d t Y (2.40) m a x a n K Cnn M I M a n M I M n n b ? M M bdCR ( X + X ) - ( 1 M ) ? M (X + X )d t Y Y (2.41) M m a x a n K C a a t M d e a I MdX ? M M b (1 M ) K X K (1 M )Xdt (2.42) CT d C T C TdS ? 1-f K C Sdt (2.43) 52 IM m a x a n K C a d e a I M a t MdX ? M M b (1 M ) K (1 M ) X K Xdt (2.44) where f is the fraction of equilibrium sites, bKC is the microbial decay rate (day-1), Kat is the attachment coefficient (day-1), Kde is the detachment coefficient (day-1), k` is the CT reaction rate (day-1), ? is the nitrate reaction rate (day-1), ? is the kinetic desorption rate (day-1), ?max is the maximum specific growth rate (day-1), Ya, Yn and Ynb are the yield rates of acetate, nitrate and biomass respectively; CCT, Ca, Cn and SCT are the aqueous concentrations of carbon tetrachloride, acetate, nitrate and the sorbed concentration of carbon tetrachloride, respectively; XM and XIM are the concentrations of mobile and immobile bacteria, respectively. Also, Ma and Mn are the Monod terms for acetate and nitrate reactions, respectively, and given by the expressions: aa sa a CM KC and nn sn n CM KC where K sa and Ksn are the half saturation coefficients of acetate and nitrate utilization reactions, respectively. The kinetic equations (2.39) to (2.41) describe biodegradation of carbon tetrachloride, utilization of an electron donor (acetate), and an electron acceptor (nitrate). Equation (2.42) describes the growth, decay, and attachment of the mobile phase bacteria, equation (2.43) describes the sorption of carbon tetrachloride using a two- site sorption model, and equation (2.44) describes the growth, decay, and detachment of immobile-phase bacteria. The grid size used was 1 cm and time step was 0.001 day. Other model parameters are summarized in Table 2.7. Figure 2.23 compares RT1D simulation results with the published model results. Figure 2.23a shows the biodegradation patterns of carbon tetrachloride within the column after 4 days. We present the published model data as well as the experimental in Figure 2.23a because the RT1D simulated concentrations were able to match the concentrations from the experimental 53 data much better than the published model data. The total amount of computer time required for simulating this benchmark problem was 28 seconds. As expected, the model results show increased biodegradation activity near the slug injection zone which was inoculated with active bacterial cells. It can be observed from the figures that the results from the RT1D simulations match well with the published model results. 54 Table 2.7: Model parameters used for problem? 7 Parameter Value Pore Velocity (cm/day) 10 Length (cm) 200 Longitudinal dispersion coefficient (cm2/day) (D) 2 Porosity (?) 0.35 Bulk density (?) (mg/L) 1.63E6 Time (days) 4 Microbial decay rate (day-1)(bKC) 0.221 Fraction of equilibrium sites (f) 0.437 Attachment coefficient (day-1) (Kat) 0.9 Detachment coefficient (day-1) (Kde) 0.043 Distribution coefficient (Kd) (L/mg) 3.9E-7 Half saturation coefficient: (mg/L) Acetate (Ksa) Nitrate (Ksn) 1.0 12.0 CT reaction rate (day-1) (k`) 0.189 Nitrate utilization coefficient (day-1) (?) 5.730 Kinetic desorption rate (day-1) (?) 0.36 Maximum specific growth rate (day-1) (?max) 3.11 Yield: Acetate (Ya) Nitrate (Yn) Biomass (Ynb) 0.4 0.25 0.46 Initial condition (ppm): Carbon tetrachloride (CCT) Acetate (Ca) Nitrate (Cn) Mobile bacteria (XM) Immobile bacteria (XIM) Sorbed CT (mg/mg) (SCT) 0.130 0 42 0 0 2.8E-8 Boundary condition (ppm): Carbon tetrachloride (CCT) Acetate (Ca) Nitrate (Cn) Mobile bacteria (XM) Immobile bacteria (XIM) 0.130 0 42 0 0 Slug injection zone inoculation (ppm): Carbon tetrachloride (CCT) Acetate (Ca) Nitrate (Cn) Mobile bacteria (XM) Immobile bacteria (XIM) 0.1 1650 42 11.8 0 55 0 . 0 0 0 . 0 2 0 . 0 4 0 . 0 6 0 . 0 8 0 . 1 0 0 . 1 2 0 . 1 4 0 . 1 6 0 . 1 8 0 50 100 150 200 C o n c e n tr a ti o n ( p p m ) D i s t a n c e ( c m ) C T ( R T 1 D ) C T ( M o d e l ) C T ( D a t a ) (a) 0 300 600 900 0 50 100 150 200 C o n c e n tr a ti o n ( p p m ) D i s t a n c e ( c m ) A c e t a t e ( M o d e l ) A c e t a t e ( R T 1 D ) (b ) 0 5 10 15 20 25 30 35 40 45 0 50 100 150 200 C o n c e n tr a ti o n ( p p m ) D i s t a n c e ( c m ) N i t r a t e ( M o d e l ) N i t r a t e ( R T 1 D ) (c ) 1 . 0 0 E + 0 3 1 . 0 0 E + 0 4 1 . 0 0 E + 0 5 1 . 0 0 E + 0 6 1 . 0 0 E + 0 7 0 50 100 150 200 250 C o n c e n tr a ti o n ( C F U /m L ) D i s t a n c e ( c m) M o b i l e K C ( R T 1 D ) M o b i l e K C ( m o d e l ) (d ) 1 . 0 E + 0 3 1 . 0 E + 0 4 1 . 0 E + 0 5 1 . 0 E + 0 6 1 . 0 E + 0 7 0 50 100 150 200 250 C o n c e n tr a ti o n ( C F U /m L ) D i s t a n c e ( c m) M o b i l e K C ( m o d e l ) M o b i l e K C ( R T 1 D ) (e ) 1 . 0 0 E + 0 3 1 . 0 0 E + 0 4 1 . 0 0 E + 0 5 1 . 0 0 E + 0 6 1 . 0 0 E + 0 7 0 50 100 150 200 C o n c e n tr a ti o n ( C F U /m L ) D i s t a n c e ( c m) M o b i l e K C ( m o d e l ) M o b i l e K C ( R T 1 D ) (f) 56 Figure 2.23 Comparison between the published model data (dots) and the RT1D simulations (line) for benchmark problem?3: (a) carbon tetrachloride after 4 days; (b) acetate after 13 days; (c) nitrate after 4 days; (d) mobile KC after 7 days; (e) mobile KC after 11 days; (f) mobile KC after 14 days; (g) mobile KC after 4 days; and (h) mobile KC after 13 days (Note: KC is the strain of the mobile bacteria) 2.4.8 Geochemical transport involving a constant capacitance model Cederberg et al. (1985) developed a research code, TRANQL, to simulate geochemical multi- component transport in a saturated groundwater system. The TRANQL code was used to study cadmium transport in the presence of chloride and bromide ions. The one-dimensional reactive transport model considered advection, dispersion, surface complexation of cadmium ion, and sorption of free cadmium to solids in the column. They used the finite-element approach to solve the governing transport equations. The geochemical problem was defined using the tableau nomenclature similar to the method presented by Westall (1979a). The chemical equilibrium problem considered a total of 6 components and 14 species. Table 2.8 provides the reaction tableau for the problem and the log K values for all the chemical reactions. Using the RT1D code, we tracked the concentrations of the following three mobile components: cadmium, bromide, and chloride. The remaining three components in this problem including hydrogen ion 57 (H+), surface hydroxyl group (SOH) and electrostatic potential (Psi) were not tracked in the transport module for the following reasons: pH and total SOH values are fixed in this problem and hence remained constant throughout the simulation. The last component, electrostatic potential (Psi), is a hypothetical component which is only used within MICROQL calculations. In order to compute the aqueous cadmium component concentration, we subtracted the concentration of sorbed cadmium species SOCd+ (species 14 in the tableau) from the total cadmium component concentration. Similar calculations can also be made for other components of interest, as described in Cederberg et al. (1985). The grid size used in this problem was 0.03 cm and the time step was 0.06 hours. Cederberg et al. (1985) solved a total of six cases with different initial and boundary concentration levels. These six cases were divided into two groups of three cases based on the total initial and boundary concentrations of cadmium, chloride and bromide ions. Further details about each of these cases are available in Cederberg et al. (1985). In this benchmark exercise, we solved two cases namely, Case-1 and Case-5, described in Cederberg et al.?s study. These two cases were chosen because Case-1 is a base case scenario where the initial and boundary conditions chloride and bromide ion concentrations remained the same at the background levels. Case-5, on the other hand, shows the system?s response when bromide and chloride ion concentration levels were allowed to be higher than the background concentration levels. The transport parameters the initial and boundary conditions used for these two cases are summarized in Table 2.9. The amount of computer time required for solving this benchmark problem was about 30 seconds. The results predicted by the RT1D are compared against TRANQL model results, reported in Cederberg et al., in Figures 2.24a and 2.24b. Figure 6a shows that the aqueous-phase and sorbed-phase cadmium profiles predicted by RT1D are in good agreement with TRANQL 58 results. Since concentrations of chloride and bromide remained constant during Case-1, they are not presented in Figure 2.24a. However, when the concentration of chloride was increased (Case-5 problem), it interacted with sorbed-phase cadmium species and this, resulted in reduced chloride ion levels; these results are shown in Figure 2.24b. Overall, the results from RT1D simulations are in excellent agreement with published results. Figure 2.24 Comparison of the RT1D results (solid lines) with the published model results (dots) for the benchmark problem?4: (a) simulation results for Case-I and (b) simulation results for Case-V 1. 00E -06 1. 00E-05 1. 00E-04 1. 00E-03 1. 00E-02 1. 00E-01 0 2 4 6 8 10 Co ncentr at ion (M ) D i stan ce, cm A q u eo u s C d (R T 1D ) So r b ed C d (R T 1 D ) A q u eo u s C d (d ata) So r b ed C d (d ata) (a) 1. 00E -06 1. 00E-05 1. 00E-04 1. 00E-03 1. 00E-02 1. 00E-01 0 2 4 6 8 10 Co ncentr at ion (M ) D i stan ce, cm A q u eo u s C d (R T 1D ) C h l o r i d e ( R T 1D ) So r b ed C d (R T 1 D ) C h l o r i d e ( d at a) So r b ed C d (d ata) A q u eo u s C d (d ata) (b ) 59 Table 2.8: Stoichiometry of chemical reactions (tableau) for the problem?8 Cl- Br- Cd+2 SOH Psi H+ log K 1 H[+] 0 0 0 0 0 1 0 2 Cd+2 0 0 1 0 0 0 0 3 Cl- 1 0 0 0 0 0 0 4 Br- 0 1 0 0 0 0 0 5 CdCl+ 1 0 1 0 0 0 1.8 6 CdCl2 2 0 1 0 0 0 2.6 7 CdBr+ 0 1 1 0 0 0 2.2 8 CdBr2 0 2 1 0 0 0 3 9 CdOH+ 0 0 1 0 0 -1 -12.69 10 OH- 0 0 0 0 0 -1 -13.91 11 SOH 0 0 0 1 0 0 0 12 SOH2 0 0 0 1 1 1 7.4 13 SO- 0 0 0 1 -1 -1 -9.24 14 SOCd+ 0 0 1 1 -1 -1 -7 XOH 60 Table 2.9: Model parameters used for problem?8 Parameter Value Pore Velocity (cm/hr) 0.33 Length (cm) 10 Longitudinal dispersion coefficient (cm2/hr) 0.0067 Porosity 0.3 Bulk density (g/l) 2500 Time (hrs) 15 Total no. of sites (mol/l) 0.046 Ionic Strength (mol/l) 0.1 pH (constant) 7 Capacitance (F/m2) 1.06 Solid surface area (m2/g) 1 Boundary condition (Case-I): Cd2+ (M) Cl- (M) Br- (M) 1.0E-4 3.0E-4 1.0E-4 Initial condition (Case-I): Cd2+(M) Cl- (M) Br- (M) SOH (M) H+ (M) 1.0E-5 3.0E-4 1.0E-4 4.6E-2 1.0E-7 Boundary condition (Case-V): Cd2+ (M) Cl- (M) Br- (M) 1.0E-4 3.0E-2 1.0E-2 Initial condition (Case-V): Cd2+(M) Cl- (M) Br- (M) SOH (M) H+ (M) 1.0E-4 3.0E-3 1.0E-3 4.6E-2 1.0E-7 2.4.9 Multiple Sequential Batch Reactor Jeppu et al. (2012) proposed a sequential equilibration reactor (SER) system to investigate transport problems involving geochemical equilibrium reactions. They studied adsorption of As(V) on goethite-coated sand using three sequentially linked reactors, identified as a multiple 61 sequential batch reactors (MSER). We simulated the results of their MSER experiment as our fifth benchmark problem. In this experiment, as an initial step, the first reactor was filled with As(V) solution while reactors 2 and 3 were filled with clean water. This is the initial condition for the problem. After equilibrating for 24 hours, the aqueous solution from the first reactor was transferred to the second reactor and was allowed to equilibrate with the solids in the second reactor. During the same time period, new arsenic laden solution was transferred to the first reactor, the solution in the second reactor was transferred to the third reactor, and the solution in the third reactor was discharged for chemical analysis. The volume of water discharged from a single reactor was designated as the ?reactor volume.? The experiment had two distinct phases; during the first phase, a total of 14 reactor volumes were discharged from the system while simultaneously renewing the solution in the first reactor with new arsenic-laden solution. During the second phase, a total of 4 reactor volumes were discharged while simultaneously replacing the solution in the first reactor with clean water. To simulate this experiment, a hypothetical column with 4 finite-difference nodes was used. The length of the each finite difference grid is set to 1 cm, and the total distance between node-1 to node-4 was 3 cms, representing the 3 reactors. The velocity was assumed to be 1 cm/day and the time step used was 1 day, maintaining a Courant number 1. Note, although there were 4 nodes in the system, the boundary node was used as an hypothetical node to define the boundary condition; geochemical reactions are allowed to occur only in nodes 2, 3, and 4, which represented the three sequentially coupled reactors. There is no hydrodynamic dispersion in this sequential batch problem and hence the dispersion module was not used. The Courant number was set to 1 in the numerical model to allow one node explicit advection that exactly mimicked the batch transfer process without any numerical dispersion effects. The total period of simulation was 18 days. The aqueous phase 62 As(V) concentration at the inlet boundary node was set at 1.25 M for 14 days, followed by zero concentration for 4 more days. Other model parameters and the tableau for representing the chemical reactions are given in Tables 2.10 and 2.11, respectively. The total amount of computer time required for solving this benchmark problem was 2 seconds. RT1D simulation results are compared against PREEQCI (Charlton and Parkhurst, 2002) results (reported in Jeppu et al. 2012) and the experimental data (also reported in Jeppu et al. 2012) in Figure 2.25. It can be observed from the figure RT1D matched the published data well. Table 2.10: Model parameters used for problem?9 Length (cm) 4 Total time (days) 18 Pulse time (days) 14 ?x 1 ?t 1 Velocity (cm/day) 1 # Mobile Species 1 # Immobile Species 1 Boundary condition: As(V) concentration (?M) 0-14 days 15-18 days 1.25 0 pH (constant/fixed) 7 Ionic strength 0.01 Surface site density (sites/nm2) 1.04 Surface area (m2/g) 1.08 Sorbent concentration (mg/L) 1.0 63 Table 2.11: Stoichiometry of chemical reactions (tableau) used for problem?9 AsO4- SOH Psi H+ Log K H+ 0 0 0 1 0 AsO4[-3]Aq 1 0 0 0 0 FeOH 0 1 0 0 0 OH- 0 0 0 -1 -13.91 HAsO4[-2] 1 0 0 1 11.23 H2AsO4[-1] 1 0 0 2 18.01 H3AsO4 1 0 0 3 20.16 >FeH2AsO4 1 1 0 3 31.44 >FeHAsO4- 1 1 -1 2 26.18 >FeAsO4-2 1 1 -2 1 20.1 >FeOH2[+] 0 1 1 1 7.17 >FeO[-] 0 1 -1 -1 -9.32 Figure 2.25 Comparison of the RT1D results with the published model results for problem?9 0 . 0 0 . 4 0 . 8 1 . 2 1 . 6 0 5 10 15 20 C o n c e n t r a t i o n ( ? M) R e a c t o r v o l u m e R T 1 D P H R E E Q C I M o d e l J e p p u e t a l ( 2 0 1 2 ) D a t a 64 2.4.10 Ion exchange Valocchi et al (1981) presented an analytical framework that allowed the characterization of certain key concentration profiles during the transport of ion-exchange solutes based on the chromatography theory. The validity of this theory was tested by applying it to a field situation in Palo Alto Baylands in California. The field project involved the injection of advanced municipal effluent into the aquifer. The principal chemical mechanism involved is the heterovalent ion exchange of Na+, Mg2+ and Ca2+. The pore velocity was about 1.01 m/day and the dispersivity was 1 m and the total length of the column was 16 m. The stoichiometric table for this problem is presented in Table 2.12 and the results from the RT1D simulations and the comparison against the published model data are available in Figure 2.26. It can be observed from the figures that RT1D was able to simulate the results well. No. Species Name Na+ Mg2+ Ca2+ X H+ Log K 1 Na+ 1 0 0 0 0 0 2 Mg2+ 0 1 0 0 0 0 3 Ca2+ 0 0 1 0 0 0 4 Na-X 1 0 0 1 -1 0 5 Mg-X 0 1 0 2 -2 0.355 6 Ca-X 0 0 1 2 -2 0.602 7 H+ 0 0 0 0 1 0 Table 2.12 Stoichiometric matrix for the simulation of Test Problem ? 10 65 Figure 2.26. Results for example problem 10 a) Breakthrough profile for Na+ b) Breakthrough profile for Mg2+ c) Breakthrough profile for Ca2+ 2.5 Summary and Conclusions In this study, we have presented the details of a numerical modeling tool for solving a variety of biochemical and geochemical reactive transport problems. The code was developed within the EXCEL Visual Basic platform and it can be run within the standard EXCEL without any additional software. The tool is capable of solving a wide range of kinetic-limited reactive transport problems that could be defined through a reaction package. RT1D can also solve a variety of equilibrium-controlled geochemical transport problems defined through a chemical reaction matrix (also known as the tableau). The capabilities of the tool were demonstrated by solving several benchmark problems of varying level of complexity. The results show that 1.00E +02 1.00E +03 1.00E +04 1.00E +02 1.00E +03 1.00E +04 1.00E +05 Conc ( m g/L ) Cub ic M e te r In je c te d Na+ N a+ ( A n j an i) 1.0 0E +0 1 1.00E +02 1.0 0E +0 3 1.00 E+02 1.00 E+03 1.00 E+04 1.00 E+05 Con ce n ration ( m g/ L ) Cu b ic M et er In j ec t ed M g2+ M g2+ ( an j an i) 1.00E +01 1.00E +02 1.00E +03 100 1000 10000 100000 Conc ( m g/L ) Cub ic M e te r In je c te d Ca2+ Ca2+( an jani ) a b c 66 RT1D simulations closely matched previously published results. RT1D is a flexible tool that allows users to add their own routine to define any type of user-defined kinetic reactions. The geochemistry package can be used to define and solve transport problems involving a variety of surface complexation reactions. The tool is equipped with a robust TVD advection solver, an implicit dispersion solver, and an adaptive time stepping ODE solver to handle any complex problem. RT1D code is a useful tool for laboratory researchers who are interested in analyzing batch and column data within a user-friendly EXCEL platform. 67 Chapter 3 MODELING Dehalococcoides Sp. AUGMENTED BIOREMEDIATION IN A SINGLE FRACTURE SYSTEM 3.1 Introduction Extensive use of various forms of chlorinated ethenes as solvents for dry cleaning and metal degreasing efforts has resulted in widespread contamination of groundwater and soil systems. The toxicity and carcinogenicity potential of these compounds can be high hence they pose a significant threat to human and ecological health (Coleman et al., 2002). At contaminated field sites, depending on the history of the spill and the heterogeneity of the subsurface, the chlorinated solvents discharged as a dense non-aqueous phase liquid (DNAPL) might be immobilized in the form of trapped blobs or as pools. These trapped DNAPL phases continue to dissolve for a long time to form large aqueous chlorinated ethene plumes (Clement et al., 2004). The most commonly observed chlorinated ethenes in groundwater are: perchloroethene (PCE), trichloroethene (TCE), the dichloroethenes (cis-1,2- (cDCE), trans-1,2- (tDCE), and 1,1- (1,1DCE)), and vinyl chloride (VC). Microbial metabolism plays a crucial role in the degradation of these chlorinated compounds. The biodegradation process can occur under both aerobic and anaerobic conditions (Beeman and Bleckmann, 2002). During anaerobic degradation, PCE and TCE are reductively dechlorinated to form cis-DCE, VC and the end product ethane, mostly via halorespiration (Bradley, 2003; El Fantroussi et al., 1998; Smidt and de Vos, 2004). 68 Bioaugmentation remediation methods that employ Dehalococcoides sp. (DHC) have been widely tested for treating chlorinated solvent contaminant plumes (Maymo-Gatell et al., 1997; Schaefer et al., 2009b; Schaefer et al., 2010b). Several laboratory and field experiments have been conducted to study the degradation patterns of chlorinated ethenes in the presence of DHC. Cupples et al (2004) studied the reductive dechlorination of PCE using DHC and developed the Monod kinetics model to account for competition between the electron acceptors. They attributed the accumulation of chlorinated intermediate compounds DCE and VC to the lack of appropriate microorganisms, insufficient supply of donor substrate, or reaction kinetic limitations. Lee et al (2004) used glucose as a model carbohydrate to understand the effectiveness of dechlorination process using a culture obtained from a PCE-contaminated site in Victoria, TX. They developed a numerical model that simulated the batch experiments.This model included kinetic expressions to simulate the competition between fermentors, methanogens, and two separate dehalogenator groups. Their model simulations suggested that the amount of dechlorination achieved was significantly affected by the initial relative population of dehalogenators and the H2 utilizing methanogens. Yu et al (2005) modeled the reductive dechlorination reaction kinetics using two models that employed the Michaelis-Menten equation. In this study, the competitive and Haldane inhibition models were tested by fitting batch kinetic data obtained using three types of dechlorinating populations: PM culture (obtained from a chlorinated solvent contaminated site in Point Mugu, CA), EV culture (obtained from Evanite in Corvallis, OR), and BM culture (a binary mixture of PM and EV). The study demonstrated that for accurate modeling a combination of competitive and Haldane inhibition kinetics is necessary, and such models could 69 make accurate predictions over a broad range of PCE/TCE concentrations. Schaefer et al (2009b) studied the fate and transport of DCE in the presence of DHC. They performed both batch and column experiments to evaluate the transport, growth, and dechlorination activity of DHC in a bioaugmentation column experiment. The results from the column experiment showed that the dechlorination occurred over the entire length of the column. They also observed that the reaction rates of DHC in the column experiments were 200 times more efficient than those observed in batch experiments; however, this 200-fold enhancement was not observed at the field scale (Schaefer et al., 2010a). Schaefer et al. (2010b) performed laboratory experiments in discretely fractured sandstone blocks to study the use of bioaugmentation to treat residual PCE-DNAPL. Results from these experiments indicated significant dechlorination activity and growth of DHC within the fracture. The DNAPL dissolution was enhanced during bioaugmentation by a factor of 5 and the dissolved PCE concentration levels were close to the solubility level. The extent of dechlorination and DNAPL dissolution enhancement depended on the fracture characteristics, residence time, and the dissolved concentration of PCE. Although considerable amount of experimental data are available to test the feasibility of bioaugmentation process to treat chlorinated ethenes, very few reactive transport models are available that can describe these experimental data by simulating the bacterial growth and dechlorination activity coupled with transport processes. Development of such a model would help facilitate the design and operation of bioaugmentation applications. In this study, we propose a comprehensive reactive transport modeling framework for modeling bioaugmentation remediation process that employs DHC to treat a PCE-contaminated fracture system. The model 70 was calibrated and tested using the batch and transport experimental results presented in Schaefer et al (2009b) and Schaefer et al. (2010b). 3.2 Experimental method and governing equations Schaefer et al. (2010b) examined the enhanced dissolution of residual DNAPL sources in bench- scale fractured sandstone blocks during bioaugmentation. The data reported for two experiments used in this study were conducted in an Arizona buff sandstone block (29 cm ? 29 cm ? 5 cm). A discrete linear fracture of aperture size 0.054 cm (Schaefer et al., 2009a) was created along the naturally occurring bedding plane in this sandstone block. The outer edges were sealed and small holes were drilled into the rock along the influent and effluent fracture edges. Artificial groundwater was used in all the experiments. The DHC used in these bioaugmentation experiments was a commercially available microbial culture SDC-9 (Vainberg et al., 2009). The biodegradation experiments were performed after residual PCE-DNAPL saturation was attained in the rocks as described in Schaefer et al (2009a). Residual DNAPL within the fracture served as a long-term source for dissolved PCE throughout the duration of the experiment. The fractures were flushed with an anoxic solution for 2 to 7 days before bioaugmentation with DHC was performed. Two experimental datasets were reviewed and used in this study; these were identified as Experiment-1 and Experiment-3 in the Schaefer et al. (2010b) paper. As reported in Schaefer et al. (2010b), Experiment-1 was a ?high? flow rate experiment whereas Experiment-3 was a ?low? flow rate experiment. In this study, we will identify Experiment-3 and Experiment-1 as Experiment-A and Experiment-B, respectively. We used Experiment-A dataset for model calibration and Experiment-B dataset for testing the calibrated model. 71 Similar to Schaefer et al. (2009a) study, the single fracture system was conceptualized as an equivalent, one-dimensional porous media system with a constant transport velocity. Reductive dehalogenation of PCE into TCE, DCE, VC and ethene and the growth of mobile and immobile DHC were modeled using Monod kinetics (Schaefer et al., 2009b; Yu et al., 2005). The governing transport equations with appropriate kinetic biochemical reaction terms are summarized below: (3.1) 2 T CE T CE T CE T CE T CE P CE P CE T CE T CE2 P CE P CEP CE T CE T CE P CE C C C q X C q X C= - V + D - + - ?C t x x C + KCC + K 1 + I (3.2) 2 D CE D CE D CE D CE D CE T CE T CE D CE D CE2 D CE PCE D CE D CE T CE T CE D CE PCE C C C q X C q X C= - V + D - + - ?C t x x CCC + K 1 + C + K 1 + II (3.3) 2 V C V C V C V C V C D C E D C E V C V C2 D C E T C E D C E V C V C D C E D C E D C E T C E D C E C C C q X C q X C= - V + D - + - ?C t x x C C CC + K 1 + + C + K 1 + I I I (3.4) 2PC E PC E PC E PC E PC E 2 PC E PC EC C C q X C= - V + D -t x x C + K 72 2 E th E th E th V C V C E th E th2 D CE T CE V C V C D CE T CE C C C q X C= - V + D - ?C t x x CCC + K 1 + + II (3.5) - - - 2 C l C l C l PCE PCE T C E T C E 2 PCE PCE PCE T C E T C E PCE D C E D C E V C V C D C E D C E T C E D C E D C E V C V C D C E D C E T C E C C C q X C q X C = - V + D + t x x C + K C C + K 1+ I q X C q X C ++ C C C C + K 1+ C + K 1+ + I I I (3.6) 2 M M M de t I M2 P C E P C E T C E T C E P C E P C E P C E T C E T C E P C E M D C E D C E V C V C D C E D C E T C E D C E D C E V C V C D C E D C E T C E X X X =- V +D +k X t x x q C q C ++ C +K C C +K 1+ I +Y X q C q C + C C C C +K 1+ C +K 1+ + I I I (3.7) 73 P C E P C E T C E T C E P C E P C E P C E T C E T C E P C E IM de t I M I M DC E DC E VC VC DC E DC E T C E DC E DC E VC VC DC E DC E T C E q C q C ++ C +K C C +K 1+ IX =- k X +Y X t q C q C + C C C C +K 1+ C +K 1+ + I I I (3.8) where Ci (mM) is the aqueous concentration of the compound i (where i is either PCE, TCE, DCE, VC, ethene or chloride), V is the transport velocity (cm/hr), D is the hydrodynamic dispersion co-efficient (cm2/hr) Ii (mM)is the competition coefficient, qi (mmol L-1 (cell h) -1) is the DHC maximum utilization rate coefficient, Ki (mM) is the half velocity coefficient, XM (cells/L) is the mobile phase DHC concentration XIM (cells/L) is the immobile phase DHC concentration, X (cells/L) is the sum of the mobile and immobile DHC (X=XM+XIM), Y (cell/mM) is the yield and ?i (h-1) is the back-partitioning coefficient, and Kdet (h-1) is the bacteria cell detachment coefficient. The reactive transport model contains 7 mobile species: PCE, TCE, DCE, VC, ethene, chloride and mobile DHC. The only immobile species present in the system is the DHC bacterial population attached to the solid phase. Note that the immobile species equation does not have advection and dispersion terms. Equations (3.1) ? (3.5) describe biodegradation of PCE into its daughter products by both mobile and immobile DHC. Equation (3.1) describes the degradation of PCE to TCE; equation (3.2) describes formation of TCE from PCE and its subsequent degradation to DCE. Similarly, 74 equations (3.3) ? (3.5) describe the formation DCE, VC and Ethene, along with the formation of the daughter products VC and ethene, respectively. The equations (3.2) ? (3.5) also have a model coefficient (?), which was used to account for the expected loss of aqueous phase contaminant concentration due to back-partitioning of the daughter products into residual DNAPL; this back- partitioning mechanism was described in Schaefer et al. (2010b) and Ramsburg et al (2010). Note the back-partitioning process was not included in equation (3.1) because PCE cannot back- partition onto itself. Equation (3.6) describes the accumulation process for chloride where one mole of chloride is formed for every mole of the daughter product. Equations (3.7) and (3.8) describe the bacterial growth along with the attachment and detachment kinetics for mobile and immobile bacteria cells (Clement et al., 1997; Peyton et al., 1995). In all the numerical simulations we used a one-dimensional finite difference grid of size 1 cm (total of 30 nodes). Other numerical parameters used are summarized in Table 3.1. The above set of equations was solved using the operator-split strategy. The fully-implicit finite difference approximation was used for solving the advection and dispersion terms, and a Runge- Kutta procedure was used for solving the reaction terms on a node-by-node basis. Further details of the numerical scheme used for solving the coupled multi-species reactive transport problem are available in the literature (Clement et al., 2004; Clement et al., 1996; Clement et al., 1997; Clement et al., 1998; Walter et al., 1994; Zheng and Wang, 1999). 75 Table 3.1: Summary of numerical model parameters used for Experiment A Length 29 cm Velocity 6.45 if t < 13 days 1.46 if t > 13 days Dispersivity 5 cm Simulation time 2640 hrs Grid size (?x) 1 cm Time step (?t) 0.1 hr Immobile bacteria per node at t=0 1.00E+06 cells/L Detachment factor (Kdet) 6E-07 hr-1 3.3 Results and discussion 3.3.1 Testing the batch kinetic model against Schaefer et al (2009a) model predictions The biodegradation kinetic model was first tested to reproduce batch simulation results reported in Schaefer et al (2009b). The purpose of this exercise was to test the numerical code by reproducing published batch simulation results. In addition, this modeling step also tested whether a simplified version of the kinetic model presented in the earlier section can be used to reproduce literature results. The governing reaction kinetic equations used in this batch simulation are summarized in Section 2.4.4. Note that the equations shown in the appendix are a simplified version of the reaction model used in equations (3.2)-(3.5). These simplifications were required since this batch study used TCE instead of PCE (Schaefer et al., 2009b). The biodegradation parameters used in batch simulations, summarized in Table 3.2, were obtained from Schaefer et al (2009b). The results from this model comparison study are shown in Figure 3.1. The figure shows that the current model predictions were almost identical to those predicted using the Schaefer et al (2009b) model, indicating an excellent match. 76 Table 3.2: Bio-augmentation parameters regressed from batch experiments in Schaefer et al (2009a). Species K(mM) q (mmol L-1 (cell hr)-1) I (mM) Yield (Cells/mM) PCE* 0.42 1.70E-12 2.50E -01 4.4E+09 TCE 0.0032 1.30E-12 1.00E+06 DCE 0.002 7.00E-13 5.20E -03 VC 0.014 1.40E-12 1.00E+06 *Note: PCE model parameters were estimates provided by the research group identified in Schaefer et al (2009b). Figure 3.1. Batch model results: Comparison of proposel model results against the model results reported in Schaefer et al (2009b) 77 3.3.2 Calibration of the reactive transport model The successful recreation of the batch results reported in Schaefer et al (2009b) provided sufficient basis to build a simulation model for the bio-reactive transport experiment completed by Schaefer et al. (2010b). Note both Schaefer et al (2009b) and Schaefer et al. (2010b) studies have used identical microorganisms for simulating biodegradation. As discussed in the methods section, the low velocity experiment, Experiment-A, was used for the calibration effort. Experiment-A was the most comprehensive dataset since it included chloride ion concentrations. Chloride data was important information for verifying the mass-balance closure of calibrated model results. It was inferred from the experimental data that there was substantial variations in dissolved PCE concentrations in both experiments. Such variation is not uncommon for DNAPL dissolution studies in bedrock fractures systems (Dickson and Thomson, 2003). For the purpose of modeling, we conceptualized that most of the residual DNAPL was trapped near the inlet and contributed to the dissolved aqueous PCE concentrations. However, in this study, we did not explicitly model DNAPL dissolution processes; instead an equivalent input PCE concentration signal was estimated from the measured values of PCE and its daughter product concentrations. This function was used to define the average daily input PCE concentration levels at the inlet. To model the bioaugmentation step, it was assumed that the initial inoculation process equally distributed DHC among all the 30 nodes. The model was constrained to allow bacterial accumulation at a node to a maximum limiting value of 1.0E+11 cells/L; this is a common approach used for preventing unrealistic accumulation of cells within pore spaces. For example, Zysset et al (1994) studies used a parameter ?max to describe the maximum capacity for the adhering bacteria. The limiting parameter used in this study employs a similar methodology. 78 Two distinct values of DHC maximum utilization rate for DCE degradation, qDCE, were used to account for possible inhibitory effects due to presence of high concentrations of PCE (Amos et al., 2009; Yu and Semprini, 2004). When the PCE concentration was less than 0.3mM, the value of qDCE was assume to be identical to the batch value. Yu and Sempirini (2004) have indicated that the system might show possible inhibitory effects when the concentration of PCE was more than 0.3mM. Therefore, under high PCE concentration conditions, the value of qDCE was allowed go below the batch-estimated value. The lower rate was estimated as a part of the calibration process. During the calibration step, we used the Monod parameters from the batch experiments, shown in Table 3.2, as reference values and perturbed them by an order of magnitude and used a trial-and-error process to fit the observed experimental data. In addition to adjusting the Monod parameters, we also fitted the back-partitioning coefficient. Based on chloride mass balance results, Ramsburg et al. (2010) and Schaefer et al. (2010b) have postulated that certain fraction of degraded daughter products can back-partition to the original DNAPL phases, thus the original DNAPL can serve as a sink for the daughter species. It was observed from the experimental data that DCE was the major daughter product which indicated maximum back partitioning. Hence, the back-partitioning coefficient for DCE was identified as the primary fitting parameter. It is logical to assume that the value of back partitioning coefficient ? would depend on the solubility level (which is affinity of the species to remain in the aqueous phase). Based on this assumption, the values of ? for the remaining daughter products were simply scaled using their respective values of solubility. The equation used for this scaling process was: 79 i D C ED C E s o l u b i l i t y i n m M? = * ?S o l u b i l i t y o f d a u g h t e r p r o d u c t i (3.9) Appendix A2 provides various values ? and the detailed calculation procedure. Figure 3.2 shows the comparison of the final (fitted) model results against the Experiment-A data. The calibrated values of Monod parameters are summarized in Table 3.3. It can be observed from the Figure 3.2 that the results from the model simulations closely follow the trends observed in the experimental results. A qualitative assessment of the model indicated that the model results showed good overall mass balance and predicted the observed chloride ion concentrations reasonably well. Comparison of the values shown in Table 3.3 and Table 3.2 indicate that the fitted model parameters for the transport experiment are within an order of magnitude of the batch parameters. The observed differences in model parameters estimated for the transport and batch model could be due to the heterogeneities present in the fractured system, and/or the elevated PCE and DCE concentrations observed in the fracture experiments compared to the batch experiments that were used to derive the Monod parameters. Several published studies have concluded that batch and column parameters could differ due to pore-scale variations and other heterogeneities (Brusseau, 2006; Jeong-Hun Park et al., 2001). Table 3.3: Calibrated Monod parameters for the Experiment A Species K(mM) q (mmol L-1 (cell hr)-1) I (mM) Yield (Cells/mM) PCE 0.42 0.87E-12 2.50E -01 1.85E+10 TCE 0.0032 1.05E-12 1.00E+06 DCE 0.002 1.85E-13, 7.00E-13 5.20E -03 VC 0.014 1.05E-12 1.00E+06 80 Figure 3.2 Model calibration results. Continuous lines represent the model results and the symbols represent Experimental-A data 3.3.3 Testing the reactive transport model The calibrated model developed using the low flow rate data, Experiment-A, was used to simulate the high-flow experimental system, Experiment-B. It should be noted that the all the Monod parameters used in the validation simulation were identical to those used estimated in the calibration step. However, the validation step required scaling of two physical transport parameters?the back-partitioning coefficient (? for DCE) and the detachment factor (Kdet) to scale for the high velocity conditions. The back portioning coefficient is likely a function of the 81 residence time (higher velocity would allow less time for daughter products to back-partition into the DNAPL). Also, the detachment process would depend on shear forces in the system (higher velocity would induce more shear, as observed in Schaefer et al. 2010b). Therefore, we recalibrated these two values to reflect the new transport conditions. The modified values of ?DCE and kdet are 0.004 hr-1, and 0.0018 hr-1, respectively. Other parameters for this experimental study are presented in Table 3.4. The simulation results along with the experimental data for all daughter products are summarized in see Figure 3.3. A comparison of the measured PCE concentration data and the modeled PCE concentrations for Experiment-B is shown in Figure 3.4. The parent and daughter product data were reported in separate figures since the concentrations of PCE were much higher than any of the daughter products. The results from model-simulations closely follow the trends observed in the experimental data. Table 3.4: Summary of numerical model parameters used for Experiment B Length 29 cm Velocity 6.45 cm hr-1 Dispersivity ) 5 cm Simulation time 5712 hr Grid size (?x) 1 cm Time step (?t) 0.1 hr Immobile bacteria per node at t=0 1.33E+07 cells/L Detachment factor (Kdet) 0.0018 hr-1 82 Figure 3.3 Model validation results. Continuous lines represent the model results and the symbols represent Experimental-B data 83 Figure 3.4 Comparison of modeled and observed PCE concentration levels 3.3.4 Sensitivity Analysis To further understand the sensitivity of the model to variations in the parameter values, we completed a sensitivity study. During the calibration step we have identified that the yield coefficient and DHC maximum utilization rate coefficient (q) were the most sensitive model parameters. We focused on these two model parameters and perturbed them within an order of magnitude and explored its effects on model predicted results for Experiment-A. The results are summarized in the following sections. 84 3.3.4.1 Model response to variations in the yield coefficient The yield coefficient is an important parameter that governs the rate of growth of the bacteria. In this analysis, the yield coefficient was perturbed by the following two factors: 0.5 and 5. Figures 3.5 and 3.6 provide a summary of these simulation results. It can be observed from these figures that the growth of the bacteria increased with increase in yield coefficient and thereby causing significant amount of biodegradation. It can also be observed from Figure 3.6 that when the yield coefficient is lower, the amount of bacterial growth was low and the biodegradation activity is delayed until 50 days instead of the expected 20 days. However, when the yield coefficient is higher, as shown in Figure 3.6, the bacterial growth is faster and the biodegradation occurs earlier. 85 Figure 3.5 Model response to an decrease in the yield value (by 0.5) for Experiment-A 86 Figure 3.6 Model response to an increase in the yield value (by 5) for Experiment-A 3.3.4.2 Model response to variations in the DHC maximum utilization rate constant (q) The parameter q governs the rate at which the chlorinated ethenes are degraded into their daughter products. The DHC maximum utilization rate constant (q) was varied by the following two factors: 0.5 and 2. Figures 3.7 and 3.8 provide a summary of these simulation results. It can be observed from the figures that when the values of q were higher, the biodegradation of PCE and its daughter products was higher as shown in Figure 3.8. This increase in biodegradation activity was marked by a significant increase in the chloride ion concentrations. However, decrease in the q value, decreased the biodegradation rate. Therefore, the predicted aqueous 87 concentrations of PCE and its daughter products are considerably lower when compared to the experimental data. Figure 3.7 Model response a decrease in the q value (by 0.5) for Experiment-A 88 Figure 3.8 Model response to an increase in the q value (by 2) for Experiment-A 3.4 Summary and Conclusions We present a mathematical modeling framework that can be used to model bioremediation of PCE contamination using Dehalococcoides in a single fracture system. The model was calibrated to a previously published experimental dataset. The performance of the calibrated model was tested by completing another simulation where the model was used to predict a transport scenario involving higher flow rate. The model was able to predict the high-flow rate dataset using the calibrated values of biological process parameter. Only minimal changes had to be made to scale transport parameters such as back-partitioning coefficient and the detachment 89 rate, to reflect the high flow conditions. Our sensitivity studies show that the yield coefficient and the DHC maximum utilization rate coefficient are highly sensitive parameters that need to be carefully calibrated. However, the batch estimates of yield coefficient and maximum utilization rate coefficient provide good estimates to guide the calibration process, and final calibrated values were within an order of magnitude of the batch-kinetic values presented in Schaefer et al (2009b). The disparities in the process scales between the batch and fracture experiments are probably the cause of the variations in the parameter values. The sensitivity analysis studies show that some of these parameters are highly sensitive and can alter the biodegradation processes significantly. Both batch and flow datasets considered in this study were obtained from well-controlled experimental conditions that might not reflect actual field conditions. Therefore, future experiments and simulations studies should include field observations before this modeling framework could be up scaled to predict field scenarios. 90 Chapter 4 ASSESSMENT OF PARALLEL GENETIC ALGORITHMS FOR CALIBRATING ONE- DIMENSIONAL MULTI-COMPONENT REACTIVE TRANSPORT MODELS 4.1 Background Reactive transport models have been commonly used to simulate the fate and transport of contaminants in both laboratory and field-scale problems. The accuracy and reliability of these models would strongly depend on the values of model parameters, which are commonly estimated from controlled laboratory and/or field experiments. These experiments are often conducted by isolating certain reaction steps, to fully understand the complex bio-geochemical interactions occurring in the subsurface. The experimental data obtained from the laboratory experiments are then used to formulate more general bio-kinetic or geochemical models that can describe contaminant transformation processes. Once the process model is formulated, several unknown parameters in the overall model are normally estimated by a trial and error process to minimize the sum-squared errors between the experimental data and the model fitted data (Engesgaard and Kipp, 1992; Gramling et al., 2002; Schaefer et al., 2009b; Torlapati et al., 2012). However, the trial-and-error process could become inefficient as the number of parameters in the model increases. Therefore, some type of numerical inverse routines are employed (e.g., CXTFIT (Toride et al., 1995)) to automatically estimate the model parameters. However, several of these inverse methods might converge to a local minima and their overall performance would depend on the robustness of the search algorithm and the choice of the 91 initial parameters supplied by the user (Toride et al., 1995). Doherty and Hunt (2010) developed a more robust parameter estimator, PEST, for solving highly parameterized groundwater problems using regularized inversion schemes. Baginska et al. (2003) applied the Annualized Agricultural Nonpoint Source Model (AnnAGNPS) for the prediction of export of nitrogen and phosphorous in Currency Creek of the Sydney Region. In addition, they have also used PEST to determine the sensitivity and importance of the key parameters of the model. Yabusaki et al. (2007) used PEST by coupling with BIOGEOCHEM to automate the calibration procedure in understanding the transport and bioreduction of Uranium. More recently, genetic algorithms have been employed in parameter estimation of column and batch reactive transport experiments (Majdalani et al., 2009; Massoudieh et al., 2008). Genetic Algorithms (GAs) are a branch of evolutionary algorithms which are primarily used to optimize nonlinear problems in various fields (Massoudieh et al., 2008). The development of GAs are based on the concept of natural selection and the rearrangement of genetic material (Holland, 1975). In the field of groundwater hydrology, the GAs have been used in the optimization of the pumping problem and for the estimation of system parameters in heterogeneous aquifers (El Harrouni et al., 1996). Wagner (1992) applied GAs to estimate the transport parameters in kinetically controlled one- or two-site sorption models. Wang (1997) studied the usefulness of GAs for calibrating rainfall-runoff models with nine parameters and found that the GA was able to attain the global minimum for a hypothetical catchment. Wang and Zheng (1997) have coupled MODFLOW and MT3D with a GA routine to find the optimal pumping and injection rates for a remediation process. They applied the model to a three dimension field problem and demonstrated the superiority of their GA solution to an existing solution obtained using a trial and error approach. Mulligan and Brown (1998) used a GA to 92 optimize the water quality model parameters and found that GA was a useful calibration tool to estimate least squares parameters by accumulating useful information about the response surface. Reed et al. (2000) studied GAs to find a theoretical relationship for population size and number of generations required for convergence in groundwater well monitoring design applications. Giacobbo et al. (2002) investigated the feasibility of using GAs for estimating groundwater contaminant transport parameters for a three-layered one-dimensional saturated flow and transport problem. Singh et al. (2005) presented an interactive GA to solve an inverse problem that estimated the conductivity of a heterogeneous hypothetical aquifer whose value was known a priori. B?ranger et al. (2005) coupled a GA with an analytical, one-dimensional, multi- component, reactive transport model to estimate the first-order decay coefficients and enrichment factors. Singh et al. (2008) developed a novel interactive framework, called the ?Interactive Multi-Objective Genetic Algorithm? (IMOGA), to solve the groundwater inverse problem considering different sources of quantitative data and qualitative expert knowledge. Massoudieh et al. (2008) used GA to minimize the error between measured and modeled breakthrough data for reactive transport involving Cd and tributyltin, and estimated the equilibrium constants. Lee and Heber (2010) combined GA with biofiltration models to estimate unknown model parameters, and the model was subsequently used to predict ethylene removal efficiencies. Madsen and Perry (2010) coupled a simple GA with MODFLOW to optimize the net groundwater flow into a river by optimizing the following four input parameters: recharge rate, river conductance, and water levels at two general head boundaries. While GAs are useful tools, they are also computationally intensive routines since they search through a large population to find the optimal solution. The search process could take a substantial amount of time, if it is not optimized. Parallel computing techniques can be used to 93 improve the efficiency of GAs by exploiting the concurrency of calculations performed in genetic algorithms. Depending on their architecture, the computers capable of running parallel codes can be classified as either distributed memory computers or shared memory computers (Cant?-Paz, 1998). Most of the earlier work on parallel computing efforts focused on distributed memory computers, where several computers are connected using a fast network to reduce communication time between the processors to implement parallel genetic algorithms (Abramson et al., 1994; Baluja, 1992; Fogarty and Huang, 1991; Tanese, 1989). In the field of groundwater, McKinney and Lin (1994) used parallel genetic algorithms to solve three groundwater management problems involving maximization of pumping from an aquifer, minimization of cost for a water supply problem, and minimization of cost for an aquifer remediation problem. They observed that the genetic algorithms performed effectively to obtain globally optimal solutions and the speedup of the parallel genetic algorithm was almost linear. Tsai et al. (2009) developed a production well management model for the water resource management in semi-arid areas by integrating a large-scale pressurized water distribution system management model, EPANET, and a three-dimensional groundwater model, MODFLOW, under a unified optimization framework. They used a 64 processor cluster to run the computer code in a parallel mode. The speedup on distributed memory computers is hindered by the communication time between the processors because each processor has its own local memory, which is not available to the other processors; hence, the programmer has to manually sync the variables after each generation. However, in shared memory computers, all the processors have access to the same memory and the synchronization step can be avoided (Abramson and Abela, 1991). Sarma and 94 Adeli (2002) used parallel fuzzy genetic algorithms for optimizing steel structures using two different schemes. The authors also presented two bi-level parallel genetic algorithms that combine message passing interface (MPI) and OpenMP programming languages for optimization. They observed almost linear speedup for 16 processors. Fredrickson et al. (2003) evaluated the performance of parallel genetic algorithm (PGA) using OpenMP constructs, kernels and application benchmarks on large-scale SMP systems using a 72 node Sun Fire 15k SMP node. They reported the basic timings, scalability and run times for different parallel regions. GAs are robust algorithms that have been proven to be suitable for solving different types of parameter estimation problems using an appropriate encoding method. The process by which a population is coded into a suitable form that enables genetic recombination is called encoding. The early studies of GA in reactive transport problems are limited by their usage of binary encoding especially when the parameters of different magnitudes are present (El Harrouni et al., 1996; Massoudieh et al., 2008). Also, most of these algorithms have been optimized to solve a single problem and their ability to run different kinds of reactive transport problems has not been explored. Moreover, none of these studies considered optimizing the implementation of parallel GA algorithms for multicore personal computers that use shared memory architecture. Shared memory, multi-core PCs have become common computational platforms in the recent years with the introduction of INTEL and other multicore processors in desktop and laptop computers. These multicore systems are powerful processors that can be used to improve the efficiency of current GA algorithms by implementing them using a shared memory, parallel computing languages such as OpenMP FORTRAN. Currently, there are no studies available in the hydrogeology literature that uses an OpenMP platform to parallelize a GA algorithm for 95 estimating model parameters in multi-component reactive transport models. The objective of this study is to develop a general parallel genetic algorithm (PGA) that is capable of estimating both transport and kinetic parameters in reactive transport models. We compare the performance of the PGA using four different benchmark problems. Speedup data for the PGA are also presented. 4.2 Methodology ? General Steps in Genetic Algorithm The six key steps involved in a traditional GA are: encoding, population generation, selection, crossover, mutation, and termination (Holland, 1975). The GA starts with a randomly-generated initial set of solutions (also known as chromosomes) and this is called the initial population. The fitness of this population is calculated using the objective function. The fitness of each chromosome in the population is used to assess its ability to survive the current generation. For a minimization optimization problem, a lower value of fitness is desirable. Based on this fitness value, two parents are selected using a selection process. The selected parents undergo a crossover, where the genetic information is exchanged between the parents using a crossover function. Since, the genetic information is transferred to the subsequent generation of children it is always preferable to choose individuals with better fitness in the selection process. It is also possible that an offspring generated from the crossover of the parents could undergo a mutation operation governed by a mutation probability. The fitness of the offspring is calculated and is combined with the entire population. The individuals with poor fitness are removed from the population (death) at the end of the generation. There are several strategies available for the discarding bad solutions, and for implementing the process of encoding, selection, crossover, and mutation. The specific methods used in this study are discussed below. 96 4.3 Genetic Algorithm Implemented in this Study In this study, a real value encoding is used because each parameter value in our problem could have a different magnitude. It has been observed that for engineering applications, which are sensitive to parameter variations, real value encoding method performs better than the binary encoding method (Gaffney et al., 2010; Michalewicz, 1996; Reed et al., 2000) . We generated an initial population of 32 solutions within a specified range given by the user. The parameter values were then transformed to log (of base 10) scale and a uniform random number (distributed between 0 and 1) was used to generate various random parameter values using the formula: log(low)+r*[log(high) - log(low)], where r is the random number. The random parameter values were then raised to the power of 10 (to transform back to real number scale) and the value was used for populating the chromosome. A one-dimensional multi-component reactive transport model was used to simulate the concentrations. The concentrations generated from these parameters were used to calculate the sum square of errors (SSE) between model predicted concentrations and concentrations obtained from experimental data. This calculated SSE value was assigned as the fitness parameter for that particular chromosome. The selection of parents was done using a tournament selection method (Koza, 1992). In this method, the algorithm randomly selected 5 possible candidates for the parents from the population and the individual with the best fitness is chosen as the parent. The process was repeated to find the second parent. Tournament selection allows the selection of individuals with best fitness so that their genetic material can be passed on to the next generation (Koza, 1992). The selected parents underwent a crossover using a weighted average. The weights between the parents are chosen by randomly generating a real number between 1 and 0.5. If the chosen random number is r, then the new parameters are calculated by the formula: r* p a re n t1 + (1 -r)* p a re n t2(Cant?-Paz, 1998). This 97 weightage within the bounds of 1 and 0.5 allows us to keep the offspring within the boundaries specified at the beginning of the program. If the random number generated is close to 0.5, then we have an average of both the parents whereas if the random number generated was close to the higher bound (of 1), then the value of the offspring will be closed to the first parent. A total of 8 children were generated by performing the crossover 8 times. The total number of children and the initial population were ensured to be multiples of 4 so that the total load distributed on each processor (we used 4 Pentium processors, details given below) during parallelization is equal. These children could undergo a mutation step if a randomly generated number is less the probability of mutation (Pm). The mutation operator used in this algorithm multiplies the parameter by 0.5 before ensuring that it does not cross the bounds set at the beginning of the program. The fitness of the offspring is calculated and is combined with the initial population. The population is then sorted according to its fitness and the best 32 solutions are preserved for the next generation. The best solution is always preserved in this fashion and hence this algorithm can be classified as an elitist approach. The process of selection, crossover, mutation and death were repeated for about 100 generations, and it was observed from our sensitivity analysis studies that the GA solution does not improve after about 100 generations. 4.4 Parallelization of the Genetic Algorithm GA provides a natural and easy approach for parallelization within each generation since most of the loops within a generation contain variables that are not dependent on its value at the previous iteration. This allows for little to virtually no communication time between the processors for synchronization. The parallelization of the GA was achieved by using the shared memory programming procedure OpenMP available within the Intel FORTRAN90 compiler. The desktop 98 computer used for performing simulations used an Intel Xeon processor with two dual core processors, with a total of four processors available for parallelization. The parallelization was accomplished by placing OpenMP constructs at the beginning and the end of the loop that is desired to be run in parallel mode. The OpenMP constructs are also used to specify the number of processors to be used for parallelization and the variables that are private or public to each processor and the kind of schedule to be used to distribute the load among the processors. A guided schedule was used in this study. The loops that were parallelized include the fitness calculation of the initial population since this was the most time consuming part of the program. Also, the fitness calculations of the offspring were completed in a parallel mode. Figure 4.1 illustrates the computational steps involved in implementing the PGA algorithm for a four processor system. Although the selection, crossover and mutation processes can be performed in parallel, these are not computationally intensive tasks; we found the performance gains to be marginal when these loops were optimized. 99 Figure 4.1 Illustration showing the flow of a parallel genetic algorithm 4.5 Details of the Numerical Model used for Fitness Calculation To calculate the fitness of the chromosomes, a multi-component one-dimensional reactive transport model was used. The simulation model is a Fortran version of a previously published Visual Basic software RT1D (Torlapati and Clement, 2012b). The model solves a set of advection-dispersion-reaction equations that describe the transport of ?m? mobile components 100 and ?n? immobile components. The general form of the transport equations are (Torlapati and Clement, 2012b): 2i i i ii2C C CR = - V + D + ?t x x where i = 1, 2,3? m (4.1) jjjSR=?t where j = (m+1), (m+2), (m+3)+? (m+n) (4.2) Where V is the velocity (m/day), D is the hydrodynamic dispersion coefficient (m2/day), Ci is the aqueous phase concentration (mg/L) of mobile component ?i,? where i = 1, 2...m; Sj is the solid phase concentration (mg/mg) of immobile component ?j,? where j = m+1, m+2... m+n; and ?i & ?j are the reaction terms for the mobile and immobile components, respectively. Note the immobile component equations do not have the advection dispersion terms but will have a reaction term that might be coupled to other reaction terms in mobile components. The equations (4.1) and (4.2) are numerically solved using the operator split strategy (Clement et al., 1998; Torlapati and Clement, 2012b). An implicit finite difference scheme was used to solve the advection-dispersion part of the equation and the reaction part, which reduces to a set of ordinary differential equations, is solved using a Runge-Kutta-Felhberg with an adaptive time stepping (Chapra and Canale, 1998). The model concentrations obtained for each chromosome were used to calculate the absolute error using the given experimental dataset. This error was squared and a sum of all these errors was calculated and was designated as the fitness for the chromosome. The objective of the PGA was to minimize this sum square error. 101 4.6 Results and Discussion To test the performance of PGA, four test problems that have different levels of complexity were selected. All of these problems either have published experimental data, or analytical solution that can be used to verify the results. After completing the parameter identification step, we performed speedup tests by varying the number of threads and quantified the advantages of adding additional processors. We also performed sensitivity analysis tests to quantify the sensitivity of the solution to changes in the initial population size and total number of generations. 4.6.1 Benchmark Problem 1 ? Parameter estimation in a rate-limited sorption problem In this benchmark problem, we solve a rate-limited sorption process where non-equilibrium conditions exist. Clement et al. (1998) and Torlapati and Clement (2012b) modeled these kinetic processes using the following governing equations. 2 2 d C C C SVD ? C -t x x K (4.3) d dS ?? SC-dt ?K (4.4) where C is the concentration in the component in the aqueous phase (mg/L), S is the concentration of the component in the solid phase (mg/mg), ? is the bulk density (mg/L), ? is the porosity, Kd is the linear sorption constant (L/mg), k is the first order decay constant (day-1), and ? is the mass transfer coefficient (day-1). 102 The unknown parameters in this problem are: D, ? and Kd. The model was run in the forward simulation model using known parameter values, and the aqueous phase concentrations predicted after 50 days was used as the data. The pore velocity used in this problem was 0.53 cm/day. The concentration values were made available at every 2 cm over 30 cm long column. The temporal and spatial time steps used are 0.4 cm and 0.01 days respectively. The lower and higher bounds for each unknown model parameters were perturbed by two orders of magnitude (one in each direction), as Table 4.1. The parameters estimated by the PGA after 100 generations are given in Table 4.1 along with their true values. The minimum value of fitness obtained at the end of PGA simulations was 2.5E-05. It can be seen from Table 4.1 that the parameter values estimated by the code are close to the original values. Comparison of the concentration profiles simulated using PGA-estimated parameter values and the ?true? parameter values are shown in Figure 4.2. It can be seen from the figure that PGA-estimated model predictions fit this synthetic dataset well. Table 1 also shows that the parameter values estimated by PGA only have about 3 to 5% difference from the original estimates. Table 4.1: Comparison of the PGA estimated parameters with true solutions along with their bounds and percentage error for benchmark problem -1 Parameter True Value Low High PGA Estimate Error % Longitudinal dispersion coefficient, D (cm2/day) 8.00E-02 1.00E-02 1.00E-01 7.80E-02 2.50% Mass transfer coefficient, ? 1.50E-02 1.00E-02 1.00E-01 1.43E-02 4.67% Linear sorption constant, Kd (L/mg) 1.84E-04 1.00E-04 1.00E-03 1.90E-04 3.26% 103 Figure 4.2 Comparison of results from PGA estimated parameters and the true parameters for benchmark problem ? 1 4.6.2 Benchmark Problem 2 ? Parameter estimation in a sequential decay problem Quezada et al. (2004) presented analytical solutions for solving coupled multi-dimensional multi- species transport equations involving first-order kinetic interactions. This is a four-component problem with four kinetic parameters. The governing transport equations are: 21 1 1 1 1 12C C CR = - V + D - k Ct x x (4.5) 22 2 2 2 c 2 / c 1 c 2 / c 1 1 1 2 2 c 2 / c 3 c 2 / c 3 3 32C C CR = - V + D + F Y k C - k C + F Y k Ct x x (4.6) 23 3 3 3 c 3 / c 1 c 3 / c 1 1 1 c 3 / c 2 c 3 / c 2 2 2 3 32C C CR = - V + D + F Y k C + F Y k C - k Ct x x (4.7) 24 4 4 4 c 4 / c 2 c 4 / c 2 2 2 c 4 / c 3 c 4 / c 3 3 3 4 42C C CR = - V + D + F Y k C + F Y k C - k Ct x x (4.8) 104 where Ci is the aqueous concentration of the component i (i=1, 2, 3 or 4) (mg/L), ki is the first order degradation constants for the component i (day-1), Y is the yield coefficient between two components, F is the fraction of yield between two components. The unknown parameters in this benchmark problem are: D, k1, k2, k3, and k4. The concentrations of all the components at an interval of 1 cm along the 30 cm long column predicted after 50 days of transport were made available for the fitness calculation. The pore velocity used in this problem was 0.4 cm/day. The time step and the grid size used for the simulations were 0.1 days and 0.1 cm respectively. The yield values were set to 1. The simulation was run for 100 generations and the minimum SSE observed was about 1.5E-2. The comparison of the PGA estimated parameters along with their low and high bounds used assumed in the simulation are given in Table 4.2. The concentration profiles generated using PGA-estimated parameters and the ?true? parameters are shown in Figure 4.3. It can be observed from the figure that both solutions match well. The error in the parameter values vary by 4 to 38%. Table 4.2 Comparison of the GA estimated parameters with true solutions along with their bounds for benchmark problem - 2 Parameter True Value Low High PGA Estimate Error % Longitudinal dispersion coefficient, D (cm2/day) 8.00E-02 1.0E-03 1.0E-01 8.40E-02 5.00% Decay constant for Component 1, k1 (day-1) 7.50E-02 1.0E-03 1.0E-01 7.80E-02 4.00% Decay constant for Component 2, k2 (day-1) 5.00E-02 1.0E-03 1.0E-01 5.80E-02 16.00% Decay constant for Component 3, k3 (day-1) 2.00E-02 1.0E-03 1.0E-01 2.77E-02 38.50% Decay constant for Component 4, k4 (day-1) 4.50E-02 1.0E-03 1.0E-01 4.00E-02 11.11% 105 Figure 4.3 Comparison of results from PGA estimated parameters and the true parameters for benchmark problem ? 2 4.6.3 Benchmark Problem ? 3: Parameter estimation for a TCE Biodegradation Model Schaefer et al. (2009b) conducted batch experiments to study the degradation of TCE in the presence of Dehalococcoides Sp. They used modified Monod kinetics to model the bioaugmentation process and the associated biochemical reactions. The kinetic equations are: TC E TC E TC E TC E TC E TC E d C q X C1=-d t R C + K (4.9) D CE D CE D CE T CE T CE D CE T CE T CE T CET CE D CE D CE T CE d C q X C q X C11=- d t R R C + KCC + K 1 + I (4.10) 106 V C V C V C D C E D C E V C D C ET C E D C E T C E V C V C D C E D C E T C E D C E T C E d C q X C q X C11= - + d t R RC C CC + K 1 + + C + K 1 + I I I (4.11) T CE T CE D CE D CE V C V C T CE T CE T CE D CE V CT CE T CE D CE D CE D CE V C V C T CE T CE D CE q C q C q Cd X 1 1 1= Y X + + d t R C + K R RC C CC + K 1 + C + K 1 + + I I I (4.12) Where Ci (mM) and X (cells/L) are the concentration of ith component and biomass respectively; i can be either TCE, DCE and VC; qi is the maximum biomass utilization rate (mmol/L/(cell h)), Ki is the half velocity coefficient of the compound (mM), I is the competition coefficient (mM), Ri is a retardation term that accounts for the presence of air in the system (Schaefer et al., 2009b). The unknown model parameters are: qTCE, qDCE, qVC, KTCE, KDCE, KVC and IDCE. The batch simulation experiments were performed for a total of 12 days and the time step used was about 0.01 days. The initial concentration of TCE, ethene and biomass were 0.08 mM, 0.003 mM and 2.8E+10 cells/L respectively. The biomass yield coefficient used in this problem was 4.4E+09. The concentrations of TCE, DCE, VC and ethene after every hour were provided for the calculation of the fitness of the PGA. The PGA was run for 100 generations and the PGA estimated model parameters are compared with the ?true values? and the corresponding relative errors are given in Table 4.3. The table also provides the low and upper bounds used in this PGA search. The comparison of the concentration profiles predicted using the PGA-estimated parameters and the ?true? parameters are shown in Figure 4.4. The figure also shows the actual 107 experimental measurements. It can be observed from the figure that the parameters estimated from the PGA are able to predict the experimental results reasonably well. Table 4.3: Comparison of the PGA estimated parameters with true solutions along with their bounds for benchmark problem - 3 Parameter True Value Low High PGA Estimate Error % Biomass utilization rate for TCE, qTCE (mmol/L/cells h) 3.2E-03 1E-04 1E-02 5.47E-03 70.94% Biomass utilization rate for DCE, qDCE (mmol/L/cells h) 2.0E-03 1E-04 1E-02 6.12E-03 206.00% Biomass utilization rate for VC, qVC (mmol/L/cells h) 1.4E-02 1E-03 1E-01 7.10E-02 407.14% Half velocity constant for TCE, KTCE (mM) 1.3E-12 1E-13 1E-11 1.42E-12 9.23% Half velocity constant for DCE, KDCE (mM) 7.0E-13 1E-14 1E-12 7.46E-13 6.57% Half velocity constant for VC, KVC (mM) 1.4E-12 1E-13 1E-11 5.10E-12 264.29% Competition coefficient for DCE, IDCE (mM) 5.2E-03 1E-04 1E-02 4.20E-03 19.23% 108 Figure 4.4 Comparison of results from PGA estimated parameters and the true parameters for benchmark problem ? 3 4.6.4 Benchmark Problem ? 4: Parameter estimation for a carbon tetrachloride bioremediation problem Phanikumar et al. (2002) conducted experiments to study the bioremediation of carbon tetrachloride (CT) contaminated column which was intermittently fed with nutrient such as acetate and nitrate. They developed a reactive transport model for the system and used a modified version of the RT3D code to simulate their experimental data. The model included a total of 4 mobile and 2 immobile components. The mobile components in the system were CT, acetate, nitrate and the mobile bacteria whereas the immobile components were the sorbed CT and immobile phase bacteria. The kinetic reaction equations used in the reactive transport model are: 109 'd C T C T M I M d C T C T? fK d C ??1 - k C ( X + X ) - 1 - f K C S? d t ? (4.13) a m a x a na M I M a dC ? M MR - ( X X )d t Y (4.14) m a x a n K Cnn M I M a n M I M n n b ? M M bdCR ( X + X ) - ( 1 M ) ? M (X + X )d t Y Y (4.15) M m a x a n K C a a t M d e a I MdX ? M M b (1 M ) K X K (1 M )Xdt (4.16) CT d C T C TdS ? 1-f K C Sdt (4.17) IM m a x a n K C a d e a I M a t MdX ? M M b (1 M ) K (1 M ) X K Xdt (4.18) where f is the fraction of equilibrium sites, bKC is the microbial decay rate (day-1), Kat is the attachment coefficient (day-1), Kde is the detachment coefficient (day-1), k` is the CT reaction rate (day-1), ? is the nitrate reaction rate (day-1), ? is the kinetic desorption rate (day-1), ?max is the maximum specific growth rate (day-1), Ya, Yn and Ynb are the yield rates of acetate, nitrate and biomass respectively; CCT, Ca, Cn and SCT are the aqueous concentrations of carbon tetrachloride, acetate, nitrate and the sorbed concentration of carbon tetrachloride, respectively; XM and XIM are the concentrations of mobile and immobile bacteria, respectively. Also, Ma and Mn are the Monod terms for acetate and nitrate reactions, respectively, and given by the expressions: aa sa a CM KC and nn sn n CM KC where Ksa and Ksn are the half saturation coefficients of 110 acetate and nitrate utilization reactions, respectively. Further details of the model and the experiment are given in Phanikumar et al. (2002) and Torlapati and Clement (2012b). The unknown parameters in this model are: k`, ?, Kde, and bKC. The known parameters for this model are summarized in Table 4.4. The known concentrations values after 4 days of operation at each node point were supplied to the PGA to calculate the fitness. The PGA was run for 100 generations and the SSE after 100 generation was 666. The PGA-estimated parameters are compared against the ?true? values in Table 4.5; the table also provided estimated error for each parameter value and the lower and higher bounds values. The concentration profiles predicted by the PGA-estimated parameters and compared against the simulated data points in Figure 4.5. It can be observed from the figure that the results compare well. 111 Table 4.4: Model parameters used for benchmark problem ? 4 Parameter Value Pore Velocity (cm/day) 10 Length (cm) 200 Longitudinal dispersion coefficient (cm2/day) (D) 2 ?x (cm) 1 ?t (days) 0.001 Porosity (?) 0.35 Bulk density (?) (mg/L) 1.63E6 Time (days) 4 Fraction of equilibrium sites (f) 0.437 Attachment coefficient (day-1) (Kat) 0.9 Distribution coefficient (Kd) (L/mg) 3.9E-7 Half saturation coefficient: (mg/L) Acetate (Ksa) Nitrate (Ksn) 1.0 12.0 Kinetic desorption rate (day-1) (?) 0.36 Maximum specific growth rate (day-1) (?max) 3.11 Yield: Acetate (Ya) Nitrate (Yn) Biomass (Ynb) 0.4 0.25 0.46 Initial condition (ppm): Carbon tetrachloride (CCT) Acetate (Ca) Nitrate (Cn) Mobile bacteria (XM) Immobile bacteria (XIM) Sorbed CT (mg/mg) (SCT) 0.130 0 42 0 0 2.8E-8 Boundary condition (ppm): Carbon tetrachloride (CCT) Acetate (Ca) Nitrate (Cn) Mobile bacteria (XM) Immobile bacteria (XIM) 0.130 0 42 0 0 Slug injection zone inoculation (ppm): Carbon tetrachloride (CCT) Acetate (Ca) Nitrate (Cn) Mobile bacteria (XM) Immobile bacteria (XIM) 0.1 1650 42 11.8 0 Note: Carbon tetrachloride is abbreviated as CT. 112 Table 4.5: Comparison of the PGA estimated parameters with true solutions along with their bounds for benchmark problem - 4 Parameter True Value Low High PGA Estimate Error % CT reaction rate (day-1) (k`) 0.189 0.1 0.5 0.282 49.21% Nitrate utilization coefficient (day-1) (?) 5.73 1 10 5.89 2.79% Detachment coefficient (day-1) (Kde) 0.043 0.01 0.1 0.063 46.51% Microbial decay rate (day-1)(bKC) 0.221 0.1 1 0.219 0.90% Note: Carbon tetrachloride is abbreviated as CT. Figure 4.5 Comparison of results from PGA estimated parameters and the true parameters for benchmark problem ? 4 4.6.5 Scalability of Parallel GA The computational performance of the PGA was tested on a computer with four processors. All four benchmark problems were run using 1, 2, 3 and 4 processors and the total program run time 113 was calculated using the omp_wall_time() function. This internal clock time function was called at the beginning of the program and subsequently at the end of the program. The difference between these two times gave an estimate of total program run time. The speedup of the parallel program was calculated using the formula: s e q u e n t i a l t i m es p e e d u p = p a r a lle l t i m e (4.19) The sequential time in equation (4.19) was obtained by solving the problems using a single processor. The performance of the PGA for different benchmark problems and its comparison against an idealistic linear speedup function are shown in Figure 4.6. The details of simulation times are summarized in Table 4.6. It can be observed from the figure that the performance of the PGA for all the benchmark problems is close to the ideal linear speedup function, except for the third problem. This is because, the total simulation time taken for Problem-3 was extremely small and hence the calculation of fitness was not as time intensive compared to other benchmark problems. In a problem where the parallelized parts of the problem are not as computationally intensive, as expected, the time taken to perform the non-parallelized tasks become a limiting factor and thereby reduce the total efficiency of the parallel operations. Table 4.6: Simulation times for all the benchmark problems for different number of processors in seconds # Processors BP-1 BP-2 BP-3 BP-4 1 885.00 132.00 7.00 3920.00 2 446.00 66.00 3.92 1976.00 3 339.73 50.00 3.12 1479.00 4 249.63 39.00 2.71 1100.00 114 Figure 4.6 Speedup for all benchmark problems for different number of processors 4.6.6 Sensitivity to Initial Population and Generations The final quality of the solution that GA would find would depend on the size of the initial population generated by the GA and the number of generations. Ideally, a larger number of initial solutions would allow the GA to search a larger solution space, but this also increases the computational burden as the fitness has to be calculated for all the initial solutions. On the other hand, having a smaller initial population would limit the solution space of the GA and this could cause the GA to be trapped in a local minimum. To check the sensitivity of GA to the size of the initial population and the number of generations, we ran all four benchmark problems using four different initial population sizes, 8, 16, 32 and 64; in addition, we also increased the total number of generation to 300. The best solution for each generation was stored to compare the general convergence pattern for different initial population sizes. The results from the sensitivity analysis for all four benchmark problems are shown in Figure 4.7 (a)-(d). It can be observed from the results that the convergence rate was faster when the initial population size was 115 increased; however, increasing the number of generations did not affect the quality of the solution found by the algorithm. Also, in the case of simulations with fewer population sizes, the minimum value reached was away from global minimum value. It is necessary to find an optimal initial population size and this could vary based on the number of parameters being estimated. Increasing the population size might not always result in a better solution as it can be seen that an increase in population size from 32 to 64 did not improve in the solutions. For problems involving high levels of computational complexity, evaluation of another 32 candidates for fitness could increase the overall simulation time drastically without improving the solution. In this study, we found a population size of 32 to be the optimal number in all our simulations. Figure 4.7 Fitness of the best solution for each generation by varying initial population for each benchmark problem 116 4.6.7 Sensitivity to the Bounds of the Unknown Parameters The rate of convergence and the quality of the final solution is also a function of level of uncertainty (characterized by bounds used to define the minimum and maximum values) associated with the unknown parameter. To understand the sensitivity of GA to these bounds, we ran the test Problem 1 multiple times using different bounds. This problem was selected because some parameters in this problem were highly sensitive, and even a minor change would cause considerable fluctuations in the concentration profiles. We developed four different scenarios to test the sensitivity of GA to parameter bounds. In the first scenario, the bounds were kept within an order of magnitude of the true parameter; in the second scenario, the lower bound was reduced by an order of magnitude; in the third scenario, the lower bound is kept the same as the first scenario, but the higher bound was increased by an order of magnitude, and in the fourth scenario both lower and higher bounds were increased by an order of magnitude. Table 4.7 summarizes the lower and higher bounds used in all four the scenarios. The concentration profiles predicted under different scenarios are shown in Figure 4.8. It was observed from the figure that the quality of the solution deteriorated as the bounds for the unknown parameters were increased above an order of magnitude in either direction. The results were far away from the original solution when both upper and lower bounds were increased by an order of magnitude. This shows that for highly sensitive problems one would need reasonable guesses to get meaningful solutions. For problems like this, GA should perhaps be viewed as a polishing algorithm that can refine the solution within a known domain, rather than a true search algorithm that can operate on totally unknown domain. 117 Table 4.7: High and low bounds for the four difference scenarios along with the GA estimated value Parameter Lower bound Higher Bound GA Estimated Value Scenario-1 D 1.00E-02 1.00E-01 8.17E-02 ? 1.00E-02 1.00E-01 2.25E-02 Kd 1.00E-04 1.00E-03 1.37E-04 Scenario-2 D 1.00E-03 1.00E-01 3.85E-02 ? 1.00E-03 1.00E-01 4.44E-02 Kd 1.00E-05 1.00E-03 8.53E-05 Scenario-3 D 1.00E-02 1.00E+00 2.08E-01 ? 1.00E-02 1.00E+00 4.04E-02 Kd 1.00E-04 1.00E-02 1.00E-04 Scenario-4 D 1.00E-03 1.00E+00 8.79E-01 ? 1.00E-03 1.00E+00 6.93E-01 Kd 1.00E-05 1.00E-02 6.10E-05 Figure 4.8 Comparison of different scenarios of the high and low bounds with the true solution for the benchmark problem ? 1 118 The GA used in this study works by systematically searching through a given solution space to find the best solution. Therefore, it is important that the initial population contains good quality solutions for the GA to attain the global optimum. The quality of the initial population is largely dependent on the bounds set by the user at the beginning of the simulation. Most reactive transport problems contain a few highly sensitive parameters and small changes in these parameter values would significantly alter the solution. Therefore, depending on the type of problem, the user might have to restrict the GA within a narrow boundary, if the solution had convergence problems. Future studies explore options to modify the algorithm to automatically identify a set of sensitive parameters and automatically narrow the bounds. 4.7 Conclusions In this study, we have tested the performance of a parallel version of GA by solving four benchmark problems that simulated both batch and reactive-transport scenarios. The PGA algorithm was able to successfully estimate the model parameters for different types of reaction models. In all the cases, the PGA estimated parameters were close to the original parameters; furthermore, the simulation results from these parameters were able to match the original experimental or analytical/numerical data well. The PGA routines were general enough to solve all four benchmark problems without the need for any problem-specific modifications to the routines. One of the limitations, however, is that the user must provide a good initial guess (which should be known within an order of magnitude) for all the model parameters to obtain good quality, convergent final results. In most real cases, the order of magnitude estimates for parameter values can be obtained from literature data and PGA can be used to refine these estimates. The user should, however, be careful not to over-constrain the problem and ensure 119 that the upper and lower parameter bounds are adequately large to generate sufficient number of solutions to avoid convergence to a local optimum. The PGA was optimized to run in parallel mode using OpenMP framework available within Intel FORTRAN v9.0 compiler on a shared memory system. The speedup was quantified for four benchmark problems and the results indicate close to linear speedup for three benchmark problems. The fourth benchmark (designated as Problem 3) was a much simpler problem that required very little computational effort and hence the parallel computing steps did not reduce the overall computational time. These results show that the use of PGA was more appropriate for solving computationally intensive reactive transport problems. The PGA used in this study was successfully demonstrated to run on a standard multi-core desktop Pentium PC platform. The overall computation gain obtained using this hardware was significant. Since most modern desktop PCs are now equipped with multi-core processors, the methods used in this study can be easily adapted to take advantage of these platforms. The proposed optimization framework, which was used for estimating unknown kinetic and transport parameters in our multi-component reactive transport problems, is a generic procedure that can be extended to solve a variety of environmental problems. 120 Chapter 5 SUMMARY AND CONCLUSIONS 5.1 Summary and Conclusion The overall goal of this dissertation was to develop a comprehensive set of tools that can be used to simulate the fate and transport of reactive contaminants in one-dimensional groundwater systems, as well as estimate the model parameters relevant to these systems by solving the inverse problem. A comprehensive one-dimensional model, RT1D, was developed and the capabilities of the tool were demonstrated by simulating a variety of kinetic and geochemistry problems. The mathematical model was then used to simulate a bioaugmentation experiment completed to remediate PCE-DNAPL in a single fracture system. We developed a mathematical framework to simulate the bioaugmentation of PCE-DNAPL in single fracture system. The mathematical framework describes multi-species bioreactive transport processes that include bacterial growth and detachment dynamics, biodegradation of chlorinated species, competitive inhibition of various reactive species, and the loss of daughter products due to back-partitioning effects. The kinetic Monod parameters evaluated from batch experiments were scaled to estimate the parameters for the fracture system using a trial and error method. We were able to verify and validate these parameters by simulating two different experimental datasets conducted using high and low flow velocities. The model was able to predict the data well. During the calibration process, we realized the limitations of trial-and-error methods and the need for a flexible parameter identification tool for assisting in the calibration process. Therefore, in the final phase 121 of this study, we explored the use of parallel genetic algorithms (PGA) for automatic parameter estimation in reactive transport models. We developed a PGA to estimate the kinetic parameters in a reactive transport model for a given experimental dataset. We demonstrated the flexibility and usefulness of the PGA code by solving four different benchmark problems that have published model results or analytical solutions. The PGA was able to estimate the parameters that were close to the true parameters. In some cases, the percentage error between the estimated parameters and the true parameters were high (over 400% in the case of benchmark problem - 3). It should be noted that the goal of the PGA was to minimize the sum square of errors between the given experimental data and the simulated concentrations from the PGA estimated parameters. The PGA was able to accomplish this goal as the predicted concentrations from the PGA estimated parameters and the true parameters are very close; however, in some cases, the error percentages are high. This suggests that the underlying parameter estimation problem is non-unique. As the number of unknown parameters increases, the non-uniqueness of the problem could become a major issue. In such cases, the PGA estimated parameters have to be validated using a different experimental dataset to ensure that the PGA estimated parameters are unique to that dataset. It was also observed that the GAs are computationally intensive search algorithms; however, we made it computationally efficient by running the GA in a parallel mode using the shared memory parallel computing platform and the OpenMP FORTRAN complier. The efficiency of the parallel code showed close to linear speedup for a desktop computer with 4 processors. 5.2 Recommendations for future work The modeling tool developed as a part of the first objective is an advanced numerical tool equipped with a variety of solvers. Further improvements can be made to the tool by either 122 adding better solvers or scaling it for multi-dimensional framework. In order to improve the solvers, accurate advection solvers that are being developed in the fluid dynamics literature can be incorporated and the ODE solver could also be upgraded to provide the user additional options for solving stiff kinetic problems. The EXCEL-VB modeling tool can also be further extended for solving multi-dimensional problems. The necessary matrix and ODE solvers can be written in EXCEL-VB and could be integrated within a multi-dimensional modeling framework. However, the user-friendliness of the EXCEL-VB platform would always be undermined by the lack of parallel processing capabilities. Also, the overall speed of EXCEL- VB could be slightly slower compared than other computer languages such as FORTRAN and C. Therefore, careful trade-off assessments should be made before developing multi-dimensional tools within the EXCEL framework. The mathematical modeling framework developed as a part of the second objective for the bioaugmentation of PCE-DNAPL would be more robust if a dissolution model for DNAPL could be included. The dissolution of DNAPL in fractures is fraught with statistical variations and a dissolution model capable of accounting for these stochastic variations would be a significant addition to the existing mathematical framework. The PGA developed as a part of the third objective is a generic genetic algorithm. Although, this algorithm is capable of estimating variety of parameters without modifying the algorithm, it is dependent on the quality of initial guesses to provide a good solution space for its convergence. Several strategies can be employed to improve the algorithm to overcome the limitations caused by bad initial guess values. A preliminary analysis could be performed in the algorithm to identify the sensitive parameters and their sensitivity in a particular range and the GA could limit the variation of these parameters to this specific range. Another strategy would 123 be to use a multi-objective genetic algorithm (MOGA) by using the niche Pareto curve to find the ideal solution space (Singh et al., 2008). In this strategy, the MOGA can identify the ideal solution space by optimizing the parameters for each component and this solution space can be used to identify the global minimum for that problem. Genetic algorithms are a vastly improving modern field and there are several ways to apply some of these modified algorithms to improve the parameter estimation schemes. 124 REFERENCES Abramson, D., Abela, J., 1991. A parallel genetic algorithm for solving the school timetabling problem. Division of Information Technology, CSIRO. Abramson, D., Mills, G., Perkins, S., 1994. Parallelisation of a genetic algorithm for the computation of efficient train schedules. Parallel Computing and Transputers 37 139-149. Allison, J.D., Brown, D.S., Novo-Gradac, K.J., 1990. MINTEQA2/PRODEFA2, A geochemical assessment model for environmental systems: Version 3.0. Environmental Research Laboratory, Office of Research and Development, U.S Environmental Protection Agency, Athens, Georgia, p. 106. Amos, B.K., Suchomel, E.J., Pennell, K.D., Loffler, F.E., 2009. Spatial and Temporal Distributions of Geobacter lovleyi and Dehalococcoides spp. during Bioenhanced PCE-NAPL Dissolution. Environmental Science and Technology 43(6) 1977-1985. Atlas, R.M., Cerniglia, C.E., 1995. Bioremediation of petroleum pollutants: diversity and environmental aspects of hydrocarbon biodegradation. Bioscience 45(5) 332-338. ATSDR, March 3, 2011, Toxic Substances Portal - Vinyl Chloride http://www.atsdr.cdc.gov/mmg/mmg.asp?id=278&tid=51, Last accessed on (January 13, 2012) Aziz, C.E., 2000. BIOCHLOR: Natural Attenuation Decision Support System: User's Manual Version 1.0. National Risk Management Research Laboratory, Office of Research and Development, US Environmental Protection Agency. 125 Babbar, M., Minsker, B.S., 2006. Groundwater Remediation Design Using Multiscale Genetic Algorithms Journal of Water Resources Planning and Management 132(5) 341-350. Baginska, B., Milne-Home, W., Cornish, P., 2003. Modelling nutrient transport in Currency Creek, NSW with AnnAGNPS and PEST. Environmental Modelling & Software 18(8) 801-808. Baluja, S., 1992. A massively distributed parallel genetic algorithm. DTIC Document. Bauer, P., Attinger, S., Kinzelbach, W., 2001. Transport of a decay chain in homogenous porous media: analytical solutions. Journal of Contaminant Hydrology 49(3) 217-239. Beeman, R.E., Bleckmann, C.A., 2002. Sequential anaerobic-aerobic treatment of an aquifer contaminated by halogenated organics: field results. Journal of Contaminant Hydrology 57(3-4) 147-159. B?ranger, S.C., Sleep, B.E., Lollar, B.S., Monteagudo, F.P., 2005. Transport, biodegradation and isotopic fractionation of chlorinated ethenes: Modeling and parameter estimation methods. Advances in Water Resources 28(1) 87-98. Bradley, P., 2003. History and ecology of chloroethene biodegradation: a review. Bioremediation Journal 7(2) 81-109. Brusseau, M.L., S. K. Sandrin, L. Li, I. Yolcubal, F. L. Jordan, R. M. Maier 2006. Biodegradation during contaminant transport in porous media: 8. The influence of microbial system variability on transport behavior and parameter determination Water Resources Research 42. 126 Canada, E., 2011, Groundwater, http://www.ec.gc.ca/eau- water/default.asp?lang=En&n=300688DC-1, Last accessed on (September 30th, 2011) Cant?-Paz, E., 1998. A survey of parallel genetic algorithms. Calculateurs paralleles, reseaux et systems repartis 10(2) 141-171. Cederberg, G.A., Street, R.L., Leckie, J.O., 1985. A groundwater mass transport and equilibrium chemistry model for multicomponent systems. Water Resources Research 21(8) 1095-1104. Chapra, S.C., Canale, R.P., 1998. Numerical Methods for Engineers with Programming and Software Applications, Third ed. WCB/McGraw-Hill. Charlton, S.R., Parkhurst, D.L., 2002. PHREEQCI: A graphical user interface to the geochemical model PHREEQC. US Department of the Interior, US Geological Survey. Cleanup-Information, April 15, 2011, Dense Nonaqueous Phase Liquids (DNAPLs). Chemistry and Behavior. cis 1,2-Dichloroethene, http://www.clu- in.org/contaminantfocus/default.focus/sec/Dense_Nonaqueous_Phase_Liquids_%28DNAPLs%2 9/cat/Chemistry_and_Behavior/p/3/n/2, Last accessed on (January 13, 2012) Clement, T.P., 2001. Generalized solution to multispecies transport equations coupled with a first-order reaction network. Water Resources Research 37(1) 157-163. Clement, T.P., 2011. Bioremediation of contaminated groundwater systems, In: Aral, M.M., Taylor, S. (Eds.), Groundwater Quantity and Quality Management. American Society of Civil Engineers (ASCE), pp. 522-559. 127 Clement, T.P., Gautam, T.R., Lee, K.K., Truex, M.J., Davis, G.B., 2004. Modeling coupled NAPL-dissolution and rate-limited sorption reactions in biologically active porous media. Bioremediation Journal 8(1-2) 47-64. Clement, T.P., Hooker, B.S., Skeen, R.S., 1996. Numerical modeling of biologically reactive transport near nutrient injection well. Journal of Environmental Engineering 122 833-839. Clement, T.P., Johnson, C.D., Sun, Y., Klecka, G.M., Bartlett, C., 2000. Natural attenuation of chlorinated ethene compounds: model development and field-scale application at the Dover site. Journal of Contaminant Hydrology 42(2) 113-140. Clement, T.P., Peyton, B.M., Skeen, R.S., Jennings, D.A., Petersen, J.N., 1997. Microbial growth and transport in porous media under denitrification conditions: experiments and simulations. Journal of Contaminant Hydrology 24(3-4) 269-285. Clement, T.P., Sun, Y., Hooker, B.S., Petersen, J.N., 1998. Modeling multispecies reactive transport in ground water. Ground Water Monitoring & Remediation 18(2) 79-92. Coleman, N.V., Mattes, T.E., Gossett, J.M., Spain, J.C., 2002. Biodegradation of cis- Dichloroethene as the Sole Carbon Source by a {beta}-Proteobacterium. Applied and environmental microbiology 68(6) 2726. Cupples, A., Spormann, A., McCarty, P., 2004. Vinyl chloride and cis-dichloroethene dechlorination kinetics and microorganism growth under substrate limiting conditions. Environmental Science and Technology 38(4) 1102-1107. 128 de Blanc, P., McKinney, D., Speitel, G., Delshad, M., Pope, G.A., Sepehrnoori, K., 1996. UTCHEM biodegradation model description and capabilities. Center for Petroleum and Geosystems Engineering 1-17. Dickson, S.E., Thomson, N.R., 2003. Dissolution of entrapped DNAPLs in variable aperture fractures: Experimental data and empirical model. Environmental Science and Technology 37(18) 4128--4137. Doherty, J.E., Hunt, R.J., 2010. Approaches to highly parameterized inversion: a guide to using PEST for groundwater-model calibration. US Department of the Interior, US Geological Survey. El Fantroussi, S., Naveau, H., Agathos, S., 1998. Anaerobic dechlorinating bacteria. Biotechnology progress 14(2) 167-188. El Harrouni, K., Ouazar, D., Walters, G.A., Cheng, A.H.D., 1996. Groundwater optimization and parameter estimation by genetic algorithm and dual reciprocity boundary element method. Engineering Analysis with Boundary Elements 18(4) 287-296. Engesgaard, P., Kipp, K.L., 1992. A geochemical transport model for redox-controlled movement of mineral fronts in groundwater flow systems: A case of nitrate removal by oxidation of pyrite. Water Resources Research 28(10) 2829-2843. Farthing, M.W., Miller, C., 2001. A comparison of high-resolution, finite-volume, adapative- stencil schemes for simulating advective-dispersive transport. Advances in Water Resources 24(1) 29-48. 129 Fogarty, T., Huang, R., 1991. Implementing the genetic algorithm on transputer based parallel processing systems. Parallel Problem Solving from Nature 145-149. Fredrickson, N.R., Afsahi, A., Qian, Y., 2003. Performance characteristics of OpenMP constructs, and application benchmarks on a large symmetric multiprocessor. ACM, pp. 140-149. Gaffney, J., Pearce, C., Green, D., 2010. Binary versus real coding for genetic algorithms: A false dichotomy? ANZIAM Journal 51 C347--C359. Giacobbo, F., Marseguerra, M., Zio, E., 2002. Solving the inverse problem of parameter estimation by genetic algorithms: the case of a groundwater contaminant transport model. Annals of Nuclear Energy 29(8) 967-981. Gramling, C.M., Harvey, C.F., Meigs, L.C., 2002. Reactive transport in porous media: A comparison of model prediction with laboratory visualization. Environmental science & technology 36(11) 2508-2514. He, K., Zheng, L., Dong, S., Tang, L., Wu, J., Zheng, C., 2007. PGO: A parallel computing platform for global optimization based on genetic algorithm. Computers & geosciences 33(3) 357-366. Holland, J.H., 1975. Adaptation in natural and artificial systems. University of Michigan press. Inchem, 2011, SIDS Initial assessment profile, http://www.inchem.org/documents/sids/sids/74851.pdf, Last accessed on (January 13, 2012) Information, C., May 18, 2011, Trichloroethylene (TCE). Chemistry and Behavior, http://www.clu- 130 in.org/contaminantfocus/default.focus/sec/Trichloroethylene_%28TCE%29/cat/Chemistry_and_ Behavior/, Last accessed on (January 13, 2012) Islam, J., Singhal, N., 2002. A one-dimensional reactive multi-component landfill leachate transport model. Environmental Modelling & Software 17(6) 531-543. Jeong-Hun Park, Xianda Zhao, Voice, T.C., 2001. Comparison of Biodegradation Kinetic Parameters for Naphthalene in Batch and Sand Column Systems by Pseudomonas Putida. Environmental Progress 20(2) 93-102. Jeppu, G.P., Clement, T.P., Barnett, M.O., Lee, K.-K., 2012. A modified batch reactor system to study equilibrium-reactive transport problems. Journal of Contaminant Hydrology 129-130(0) 2- 9. Kipp, K.L., 1987. HST3D: A computer code for simulation of heat and solute transport in three- dimensional ground-water flow systems. US Geological Survey. Koza, J., 1992. Genetic programming. MIT Press, Cambrige, MA. Kumar, A., 2006. Coupling Transoport Codes with Geochemical Models, Department of Civil Engineering. Auburn University: Auburn, AL, p. 93. Lee, I., Bae, J., Yang, Y., McCarty, P., 2004. Simulated and experimental evaluation of factors affecting the rate and extent of reductive dehalogenation of chloroethenes with glucose. Journal of Contaminant Hydrology 74(1-4) 313-331. Lee, S., Heber, A.J., 2010. Ethylene removal using biotrickling filters: part II. Parameter estimation and mathematical simulation. Chemical Engineering Journal 158(2) 89-99. 131 Lu, G., Clement, T.P., Zheng, C., Wiedemeier, T.H., 1999. Natural attenuation of BTEX compounds: Model development and field-scale application. Ground Water 37(5) 707-717. Madsen, K.M., Perry, A.E., 2010. Using Genetic Algorithms on Groundwater Modeling Problems in a Consulting Setting, Proceedings of the Annual International Conference on Soils, Sediments, Water and Energy, p. 11. Majdalani, S., Fahs, M., Carrayrou, J., Ackerer, P., 2009. Reactive transport parameter estimation: Genetic algorithm vs. Monte carlo approach. American Institute of Chemical Engineers 55(8) 1959-1968. Massoudieh, A., Mathew, A., Ginn, T.R., 2008. Column and batch reactive transport experiment parameter estimation using a genetic algorithm. Computers & geosciences 34(1) 24-34. Maymo-Gatell, X., Chien, Y., Gossett, J., Zinder, S., 1997. Isolation of a bacterium that reductively dechlorinates tetrachloroethene to ethene. Science 276(5318) 1568. McKinney, D.C., Lin, M.D., 1994. Genetic algorithm solution of groundwater management models. Water Resources Research 30(6) 1897-1906. Michalewicz, Z., 1996. Genetic algorithms+ data structures = Evolution Programs, Third ed. Springer. Miller, C., Benson, L., 1983. Simulation of solute transport in a chemically reactive heterogeneous system: Model development and application. Water Resources Research 19(2) 381-391. 132 Moreau, M., 1987. Some European Perspectives on Prevention of Leaks from Underground Petroleum Storage Systems. Ground Water Monitoring & Remediation 7(1) 45-48. Mulligan, A.E., Brown, L.C., 1998. Genetic algorithms for calibrating water quality models. Journal of environmental engineering 124(3) 202-211. Parkhurst, D.L., 2004. PHAST--a Program for Simulating Ground-water Flow, Solute Transport, and Multicomponent Geochemical Reactions. US Department of the Interior, US Geological Survey. Parkhurst, D.L., Appelo, C.A.J., 1999. User's guide to PHREEQC (version 2) - A computer program for speciation, batch-reaction, one-dimensional transport, reaction path, and inverse geochemical calculations,. U.S Geological Survey Water-Resources Investigations Report, p. 312 pp. Perlman, H., 2011, Groundwater Use, The USGS Water Science School, http://ga.water.usgs.gov/edu/wugw.html, Last accessed on (May 9th, 2012) Peyton, B., Skeen, R., Hooker, B., Lundman, R., Cunningham, A., 1995. Evaluation of bacterial detachment rates in porous media. Applied Biochemistry and Biotechnology 51(1) 785-797. Phanikumar, M.S., Hyndman, D.W., Wiggert, D.C., Dybas, M.J., Witt, M.E., Criddle, C.S., 2002. Simulation of microbial transport and carbon tetrachloride biodegradation in intermittently-fed aquifer columns. Water Resources Research 38(4) 1-13. 133 Prommer, H., Barry, D.A., Davis, G.B., 1998. A one-dimensional reactive multi-component transport model for biodegradation of petroleum hydrocarbons in groundwater. Environmental Modelling and Software 14(2-3) 213-223. Quezada, C.R., Clement, T.P., Lee, K.K., 2004. Generalized solution to multi-dimensional multi- species transport equations coupled with a first-order reaction network involving distinct retardation factors. Advances in Water Resources 27(5) 507-520. Rafai, H.S., Bedient, P., Borden, R.C., Haasbeek, J.F., 1987. BIOPLUME II- Computer model of two-dimensional contaminant transport under the influence of oxygen limited biodegradation in ground water. Robert S. Kerr Environmental Research Laboratory: Ada, Oklahoma, pp. 17-19. Rafai, H.S., Newell, C., Gonzales, J., Dendrou, S., Dendrou, B., 1998. BIOPLUME III: Natural Attenuation Decision Support System, User's Manual Version 1. 0. National Risk Management Research Laboratory, Office of Research and Development, U.S Environmental Protection Agency: Cincinnati, OH. Ramsburg, C.A., Thornton, C.E., Christ, J.A., 2010. Degradation Product Partitioning in Source Zones Containing Chlorinated Ethene Dense Non-Aqueous-Phase Liquid. Environmental Science and Technology 44(23) 9105-9111. Reed, P., Minsker, B., Goldberg, D.E., 2000. Designing a competent simple genetic algorithm for search and optimization. Water Resources Research 36(12) 3757-3761. Rifai, H., Newell, C., Gonzales, J., Dendrou, S., Kennedy, L., Wilson, J., 1998. BIOPLUME III - Natural attentuation decision support system, User's Manual v1.0. National Risk Management Research Laboratory: Cincinnati, Ohio. 134 Rolle, M., Clement, T.P., Sethi, R., Di Molfetta, A., 2008. A kinetic approach for simulating redox-controlled fringe and core biodegradation processes in groundwater: model development and application to a landfill site in Piedmont, Italy. Hydrological Processes 22(25) 4905-4921. Sarma, K.C., Adeli, H., 2002. Bilevel parallel genetic algorithms for optimization of large steel structures. Computer?Aided Civil and Infrastructure Engineering 16(5) 295-304. Schaefer, C.E., Callaghan, A.V., King, J.D., McCray, J.E., 2009a. Dense Nonaqueous Phase Liquid Architecture and Dissolution in Discretely Fractured Sandstone Blocks. Environmental Science and Technology 43(6) 1877--1883. Schaefer, C.E., Condee, C.W., Vainberg, S., Steffan, R.J., 2009b. Bioaugmentation for chlorinated ethenes using Dehalococcoides sp.: Comparison between batch and column experiments. Chemosphere 75(2) 141-148. Schaefer, C.E., Lippincott, D.R., Steffan, R.J., 2010a. Field-Scale Evaluation of Bioaugmentation Dosage for Treating Chlorinated Ethenes. Ground Water Monitoring & Remediation 30(3) 113-124. Schaefer, C.E., Towne, R.M., Vainberg, S., McCray, J.E., Steffan, R.J., 2010b. Bioaugmentation for treatment of dense non-aqueous phase liquid in fractured sandstone blocks. Environmental science & technology 44(13) 4958-4964. Singh, A., Minsker, B., Takagi, H., 2005. Interactive Genetic Algorithms for Inverse Groundwater Modeling: Issues with Human Fatigue and Prediction Models. ASCE. 135 Singh, A., Minsker, B.S., Valocchi, A.J., 2008. An interactive multi-objective optimization framework for groundwater inverse modeling. Advances in Water Resources 31(10) 1269-1283. Sinha, E., Minsker, B.S., 2007. Multiscale island injection genetic algorithms for groundwater remediation. Advances in Water Resources 30(9) 1933-1942. Smidt, H., de Vos, W., 2004. Anaerobic microbial dehalogenation. Annual Review of Microbiology 58 43-73. Srinivasan, V., Clement, T., 2008. Analytical solutions for sequentially coupled one-dimensional reactive transport problems-Part I: Mathematical derivations. Advances in Water Resources 31(2) 203-218. Srinivasan, V., Clement, T., Lee, K., 2007. Domenico Solution?Is It Valid? Ground Water 45(2) 136-146. Sun, Y., Clement, T.P., 1999. A Decomposition Method for Solving Coupled Multi?Species Reactive Transport Problems. Transport in porous media 37(3) 327-346. Sun, Y., Petersen, J., Clement, T., 1999. Analytical solutions for multiple species reactive transport in multiple dimensions. Journal of Contaminant Hydrology 35(4) 429-440. Tanese, R., 1989. Distributed genetic algorithms for function optimization, Proceedings of the Third International Conference on Genetic Algorithms and their Applications: San Mateo, CA, pp. 434-439. Todd, D.K., 1980. Groundwater Hydrology, Second Edition ed. John Wiley & Sons. 136 Toride, N., Leij, F.J., Genuchten, M.T., 1995. The CXTFIT Code for Estimating Transport Parameters from Laboratory or Field Tracer Experiments (ver 2.1). US Salinity Laboratory, Agricultural Research Service, US Department of Agriculture: Riverside, CA. Toride, N., Leij, F.J., Van Genuchten, M.T., 1993. A comprehensive set of analytical solutions for nonequilibrium solute transport with first-order decay and zero-order production. Water Resources Research 29(7) 2167-2182. Torlapati, J., Clement, P.T., 2012a, Assessment of Parallel Genetic Algorithms for Calibrating One-Dimensional Multi-species Reactive Transport Models, Environmental Modelling and Software, Submitted Manscript Torlapati, J., Clement, T.P., 2012b, Benchmarking a Visual-Basic based Multi-Component One- Dimensional Reactive Transport Modeling Tool, Computers & geosciences, 10.1016/j.cageo.2012.08.009 Torlapati, J., Clement, T.P., Schaefer, C.E., Lee, K.-K., 2012. Modeling Dehalococcoides sp. Augmented Bioremediation in a Single Fracture System. Ground Water Monitoring & Remediation 32(3) 75-83. Truesdell, A.H., Jones, B.F., 1974. WATEQ, A computer program for calculating chemical equilibria of natural water. U.S. Geological Survey, Journal of Research, pp. 233-248. Tsai, F.T.C., Katiyar, V., Toy, D., Goff, R.A., 2009. Conjunctive management of large-scale pressurized water distribution and groundwater systems in semi-arid area with parallel genetic algorithm. Water resources management 23(8) 1497-1517. 137 Vainberg, S., Condee, C.W., Steffan, R.J., 2009. Large-scale production of bacterial consortia for remediation of chlorinated solvent-contaminated groundwater. Journal of Industrial Microbiology and Biotechnology 36(9) 1189-1197. Valocchi, A.J., Werth, C.J., 2004. Web-based interactive simulation of groundwater pollutant fate and transport. Computer Applications in Engineering Education 12(2) 75-83. van Genuchten, M.T., Alves, W., 1982. Analytical solutions of the one-dimensional convective- dispersive solute transport equation. Technical Bulletin(1661). Wagner, B.J., 1992. Simultaneous parameter estimation and contaminant source characterization for coupled groundwater flow and contaminant transport modelling. Journal of Hydrology 135(1- 4) 275-303. Walter, A., Frind, E., Blowes, D., Ptacek, C., Molson, J., 1994. Modeling of multicomponent reactive transport in groundwater: 1. Model development and evaluation. Water Resources Research 30(11) 3137-3148. Wang, M., Zheng, C., 1997. Optimal remediation policy selection under general conditions. Ground Water 35(5) 757-764. Wang, Q., 1997. Using genetic algorithms to optimise model parameters. Environmental Modelling & Software 12(1) 27-34. Westall, J.C., 1976. MINEQL: A computer program for the calculation of the chemical equilibrium composition of aqueous system. Department of Civil Engineering, Massachussets Insitute of Technology: Cambridge, p. 91. 138 Westall, J.C., 1979a. MICROQL. I. A Chemical Equilibrium Program in BASIC. Swiss Federal Institute of Technology, EAWAG: Dubendorf, Switzerland. Westall, J.C., 1979b. MICROQL. II. Computation of Adsorption Equilibria in BASIC. Swiss Federal Institute of Technology: Dubendorf, Switzerland. Yabusaki, S.B., Fang, Y., Long, P.E., Resch, C.T., Peacock, A.D., Komlos, J., Jaffe, P.R., Morrison, S.J., Dayvault, R.D., White, D.C., 2007. Uranium removal from groundwater via in situ biostimulation: Field-scale modeling of transport and biological processes. Journal of Contaminant Hydrology 93(1) 216-235. Yu, S., Dolan, M., Semprini, L., 2005. Kinetics and inhibition of reductive dechlorination of chlorinated ethylenes by two different mixed cultures. Environmental Science and Technology 39(1) 195-205. Yu, S., Semprini, L., 2004. Kinetics and modeling of reductive dechlorination at high PCE and TCE concentrations. Biotechnology and bioengineering 88(4) 451-464. Zheng, C., Wang, P.P., 1999. MT3DMS, A modular three-dimensional multispecies transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; documentation and user?s guide. U.S. Army Engineer Research and Development Center Contract Report SERDP-99-1: Vicsburg, MS, p. 202. Zysset, A., Stauffer, F., Dracos, T., 1994. Modeling of reactive groundwater transport governed by biodegradation. Water Resources Research 30(8) 2423-2434. 139 APPENDIX-A1 RT1D PROGRAM DETAILS A1.1 Understanding the spreadsheet In this section, we will briefly discuss the different sections of the spreadsheet and illustrate their use to set up various problems. As explained in the earlier section, RT1D has four different simulation options. The parameters for the transport module can be set in Section-1 as shown in Figure A1.1. These parameters are used by the model when a simulation option involving transport module is selection (2 & 4). The left side of Section-2 is dedicated for kinetic type problem whereas the right side is dedicated for the geochemistry equilibrium type problem. Depending on the type of problem, only one of the sections is used. The program automatically sets the unused type of problem to ?N/A?. The Section-3 in Figure A1.1 is a ?Generate Input Template? button that will automatically generate a spreadsheet input table depending on the simulation option and the reaction-type parameters set in the Section-2. It can be observed from this figure that the input for kinetic-type reaction has been set to ?N/A? by the program after pressing the button in Section-3. An example for a kinetic-type input sheet is shown in Figure A1.2 and a geochemistry-type input sheet is shown in Figure A1.3 respectively. 140 Figure A1.1. Spreadsheet layout for RT1D parameters to generate reaction-specific input template In the following section, we will briefly discuss the various steps involved in running a simulation using RT1D. The parameters and the required short codes for each module and options is made available in a textbox in the spreadsheet for the convenience of the user. Step 1: Transport module - To perform simulations that involve the use of a transport module, the advection-dispersion parameters along with the length of the column need to be input in their respectively columns in Section-1. The total simulation time and the amount of time for which the concentration pulse has been supplied should also be provided by the user. For a continuous pulse, the pulse time should be equal to the total simulation time. The user should also provide the type of advection-dispersion solver to be used for the transport module. The simulation option (SO) describes the different types of problems that RT1D can solve. The options are: 141 batch kinetics (1), kinetic reactive transport (2), geochemistry (3), and geochemistry coupled with transport (4) Step 2: Simulation option parameters ? Depending on the simulation option chosen in the Section-1, the kinetic or the geochemistry parameters should be input in Section-2. For a kinetic- type problem, the user should input number of mobile and immobile components, the reaction reaction package number and type of ODE solver module used to solve this reaction package. The user can also set the total number of user-defined reaction parameters that will be used with the reaction package. For a geochemistry problem, the user should input the number of components, species, components whose concentrations are fixed (example: fixed pH), number of aqueous component concentrations that need to be tracked, number of sorbed concentrations that are required for correcting the aqueous concentrations during transport and the type of surface complexation model. Step 3: Generate the input template ? The ?Generate Input Template? button shown in Section-3 of Figure A1.1 will generate a table to enter the parameters specific to the kind of simulation option selected by the user. Example input template for a kinetic-type and geochemistry equilibrium-type problem are shown in Figures A1.2 and A1.3 respectively. The Section-4 shown in these figures changes depending on the parameters set in Section-2 of the spreadsheet. For the simulation options 1 and 2 (kinetic), the button generates a template to input retardation factor, initial and boundary conditions based on the number of mobile and immobile species entered in Section-2. The input template also provides an area for the user to input reaction parameters specific to that reaction package. 142 For the simulation options 3 and 4 (geochemistry), the button (3) will generate an input layout similar to Figure A1.3 to enter the tableau, surface complexation parameters, initial and boundary conditions for the aqueous component concentrations, initial condition for the sorbed concentrations, the composition of the sorbed concentrations and the surface complexation parameters. Figure A1.2. Example input template for a kinetic-type reactive transport 143 Figure A1.3. Example input template for a geochemistry equilibrium coupled with transport Step 4: Solve ? Prior to running a simulation, the user should ensure that the following list of things is completed a) Kinetic problem - Simulation options (1 & 2) i. Pick an existing reaction package from the several built-in packages using the short code or go to step (b) ii. Program your own kinetic reaction package by opening the code editor. 144 iii. Input the reaction parameters as necessary. iv. Click the solve button shown in Section-5 to perform the simulation b) Geochemistry equilibrium problem ? Simulation options (3 & 4) i. Input the total component concentrations at the boundary and the initial concentrations at the nodes, if any. ii. Input the guess concentrations, this is used by the solver for the starting solution to converge to the correct solution. A good starting point would be the total concentration. Please do not input ?0? for the guess values as the logarithm is calculated for these guess concentrations as a part of the solution process. iii. Input the tableau and the surface complexation parameters, if any. iv. Input the initial sorbed concentration present in the column and the species that combine to form the sorbed concentration. This sorbed concentration is used to calculate the aqueous species concentrations at each time step. v. Click the solve button shown in Section-5 to perform the simulation Step 5: Viewing the solutions: The results are presented in Sheets 2 & 3 of the model spreadsheet. Sheet 2 provides the spatial variation of the component concentrations after the completion of the simulation time. Sheet 3 provides the breakthrough component concentrations at the outlet after each time step. By default, the breakthrough component concentrations are printed at the last node (end of the column). The program can be modified to print the breakthrough component concentrations at any node in the column. WARNING!! - The spreadsheet erases the data below the buttons (after cells A27) when either button 3 or 6 are pressed. The output data in Sheets 2 & 3 is cleared every time the button 5 is 145 pressed. The user is advised to use a different a spreadsheet for post-processing of the results and other intermediate calculations. A1.2 Input spreadsheet label details A1.2.1 Transport module i. Length (L) ? Length of the column for the transport module ii. Total time (T) ? Total simulation time iii. Pulse time (T) ? Total time for which the concentration is pumped through the column during the experiment at the boundary. This is less than or equal to the Total time set in (ii). iv. delx (L) ? Grid size. This is used to calculate the total number of nodes using the formula (length/delx)+1. v. delt (T) ? Time step. The total number of iterations is calculated using the formula (total time/delt). vi. Velocity (LT-1) ? This is the pore velocity of the liquid flowing through the column. This is calculated from the flowrate (Q) as follows (Q/cross-sectional area)/porosity. Please make sure that units are same as the length/total time vii. Dispersion coefficient (LT-2) ? This is the hydrodynamic dispersion coefficient calculated by multiplying the value of dispersivity with the velocity. Please provide the hydrodynamic dispersion coefficient value and not the dispersivity value. The code does not multiply with the velocity. viii. Adv-Disp type ? This takes a short code for the type of Advection-Dispersion solver. The choices are 0, 1, and 2 for explicit advection method, fully implicit 146 advection-dispersion or TVD advection respectively. For choices 0 and 2, dispersion is solved by fully implicit dispersion module. ix. Simulation option ? Type of problem 1 - Batch Kinetics 2 - Reactive Transport 3 - Batch Geochemistry 4 - Geochemistry coupled with Transport A1.2.2 Kinetic reaction module i. # Mobile components ? This is the number of aqueous components inside the system. The mobile components will have an advection/dispersion and reaction term. There should be at least 1 mobile species for the program to work. ii. # Immobile components ? This is the number of solid phase components in the system that do not undergo advection-dispersion. They only undergo reaction. The immobile components will be placed at the bottom after all the mobile components. In the reaction package, the immobile species number will start after the mobile species. That is, if there are 6 mobile species and 2 immobile species, we will count the immobile species as 7 & 8. Please look at the example problem for a detailed explanation. iii. Reaction package # ? This is a short code for the type of reaction kinetics. The different kinds of reaction kinetics that have been pre-programmed into the RT1D model. We are adding more packages as we continue developing the model. The user can use their reaction kinetics by picking the option 10. 1 ? First order sequential degradation 2 ? Four-component coupled sequential degradation 147 3 ? Four-component decay chain (Bauer et al., 2001) 4 ? Modified Monod kinetics of TCE bioaugmentation 5 ? Rate-limited sorption (Benchmark problem ? 1 Torlapati et al., 2012) 6 ? Denitrification (Benchmark problem ? 2 Torlapati et al., 2012) 7 ? Biodegradation of Carbon Tetrachloride (Benchmark problem ? 3 Torlapati et al., 2012) 8 ? Open 9 ? Open 10 ? User-defined reaction package iv. ODE Solver type ? Short code for the type of reaction solver 0 ? Adaptive Runge-Kutta-Fehlberg solver 1 ? Fourth order Runge-Kutta solver v. # Reaction parameters ? Set the number of user-defined reaction terms needed A1.2.3 Geochemistry equilibrium module These variables come into the picture only when you?re using the geochemistry module of the code. This is when you enter either 3 or 4 short codes for the reaction type i. # Components ? Number of components ii. # Species ? Number of species iii. # Fixed component concentrations ? set the number of components whose concentrations are fixed (Example: fixed pH) iv. # Aqueous components ? number of aqueous components that we are tracking v. # Sorbed concentrations ? the number of solid phase concentrations that are necessary for correcting the aqueous phase component concentrations 148 vi. SCM TYPE ? Type of surface complexation method. 0 ? No surface complexation 1 ? Constant capacitance 2 ? Diffuse layer 3 ? Stern Layer 4 ? Triple Layer 5 ? Generalized Two Layer Modem (Dzombak & Morel, 1990) A1.2.4 Kinetic reaction parameters These labels are generated in the Section-4 of Figure 3 when either simulation option 1 or 2 are chosen. The spreadsheet displays four different columns requiring the input for the following parameters: i. Component name: The input template automatically populates it with a default component name. This can be changed to a more appropriate name by the user. The program reads this name and uses for output in Sheets 2 and 3. ii. R (Retardation factor): In case of linear sorption, we have a retardation factor. This retardation factor is given by d ?K1 where ? is bulk density of the soil (mg/L), Kd is the linear sorption constant (L/mg) and ? is the porosity. This is equal to 1 when there is no sorption. Cannot be less than 1. iii. Initial: This is the initial concentration of the species in the column. It is possible that there is a residual component concentration present in the column before the simulation time has begun. This option sets a constant initial concentration, as specified in the 149 spreadsheet, across the all the nodes in the column. It is possible to set a variable initial concentration by modifying the setinit() subroutine in the code. iv. Boundary: This is the inlet concentration of the species at time = 0. This boundary concentration will be supplied to the column until the end of the pulse time. v. Reaction terms: Generates a table for entering the reaction terms and labels based on the total reaction terms. A1.2.5 Geochemistry equilibrium parameters (without transport) The input is different geochemistry package with advection and without advection. The input template also changes based on the kind of surface complexation reaction chosen in the geochemistry input sheet. i. Mobile species: This presents with 4 columns of different parameters for the mobile species. a. Index # is the serial number of the component name. This is automatically populated and should not be changed b. Comp. Name: A default component name is generated automatically. It is advised that this component name be changed to something more suitable. The components with fixed concentration are input after all the variable components have been entered into the spreadsheet cells. c. Total concentration: If the total concentration of the component is known, please enter the value here. d. Guess concentration: If the values are unknown, please enter a guess value so the program has a starting value for the iteration process. The iterative process converges faster if suitable starting concentrations are chosen. A good guess for 150 the guess concentration is the initial concentration itself unless it is zero. Since the logarithm of the guess concentrations are calculated during the solution procedure, please do not use 0?s as guess concentrations. If the solution does not converge, try a different guess value or check the problem inputs. ii. Tableau: This is the tableau where you fill in the stoichiometric matrix. Make sure the components are in the same order as the Total concentrations in the column. Also enter the Log K values in the end for each species for the mass action matrix. iii. Additional parameters are required based on the surface complexation model chosen. a. LSIG0: Index for PSI0 in the component list [L0 in Westall (1979)] b. LSIG1: Index for PSI1 in the component list [L1 in Westall (1979)] c. LSIG2: Index for PSI2 in the component list [L2 in Westall (1979)] d. LSIG3: Index for SOH in the component list [L3 in Westall (1979)] e. SSD: Surface site density (sites/m2) [C1 in Westall (1979)] f. SURFA: Specific surface area (m2/g) [C2 in Westall (1979)] g. CONCS: Concentration of sorbing solid (g/L) [C3 in Westall (1979)] h. XMU: Ionic strength (moles/L) [C4 in Westall (1979)] i. CAP1: Inner capacitance (F/m2) [C5 in Westall (1979)] j. CAP2: Outer Capacitance (F/m2) [C6 in Westall (1979)] iv. Sorbed Species ID: When generating isotherms, the program requires Species index numbers of the sorbed species in the tableau. The number of sorbed species displayed here is based on the ?Sorbed concentrations? set before pressing the ?Generate Input Template? button. Please input all the ID?s of the sorbed species here. This is used to calculate the amount of total sorbed species in moles/L. The program reads the 151 concentrations of all the sorbed species entered here and displays the total sum of their concentrations as total sorbed concentration. In case of a speciation problem where the sorbed concentrations are not required, this should be left blank. v. Component # of the aqueous species: This is only required when generating isotherms. Enter the component number of the aqueous species. The aqueous concentration of the component is calculated at the end of the simulation by subtracting the sorbed concentration from the total concentration of the component in moles/L. A1.2.6 Geochemistry equilibrium parameters (with transport) The input template is similar to the above except we have a few additional parameters. This section discusses about the additional parameters to avoid repetition. i. Boundary: The total concentrations of all the components at the inlet have to be entered here. The initial guess value for the concentrations is also entered here. ii. Initial: The total residual concentration of the existing in the column before the beginning of the simulation. iii. The tableau information is similar to the batch geochemistry problem iv. Sorbed phase concentrations: The initial concentration of the sorbed species is input in the cells. The number of sorbed concentrations is dependent on the input parameter set in Section-2. The user needs to enter the sorbed species index and the initial concentration for the sorbed phase. v. Sorbed phase species composition: In this section, we define the composition of sorbed phase concentrations. We have to provide the index of the sorbed phase concentration, and the number of species it is composed of and the species index # of all the species in the same row. 152 APPENDIX-A2 CALCULATION OF BACK-PARTITIONING COEFFICIENT Solubility of TCE = 1100 mg/L or 8.37 mM (Molecular weight = 131 g/mol) (Information, May 18, 2011) Solubility of cis-DCE = 3500 mg/L or 36 mM (Molecular weight = 96.95 g/mol) (Cleanup- Information, April 15, 2011) Solubility of VC = 2763 mg/L or 44.20 mM (Molecular weight = 62.5 g/mol) (ATSDR, March 3, 2011) Solubility of Ethene = 131 mg/L or 4.67 mM (Molecular weight = 28.05 g/mol) (Inchem, 2011) Regressed value of DCE back-partitioning coefficient for Experiment-A = 0.04 hr-1 Regressed value of DCE back-partitioning coefficient for Experiment-B = 0.004 hr-1 Similar calculations can be performed and the back-partitioning coefficients for the other daughter products of Experiment-B could be determined. 136 * 0 .0 4 0 .1 7 8 .3 7T C E hr 136 * 0 .0 4 0 .0 3 44VC hr 136 * 0 .0 4 0 .3 4 .6 7E th hr