Citation
Relationships between Mass and Flux in Dense Nonaqueous Phase Liquid Source Zones

Material Information

Title:
Relationships between Mass and Flux in Dense Nonaqueous Phase Liquid Source Zones
Creator:
FURE, ADRIAN D. ( Author, Primary )
Copyright Date:
2008

Subjects

Subjects / Keywords:
Architectural models ( jstor )
Average linear density ( jstor )
Groundwater ( jstor )
Modeling ( jstor )
Parametric models ( jstor )
Porosity ( jstor )
Simulations ( jstor )
Statistical discrepancies ( jstor )
Tracer bullets ( jstor )
Travel time ( jstor )

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright Adrian D. Fure. Permission granted to University of Florida to digitize and display this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Embargo Date:
7/30/2007
Resource Identifier:
71376109 ( OCLC )

Downloads

This item has the following downloads:


Full Text

PAGE 1

RELATIONSHIPS BETWEEN MASS AND FLUX IN DENSE NONAQUEOUS PHASE LIQUID SOURCE ZONES By ADRIAN D. FURE A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2005

PAGE 2

Copyright 2005 by Adrian D. Fure

PAGE 3

ACKNOWLEDGMENTS Primary thanks goes to my committee members, Drs. Mike Annable, Jim Jawitz, Kirk Hatfield, and Suresh Rao. I would especially like to thank my cochairmen Drs. Mike Annable and Jim Jawitz for their support and encouragement and in general just being a lot of fun to work with. I want to thank Dr. Annable for being supportive when I decided to switch research topics after my first year and always encouraging me to pursue avenues of inquiry which I found most interesting. I want to thank Dr. Jawitz for getting me fired up about streamtube models and turning me on to some alternative ways to approach contaminant transport problems. iii

PAGE 4

TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................iii LIST OF TABLES ............................................................................................................vii LIST OF FIGURES .........................................................................................................viii ABSTRACT .....................................................................................................................xiv 1 PARADIGM SHIFTS: THE CHALLENGE OF MANAGING DENSE NONAQUEOUS PHASE LIQUID CONTAMINATED SITES.................................1 Thomas Kuhn and Paradigm shifts ...............................................................................1 A Crisis of Technological Limitations .........................................................................3 Definition of Risk .........................................................................................................4 Alternative Performance Metrics for Source Zone Remediation .................................8 Technical Considerations ............................................................................................10 2 SIMPLIFIED SOURCE DEPLETION MODELS.....................................................12 Introduction .................................................................................................................12 General Mathematical Framework .............................................................................13 Streamtube Models .....................................................................................................14 Lagrangian Definition of NAPL Architecture.....................................................14 The Lognormal Distribution ................................................................................15 Equilibrium Streamtube Model ...........................................................................17 Concept of contribution time ........................................................................19 Derivation of contribution time ....................................................................19 Expected value for streamtube ensemble .....................................................20 Breakthrough curves and mass reduction/flux reduction relationships .......21 Rate Limited Streamtube Model .........................................................................25 Advection Dispersion Model ...............................................................................26 BTCs and Mass Reduction/Flux Reduction ........................................................27 Power Function Model ........................................................................................27 BTCs and Mass Reduction/Flux Reduction ........................................................31 Comparison of Source Depletion Models with Numerical Simulations ....................33 Constitutive Relationships ...................................................................................33 Simulation Setup .................................................................................................37 Simulation Set 1 ..................................................................................................38 iv

PAGE 5

Simulation domain and random field generation .........................................38 NAPL architectures ......................................................................................40 Dissolution profiles, model fits, and mass reduction/flux reduction relationships ..............................................................................................40 Simulation Domain and Grid Setup for Simulation Sets 2 through 4 .................44 Simulation Set 2 ..................................................................................................45 NAPL architectures ......................................................................................45 Dissolution profiles, model fits, and mass reduction/flux reduction relationships ..............................................................................................46 Simulation Set 3 ..................................................................................................48 NAPL architectures ......................................................................................48 Dissolution profiles, model fits, and mass reduction/flux reduction relationships ..............................................................................................49 Simulation Set 4 ..................................................................................................54 NAPL architectures ......................................................................................54 Dissolution profiles, model fits, and mass reduction/flux reduction relationships ..............................................................................................54 Conclusions and Recommendations for Further Research .........................................58 3 LABORATORY STUDIES OF RELATIONSHIPS BETWEEN ARCHITECTURE AND FLUX IN MODEL SOURCE ZONES..............................61 Introduction .................................................................................................................61 Modeling Approach ....................................................................................................64 Experimental Methods ................................................................................................64 Flow Cell Design .................................................................................................65 Porous Media Packing and NAPL Injection .......................................................67 Image Analysis ....................................................................................................69 Experimental Results, Modeling, and Discussion ......................................................70 NAPL Migration and General Dissolution Behavior ..........................................70 Flux Response to Changes in NAPL Architecture ..............................................72 Source Depletion Models ....................................................................................77 Additional NAPL Architectures ..........................................................................83 Conclusions .................................................................................................................87 4 AN EVALUATION OF PARTITIONING TRACERS FOR PARAMETERIZING SOURCE DEPLETION MODELS............................................................................89 Overview of Partitioning Tracers................................................................................89 Moment Based Parameter Estimation.........................................................................91 Measured moments..............................................................................................93 First and second moments of nonreactive tracer..................................................94 Domain average trajectory integrated NAPL content..........................................94 Excess spreading..................................................................................................96 Moment equations for quantifying Lagrangian NAPL architecture....................97 Inverse Modeling Based Parameterization (Lognormal Distribution)........................99 Additional Contribution Time Parameterization Methods.......................................101 v

PAGE 6

ADE based paramterization...............................................................................101 Inverse modeling with no travel time variability...............................................102 Simulation Setup.......................................................................................................103 Results.......................................................................................................................103 Predicted domain averaged NAPL saturation....................................................103 Method of moments dissolution predictions......................................................104 Inverse lognormal dissolution predictions.........................................................113 ADE dissolution predictions...............................................................................123 Comparison of techniques...................................................................................125 5 EXTENDING SOURCE DEPLETION MODELS TO FLUSHING BASED REMEDIATION.......................................................................................................139 Streamtube Models .................................................................................................139 Advection Dispersion Model ..................................................................................142 Extensions for Dispersive Mixing ..........................................................................143 Experimental Setup .................................................................................................146 Initial Architectures ................................................................................................148 Modeling of Effluent BTCs ....................................................................................149 Remedial Endpoints ................................................................................................152 6 SYNTHESIS: A FLUX BASED FRAMEWORK FOR MANAGING DNAPL CONTAMINATED SITES......................................................................................156 Three Definitions of Efficiency ..............................................................................156 Combining Three Efficiency Definitions ................................................................164 Source Identification (Pre-Reactor Based Characterization) ..................................165 Transitional Reactor-Based Characterization: Mass Flux as a Remedial Screening Tool and Qualitative Estimator of Aging ......................................166 Reactor-Based Characterization ..............................................................................168 LIST OF REFERENCES .................................................................................................175 BIOGRAPHICAL SKETCH ...........................................................................................183 vi

PAGE 7

LIST OF TABLES Table page 2-1. Parameters used in UTCHEM simulations. ................................................................38 2-2. Grid and permeability field parameters. .....................................................................39 2-3. Parameters assigned to each indicator class. ..............................................................40 2-4. Grid/random field parameters for simulation sets 2-4. ...............................................44 3-1. Summary of experimental conditions for conducted experiments. ............................65 4-1. Domain averaged NAPL saturations estimated with the different techniques. ........105 4-2. Tracer derived parameters using method of moments. .............................................107 4-3. Moment derived contribution time mean and variance. ...........................................110 4-4. Inverse lognormal derived mean and variance of contribution time. .......................118 4-5. Tracer derived parameters for ADE method. ...........................................................124 4-6. Tracer determined source depletion parameters for ADE model. ............................128 4-7. Moment derived parameters for dissolution experiment with tracers. .....................138 5-1. Fitted model parameters from UTCHEM surfactant simulations. ...........................154 vii

PAGE 8

LIST OF FIGURES Figure page 1-1. Components to consider in remedial design. ................................................................5 1-2. General approaches for managing risk at contaminated sites. ......................................7 1-3. Source zone TMDL approach of Einarson and Mackay (2001) ................................10 2-1. Eulerian and Lagrangian definitions of NAPL architecture. ......................................15 2-2. Illustration of the tau concept. ....................................................................................20 2-3. BTCs for increasing variability in contribution time (m2). ......................................23 2-4. Mass reduction/flux reduction relationships for BTCs in figure 2-3. .........................24 2-5. BTCs for increasing Peclet number. ...........................................................................28 2-6. Mass reduction/flux reduction relationships for BTCs in figure 2-5 ..........................28 2-7. BTCs for the Da model. ..............................................................................................32 2-8. Mass reduction/flux reduction relationships for Da model ........................................32 2-9. Comparison of NAPL spill using Brooks-Corey and Van Genuchten models. .........35 2-10. Raining effect observed when running UTCHEM in restart mode. .........................35 2-11. NAPL spills for set 1 (contour range: 0 to 25% saturation). ....................................40 2-12. Streamtube model fit of UTCHEM dissolution data (set 1). ....................................42 2-13. ADE model fit of UTCHEM dissolution data (set 1). ..............................................42 2-14. Mass reduction/flux reduction fit for streamtube model (set 1). ..............................43 2-15. Mass reduction/flux reduction fit for ADE model (set 1). .......................................43 2-16. NAPL spills for set 2 (contour range 0 to 25% saturation). .....................................45 2-17. Streamtube model fit of UTCHEM dissolution data (set 2). ....................................46 viii

PAGE 9

2-18. ADE model fit of UTCHEM dissolution data (set 2). ..............................................47 2-19. Mass reduction/flux reduction fit for streamtube model (set 2). ..............................47 2-20. Mass reduction/flux reduction fit for ADE model (set 2). .......................................48 2-21. NAPL spills for set 3 (contour range: 0-25% saturation). ........................................49 2-22. Streamtube model fit of UTCHEM dissolution data (Set 3). ...................................51 2-23. ADE model fit of UTCHEM dissolution data (Set 3). .............................................51 2-24. Discontinuous fit of simulation 4 (set 3). .................................................................52 2-25. Streamtube model fit of mass reduction/flux reduction relationships. .....................53 2-26. ADE model fit of mass reduction/flux reduction relationships. ...............................53 2-27. NAPL spills for set 4 (contour range: 0 to 25% saturation) .....................................55 2-28. Streamtube model fit of UTCHEM dissolution data (set 4) .....................................56 2-29. ADE model fit of UTCHEM dissolution data. .........................................................57 2-30. Streamtube model fit mass reduction/flux reduction relationships (set 4). ..............57 2-31. ADE model fit of mass reduction/flux reduction relationships (set 4). ....................58 3-1. Schematic representation of segmented box. .............................................................66 3-2. Flux plane response to changes in NAPL architecture for experiment DCA-1. Note that MF is the mass fraction ((mass at time t)/(initial mass at t = 0)) and FF is the flux fraction ((flux at time t) /(initial flux at t = 0)). .......................................73 3-3. Flux plane response to changes in NAPL architecture for experiment TCE-1. Note that MF is the mass fraction ((mass at time t)/(initial mass at t = 0)) and FF is the flux fraction ((flux at time t) /(initial flux at t = 0)). ................................................74 3-4. Flux plane response to changes in NAPL architecture for experiment DCA-2. Note that MF is the mass fraction ((mass at time t)/(initial mass at t = 0)) and FF is the flux fraction ((flux at time t) /(initial flux at t = 0)). .......................................75 3-5. Cumulative reductions in mass fraction and flux fraction with increasing dissolution time for selected ports from experiments DCA-1, TCE-1, and DCA-2. Note that the cumulative reductions are expressed in relation to the total mass reduction and the total flux reduction (e.g. 27% of the total reduction in flux for experiment DCA-2 comes from port 5). ..................................................................77 ix

PAGE 10

3-6. Mass reduction/flux reduction relationships for the entire system (upper three graphs) and for selected well segments (bottom three graphs). ...............................78 3-7. Streamtube mass fractions estimated from the image analysis technique. .................79 3-8. Comparison of the rate-limited streamtube, equilibrium streamtube, and effective Damkohler approaches for modeling source depletion from experiments DCA-1, TCE-1, and DCA-2. .................................................................................................80 3-9. ADE model fit for source depletion experiments DCA-1, TCE-1, and DCA-2. ........82 3-10. Initial NAPL distribution and effluent BTC for experiment DCA-3. ......................83 3-11. Mass reduction flux reduction relationship for experiment DCA-3 and hypothetical system developed by the superposition of selected segmented BTCs. ........................................................................................................................85 3-12. Comparison of the rate-limited streamtube, equilibrium streamtube, effective Damkohler, and ADE approaches for modeling effluent BTCs from experiments DCA-3. .....................................................................................................................86 4-1. UTCHEM generated reactive tracer fit with product of two lognormals. ..................92 4-2. RMS matrix for tracer data in figure 4-1. ...................................................................93 4-3. BTCs representing transport of a nonreactive tracer, a reactive tracer with homogenously distributed NAPL, and a reactive tracer with heterogeneously distributed NAPL. ....................................................................................................97 4-4. Summary of process used to predict dissolution using partitioning tracers and moment analysis. ......................................................................................................98 4-5. Cubic spline interpolation of arbitrary UTCHEM generated tracer data. ..................99 4-6. Set 1 tracer fit (method of moments). .......................................................................108 4-7. Set 2 tracer fit (method of moments). .......................................................................108 4-8. Set 3 tracer fits (method of moments). .....................................................................109 4-9. Set 4 tracer fits (method of moments). .....................................................................109 4-10. Set 1 predicted source depletion (method of moments ...........................................110 4-11. Set 2 predicted source depletion (method of moments). ........................................111 4-12. Set 3 predicted source depletion (method of moments). ........................................112 4-13. Set 4 predicted source depletion (method of moments). ........................................112 x

PAGE 11

4-14. Set 1 tracer fit (inverse lognormal). ........................................................................116 4-15. Set 2 tracer fits (inverse lognormal). ......................................................................116 4-16. Set 3 tracer fit (inverse lognormal). ........................................................................117 4-17. Set 4 tracer fit (inverse lognormal). ........................................................................117 4-18. Set 1 predicted source depletion (inverse lognormal). ...........................................119 4-19. Set 1 predicted source depletion (inverse lognormal-no travel time variability). ..119 4-20. Set 2 predicted source depletion (inverse lognormal). ...........................................120 4-21. Set 2 predicted source depletion (inverse lognormal w/ no travel time variability). .............................................................................................................120 4-22. Set 3 predicted source depletion (inverse lognormal). ...........................................121 4-23. Set 3 predicted source depletion (inverse lognormal w/ no travel time variability). .............................................................................................................121 4-24. Set 4 predicted source depletion (inverse lognormal). ...........................................122 4-25. Set 4 predicted source depletion (inverse lognormal w/ no travel time variability). .............................................................................................................122 4-26. Set 1 tracer fit (ADE). .............................................................................................126 4-27. Set 2 tracer fit (ADE). .............................................................................................126 4-28. Set 3 tracer fit (ADE). .............................................................................................127 4-29. Set 4 tracer fits. (method of moments). ..................................................................127 4-30. Set 1 predicted source depletion (ADE). ................................................................129 4-31. Set 2 predicted source depletion (ADE). ................................................................130 4-32. Set 3 predicted source depletion (ADE). ................................................................130 4-33. Set 4 predicted source depletion (ADE). ................................................................131 4-34. Root mean square errors for parameterization techniques. Points indicate average RMS for simulation set. Lines indicate average RMS for all simulations. 131 4-35. Impact of increased longitudinal dispersivity on tracer response and dissolution profile. ....................................................................................................................132 xi

PAGE 12

4-36. Relationship between ln(K) and contribution time parameters. .............................135 4-37. Relationship between travel time, S and contribution time variability. ...............136 4-38. Initial TCE distribution for dissolution experiment with tracers. ...........................136 4-39. Tracer BTCs (top graph) and tracer BTCs selected for moment analysis (bottom graph). ....................................................................................................................137 4-40. Experimental and predicted dissolution data. .........................................................138 5-1. Dissolution profiles for aqueous dissolution and flushing based remediation. ........140 5-2. Flushing BTCs for increasing second moment of .................................................142 5-3. Flushing BTCs for increasing second moment of .................................................144 5-4. Initial PCE distribution for experiments 1 and 2. .....................................................147 5-5. Model fit of ethanol BTCs for experiments 1 (top graph) and 2. .............................148 5-6. Model fit of experiment 1 data. ................................................................................149 5-7. Equivalent aqueous dissolution source depletion for experiment 1 .........................150 5-8. Model fit of experiment 2 data. ................................................................................151 5-9. Equivalent aqueous dissolution source depletion for box 2. ....................................151 5-10. Pre and post flux-averaged concentrations in experiments 1 and 2. ......................153 5-11. Model fits of UTCHEM surfactant floods. .............................................................155 6-1. Comparison of source zones A, B, & C using mass reduction efficiency. ...............157 6-2. Comparison of source zones A, B, & C using reduction in source discharge per reduction in mass. ...................................................................................................158 6-3. Comparison of source zones A, B, & C using reduction in source discharge efficiency. ...............................................................................................................159 6-4. Aging effects on relationships between reduction in mass and reductions in source discharge for source zone A, B, & C. .....................................................................161 6-5. Aging effects on source discharge efficiency for source zone A, B, & C. ...............163 6-6. Source discharge profiles for a fixed initial mass and alternative initial source discharge. ................................................................................................................164 xii

PAGE 13

6-7. Figure 6-2 modified to include net work input and site specific flux-based performance metric. ...............................................................................................165 6-8. Source discharge profile and down-gradient flux planes at selected time steps. .....168 6-9. Key assumptions required when using spatially resolved flux information for flux mapping. .................................................................................................................170 6-10. Use of unique nonreactive tracers for flux mapping. .............................................171 6-11. Conceptual example of using figure 6-5 in the remedial design process. ..............174 xiii

PAGE 14

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy RELATIONSHIPS BETWEEN MASS AND FLUX IN DENSE NONAQUEOUS PHASE LIQUID SOURCE ZONES By Adrian D. Fure August, 2005 Chair: Michael D. Annable Cochair: James W. Jawitz Major Department: Environmental Engineering Sciences The relationship between reductions in dense nonaqueous phase liquid (DNAPL) mass and the corresponding reduction in contaminant source discharge is evaluated in experimental and numerical studies. Analytical solutions for modeling the relationships between mass and source discharge under natural gradient conditions are developed. The models are extended to also consider flushing based source zone remediation. The models are compared with numerical and experimental studies of natural gradient dissolution and flushing based remediation. The use of reactive and nonreactive tracers for parameterizing the analytical models is also evaluated. The utilization of the models in a flux based framework for addressing DNAPL contaminated sites is discussed. xiv

PAGE 15

CHAPTER 1 PARADIGM SHIFTS: THE CHALLENGE OF MANAGING DENSE NONAQUEOUS PHASE LIQUID CONTAMINATED SITES In this chapter I will provide a brief historical sketch of the approaches utilized to manage sites contaminated with dense nonaqueous phase liquids (DNAPLs). I have attempted to frame this discussion in terms of Thomas Kuhns The Structure of Scientific Revolutions as I believe the transition from mass-based to flux-based performance metrics for DNAPL source zones is a paradigm shift, marked by many of the characteristics outlined by Kuhn. Another reason I have elected to incorporate Kuhn is his notion of normal science or paradigm articulation. I would consider most of my own work as fitting into this classification, i.e. I feel that much of my work has been directed towards advancing and refining an already introduced flux-based paradigm (Rao et al., 2001). Thomas Kuhn and Paradigm shifts The notion of paradigmsthe idea that the direction of scientific inquiry and the interpretation of results is significantly influenced by background contextscan be to traced to Thomas Kuhns influential book The Structure of Scientific Revolutions. Kuhn was a graduate student in theoretical physics at Harvard within sight of the end of [his] dissertation when he became involved with an experimental college course treating the physical sciences for the non-scientist. This course would drastically alter his career path, taking him from physics to the history of science and eventually back to the 1

PAGE 16

2 philosophical underpinnings of science. He would publish The Structure of Scientific Revolutions 15 years later in 1962. Kuhn suggests that paradigm shifts share two essential characteristics, they are sufficiently unprecedented to attract a group of adherents away from competing modes of scientific inquiry and they are open-ended enough to leave all sorts of problems for the redefined group of practitioners to resolve (Kuhn refers to the resolution of these problems as normal science or paradigm articulation). Although Kuhn dealt only with truly revolutionary paradigm shifts, such as the shift from the Ptolemaic system of astronomy (which placed the earth at the center of the solar system) to Copernican astronomy (which placed the sun at the center of the solar system), the same patterns of paradigm shifts have been applied extensively and perhaps often inappropriately to other less revolutionary, non-scientific modes of inquiry that have undergone fundamental change in the background context within which they operate. In my opinion, a shift from mass-based to flux-based performance metrics for dense nonaqueous phase liquid (DNAPL) source zones constitutes a paradigm shift which, to borrow again from Kuhn, requires the reconstruction of prior theory and the re-evaluation of prior fact. This is especially true with regards to surfactant-based flushing technologies, a technology adapted from the oil industry which, by definition, has been previously designed and optimized to most effectively remove mass from the system, not to yield the largest reduction flux. The mass reduction/flux reduction relationships discussed throughout serve as a sort of paradigm transfer function seeking to take at least some of the large body of work that has operated within the mass reduction paradigm and cast it into the new flux based paradigm.

PAGE 17

3 A Crisis of Technological Limitations Kuhn suggests that paradigm shifts occur as the result of crisis or the inability of an existing paradigm to describe certain observed behavior. In the case of Ptolemaic astronomy, the crisis arose from complications describing astronomical observations outside of planetary positions, such as equinoxes. Kuhn suggests that out of crises arise any number of competing paradigms, which vie for preeminence within the scientific community. The initial stages of a paradigm shift are invariably characterized by elevated tension between members belonging to the previous paradigm and those members aligning themselves with any one of the competing paradigms. Eventually, a dominant paradigm arises (e.g. that of Copernican astronomy) leading to the gradual disappearance of both the original paradigm and any competing paradigms. When, in the development of a natural science, an individual or group first produces a synthesis able to attract most of the next generations practitioners, the older schools gradually disappear. In part their disappearance is caused by the members conversion to the new paradigm. But there are always some men who cling to one or another of the old views, and they are simply read out of the profession, which thereafter ignores their work. The new paradigm implies a new and more rigid definition of the field. Those unwilling to or unable to accommodate their work to it must proceed in isolation or attach themselves to some other group. The crisis that precipitated the proposed shift from a mass based to a flux based paradigm is not one in which the preexisting paradigm was unable to describe a given set of physical observations, but rather a crisis of technological limitations. In this case, the technological limitations are those associated with the suite of potential remediation technologies available to those charged with managing sites contaminated with DNAPLs, none of which perform as a silver bullet technology capable of efficiently removing or destroying all of the contaminant mass from the subsurface and thereby, prima facie, appearing as failures. A paradigm shift occurs by considering these same technologies

PAGE 18

4 not for their ability to remove or destroy mass, but for their ability to reduce risk, which does not require the complete removal of mass. Such a shift, as Kuhn suggests, has been characterized by initial resistance, both from those within the regulatory and scientific communities. I hope the flux paradigm will soon be characterized, also as Kuhn suggests, by a large influx of research that advances and refines the new paradigm (paradigm articulation) as it gains acceptance as the more rational paradigm. Definition of Risk Of critical importance to the flux-based paradigm is the contention that there is no direct linkage between risk and the amount of DNAPL mass present at a given site as a receptor will virtually never come into direct contact with a contaminant in the nonaqueous phase. The risk associated with a given DNAPL contaminated site is thereby derived from that particular sites ability to transfer contaminant mass from the nonaqueous phase to aqueous or gaseous phases which can migrate beyond the boundaries of the source zone and potentially interface with contaminant receptors. An evaluation of risk would thereby require an evaluation of the rate at which contaminant mass is being transferred from the nonaqueous to the aqueous and/or gaseous phases (source discharge) coupled with predictive contaminant fate and transport analysis to determine potential contaminant receptors, exposure pathways, and associated toxicological impacts. While acknowledging gaseous phase transport, the following analysis is primarily directed towards aqueous phase transport. Therefore, unless specifically stated the term risk will be used hereafter to specifically refer to receptor exposure to a dissolved contaminant plume.

PAGE 19

5 Historical Sketch of Approaches for Managing DNAPL Contaminated Sites Based upon the discussion outlined in the previous section there are three components to consider when managing the risk associated with a DNAPL contaminated site: the source, the receptor, and the pathway by which the contaminant is transferred from the source to the receptor (figure 1). Management strategies can conceptually be delineated into three groups based upon which of the three components (the source zone, the exposure pathway, or the receptor) are targeted (figure 2). Source ZoneTransport Vehicle (e.g. flowing groundwater) ReceptorExposure Pathway Nonaqueous PhaseMass Transfer Figure 1-1. Components to consider in remedial design. Institutional based approaches provide an immediate short-term reduction in risk by removing the receptor from the exposure pathway. An example of such an institutional-based approach would be taking a drinking water well offline due to the interaction of the well with a dissolved contaminant plume. Institutional approaches of course make no

PAGE 20

6 attempt to remediate or restore sites to a level where the site can again be used for industrial, residential, or recreational purposes. Plume management based approaches target the exposure pathway by either removing/stripping the contaminant from the transport vehicle (e.g. reactive barriers), removing the transport vehicle from the exposure pathway (e.g. pump and treat), or by enhancing transformation of the contaminant to a more benign derivative as it migrates along the exposure pathway (e.g. enhanced bioremediation). Pump and treat, which extracts contaminated groundwater from the subsurface and treats it above ground, historically has been the approach most employed at DNAPL contaminated sites (NRC, 1994). The major limitation associated with pump and treat technologies when applied to DNAPL sites stems from the longevity of most DNAPL source zones, which tend to function as sources of groundwater contamination for extended periods of time. Wide recognition in the early to mid 90s that pump and treat was by and large ineffective at most DNAPL sites due to the longevity of DNAPL source zones (NRC, 1994; Travis and Doty, 1990) sparked interest in approaches directed towards the removal or destruction of the source of contamination.

PAGE 21

7 Source Zone Remediation Plume Management Institutional Eliminate source discharge by complete removal/destruction (e.g. surfactant flushing; chemical oxidation)Reduce source discharge via partial removal/destructionReduce source discharge by limiting access of transport vehicle to source (e.g. slurry walls)Remove transport vehicle from exposure pathway (e.g. pump and treat)Remove contaminant from transport vehicle in-situ (e.g. permeable reactive barriers)Enhance in-situ transformation of contaminant as it travels along the exposure pathway (e.g. enhanced bioremediaton)Source ZoneExposure PathwayReceptorRemove receptor from exposure pathway (e.g. taking drinking water well offline) Figure 1-2. General approaches for managing risk at contaminated sites. It was of course initially hoped for that the research and funding directed towards source zone remediation would yield technologies capable of complete contaminant removal/destruction. Source zone remediation is inherently complex however, from the initial characterization of the site to the eventual implementation of the given technology, and has, as a result, prohibited the development of a silver bullet type technology which can efficiently remove or destroy all of the contaminant mass in the source zone (Soga et al., 2004). Several source zone remediation technologies have however shown the ability to remove/destroy a large portion of the contaminant mass which introduces questions as to the benefits of partial source removal and whether those benefits are

PAGE 22

8 commensurate to the often considerable cost required to achieve partial source removal (E.P.A., 2003). Alternative Performance Metrics for Source Zone Remediation A recent E.P.A. report (EPA, 2003) groups source zone remediation performance metrics into three categories. Type I metrics are measurements that can be reliably acquired and are commonly used such as the total mass of DNAPL removed from the subsurface and changes in resident groundwater concentrations. Type II metrics are metrics that can sometimes be measured, but are not in wide use such as remaining DNAPL mass and DNAPL architecture. Type III metrics are metrics that are theoretically possible and under development such as mass flux and mass discharge. To consider partial source removal as a viable remedial alternative requires an alternative performance metric that is not based upon reducing local resident concentrations in the source zone below regulatory maximum contaminant levels (MCLs) (type I metrics), but rather on reducing risk to down-gradient receptors below a given threshold. A risk-based analysis requires consideration of both the source strength (type III metrics) and fate and transport modeling of the dissolved contaminant plume. Rao et al. (2001) discuss attenuation based performance metrics for aggressive source zone remediation technologies based on whether the contaminant flux is reduced to a level where the following four criteria are met: (a) the spatial extent of the existing dissolved plume is stable or decreases; (b) the total contaminant mass within the plume is constant or diminishing; (c) both the average concentration and the range in concentrations is diminishing; (d) contaminant fluxes decrease at succeeding control planes along the dissolved plume. If these conditions are met, then the dissolved plume is stable (source

PAGE 23

9 strength equals attenuation capacity) or decreasing (source strength is less than the attenuation capacity) (Rao et al., 2001). When considering the ecological and human impacts of a migrating dissolved contaminant plume, one is typically concerned with communication between the plume and down-gradient surface water bodies and/or drinking water supply wells. In relation to surface water bodies, where there are potentially multiple contaminant sources contributing to the overall contaminant load, a mass discharge approach for evaluating and prioritizing source zone remediation is clearly more appropriate than an MCL approach as a mass discharge approach would fit into the preexisting Total Maximum Daily Load (TMDL) framework. Einarson and Mackay (2001) suggest a similar TMDL type approach for the protection of water supply wells with the pollutant mass balance being performed on the capture zone as opposed to the surface water body. The approach used in the source zone TMDL approach to estimate the contaminant concentration at a drinking water supply well is to divide the total pollutant discharge to the capture zone by the total extraction flow rate from the capture wells (figure 3). As noted by Einerson and Mackay (2001), the TMDL approach, if implemented, would establish a more rational framework for addressing multiple sources and coping with the reality that the resources available for addressing groundwater contamination are substantially smaller than the scope of the problem. In the case of multiple sources, a TMDL approach would also provide an opportunity for emissions trading in manner directly analogous to the current approach utilized for air and water based emissions.

PAGE 24

10 MdMd Capture ZoneDrinking Water Supply Wells wellsjjsourcesidwellQMCi11Dissolved Contaminant Plumes Q Figure 1-3. Source zone TMDL approach of Einarson and Mackay (2001) Technical Considerations The purpose of this opening chapter is to provide a socio-political context for the more technical portion that follows. I want to emphasize that a flux-based paradigm/source zone TMDL approach is not widely accepted and is somewhat controversial. Aside from mentioning this controversy in passing and simply stating that the source zone TMDL approach is, in my opinion, the most rational and appropriate way to approach DNAPL contaminated sites, I will not evaluate the socio-political components of this controversy in any sort of systematic fashion. My main focus here is to look at some of the technical ramifications of adopting a flux based paradigm with my main thesis being that many of the technical challenges associated with the management of DNAPL contaminated sites are more tractable in a flux-based paradigm.

PAGE 25

11 In the following chapters, analytical models are introduced which are useful for addressing some of the difficult questions surrounding the benefits of partial source removal. In chapter 2, the models are introduced and compared with numerical simulations conducted with the multiphase flow and transport simulator UTCHEM to determine if the analytical models are capable of reproducing the range of dissolution profiles produced by UTCHEM. Chapter 3 consists of experiments designed to qualitatively investigate linkages between source discharge and NAPL architecture. The experiments of chapter 3 were also designed to validate some of the underlying assumptions of the analytical models. In chapter 4, the use of tracers to parameterize the analytical models was investigated. The overall performance of the tracers served as an additional validation of the model assumptions. In chapter 5, the models are extended to handle flushing-based (surfactant/cosolvent) remediation. In chapter 6, after the models have been validated, a framework for utilizing the models is suggested. The models are also used to answer some important questions related to the benefits of partial source depletion.

PAGE 26

CHAPTER 2 SIMPLIFIED SOURCE DEPLETION MODELS Introduction The decision process for risk-based corrective action at sites contaminated with DNAPLs introduces several difficult questions, many of which do not have answers that are widely accepted throughout the technical and regulatory communities. Of critical importance when considering plume management vs. aggressive source zone remediation is the prediction of the long term natural-gradient dissolution of the source zone to determine the required duration of plume management. The use of numerical simulators for the prediction of long term source depletion is perhaps impractical due to the sensitivity of a migrating DNAPL to small scale changes in permeability, which leads to complex entrapment configurations (e.g., Keuper et al., 1993) which are next to impossible to characterize using currently available technologies (Parker and Park, 2004). The practical limitations of conventional numerical simulators has led researchers to propose simplified screening models which serve as analogs to more sophisticated numerical simulators (Parker and Park, 2004; Zhu and Sykes, 2004). These screening models require relatively few parameters and can be efficiently coupled to contaminant fate and transport models as source terms. Screening models are also useful for evaluating uncertainty in key input parameters (e.g., NAPL mass) and the associated effect of that uncertainty on predicted source depletion. In this chapter three simplified NAPL source depletion models are considered, a streamtube model, an advection dispersion model, and the power function model of Parker and Park (2004). 12

PAGE 27

13 There are two key questions to consider in relation to utilizing simplified source depletion models in the decision process for DNAPL contaminated sites. The first question is whether the model is capable of reproducing the range of source depletion behavior observed in field, laboratory and numerical studies of source depletion. The second question is how would the parameters for these models be determined in order to be of use to site managers? There are two approaches that can be employed. The first approach is to evaluate an expected range of model parameters based on typical values encountered in field, laboratory and numerical studies and incorporate that range/uncertainty into the site specific decision process. The second approach is to try and physically measure the parameters in the field. In order to do this, the parameters must take on some physical meaning (as opposed to being just fitting parameters) and there must be an available technology to measure the parameter. In this chapter I will investigate part of the first question: whether these simplified source depletion models can serve as mathematical analogs to more sophisticated numerical simulators. In chapter 3, I will evaluate these models with respect to experimental studies of source depletion. In chapter 4, I will investigate the second question, whether these models can be parameterized using existing technologies. General Mathematical Framework The three models introduced below take on a similar formulation where the flux-averaged concentration exiting a source zone is related to a source depletion term (S D ) as in equation 3-1 below. )(1)(TSCfTCDscf (2-1)

PAGE 28

14 Where C f is the flux averaged concentration exiting the source zone [M/L 3 ], C s is the solubility of the NAPL in groundwater (or flushing solution) [M/L 3 ], and T is the flushing duration. The contaminated fraction, f c is the fraction of streamtubes initially containing NAPL. It is expected that the flux averaged concentration exiting a source zone will be well below solubility (MacKay et al., 1985) resulting from both a physical nonequilibrium brought about by streamtubes that do not intersect with the contaminated fraction of the source zone (manifested in f c ) and are subsequently blended in the extraction well with contaminated streamtubes, and to a lesser extent by chemical nonequilibrium brought about by rate limited mass transfer from the NAPL to the flowing aqueous phase. The three model types introduced below share the same general form as equation 3-1, differing only in the way in which the source depletion term (S D ) is formulated. Streamtube Models Lagrangian Definition of NAPL Architecture Sale and McWhorter (2001) use the term source zone architecture to refer to the geometry, spatial distribution, and saturation of a collection of subzones comprising a source zone. Quantitatively there are two ways to define NAPL architecture, an Eulerian approach and a Lagrangian approach. In the Eulerian definition of source zone architecture, point values of NAPL saturation are expressed in terms of their magnitude and location in an arbitrary coordinate system. The concept behind the Lagrangian-based definition of NAPL architecture is to resolve the source zone into a collection of streamtubes with variability in NAPL saturation (and mass transfer rate coefficients) along each streamtube integrated and transformed into an effective value for each streamtube (Cvetkovic et al., 1998; Jawitz et al., 2003a). Following Jawitz et al. (2003a),

PAGE 29

15 this integrated parameter is referred to as the trajectory integrated NAPL content, wNSS (where is the trajectory integrated NAPL saturation, is the porosity [L NS 3 L -3 ], and w is the water content [L 3 L -3 ]). The paramteris utilized in place of because it exhibits a range of (0,) as opposed to the (0,1) range exhibited by which is consistent with lognormal and gamma probability distributions. The Lagrangian definition of source zone architecture would then be defined through some form of a probability distribution which would describe the distribution of trajectory integrated NAPL contents among the streamtube network. S NS NS NAPL saturation/velocity defined in terms of x,y,zlocation Eulerian Definition of NAPL Architecture Source Zone discretized into a network of streamtubes Mass FluxLagrangian Definition of Source Zone Architecture Figure 2-1. Eulerian and Lagrangian definitions of NAPL architecture. The Lognormal Distribution In the streamtube models discussed below, the travel time, t, and are considered random variables, described by a lognormal probability density function (pdf). S

PAGE 30

16 222))(ln(exp21)(xxxf (2-2) Where is the mean of the log transformed variable x and is the standard deviation. When considering process models with two random variables, a joint pdf must be utilized if the random variables are correlated. )1(22exp121),(2222yyxxyxYYYYxyyxf (2-3) xxxxY)ln( (2-4) yyyyY )ln( (2-5) Where is the correlation coefficient between x and y. If x and y are uncorrelated then (2-3) can be expressed as the product of f(y) and f(x). The expected value of x is given by the following: (2-6) dxxfxxEnn0)( The analytical solution of (2-6) is the well known moment generating function for a lognormal distribution (e.g. Jury and Roth 1990): 2exp22NNmxEXNN (2-7) Given the first two moments, any higher moments can be determined by using (2-7) and solving for and 2)ln()ln(221XxXmm (2-8)

PAGE 31

17 (2-9) )ln(2)ln(12XXXmm A final property of the lognormal distribution that has relevance to the streamtube models below is that the product of two lognormals, Z = XY, can be substituted for with another lognormal with mean and variance: YXZ (2-10) 222YYXXZ (2-11) Equilibrium Streamtube Model The first model considered is a version of a model presented in Jawitz et al. (2005) that has been modified for natural gradient source depletion. The model is built upon the contention that the Lagrangian-based NAPL architecture and its relation to the velocity field exert primary influence over dissolution dynamics, where additional mechanisms not explicitly accounted for in the streamtube model, such as rate-limited mass transfer, transverse dispersion, and relative permeability exert only secondary influence, causing perturbation about the mean behavior. The assumption that a parcel of water that traverses the length of a contaminated source zone, from injection plane to extraction plane, is at solubility by the time it reaches the extraction plane is well supported in the literature as column studies of aqueous dissolution have shown the spatio-temporal scales required to reach equilibrium concentrations to be on the cm scale for typical groundwater velocities (Miller et al., 1990; Powers et al., 1994). The fundamental assumption of the streamtube models is that dissolution dynamics are governed primarily by the NAPL architecture and its relation to the velocity field such that additional mechanisms not accounted for in the mathematical formulation such as relative

PAGE 32

18 permeability and transverse dispersion introduces negligible error. The major advantage of simplifying dissolution dynamics by neglecting mechanisms such as transverse dispersion and relative permeability is that it allows for the development of relatively compact analytical expressions that can be parameterized with tracer tests (Jawitz et al., 2000a) and are potentially of more practical value to site managers (Parker and Park, 2004). The assumption that NAPL architecture is a principal controlling mechanism of dissolution dynamics is well supported by previous numerical modeling studies of NAPL infiltration and dissolution/remediation (Mayer and Miller, 1996; Dekker and Abriola, 2000; Lemke et al., 2004b; Parker and Park, 2004). The effects of neglecting transverse mixing between streamtubes has been investigated previously by Fiori (1996) and Berglund and Fiori (1997). Berglund and Fiori (1997) determined that the advection-only approach is appropriate for sorbing solutes under typical field conditions. Other researchers have incorporated mixing effects into streamtube formulations (Cirpka and Kitanidis, 2000; Ginn 2001) in order to model reaction processes that are highly dependent on transverse mixing (Kapoor et al., 1997; Kapoor et al., 1998; Cirpka et al., 1999). Laboratory studies have documented the importance of transverse mixing in the dissolution of NAPL pools (Chrysikopoulos et al., 2000; Eberhardt and Grathwohl, 2002), yet the effect of transverse mixing when considering dissolution dynamics at the field scale are somewhat less certain. For example, in the numerical simulations of surfactant-based DNAPL source zone remediation of Lemke et al. (2004b), the mass flux and PCE recovery were found to be insensitive to transverse dispersivity values ranging from .0075m to 0.06m, values extending well beyond the upper limit reported in field and

PAGE 33

19 laboratory experiments (Klenk and Grathwohl 2002). Further research is needed to determine the effect of neglecting transverse dispersion on field scale dissolution dynamics. The effect of neglecting relative permeability is also uncertain and an area for further research. If relative permeability effects are large the flow paths (i.e. the travel time and NAPL associated with that travel time) may be temporally dynamic, changing throughout the dissolution process (Geller and Hunt, 1993; Powers et al., 1998). Concept of contribution time The model formulation below is based on the concept of contribution time, denoted throughout as The parameter can be thought of as the time at which a given streamtube is depleted of NAPL under a set of specific dissolution/flushing conditions. For an arbitrary flushing/dissolution duration (T), all streamtubes with values less than T are clean, while all streamtubes with values greater than T are still contaminated. Under the equilibrium assumption, the concentration discharging from the contaminated streamtubes is at solubility while the concentration discharging from the clean streamtubes is zero. This concept is illustrated in figure 2-2. Derivation of contribution time Utilizing the assumptions listed above it is possible to consider a simple mass balance on an individual streamtube: StKvACLASQCMfwsNnsn (2-12) where M n [M] is the NAPL mass in the individual streamtube, Q is the flow rate [L 3 /T], n is the NAPL density, v is the velocity in the streamtube [L/T], L is the length of the streamtube, A is the cross sectional area of the streamtube [L 2 ], t is the residence time or the travel time and K f is ratio of the NAPL density to the solubility ( n /C s ). Utilizing

PAGE 34

20 Figure 2-2. Illustration of the tau concept. equations 2-10 and 2-11 and substituting a lognormal random variable t S for the product of two random lognormals t and, equation 2-12 reduces to S SftK with mean and variance given by: tStS 222ttSStS (2-13) Expected value for streamtube ensemble Using the concept of contribution time, the ensemble average for a collection of streamtubes is given by equation 2-1 with S D equal to the cumulative distribution function of :

PAGE 35

21 )(1)(TSCfTCDscf with: )()( pTSD (2-14) where p() is the cumulative distribution function of the contribution time (). The mean and variance for the contribution time can be determined by employing equations 2-7 through 2-9: 2)ln()ln(221mm (2-15) (2-16) )ln(2)ln(12mm (2-17) StfmKm11 (2-18) StfmKm222 Breakthrough curves and mass reduction/flux reduction relationships The magnitude of the average contribution time (represented by m 1 ) increases with a decrease in average velocity or an increase in the domain averaged NAPL saturation. Contribution time heterogeneity (manifested in m 2 ) increases by introducing more variability to the flow field and/or the NAPL distribution. Contribution time heterogeneity also increases through an increase in positive correlation between travel time and NAPL saturation. The physical meaning of a positive correlation between NAPL saturation and travel time is that more of the NAPL mass would be present in the lower conductivity, lower velocity media. Qualitative studies of NAPL spills (e.g., Keuper et al., 1993; Schwille, 1998) have demonstrated that a migrating NAPL is extremely sensitive to small scale changes in permeability, moving preferentially through

PAGE 36

22 courser grained media, suggesting a negative correlation between travel time and NAPL saturation. A suite of breakthrough curves (BTCs) for a fixed average contribution time (fixed m 1 ; i.e., fixed NAPL mass) and increasing variability in contribution time (expressed through increasing m 2 ) is displayed below in figure 2-3. Within the streamtube framework outlined above, the shape of the BTC exiting the source zone is controlled by the distribution of NAPL saturations and velocities among the streamtubes rather than the equation that describes transport at the scale of the individual streamtube. For the case of a homogenous distribution of NAPL saturations in a uniform flow field (m 2 = 1.01 in figure 2-3), the dissolution profile resembles a step function with each streamtube depleted of mass at an identical flushing duration. This type of dissolution behavior would most resemble column-based experiments of uniformly distributed residual NAPL saturations (Imhoff et al., 1994; Powers et al., 1994), the PCE component of the relatively homogenous Borden emplaced source (Frind et al., 1999), and the theoretical analysis of Sale and McWhorter (2001). In this scenario there is limited reduction in contaminant discharge until most of the NAPL mass has been depleted (figure 2-4). With increased variability in contribution time ( distribution) the dissolution profile moves from a step function, to more of an exponentially decaying dissolution profile, and eventually to a dissolution profile at high variability that is characterized by a rapid drop in flux followed by extensive tailing. It is possible to consider the BTCs of figure 2-3 as having two distinct components: an accessible component associated with an initial flux decrease and a recalcitrant

PAGE 37

23 component associated with the late-term tailing portion. As variability in contribution time increases, more mass is shifted to the recalcitrant late-term tailing portion of the Figure 2-3. BTCs for increasing variability in contribution time (m2). BTC, leading to a larger initial reduction in flux for a given reduction in mass. A similar concept used recently by Lemke et al. (2004b) and Christ et al. (2005) is the ganglia to pool ratio (GTP), which describes the ratio of NAPL existing in a residual phase to that existing as higher saturation pools. This is somewhat analogous to describing the NAPL architecture using the Lagrangian approach outlined previously. A high variance lognormal distribution of trajectory integrated NAPL contents is one in

PAGE 38

24 0.3 PV0.6 PV1 PV2 PVcumulative reduction in masscumulative reduction in source discharge Figure 2-4. Mass reduction/flux reduction relationships for BTCs in figure 2-3. which the majority of the NAPL mass is associated with a small number of high streamtubes. In numerical simulations of DNAPL spills conducted by Christ et al. (2005), results indicated that as the variance of the permeability field increased, the GTP ratio decreased, meaning that the NAPL distribution became more skewed towards the presence of NAPL in pools as opposed to residual. The results of Christ et al. (2005) are in agreement with the previous work of Mayer et al. (1996) and Dekker et al. (2001). Numerical simulations of aqueous dissolution have indicated that that as the NAPL distribution becomes increasingly skewed towards more recalcitrant pools, there is increased tailing in the BTC (Parker and Park, 2004; Mayer et al., 1996), which is in agreement with the streamtube models outlined above. In the high resolution simulation S

PAGE 39

25 of Parker and Park (2004), the NAPL architecture was characterized by residual zones of varying lengths (responsible for the large initial decrease in flux) overlying a laterally extensive pool perched on top of a clay aquitard (responsible for long term tailing). It is instructive to notice that m 1 is equal to the contribution time for a step function dissolution profile, that is m 1 is equal to the total mass divided by the initial mass discharge. An interesting property of the streamtube model (lognormal distribution) outlined above is that for a given value of m 1 the range of possible BTCs for different values of m 2 always intersect at a flushing duration equal to m 1 and at a flux-averaged concentration equal to 0.309 (figure 2-3), meaning that the tie line for a flushing duration equal to m 1 is a parallel to the mass reduction axis and equal to a flux reduction of .691 (see figure 2-4 for a flushing duration of 1 PV). An additional property of interest is that as m 2 approaches infinity, the total mass left of the intersection point approaches 0.309. The implication of this is that the lognormal distribution cannot produce mass reduction/flux reduction relationships that extend beyond the .309 mass reduction/(1-.309) flux reduction point in the upper left portion of the mass reduction/flux reduction plot. This is entirely a property of the lognormal distribution and has no physical meaning. Cases of more extreme tailing could be accounted for with alternative probability distributions. Rate Limited Streamtube Model Berglund (1997) presented a stochastic-advective model for rate-limited nonaqueous phase liquid dissolution. tStC T CN with:

PAGE 40

26 CCSktSsN (2-19) Where k is the mass transfer rate coefficient that is linearly related toSwith units of [T -1 ]. The analytical solution of (2-19) is given below in (2-20). 1)(expexp)(exp1NsNsscftTkCtSktTkCCfC (2-20) Advection Dispersion Model Jury (1990), among others, have demonstrated the mathematical similarity of the advection dispersion equation and the lognormal travel time distribution. An alternative to utilizing the cumulative distribution function of for S D (T) would then be the dimensionless analytical solution for the advection dispersion equation (step injection; injection and detection in flux). )(1)(TSCfTCDscf with: )(4)exp(21)(421)(TRRTPerfcPTRRTPerfcTSD (2-21) Where P is the Peclet number and is used here in a manner analogous to the second moment of the contribution time, functioning as variability index to describe the combined effects of hydrodynamics and NAPL architecture. The retardation factor, R, is used to describe the total mass of NAPL in the source zone in a manner similar to the first

PAGE 41

27 moment of contribution time. A useful property of equation 2-21 is that the resulting equation integrates to R. This property is widely known throughout the chemical engineering literature as column holdup time. (2-22) 0)(1RTSD Using equation 3-22, R can be used to represent the total NAPL mass in the source zone. NfSKR (2-23) Where S N is the domain average NAPL saturation. BTCs and Mass Reduction/Flux Reduction A suite of BTCs and mass reduction/flux reduction relationships for a fixed NAPL mass (fixed R) and increased contribution time heterogeneity (increased P) is displayed below in figures 2-5 and 2-6. As shown in figure 2-5 and 2-6, the advection dispersion model produces dissolution profiles and mass reduction/flux reduction relationships that are similar to the streamtube model. One difference of note is that the ADE equation is capable of producing more severe tailing, leading to mass reduction/flux reduction relationships that extend farther into the upper-left corner of the mass reduction/flux reduction plot. Power Function Model Parker and Park (2004) presented a simplified model for estimating DNAPL source zone depletion based on the concept of an effective Damkohler number ( qLkDaeff/( ) where eff [T -1 ] is the effective mass transfer coefficient, L is the source zone length in the mean flow direction, and q [LT -1 ] is the average Darcy flux for the source zone). The

PAGE 42

28 Figure 2-5. BTCs for increasing Peclet number. 0.3 PV0.6 PV1 PV2 PVcumulative reduction in masscumulative reduction in source discharge Figure 2-6. Mass reduction/flux reduction relationships for BTCs in figure 2-5

PAGE 43

29 concept of a effective Damkohler number was also utilized previously by Mayer and Miller (1996). For brevity this approach will be referred to as the Da approach or Da model throughout. The concept of the Da model is to relate changes in the field-scale mass transfer rate coefficient to changes in the global values of NAPL mass and the average groundwater velocity as in equations 2-24 and 2-252 below (Parker and Park, 2004). qLCtCeffsoutexp1)( (2-24) 21)(osoeffMtMKq (2-25) Where C out [ML -3 ] is the flux averaged concentration exiting the source zone, L is the length of the source zone, sK [LT -1 ] is the average saturated hydraulic conductivity of the source zone, M(t) is the NAPL mass in the source zone at time t, M o is the initial NAPL mass in the source zone, and o [T -1 ] 1 and 2 are fitting parameters. Conceptually, there are similarities between the Da approach and the streamtube approaches outlined above. The first similarity noted is that of the 1 parameter. For cases where 1 = 1, as in the high resolution simulation conducted by Parker and Park (2004), the flux-averaged concentration exiting the source zone is independent of the mean groundwater velocity, meaning that changes in contact time between advecting groundwater and the NAPL phase brought about by temporal changes in groundwater velocity do not result in rate-limiting effects. In this sense, the 1 =1 case for the Da model is similar to the equilibrium streamtube model where the flux-averaged concentration exiting the source zone is insensitive to changes in the residence time of the

PAGE 44

30 advecting groundwater. The second conceptual similarity is between o and f c As noted previously, in the streamtube formulation f c is a measure of the fraction of total streamtubes intersecting the contaminated portion of the source zone. Mathematically, both o of f c act to scale the solubility of the NAPL to the initial flux-averaged concentration exiting the source zone. Although o is treated in Parker and Park (2004) as a fitting parameter, it is suggested that the physical meaning of o is similar to f c and is primarily a measure of the fraction of streamtubes intersecting the contaminated portion of the source zone. The third and most important similarity between the Da model and the streamtube models is the relationship between 2 and the way in which the velocity and NAPL distributions are described in the streamtube models. In the Da model, 2 is referred to by Parker and Park (2004) as a mass depletion exponent, which they describe as being related to the groundwater velocity distribution, the DNAPL geometry, and the correlation between the two. The main difference then, between the Da model and the streamtube models, is the specificity in which the velocity and NAPL distributions are described. In the streamtube model, variability in velocity and NAPL saturation along a single streamtube is integrated and transformed into an effective value (contribution time ) for the individual streamtube, with the variability of this distribution controlling dissolution behavior. In the Da model, like the ADE model, variability in velocity and NAPL saturation is integrated not just along the individual streamtube, but over the entire source zone domain and transformed into a single parameter, 2, which relates changes in the effluent concentration to reductions in global NAPL mass. Mathematically, 2 is therefore related to the second moment of the distribution and the Peclet number in the ADE model. Two special cases of the Da model are noted. When 2

PAGE 45

31 is equal to 1, the equation describing the dissolution profile reduces to exponential decay, C(T) = C o e -kT with C o equal to the initial flux-averaged concentration exiting the source zone and k equal to a rate constant that can be solved for in terms of the initial mass and flux. The second case of note is when 2 = 0, in which case the dissolution profile reduces to a step function where the effluent concentration remains constant until all the mass within the source zone has been depleted. For steady flow conditions the Da model can be expressed in a similar fashion to the streamtube and ADE model using equation 3-1 with the source depletion term given below in 2-26. )(1)(TSCfTCDscf with: 2)(exp)(ooDMTMTS (2-26) When o is large, the combination of 2-26 and 2-1 reduces to 2-27 below (Parker and Park, 2004). 2)()(ocsMTMfCTC (2-27) Note that 2-27 is the same formulation as that of Zhu and Sykes (2004). BTCs and Mass Reduction/Flux Reduction BTCs and mass reduction/flux reduction plots for the Da model are shown below in figures 2-7 and 2-8 for a range of 2 Like the Peclet number in the ADE model and the second moment of in the streamtube model, a larger value of 2 connotes a more heterogeneous velocity field/NAPL architecture. The Da model produces similar shaped

PAGE 46

32 Figure 2-7. BTCs for the Da model. 0.3 PV0.6 PV1 PV2 PVcumulative reduction in masscumulative reduction in source discharge Figure 2-8. Mass reduction/flux reduction relationships for Da model

PAGE 47

33 BTCs for 2 > 0.7. For 2 < 0.7, the BTCs are characterized by a single inflection point and are therefore unable to simulate both an initial stage of constant concentration and a later term tailing portion. Comparison of Source Depletion Models with Numerical Simulations In this section the streamtube model and the ADE model are evaluated as mathematical analogues to the multiphase flow and transport simulator UTCHEM (Deslshad et al., 1996). If the source depletion models are robust in terms of being able to represent the range of dissolution profiles observed in field, laboratory, and numerical studies, then an alternative to trying to physically measure relevant parameters is to evaluate a typical range of observed behavior and then incorporate that range/uncertainty into the site-specific decision process. To produce a range of initial architectures in order to evaluate the robustness of the models outlined above, four simulation sets were conducted, each with a different set of parameters describing a spatially correlated random permeability field. Each simulation set contained four realizations. Random fields for the first set were generated using indicator Markov Chain techniques. Random fields for sets 2 through 4 were generated using the turning bands code of Tompson et al. (1989). In sets 2 through 4 the variance of ln(K), where ln(K) is the log transformed hydraulic conductivity in units of cm/s, was sequentially increased from 0.1 to 0.37 to 0.7 to yield progressively more heterogeneous flow fields and NAPL architectures. Constitutive Relationships UTCHEM offers several options for the various constitutive relationships required to model multiphase flow and transport. The decision as to which relationship to utilize is nontrivial as it has a significant impact on simulation results. UTCHEM offers four

PAGE 48

34 different capillary pressure/saturation (Pc/S) models, the Brooks-Corey (BC) model, the Van Genuchten (VG) model, and two versions of a hysteretic VG model, the Parker and Lenhard model (Parker et al., 1987) and the Kalurachchi and Parker model (Kalurachchi and Parker, 1990). The approach used in UTCHEM for incorporating hysteretic effects when utilizing the BC or VG model is to model the NAPL spill event in the main drainage direction for the entire spill. After the primary spill event, the model can be run in restart mode, for any subsequent water or surfactant flooding. For water or surfactant flooding the direction of the capillary pressure/saturation curve is taken to be in the imbibition direction. Hysteretic effects are incorporated by assigning different parameters to the Pc/S model (e.g. lambda parameter for BC) in the imbibition direction than were assigned in the main drainage direction (Delshad et al., 1996). A series of recent papers by Gerhard and Keuper (Gerhard and Kueper, 2003 a,b,c) investigated spill behavior for both the BC and VG capillary pressure and relative permeability models. They concluded that for accurate simulations of NAPL spills, a hysteretic BC model should be utilized. The major limitation of the VG model, which has significant implications to the simulations here, is that the continuous nature of the VG Pc/S curve (meaning it does not have a fixed entry pressure like the BC model) allows for unrealistic penetration of NAPL into finer grained NAPL (Gerhard and Keuper, 2003a). This is clearly inconsistent with virtually all studies of NAPL spills, which have shown a migrating NAPL to be extremely sensitive to small scale changes in permeability (Illangasekare et al., 1995; Kueper et al., 1993; Oostrom et al., 1999; Schwille, 1988; Taylor et al., 2001; Walker et al., 1998). A comparison of the BC and VG models for identical permeability fields and similar shaped Pc/S curves is shown below in figure 2-9.

PAGE 49

35 As shown in figure 2-9, which model is selected has a significant impact on the resulting NAPL architecture and in turn on dissolution dynamics. Van GenuchtenBrooks Corey Figure 2-9. Comparison of NAPL spill using Brooks-Corey and Van Genuchten models. As mentioned previously, one option in UTCHEM is to model the primary spill in the main drainage direction and then to conduct water or surfactant flooding in restart mode with the Pc/S curve now in the direction of imbibition. Several preliminary Figure 2-10. Raining effect observed when running UTCHEM in restart mode.

PAGE 50

36 simulations were conducted to evaluate this option. When running the simulator in restart mode a strange raining effect occurs which was determined to be unrealistic. This raining effect is shown below in figure 2-10. Note that the BC simulation shown in the main drainage direction in figure 2-9 was run in restart mode to generate the two plots spills shown in figure 2-10. In order to generate realistic NAPL spills, the code in UTCHEM should be modified to include a hysteretic BC model (Christ et al., 2005). A linear interpolation approach similar to the Parker and Lenhard model (Parker et al., 1987) and the Kalurachchi and Parker model (Kalurachchi and Parker, 1990) for the BC model was used here. The NAPL spill was conducted in the main drainage direction. The maximum NAPL saturation achieved during main drainage (denoted by the subscript max_sim below) was then used to calculate the corresponding organic residual saturation (S r ). ltheoreticasimrrSSSSmax_max_max_ (2-28) Where the subscript r_max is used to denote the maximum residual saturation and the subscript max_theoretical is the theoretical maximum saturation achievable during main drainage. A theoretical maximum organic saturation (S max_theoretical ) of 0.8 and a maximum organic residual of 0.25 (S r_max ) were utilized in all simulations. A maximum organic residual of 0.25 is towards the upper end of the typical values used in numerical simulations of NAPL spills which incorporate hysteretic effects (Christ et al., 2005; Dekker and Abriola, 2000a; Gerhard and Kueper, 2003a; Lemke et al., 2004a; Mayer and Miller, 1996; Parker and Park, 2004; Rathfelder and Abriola, 1998).

PAGE 51

37 The entry pressure in UCTCHEM is calculated based on Leverett scaling: kChd (2-29) where C is a user defined coefficient which scales the entry pressure (h d ) to the relative permeability (k), with k in units of millidarcies (md), and is the porosity. Simulation Setup Parameters and simulation conditions coincident to all simulations are outlined in this section. The NAPL was introduced into the simulation domain through an injection well at a rate of 15 L/day over a period of 100 days for simulations sets 1 and 2 and 150 days for simulation set 3 and 4. The maximum saturations after the 100 day period were then scaled external to the simulator to the appropriate residual saturation using equation 2-28 above. For certain simulations in set 1 and 2, the permeability field was such that NAPL had reached the lateral no flow boundary at the base of the simulation domain leading to extensive pooling. In those instances the NAPL distribution from a time step less than 100 days was utilized in order to generate dissolution profiles that were representative of the permeability field and not an artificial byproduct of the lateral no flow boundaries. There are certainly cases in the field in which extensive lateral pooling occurs, such as on top of clay aquitard. While recognizing this, the focus here is on NAPL architectures governed by the statistics of the permeability field. A list of parameters utilized is listed below in Table 2-1. The values selected were designed to represent typical values, consistent with recent numerical simulations of multiphase flow and transport (Christ et al., 2005; Dekker and Abriola, 2000a; Gerhard and Kueper, 2003a; Lemke et al., 2004a; Mayer and Miller, 1996; Parker and Park, 2004; Rathfelder and Abriola, 1998).

PAGE 52

38 After the NAPL spill, the dissolution process was initiated through a pair of injection and extraction wells comprising the left and right boundaries respectively. The injection and extraction rate was 10 m 3 /day. Table 2-1. Parameters used in UTCHEM simulations. Parameter Value Notes General Media Properties Porosity 0.3 NAPL fluid properties Density 1.625 g/cm 3 Viscosity 0.89 cp Solubility 1000 mg/L increased to cut down on simulation time for aqueous dissolution Water Properties Density 1 gm/cm 3 Viscosity 1 cp Capillary Pressure C 2 (1/md) 0.5 scales entry pressure to permeability 1.7 BC pore size distribution index S max_theoretical 0.8 S r_max 0.25 Relative Permeability 2.85 exponent relating NAPL saturation to relative permeability NAPL injection well Flow rate 15 L/day Duration of Injection 100 days Injection and Extraction Wells Flow rate 10 m 3 /day Dispersivity Longitudinal 0.15 m Transverse .0005 m Simulation Set 1 Simulation domain and random field generation For simulation set 1, a model domain of 20m x 10m x 1m was utilized. The grid discretization selected was 10cm in the z direction and 20cm in the x direction. Correlation lengths were assumed to be 4m is the x direction and 0.25 m in the z direction. Grid parameters are summarized in table 2-2.

PAGE 53

39 Table 2-2. Grid and permeability field parameters. Value x dimension 20 m z dimension 10 m y dimension 1 m discretization x 20 cm discretization z 10 cm discretization y 1 m correlation length x 4 m correlation length z 0.25 m Random fields for the four simulations of set 1 were generated using Transition Probability Geostatistical Software (T-Progs). T-Progs is a Markov chain based geostatistical package that is designed to simplify the process of interpreting geostatistical data. For actual sites, the user would be required to delineate the porous media into five or fewer classes, inputting the sampling location and class into the simulator. T-Progs then calculates the volume percentage of each class, transition probabilities, correlation lengths and other relevant parameters and generates a user specified number of realizations that are conditioned to the measured data. For theoretical simulations, such as those presented here, the only input data required is the volume percentage of each class and the correlations lengths. One advantage of the indicator approach is that specific values for each class can be specified (e.g. entry pressure, pore size distribution, residual saturation), as opposed to having to come up with scaling relationships. This approach was employed by Lemke and collaborators (Lemke and Abriola, 2003; Lemke et al., 2004a) at a DNAPL contaminated site in Michigan. Parameters for each class are listed below in table 2-3.

PAGE 54

40 Table 2-3. Parameters assigned to each indicator class. Class 1 Class 2 Class 3 Class 4 Class 5 Volume Fraction 0.1 0.5 0.25 0.1 0.05 Permeability (md) 8000 5000 3000 1000 100 Hydraulic Conductivity (m/day) 7.7 4.8 2.9 1 .01 BC Lambda Parameter 2.1 1.9 1.4 1.1 1 Max Residual Oil Saturation .21 .22 .23 .24 .25 NAPL architectures NAPL architectures for the four set 1 realizations are displayed below in figure 2-11. The plots display NAPL saturations contoured over a range of 0 to 25% saturation. Simulation 3Simulation 4 Simulation 1 Simulation 2 Figure 2-11. NAPL spills for set 1 (contour range: 0 to 25% saturation). Dissolution profiles, model fits, and mass reduction/flux reduction relationships The UTCHEM generated dissolution profiles for the four simulations along with model fits for the streamtube and ADE model are shown below in figure 2-12 and 2-13. The R 2 value and the fitted parameters (P and 2 ) are also shown. The objective function

PAGE 55

41 selected for curve fitting was a simple root mean square function given by equation 2-30 below. 21UTCHEMSDUTCHEMCCCRMS (2-30) Where the subscript UTCHEM denotes the UTCHEM model output and SD denotes the fit value from the source depletion models (ADE and streamtube). A value of 1 was added to the denominator to avoid dividing by zero. As evidenced in figure 2-12 and 2-13, both models provide fits of the UTCHEM generated dissolution data with R 2 values > 99%, suggesting that these simplified source depletion models can serve as mathematical analogs for more sophisticated numerical simulators such as UTCHEM. Despite differences in NAPL architectures (figure 2-11), the dissolution profiles are quite similar, suggesting the variance of the distribution and the Peclet number may be highly correlated to the statistics of the porous media. There is a strong conceptual foundation for this as the permeability field is the primary governing mechanism of both the NAPL architecture (e.g., Dekker and Abriola, 2001a; Mayer and Miller, 1996) and velocity field (e.g. Dagan, 1989; Gelhar, 1993). The mass reduction/flux reduction plots for the simulations of set 1 along with the model fits from the streamtube and ADE model are shown below in figures 2-14 and 2-15 (Note that the dissolution profiles and not the mass reduction/flux reduction plots are being fit). The perturbation about the mean dissolution behavior, which is somewhat evident in figures 2-12 and 2-13 is more apparent in the mass reduction/flux reduction plots below as the perturbation occurs over a larger mass reduction interval than a flushing duration interval. This is especially true for simulation 4.

PAGE 56

42 2= 0.58R2 = 0.999 2= 0.49R2 = 0.999Simulation 3Simulation 4Simulation 1Simulation 22= 0.54R2 = 0.9982= 0.72R2 = 0.990Cumulative Pore VolumesC/Co Source Depletion ModelUTCHEM Simulation Figure 2-12. Streamtube model fit of UTCHEM dissolution data (set 1). P = 5.3R2 = 0.999 P = 7.7R2 = 0.999 P = 7.3R2 = 0.998Simulation 3Simulation 4Simulation 1Simulation 2P = 4.9R2 = 0.990Cumulative Pore VolumesC/Co Source Depletion ModelUTCHEM Simulation Figure 2-13. ADE model fit of UTCHEM dissolution data (set 1).

PAGE 57

43 Simulation 3Simulation 4Simulation 1Simulation 2 Source Depletion ModelUTCHEM SimulationCumulative Reduction in MassCumulative Reduction in Source Discharge Figure 2-14. Mass reduction/flux reduction fit for streamtube model (set 1). Simulation 3Simulation 4Simulation 1Simulation 2 Source Depletion ModelUTCHEM SimulationCumulative Reduction in MassCumulative Reduction in Source Discharge Figure 2-15. Mass reduction/flux reduction fit for ADE model (set 1).

PAGE 58

44 Simulation Domain and Grid Setup for Simulation Sets 2 through 4 For simulation sets 2 through 4, a model domain of 30 m x 10 m x 1 m was utilized. For simulation set 3 (which was the first simulation set run), the grid discretization selected was 5cm in the z direction and 50 cm in the x direction. Grid discretization was determined based on the recommendation of Russo (1991) (cited in: Mayer and Miller, 1996) who suggest at least 4 nodes per correlation length and 15 correlation lengths over the entire domain to maintain ergodicity. The same sized domain used for simulation set 3 was also used for simulation set 2 and simulation set 4. The only difference between simulation set 3 and sets 2 and 4 was that the grid discretization was increased from 5cm in the z direction to 8 cm. This resulted in a factor of 5 decrease in computation time. For example, the dissolution simulations, which took approximately 30 hrs of run time for simulation set 3, were decreased to around 6 hrs for simulation set 4. Random fields were generated using the turning bands code of Tompson (1989). Table 2-4. Grid/random field parameters for simulation sets 2-4. Set 2 Set 3 Set 4 Units X dimension 30 30 30 m Y dimension 1 1 1 m Z dimension 10 10 10 m grid cells X 60 60 60 grid cells y 1 1 1 grid cells Z 125 200 200 dx 50 50 50 cm dy 100 100 100 cm dz 8 5 5 cm variance ln(K) 0.1 0.36 0.7 cm/s mean ln(K) -5.62 -5.62 -5.62 cm/s correlation length x 2 2 2 m correlation length z 0.2 0.2 0.2 m Geostatistical parameters for set 3, which functioned as the base case, were adapted from Jensen et al. (1993). For simulation set 2 the variance of ln(K) was decreased from the

PAGE 59

45 base case value of 0.36 to 0.1. For simulation set 4, the variance of ln(K) was increased from the base case value of 0.36 to 0.7. A summary of the grid/random field parameters used for simulation sets 2 through 4 is displayed below in table 2-4. Simulation Set 2 NAPL architectures NAPL architectures for the four set 2 realizations are shown below in figure 2-16. As shown in figure 2-16, the low variability in the permeability field leads to a relatively uniform NAPL architecture. Simulation 3Simulation 4Simulation 1Simulation 2 Figure 2-16. NAPL spills for set 2 (contour range 0 to 25% saturation).

PAGE 60

46 Dissolution profiles, model fits, and mass reduction/flux reduction relationships The UTCHEM generated dissolution profiles for set 2 simulations along with model fits of the dissolution profile and mass reduction/flux reduction relationship are displayed below in figures 2-17 through 2-20. As in simulation set 1, the models fits of the UTCHEM generated dissolution data were > 99%. Similar to simulation set 1, the dissolution profiles for all simulations in set 2 are quite similar, suggesting the Peclet number and variance of contribution time may be highly correlated to the permeability field statistics. 2= 0.43R2 = 0.9992= 0.4R2 = 0.998Simulation 3Simulation 4Simulation 1Simulation 22= 0.5R2 = 0.9992= 0.46R2 = 0.999Cumulative Pore VolumesC/Co Source Depletion ModelUTCHEM Simulation Figure 2-17. Streamtube model fit of UTCHEM dissolution data (set 2).

PAGE 61

47 P = 9.9R2 = 0.999P = 11.9R2 = 0.998P = 7.2R2 = 0.999Simulation 3Simulation 4Simulation 1Simulation 2P = 8.6R2 = 0.999Cumulative Pore VolumesC/Co Source Depletion ModelUTCHEM Simulation Figure 2-18. ADE model fit of UTCHEM dissolution data (set 2). Simulation 3Simulation 4Simulation 1Simulation 2 Source Depletion ModelUTCHEM SimulationCumulative Reduction in MassCumulative Reduction in Source Discharge Figure 2-19. Mass reduction/flux reduction fit for streamtube model (set 2).

PAGE 62

48 Simulation 3Simulation 4Simulation 1Simulation 2 Source Depletion ModelUTCHEM SimulationCumulative Reduction in Source DischargeCumulative Reduction in Mass Figure 2-20. Mass reduction/flux reduction fit for ADE model (set 2). Simulation Set 3 NAPL architectures NAPL architectures for the four set 3 realizations are displayed below in figure 2-21. The plots display NAPL saturations contoured over a range of 0 to 25%. Consistent with previous numerical simulations of NAPL migration, figure 2-21 illustrates that an increase in the variability of the permeability field leads to a decrease in the total penetration depth and a distribution more skewed towards the presence of NAPL in pools (Christ et al., 2005; Dekker and Abriola, 2000a; Mayer and Miller, 1996).

PAGE 63

49 Simulation 3Simulation 4Simulation 1Simulation 2 Figure 2-21. NAPL spills for set 3 (contour range: 0-25% saturation). Dissolution profiles, model fits, and mass reduction/flux reduction relationships The UTCHEM generated dissolution profiles for the set 3 simulations along with model fits for the streamtube and ADE model are shown below in figure 2-22 and 2-23. The models perform well, with fits characterized by R 2 values >= 99%. The trend of decreased penetration depth and an initial architecture more skewed towards NAPL

PAGE 64

50 existing in pools effects the flux averaged BTC in two important ways, both of which lead to increased longevity of the source zone. First, with a decrease in penetration depth, the fraction of streamtubes contaminated with NAPL, f c decreases. For a given amount of NAPL mass, as f c decreases, the amount of mass in each streamtube increases which of course leads to an increase in average contribution time for each streamtube. Second, a NAPL architecture more skewed towards pools is equivalent to an increase in the variability in the trajectory averaged NAPL content,. As noted previously this highly variable architecture causes substantial tailing in the BTC as the lowstreamtubes are depleted of mass relatively early with the laterally extensive, high streamtubes continuing to produce a small portion of the initial flux for several additional pore volumes. When comparing the dissolution profiles of simulation set 3 to simulation set 2, the general trend predicted by the streamtube models of increased spreading of the BTC with increased heterogeneity in the combined effects of the Lagrangian NAPL architecture and velocity field (lower Peclet number/higher variance) is evident (see figures 2-3 and 2-5). S S S The dissolution profile of simulation 4 proved somewhat more difficult to fit with the simplified source depletion models. The NAPL architecture for simulation 4 was characterized by an extensive pool in the vicinity of the injection point (figure 2-24 below), leading to a discontinuous distribution with the majority of the streamtubes associated with mixed length residual components (comparatively smaller) and a small fraction of streamtubes associated with the laterally extensive pool component (much higher ). The discontinuity arises from a large gap in values from the residual components and the pooled components, i.e. the distribution is relatively smooth and S S S S

PAGE 65

51 2= 0.98R2 = 0.9982= 1.03R2 = 0.992Simulation 3Simulation 4Simulation 1Simulation 22= 0.86R2 = 0.9942= 0.924R2 = 0.989Cumulative Pore VolumesC/Co Source Depletion ModelUTCHEM Simulation Figure 2-22. Streamtube model fit of UTCHEM dissolution data (Set 3). P = 1.5R2 = 0.998P = 1.3R2 = 0.994P = 2.12R2 = 0.994Simulation 3Simulation 4Simulation 1Simulation 2P = 1.8R2 = 0.989Cumulative Pore VolumesC/Co Source Depletion ModelUTCHEM Simulation Figure 2-23. ADE model fit of UTCHEM dissolution data (Set 3).

PAGE 66

52 continuous until all the residual mass has been depleted upon which there is a large jump to the next highest streamtube. As noted previously this highly variable architecture causes substantial tailing in the BTC as the lowstreamtubes are depleted of mass relatively early with the laterally extensive, high streamtubes continuing to produce a small portion of the initial flux for several additional pore volumes. To illustrate the discontinuous nature of the BTC, the fit in figure 2-24 was obtained by fitting the first portion of the BTC associated with the residual component and the portion associated with the pooled component each with a separate Peclet number. The respective architectures for each component of the BTC are also inset into figure 2-24. This type of dissolution behavior, with two distinct components, one associated with the mixed length residual component and the other with a laterally extensive pool component, was also evident in simulation 1 of set 4. S S S Pore VolumesC/Cs UTCHEM dissolution data distribution fit to mixed length residual componentdistribution fit to pool component Figure 2-24. Discontinuous fit of simulation 4 (set 3).

PAGE 67

53 Simulation 3Simulation 4Simulation 1Simulation 2 Source Depletion ModelUTCHEM SimulationCumulative Reduction in Source DischargeCumulative Reduction in Mass Figure 2-25. Streamtube model fit of mass reduction/flux reduction relationships. Simulation 3Simulation 4Simulation 1Simulation 2 Source Depletion ModelUTCHEM SimulationCumulative Reduction in Source DischargeCumulative Reduction in Mass Figure 2-26. ADE model fit of mass reduction/flux reduction relationships.

PAGE 68

54 The difficulty in fitting the extensive late term tailing behavior is further magnified in the mass reduction flux reduction plots (figure 2-24; 2-25). The discontinuous nature of the simulation 4 NAPL architecture is more enhanced in the mass reduction/flux reduction plots where there is a sharp jump to the late term tailing behavior once the residual component has been depleted. Simulation Set 4 NAPL architectures NAPL architectures for the set 4 realizations are displayed below in figure 2-27. The plots display NAPL saturations contoured over a range of 0 to 25%. When compared with the simulation sets 2 and 3, simulation set 4 architectures are characterized by a decrease in the total penetration depth and a distribution more skewed towards the presence of NAPL in pools. Dissolution profiles, model fits, and mass reduction/flux reduction relationships The UTCHEM generated dissolution profiles for the set 4 simulations, along with model fits for the streamtube and ADE model, are shown below in figure 2-28 and 2-29. The models provide reasonably good fits of the UTCHEM simulations, yet not as accurate as the more homogenous systems (sets 1 through 3). Comparing the dissolution profiles of simulation set 4 to the previous simulation sets, the trend of decreased penetration depth and enhanced pooling increasing source longevity is again evident. The general trend predicted by the streamtube models of increased spreading of the BTC with increased heterogeneity in the combined effects of the Lagrangian NAPL architecture and velocity field (lower Peclet number/higher variance) is also evident. Although there is the expected trend of increased spreading of the BTC, the variability between individual

PAGE 69

55 realizations also increases (e.g. compared with set 2 where all of the BTCs were quite similar). This is consistent with general probability theory that suggests that as the Simulation 3Simulation 4Simulation 1Simulation 2 Figure 2-27. NAPL spills for set 4 (contour range: 0 to 25% saturation) variance of the permeability field increases, more simulations would be required to converge to the expected value, i.e. mean BTC for the particular permeability field

PAGE 70

56 statistics. An exception to the increased variability trend was simulation 3, where unlike previous simulations where the late term tailing was dominated by a single large pool in the vicinity of the injection point, the permeability field was such that three pools of comparable size formed. The effect of these three pools of comparable size was to dampen variability in the dissolution profile, such that the BTC more resembled the simulations of set 1 and 2. The initial NAPL architecture and dissolution profile of simulation 1 was quite similar to simulation 4 of set 3, with two distinct BTC components, one associated with a residual component, and another associated with a late term tailing portion as the large pool was gradually depleted of mass. 2= 0.97R2 = 0.9822= 1.16R2 = 0.994Simulation 3Simulation 4Simulation 1Simulation 22= 0.31R2 = 0.9942= 1.3R2 = 0.937Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 2-28. Streamtube model fit of UTCHEM dissolution data (set 4) Mass reduction flux reduction plots for simulation set 4 are shown below in figures 2-30 and 2-31. Difficulty in fitting the late term tailing portion of the BTC is evident in

PAGE 71

57 P = 1.5R2 = 0.981P = 1.0R2 = 0.996P = 20.67R2 = 0.995Simulation 3Simulation 4Simulation 1Simulation 2P = .754R2 = 0.937Cumulative Pore VolumesC/Co Source Depletion ModelUTCHEM Simulation Figure 2-29. ADE model fit of UTCHEM dissolution data. Simulation 3Simulation 4Simulation 1Simulation 2 Source Depletion ModelUTCHEM SimulationCumulative Reduction in Source DischargeCumulative Reduction in Mass Figure 2-30. Streamtube model fit mass reduction/flux reduction relationships (set 4).

PAGE 72

58 simulation 1 and simulation 3 where there is a large jump to the late term tailing behavior as in simulation 4 of set 3. Simulation 3Simulation 4Simulation 1Simulation 2 Source Depletion ModelUTCHEM Simulation Cumulative Reduction in Source DischargeCumulative Reduction in Mass Figure 2-31. ADE model fit of mass reduction/flux reduction relationships (set 4). Conclusions and Recommendations for Further Research In this chapter simplified source depletion models which may serve as mathematical analogs to more sophisticated multiphase flow and transport simulators where developed. The models are conceptually appealing as a single parameter defines the total NAPL mass and a single parameter defines the combined variability of the Lagrangian NAPL architecture and the velocity field. The fundamental assumption is that dissolution dynamics are primarily governed by the Lagrangian NAPL architecture and the velocity field, whereby processes such as relatively permeability, transverse dispersion and mass transfer rate coefficients, when neglected, introduce negligible error. The general trend predicted by the models is with increased variability in the NAPL

PAGE 73

59 architecture and/or the velocity field, the dissolution profile exhibits more spreading and more severe tailing, leading to increased source longevity and a more favorable relationship between reductions in mass and reductions in source discharge. The source depletion models were compared with UTCHEM simulations, comparing favorably, especially when considering a sophisticated numerical simulator with run times in the vicinity of 30 hrs has been replaced with a relatively compact analytical solution. Consistent with previous numerical investigations, increased variability in the permeability field led to increased pooling, which in turn led to enhanced tailing of the dissolution profile, suggesting that the Peclet number/variance of contribution time (the parameters that control spreading in the simplified source depletion models) are correlated with the variance of ln(K). The higher variance ln(K) simulations would benefit from further simulations to come closer to arriving at a mean dissolution profile for a given set of permeability field statistics. More simulations would perhaps yield an empirical estimator of the Peclet number or variance of contribution time based on the statistics of the permeability field. A more in depth investigation of the impact of neglecting relative permeability, transverse mixing between streamtubes, and mass transfer rate coefficients is warranted. The impact of transverse mixing in relation to bioreactive transport has been investigated by Cirpka and Kitanidis (2000) and Ginn (2001). Berglund and Fiori (1997) studied the impact of transverse mixing on a sorbing solute. In general, the limited amount of work directed towards streamtube modeling of NAPL dissolution processes (Berglund, 1997; Jawitz et al., 2005) has been of the theoretical modeling type and has not investigated the error associated with the underlying assumptions. I have attempted to address the impact

PAGE 74

60 of the underlying assumptions somewhat here but this is a potentially large area for further research. Investigating bimodal distributions and other types of probability distributions, especially for the higher variance ln(K) simulations is also an area for further research.

PAGE 75

CHAPTER 3 LABORATORY STUDIES OF RELATIONSHIPS BETWEEN ARCHITECTURE AND FLUX IN MODEL SOURCE ZONES Introduction Of critical importance when considering plume management vs. source zone remediation is the expected performance of various source zone remediation technologies under consideration, both in terms of the expected reduction in mass (and the effect on long term plume management costs) and the expected reduction in flux (and the effect on near term risk reduction). Remedial forecasting in terms of flux reduction performance is difficult in that virtually all field scale demonstrations and laboratory studies of aggressive source zone remediation to date have focused on a given technologys ability to remove mass from the source zone and have not quantified the associated reduction in contaminant mass flux. A flux-based evaluation of previous source zone remediation technologies is further complicated by the fact that relationships between reductions in mass and flux are potentially highly variable depending on both the site characteristics and the technology selected. In the case of technologies designed to remove contaminant mass from the subsurface (e.g. surfactant or cosolvent flushing), continuous measurements of mass reduction are available via sampling of the remedial fluid for the desired constituent. In contrast to mass reduction, there are typically only two points where measurements of flux are available (again note that these measurements are typically not taken); an initial measurement before the given technology is implemented and another measurement after 61

PAGE 76

62 the completion of remediation when the system has returned to a condition thought to be representative of natural flowing groundwater. Additional flux reduction data points would require the iterative process of removing a portion of the source zone mass via the given technology and then allowing the system to return to natural flowing groundwater conditions where a flux measurement could be taken. Aqueous dissolution experiments on the other hand, because they do not require switching back and forth between remediation conditions and natural flowing groundwater, have the capability of yielding nearly continuous mass reduction /flux reduction data. The experiments presented here can serve as a pump and treat baseline for comparison with bench scale experiments of aggressive remediation technologies. More challenging is directly interpreting the near continuous mass reduction/flux reduction relationships of aqueous dissolution in order to lend insight to the mass reduction/flux reduction relationships of systems in which aggressive source zone remediation has resulted in partial mass removal. The limited body of work that has addressed the issue of reductions in flux brought about by partial reductions in contaminant mass from a NAPL source zone has been of the theoretical modeling type and directed towards systems undergoing aqueous dissolution (Sale and McWhorter, 2001; Falta, 2003; Rao and Jawitz, 2003; Parker and Park, 2004), sometimes with the implicit assumption that the conclusions reached regarding mass reduction/flux reduction relationships are transferable to other types of aggressive source zone remediation (note that Rao and Jawitz (2003) consider in situ flushing with a mathematical framework that would be identical for aqueous dissolution; i.e. the flushing solution travels along identical flow paths to that of natural flowing groundwater). Based on a steady-state mathematical formulation, Sale and McWhorter

PAGE 77

63 (2001) conclude that the vast majority of the contaminant mass must be removed to achieve significant reductions in flux. Rao and Jawitz (2003) suggest that a more favorable relationship between reductions in mass and flux may be realized with increased heterogeneity of the velocity field. Falta (2003) conducted a high resolution transient simulation of the NAPL architecture considered by Sale and McWhorter (2001) and found reductions in contaminant flux were realized for partial reductions in DNAPL mass. Parker and Park (2004) suggest a continuum of relationships between mass and flux are possible depending on the site specific NAPL distribution, groundwater velocity field, and correlation between the two. It is noted here that depending on the given technology and the way in which the technology is implemented (e.g. vertical circulation vs. line drive pumping), the order and way in which mass is removed from the source zone and the associated reduction in flux could be significantly different when compared to a system which is undergoing aqueous dissolution. It is therefore inappropriate to extrapolate mass reduction/flux reduction relationships directly from aqueous dissolution work to all types of aggressive source zone technologies. Aqueous dissolution work can however be used to gain further understanding as to the way in which the NAPL architecture must be altered in order to bring about reductions in flux. That is, aqueous dissolution experiments can be constructed in a way in order to investigate which components of the initial source zone mass, if efficiently targeted and removed via a given technology, would yield the largest (or any) reduction in flux. Accordingly, the experiments presented here are designed for three major purposes: first, for comparison to the theoretical work briefly discussed above with a specific focus on the role that NAPL architecture plays in controlling

PAGE 78

64 dissolution behavior; second, to serve as a baseline of source depletion under natural flowing groundwater (or pump and treat) from which to compare flushing-based remediation experiments currently being conducted in similar systems; and third, to evaluate simplified modeling approaches for simulating long-term source depletion. Modeling Approach In this chapter, the streamtube models discussed in the previous chapter are parameterized in a deterministic manner by dividing the source zone into n streamtubes and assigning an estimated value of to each streamtube based on image analysis of the NAPL distribution. The superposition of n solutions of the model for an individual streamtube is then used to model the dissolution profile for the entire system. The dissolution profile for the equilibrium streamtube model was estimated using: S (3-1) )()(1iniscTHCfTC where H is the heaviside step function. The dissolution profile for the rate-limited streamtube model was estimated in a similar fashion. niNisiiNissctTkCtSktTkCCfTC11)(expexp)(exp1)( (3-2) For comparison purposes, the experimental results were also modeled using the Da model and the ADE model using traditional curve fitting techniques. Experimental Methods A total of four aqueous dissolution experiments were conducted in two-dimensional flow cells. Three experiments were conducted with 1,2-Dichloroethane (DCA) and one experiment with trichloroethene (TCE). A summary of the experimental

PAGE 79

65 conditions can be found in table 3-1. For brevity the four experiments will hereafter be referred to as DCA-1, DCA-2, DCA-3, and TCE-1. Tabl e 3-1. Summary of experimental conditions for conducted experiments. Experiment 1 Experiment 2 Experiment 3 Experiment 4 Name DCA-1 TCE-1 DCA-2 DCA-3 Volume injected (ml) 10 10 12 10 Injections single single double single Box design segmented segmented segmented in-line Packing discrete layers discrete layers continuous layering fine gradation Pore water velocity (cm/hr) 22.9 18.7 18.2 17.3 Flow Cell Design Two flow cell designs were utilized for the following experiments. The first design featured a segmented extraction well that was designed to provide spatially resolved flux information about temporally changing up-gradient NAPL distributions. The segmented flux measurements, coupled with visualization of the up-gradient NAPL distribution provided a means of investigating the relationship between temporal changes in NAPL architecture and reductions in flux. The design for the segmented box was similar to that of Jawitz et al. (1998a) where a frame of 1.5cm square aluminum tubing was enclosed by two sheets of glass. The 32 cm high by 40 cm long frame, which also served as an injection and extraction well, was slotted at a width of 0.03 cm and a frequency of 4 slots per cm. The integrated injection well was connected via Teflon tubing and a switch valve to either a constant head reservoir or a syringe pump, which was used for tracer injections. The extraction well was segmented at intervals of 4 cm. Teflon tubing originating from each individual well segment was terminated at a constant elevation sampling apparatus. A tray of gas chromatography vials was inserted into the apparatus to

PAGE 80

66 collect simultaneous samples from each segment. When not sampling, effluent was collected in a fraction collector and routed to a series of waste bottles, one for each segment. The individual waste bottles were weighed periodically for an additional flow measurement and for mass balance purposes. The well segments were numbered sequentially from the bottom of the box upwards (i.e. the bottom segment is referred to throughout as port 1 where the top port is referred to as port 7). Segmented WellMedia heterogeneities Constant head reservoir Syringe pump for tracer injectionsValve x Syringe pump injection Integrated well Constant elevation sampling apparatus for segmented water and solute flux measurements Waste routed to collection bottles. One bottle per segment Figure 3-1. Schematic representation of segmented box.

PAGE 81

67 An additional experiment was conducted in a 61 cm long x 41 cm high flow cell with an integrated extraction well and inline sampling system for gas chromatography (GC) analysis (Jawitz et al., 2002). A piston pump was used to continuously deliver effluent to a modified autosampler tray. Sub-samples from the flowing effluent were injected into the GC at a user defined frequency, eliminating the need for manual sampling of the effluent. The temporal discretization of the sampling could be as small as the time required for the GC to analyze an individual sample. The inline system enabled automation at the expense of spatially resolved effluent information and was therefore used only for experiment DCA-3, where due to the porous media packing, the NAPL was expected to be highly concentrated at a particular elevation. Porous Media Packing and NAPL Injection The well-characterized Accusand (Unimin Corp.) was used in all experiments (Schroth et al., 1996). Sand was introduced into the box through a funnel. A small depth of standing water and continuous vibration were maintained throughout the packing process to ensure a saturated system. Three different types of packs were used in the four experiments. In experiments DCA-1 and TCE-1 discrete low permeability lenses of 40/50 and 50/70 sands were inserted into otherwise homogenous 20/30 sand. In experiment DCA-2, larger, more continuous layers of low permeability (40/50 sand) were developed. In experiment DCA-3, four different sand types (20/30, 30/40, 40/50, and 50/70) were mixed together to create a homogenous mixture. The mixture was introduced into the box under standing water while the funnel was slowly moved back and forth across the box. The different settling velocities of the particles produced fine gradations of alternating grain sizes that spanned the width of the box. In all four experiments a syringe pump was used to inject 10 mL of NAPL into the upper portion of the domain at a rate of

PAGE 82

68 0.5mL/min. In addition to the initial 10mL NAPL injection, 2mL of NAPL was injected into one of the low permeability layers in experiment DCA-2. It was expected based on experiments conducted in similar systems that NAPL migration and eventual distribution would be exclusive to the higher permeability 20/30 sand (Oostrom et al., 1999; Taylor et al., 2001) meaning that despite the existence of heterogeneity in the porous media, the NAPL would be present in a relatively uniform velocity field. The presence of all of the NAPL in the higher permeability background media introduces two questions. The first question is whether there are mechanisms for NAPL entry into lower permeability media that may be prevalent at a given field site that are not represented in the experimental system used here and secondly if there are mechanisms that are not represented in the experimental system, should and how could those mechanisms be best represented in the experimental system. One possible explanation for DNAPL entry into lower permeability media is a system that transitions from unsaturated to saturated conditions with the initial NAPL migration occurring under unsaturated conditions and thus promoting retention of the NAPL (which is assumed to be the wetting fluid with respect to air) in lower permeability media (Illangasekare et al., 1995). This was the justification used by Brusseau et al. (2000) who utilized an emplaced source technique where the NAPL and the lower permeability media were mixed external to the experimental system prior to emplacement of the low permeability layers in the flow cell. Another possible explanation for DNAPL entry into a lower permeability layer is a system in which portions of the domain exhibit wettability properties such that NAPL is wetting with respect to water instead of vice versa (Bradford et al., 2003). In such cases the experimental system may be modified by altering the wettability properties of

PAGE 83

69 selected media (Bradford et al., 1999). A third possible explanation for DNAPL entry into lower permeability layers is perhaps tenable due to scaling issues resulting from a possibly large discrepancy in spill rates between a field site and laboratory prototype. Numerical simulations conducted by Kueper and Gerhard (1995) and Dekker and Abriola (2000) suggest that increases in spill rates lead to larger capillary pressures achieved during migration of the NAPL, which in turn allow for entry into increasingly less permeable media. It is therefore reasonable to speculate that for a field site with a similar grain size distribution and a several order of magnitude larger NAPL spill rate that NAPL migration and eventual distribution would not be exclusive to the highest permeability sand. To circumvent the exclusive migration of the NAPL in the highest permeability sand, an additional 2ml of NAPL were injected into one of the lower permeability layers in experiment DCA-2 in order to create flow variability in the NAPL contaminated portion of the flow cell. Mass balances for all experiments were in the range of 80-90%. Possible sources of mass balance errors include volatilization during sampling, inaccuracies in the measurement of flow, and residual NAPL remaining in the flow cell upon completion of the experiment. Mass removal as a function of time and total initial mass were determined by integrating the breakthrough curve without attempting to correct or add back in the various experimental mass losses. Image Analysis A simple image analysis technique was utilized to provide semi-quantitative information of NAPL architecture. Digital images of the system were taken at each sampling acquisition. Note that images were taken from only one side of the flow cell but compared favorably (in terms of similar NAPL distributions) with the side of the flow

PAGE 84

70 cell that was not imaged. The images were imported into an image processing program and overlain with a finely discretized grid. A binary map of the NAPL distribution was developed based on the presence (in which case the grid cell was filled black) or absence (in which case the grid cell was filled white) of NAPL in each individual grid cell. The black and white image was then imported into Mathcad using an image processing function which returns a matrix containing one of two possible numeric values, one value corresponding to a white grid cell (no NAPL) and another value corresponding to a black grid cell (NAPL). Trajectory integrated NAPL contents were estimated by assuming that flow was completely horizontal, such that a trajectory could be represented by one row in the binary matrix. A simple algorithm was used to count the number of elements in each row of the matrix that had a numeric value corresponding to the positive identification of NAPL. This process returns an array containing the total number of NAPL containing grid cells in each trajectory. By assuming a uniform saturation, the mass fraction in each trajectory was then estimated by dividing the number of NAPL-containing grid cells in each row by the total number of NAPL-containing grid cells. The mass fractions were converted to values of and inserted into equation 3-1 and 3-2 to model source depletion. The contaminated fraction, f S c was also estimated from this procedure. Experimental Results, Modeling, and Discussion NAPL Migration and General Dissolution Behavior The migration of the nonaqueous phase and the eventual distribution for experiments DCA-1, DCA-2, and TCE-1 (see figures 3-2a, 3-3a, and 3-4a) was as expected: the NAPL migrated vertically until a low permeability layer was encountered

PAGE 85

71 upon which the NAPL pooled and moved horizontally, eventually spilling over the side where it again began to migrate vertically downward until encountering another low permeability unit. This migration behavior led to an eventual distribution where the NAPL resided entirely in the 20/30 sand. Based on the observed dissolution behavior, the saturation of the pooled NAPL perched on top of the low permeability layers and at the base of the flow cell was concluded to not be high enough to significantly restrict flow through the interior portion of the pool due to reductions in relative permeability. An exception was the pool that formed up-gradient of well segments 5 and 6 in the vicinity of the injection point for experiment TCE-1. Again based on the dissolution behavior, advecting water was primarily restricted to the perimeter of this pool resulting in the thinning of the pool in the vertical direction and shortening of the pool in the mean direction of flow with increased dissolution time. While flow of the aqueous phase was not restricted to the exterior of the NAPL perched on top of the low permeability layers (with the exception of the one pool in the TCE experiment), this is not to suggest that relative permeability effects were insignificant. In fact, several of the segments experienced temporary flux increases as the dissolution process progressed in time, behavior likely attributable to increases in relative permeability and the exposure of additional NAPL interfaces to the flowing aqueous phase (Geller and Hunt, 1993; Powers et al., 1998; Nambi and Powers, 2000). These flux increases were somewhat masked in the integrated BTCs as other portions of the source zone were simultaneously experiencing flux decreases in response to depleting mass. In experiment DCA-2, the lower permeability layers into which the NAPL was injected were intentionally aligned with selected effluent well segments. The flow rates

PAGE 86

72 out of these segments were approximately 0.7 ml/min where the flow rates out of the well segments aligned with the higher permeability background media (20/30 sand) were approximately 0.9 ml/min. The discrepancy in flow rates between the respective ports was smaller than expected based on hydraulic conductivities of 15.0 cm/min for the 20/30 sand and 4.3 cm/min for the 40/50 sand reported by Schroth et al. (1996). It was expected that there would be detectable flow bypassing around the low permeability zone injected with NAPL, causing the NAPL to persist within the low permeability unit for extended periods of time. The small discrepancy in flow rates did not indicate a large degree of bypass flow around this zone possibly due to flow convergence at the extraction well. The contaminant fluxes out of the ports associated with the low permeability zone were characterized by a delayed flux increase, presumably due to an increase in the relative permeability as the up-gradient NAPL mass was depleted, enhancing dissolution of the down-gradient NAPL in the low permeability zone and leading to a gradual increase in the relative permeability. Flux Response to Changes in NAPL Architecture The flux plane response to temporal changes in the NAPL architecture is shown in figures 3-2 to 3-4. The bar graphs plot the fraction of the total initial flux discharged from each extraction well segment. Results are shown for six selected time steps. The fraction of the total mass remaining (MF) and the fraction of the total initial flux (FF) are also displayed for each time step. The segmented flux data is intended to be coupled with the up-gradient NAPL distribution to investigate the linkages between NAPL architecture and flux. For discussion purposes it is useful to consider each well segment (or port) as corresponding with a bundle of streamtubes with a given distribution of trajectory integrated NAPL contents. It is also useful to recall that based on the image analysis

PAGE 87

73 procedure and the definition of, the more elongated the NAPL morphology in the mean flow direction, the larger the estimated value of the trajectory integrated NAPL content,(e.g., a NAPL pool would have a much larger estimated value of than a vertical finger). S S S Figure 3-2. Flux plane response to changes in NAPL architecture for experiment DCA-1. Note that MF is the mass fraction ((mass at time t)/(initial mass at t = 0)) and FF is the flux fraction ((flux at time t) /(initial flux at t = 0)).

PAGE 88

74 Figure 3-3. Flux plane response to changes in NAPL architecture for experiment TCE-1. Note that MF is the mass fraction ((mass at time t)/(initial mass at t = 0)) and FF is the flux fraction ((flux at time t) /(initial flux at t = 0)). Figures 3-2 to 3-4 indicate that reductions in flux were realized for partial reductions in mass as certain well segments showed flux decreases at earlier dissolution times than other well segments. The well segments that continued to produce flux at later dissolution times were predominantly associated with streamtube bundles characterized by larger Svalues (NAPL morphologies that were elongated in the mean flow direction). The well segments that were depleted of mass at earlier dissolution times were predominantly associated with streamtube bundles characterized by smaller values (NAPL morphologies that were less elongated in the mean flow direction). Figures 3-2 to 3-4 indicate that variability inyields reductions in flux for partial reductions in mass as the smaller S streamtubes are depleted of mass at earlier dissolution times. S S

PAGE 89

75 Figure 3-4. Flux plane response to changes in NAPL architecture for experiment DCA-2. Note that MF is the mass fraction ((mass at time t)/(initial mass at t = 0)) and FF is the flux fraction ((flux at time t) /(initial flux at t = 0)). Variability in the distribution associated with each well segment resulted in ports where the reduction in the mass fraction (expressed as the mass removed from the individual well segment divided by the total initial mass in the model source zone) was larger than the reduction in the flux fraction (expressed as the flux discharging from the individual well segment divided by the total initial flux discharging from the entire model source zone), ports where the reduction in the mass fraction and the flux fraction were approximately equal, and ports where the reduction in the mass fraction was less than the reduction in the flux fraction. These effects are illustrated in figure 3-5, which displays cumulative reductions in the mass fraction and the flux fraction for selected ports. For all S

PAGE 90

76 three experiments, port 1 was selected as the well segment in which the reduction in flux was smaller than the reduction in mass (figure 3-5a, 3-5d, 3-5g). The streamtube bundles associated with these well segments were characterized by larger trajectory integrated NAPL contents resulting from lateral migration of the NAPL along the base of the flow cell. The higher streamtubes required the depletion of more up-gradient NAPL mass to achieve reductions in flux when compared to the smaller streamtubes. In accordance, when considered in relation to the other well segments, the reduction in the flux fraction from port 1 is smaller than the reduction in the mass fraction (figure 3-5a, 3-5d, 3-5g) where the opposite is true (the reduction in the flux fraction is larger than the reduction in the mass fraction) for ports associated with streamtube bundles characterized by smaller values which require the removal of less mass to achieve reductions in flux (figure 3-5c, 3-5f, 3-5i). The ports where the reductions in mass and flux were approximately equal (figure 3b, 3e, 3h) were not necessarily associated with a medium length NAPL morphology or a bundle of medium S streamtubes, but rather more of a triangular distribution overlaying an elongated portion on top of a low permeability layer. Put in terms of the streamtube bundle analogy, the distribution of initial in the streamtubes was such that the cumulative reduction in flux and mass were approximately equal. S S S S This variability among different portions of the source zone is further illustrated in figure 3-6 by plotting the nine ports displayed in figure 3-5 in a mass reduction/flux reduction framework along with the integrated mass reduction/flux reduction relationship for the entire system. For all three experiments, the integrated mass reduction/flux reduction behavior for the entire system was approximately linear, with a 1:1 relationship between reductions in mass and reductions in flux. This 1:1 behavior however was the

PAGE 91

77 Figure 3-5. Cumulative reductions in mass fraction and flux fraction with increasing dissolution time for selected ports from experiments DCA-1, TCE-1, and DCA-2. Note that the cumulative reductions are expressed in relation to the total mass reduction and the total flux reduction (e.g. 27% of the total reduction in flux for experiment DCA-2 comes from port 5). result of a combination of segmented mass reduction/flux reduction relationships that differed significantly from the 1:1 behavior of the entire system, with segments that were both less than (elongated or highstreamtubes) and greater than (shorter or small streamtubes) the 1:1 mass reduction/flux reduction relationship of the entire system. S S Source Depletion Models The estimated streamtube mass fractions (see image analysis section) for experiments DCA-1, TCE-1, and DCA-2 are displayed in figure 3-7. In order to simulate

PAGE 92

78 dissolution behavior, the mass fractions displayed in figure 3-7 were converted to values of and inserted into equation 3-1 and 3-2. The data was modeled under the assumption S Figure 3-6. Mass reduction/flux reduction relationships for the entire system (upper three graphs) and for selected well segments (bottom three graphs). that the dissolution dynamics were primarily controlled by the distribution. A nonreactive travel time, t, equal to the mean nonreactive travel time for the entire system was assumed for each streamtube (measured mean travel times for the three experiments were 1.2 hr for DCA-1, 1.7 hr for TCE-1, and 1.3 hr for DCA-2). The superposition of n solutions, with n being the number of trajectories from the image analysis procedure, was then used to develop the plots in figure 3-8. S The nonequilibrium model required fitting one parameter and the Da model required fitting 2 parameters. Two approaches were available for estimating f c : if equilibrium conditions were assumed, f c could be estimated from the initial flux-averaged concentration, otherwise f c could be estimated using figure 3-7. For the equilibrium streamtube model, estimating f c from figure 3-7, based on the number of trajectories

PAGE 93

79 Figure 3-7. Streamtube mass fractions estimated from the image analysis technique. containing NAPL, resulted in an overestimation of C o Explanations for the overestimation of C o when using an f c obtained from the image analysis include velocities high enough to induce rate-limiting effects in the smaller trajectories and incomplete span of the NAPL across the entire thickness of the box in certain locations, which would S

PAGE 94

80 result in bypass flow that was not captured in the two-dimensional image analysis technique. After estimating f c from C o and using NAPL distribution information obtained from the image analysis, the equilibrium streamtube models matched the observed dissolution behavior closely (figure 3-8), which supports the argument that the NAPL architecture was the primary factor controlling dissolution behavior. Additional processes not explicitly accounted for in the model, such as relative permeability, transverse dispersivity, and velocity variability are thought to exert a secondary influence on dissolution behavior. Figure 3-8. Comparison of the rate-limited streamtube, equilibrium streamtube, and effective Damkohler approaches for modeling source depletion from experiments DCA-1, TCE-1, and DCA-2. The rate-limited streamtube model required the estimation of the mass transfer rate coefficient, k, to obtain the fits shown in figure 3-8. With the incorporation of rate

PAGE 95

81 limiting effects into the model, the value for f c measured from the image analysis procedure could be utilized directly. As a comparison of the values of k used here, TCE dissolution data from a one-dimensional column study of uniform residual saturation presented in Imhoff et al. (1994) (experiment 6 in their paper) was modeled using equation 3-2 above. Using an estimated travel time from their experimental conditions, equation 3-2 provides an excellent fit of their data with a value of k equal to 65/hr, which is in the range of the values used here. It is noted that based on equation 3-2, increased accuracy from incorporating rate-limiting effects at the laboratory scale would likely be insignificant at the field-scale for typical groundwater velocities and source zone dimensions. For example, for a travel time of two days and a value of k = 42/hr, equation 3-2 does not predict rate-limiting effects until the trajectory average NAPL content,, is reduced to a value 0.002 upon which the concentration drops rapidly, closely approximating a step function. S The Da model also provided a good fit to the data with 2 = 1 and o scaled to reflect C o Again note that for the case of 2 = 1, the Da model is reduced to an exponentially-decaying function. The case of 2 = 1 is attractive because there is possibly a wide range of long-term dissolution behavior that can be reasonably approximated with an exponentially-decaying type model. The parameterization of the exponential model is also attractive as k can be expressed in terms of the initial mass, meaning that measurements of the mean groundwater velocity, the initial flux averaged concentration, and the initial mass would be sufficient to parameterize the model. For cases where 12 there is little guidance as to how to go about estimating 2 a priori. However, based on the effectiveness of the streamtube models, the parameter 2 could be

PAGE 96

82 estimated empirically using reactive travel time statistics in a manner similar to the techniques discussed in the following chapter. ADE model fits for DCA-1, DCA-2, and TCE-1 are shown below in figure 3-9. The parameter optimized in the curve fitting procedure was the Peclet number, which again note controls the shape of the BTC in a similar fashion to the second moment of the contribution time distribution. Like the Da model, the effectiveness of the streamtube models suggests that the Peclet number can perhaps be estimated using tracer techniques which are discussed in the following chapter. Figure 3-9. ADE model fit for source depletion experiments DCA-1, TCE-1, and DCA-2.

PAGE 97

83 Additional NAPL Architectures As a demonstration of NAPL architectures that would differ significantly from the 1:1 mass reduction/flux reduction behavior of experiments DCA-1, TCE-1, and DCA-2, a fourth experiment, DCA-3 is introduced. The packing procedure in this experiment was designed to create small alternating sand layers spanning the entire width of the box. During injection, the NAPL immediately migrated laterally along two 20/30 (the highest permeability sand used) layers in the vicinity of the injection point resulting in the highly elongated distribution shown in figure 3-10. The dissolution process resulted in the continual shortening of these two elongated portions with effluent concentrations remaining constant until a large fraction of the mass had been removed (figure 3-10). The mass reduction/flux reduction relationship for this system is displayed in figure 3-11. Note the similarities between the NAPL distribution of DCA-3 and the NAPL distribution up-gradient of segmented port 1 in the previous three experiments. As expected, the dissolution behavior and resulting mass reduction/flux reduction relationship of DCA-3 (figure 3-10; 3-11) was similar to port 1 in DCA-1, DCA-2, and TCE-1 (figure 3-6). Figure 3-10. Initial NAPL distribution and effluent BTC for experiment DCA-3.

PAGE 98

84 In contrast to DCA-3, it is possible to conceptualize a system where the relationship between mass and flux is quite favorable. Consider a source zone similar to the one simulated by Parker and Park (2004) where residual zones of various length overlie a large, laterally extensive pool. In this type of system a large portion of the initial flux is associated with the overlying residual zones while a large portion of the initial mass is associated with the laterally extensive pool at the base of the simulation domain. As discussed previously, the dissolution profile for such a system is similar to the high-variance lognormal distribution of case in the streamtube formulation where there is a large initial decrease in flux as the overlying residual zones are depleted of mass followed by extensive tailing resulting from the laterally extensive pool. The resulting mass reduction/flux reduction for such a system is quite favorable where there is a large decrease in flux for a relatively small decrease in mass followed by limited reduction in flux as the laterally extensive pool is slowly depleted of mass. Experimentally, this case is difficult to construct as a large laterally extensive pool invariably migrates into the injection and extraction wells. As a substitute, consider a hypothetical experiment where 75% of the well segments are characterized by BTCs identical to the BTC of port 5 in DCA-2 and 25% of the well segments are characterized by a BTC identical to that of port 1 in experiment DCA-2 (see figure 3-5). The resulting BTC is similar to the high-variance lognormal streamtube case with a large initial drop in flux followed by extensive tailing. The mass reduction/flux reduction relationship from the superposition of these data sets is displayed alongside that of DCA-3 in figure 3-11. S Experiment DCA-3 did not conform to the model assumptions as well as the other experiments because the flow paths were likely not parallel with the bottom of the flow

PAGE 99

85 cell but rather followed the curvilinear direction of the alternating sand layers. In order to estimate the streamtube mass fractions for modeling purposes, the algorithm used to process the binary image was abandoned in favor of manual counting of the NAPL containing grid cells along the curvilinear layers of alternating grain size. In contrast to the previous three experiments, f c did not need to be measured using C o for the equilibrium streamtube model, suggesting that rate-limiting effects were less important Figure 3-11. Mass reduction flux reduction relationship for experiment DCA-3 and hypothetical system developed by the superposition of selected segmented BTCs. due to the elongated distribution of the NAPL. The mass transfer rate coefficient, k, in the rate-limited streamtube model therefore acted not to scale C o to the appropriate initial value, but rather to provide a smoother fit of the data. As in the previous DCA experiments, a value of k = 42/hr was selected for the mass transfer rate coefficient. Both the equilibrium and rate-limited streamtube models provide reasonably good fits of the data, supporting the argument that the NAPL architecture is the primary mechanism controlling dissolution behavior in the systems studied. The data/modeling of DCA-3

PAGE 100

86 supports the general trend that increased variability in the distribution leads to more favorable relationships between reductions in mass and flux. In DCA-3, there is limited variability in the distribution due to most of the NAPL mass concentrated in a small number of streamtubes, which results in a less favorable relationship between mass and flux when compared to experiments DCA-1, DCA-2, and TCE-1, which exhibit more variability in terms of their distributions. S S S Figure 3-12. Comparison of the rate-limited streamtube, equilibrium streamtube, effective Damkohler, and ADE approaches for modeling effluent BTCs from experiments DCA-3. The Da model again required scaling o to the initial concentration. Unlike the previous experiments, 2 1, meaning that 2 had to be estimated by fitting the dissolution data, as opposed to solving for the rate constant of an exponential decay model in terms of the initial flux and the initial mass. A value of 2 < 1 for laterally extensive portions of the source zone was consistent with Parker and Park (2004) who fit their model to both the upper portion of their simulation domain, which consisted of

PAGE 101

87 residual zones of variable lengths, with a 2 = 1.1, and to the lower zone, which consisted of a laterally extensive pool, with a 2 = 0.4. A value of 0.1 for 2 was used in this work to obtain the fit displayed in figure 3-12. As shown in figure 3-12, the model is not capable of simulating both the initial stage of constant concentration and the late-term tailing behavior, which is a limitation of all power function-type models. Aside from the inability to simulate the late-term tailing behavior, the Da model provides a reasonable fit of the data and outlines a trend of more favorable relationships between mass and flux for larger values of 2 The general relationship between 2 and is that S 2 increases with increased variability in the distribution. S The ADE model, unlike the Da model is capable of producing the type of dissolution profile characteristic of DCA-3 and therefore provides an excellent fit of the observed data. As discussed in the previous chapter, the streamtube model (with a lognormal distribution) and the ADE model produce similarly shaped BTCs. The streamtube model therefore would also provide an excellent fit of the data by fitting the second moment of the tau distribution. Conclusions Four aqueous dissolution experiments in two-dimensional heterogeneous systems were presented with an emphasis on the dissolved contaminant flux response to reductions in DNAPL mass. The relationship between mass and flux was found to be primarily controlled by the NAPL architecture, with mechanisms such as relative permeability, velocity variability, transverse dispersion, and mass transfer rate limitations exerting only secondary influences. A Lagrangian-type definition of NAPL architecture was utilized which incorporates the concept of a trajectory integrated NAPL content,. S

PAGE 102

88 The distribution of S amongst streamtubes governed the relationship between mass and flux. For experiments DCA-1, DCA-2, and TCE-1 the relationship between reductions in mass and flux was found to be approximately 1:1. In addition to experiments DCA-1, DCA-2, and TCE-1, a fourth experiment, DCA-3, and a theoretical case, developed from the superposition of BTCs from selected well segments, were used to illustrate cases where mass reduction/flux reduction would deviate significantly from the 1:1 behavior of DCA-1, DCA-2, and TCE-1. The dissolution profiles from the four experiments were modeled using three simplified source depletion models; a streamtube model, an advection dispersion model and an effective Damkohler number model previously presented in Parker and Park (2004). The distribution of amongst the streamtubes was estimated using a simplified image analysis procedure. Based on the effectiveness of the streamtube models, dissolution dynamics for the systems studied were determined to be primarily controlled by the distribution with mechanisms such as relative permeability and velocity variability causing perturbation about the mean behavior. In the next chapter a combination tracer techniques are evaluated for there ability to parameterize these models for prediction of long-term source depletion. S S

PAGE 103

CHAPTER 4 AN EVALUATION OF PARTITIONING TRACERS FOR PARAMETERIZING SOURCE DEPLETION MODELS In this chapter partitioning tracers are investigated as a potential parameterization technique for the models introduced in chapter 2. Two common parameter estimation techniques are evaluated: the method of moments where measured moments from an extraction well BTC are equated with derived moment equations, and inverse modeling where BTCs are fit to a transport model using specified objective functions. In addition, the ADE source depletion model of chapter 2 is parameterized using a hybrid approach where the total mass is determined using moment analysis and the Peclet number is determined using traditional curve fitting techniques. For brevity these approaches will be referred to as the moment, inverse lognormal, and ADE approach respectively. Theory Overview of Partitioning Tracers Partitioning tracers are a characterization tool for estimating NAPL saturation between pairs of injection and extraction wells (Annable et al., 1998; Brooks et al., 2002; Cain et al., 2000; Falta et al., 1999; Jawitz et al., 1998b; Jawitz et al., 2000; Meinardus et al., 2002; Nelson and Brusseau, 1996). Reactive tracers that partition into the NAPL phase are injected with one or more nonreactive tracers. The retardation of the reactive/partitioning tracer in relation to the nonreactive tracer is then used to estimate the integral NAPL saturation in the swept volume between the injection well and the extraction well. An integrated measurement of NAPL saturation is conceptually 89

PAGE 104

90 appealing as it circumvents the inherent difficulty of interpolating between highly spatially variable point measurements of NAPL saturation, which is required when employing conventional point measurement tools such as soil coring. Although partitioning tracers are conceptually appealing because they provide an integrated measurement of an often highly heterogeneous saturation field, there are some questions as to the reliability of the measurements obtained from partitioning tracers. Rao et al. (2000) and Brooks et al. (2002) outline most of the technical issues associated with the use of the partitioning tracers. Some of the issues, such as partitioning nonlinearity, can be overcome through modifications to the tracer test itself (e.g., injection of multiple tracers) while other issues may be minimized through data analysis techniques (e.g., extrapolation of BTC tails). Other issues, such as hydrodynamic accessibility of the tracer to low permeability zones and interference from organic matter and mineral components are inherent limitations of the technique (Rao et al., 2000). Practical issues of concern associated with partitioning tracers are its high cost and possible regulatory issues (EPA, 2003). A comparison of the reliability of NAPL saturation estimates obtained from tracers with more conventional point measurement techniques should not simply consider the accuracy of the point measurement vs. the accuracy of the integrated measurement, but should also include the error associated with interpolating between point measurements. The centimeter scale sensitivity of a migrating NAPL to very slight changes in permeability creates highly irregular entrapment configurations (Keuper et al., 1993), which poses questions as to the reliability of conventional techniques, regardless of how accurate the individual point measurements are.

PAGE 105

91 Moment Based Parameter Estimation The moment equations developed below are based on equations presented in Jawitz et al. (2003a). The equations of Jawitz et al. (2003a) can be simplified considerably by assuming no correlation between the travel time and NAPL saturation and the contaminated fraction parameter, (f c ), is equal to the flux-averaged contaminant concentration measured in the extraction well during the tracer test. For example, if a partitioning tracer test was conducted in a PCE contaminated source zone and the measured flux-averaged PCE concentration in the extraction well during the partitioning tracer test was 60% of solubility, it is assumed that 60% of the streamtubes are contaminated. The assumption that the travel time and trajectory integrated NAPL content, are uncorrelated is somewhat supported by the recent work of Lemke et al. (2004b) who found that the extent of correlation between organic liquid saturation and permeability that has been postulated by others (e.g. Anderson et al., 1992; Berglund, 1997; and Rao et al., 2000) was not supported by their simulation results. It should be noted however that the analysis of Lemke et al. (2004b) was based on an Eulerian type comparison, plotting point values of NAPL saturations against the corresponding hydraulic conductivity for that particular point. While investigation of a correlation structure is useful for theoretical simulations there is a practical concern when dealing with mathematical formulations that include the product of two correlated lognormals. From inspection of equation 2-11 (also 4-18 below) it is clear that there are an infinite number of and combinations of the individual lognormals that lead to the same aggregate for the product of two lognormals. To illustrate this consider fitting the reactive tracer data shown below in figure 4-1. S

PAGE 106

92 Cumulative Pore VolumesFlux Averaged Concentration UTCHEM generated reactive tracer data Model fit using correlated lognormals Figure 4-1. UTCHEM generated reactive tracer fit with product of two lognormals. The 3-D plot in figure 4-2 below shows the root mean square for a matrix of variances and correlation coefficients. Evident in figure 4-2 is that for a function which is a product of two correlated lognormals, there are an infinite number of individual combinations of correlation and variance that produce the same aggregate variance for the product of two lognormals. Therefore, from a practical standpoint, assuming no correlation does not limit in any way the range of BTCs the model can produce and dramatically simplifies the equations. An exception to this statement is that if the nonreactive tracer exhibits more spreading than the partitioning tracer, a negative correlation would be needed in order to obtain a reasonable fit of the partitioning tracer (this will be discussed further below).

PAGE 107

93 correlationvariance RRoot Mean Square combinations of correlation and variance yielding same solution Figure 4-2. RMS matrix for tracer data in figure 4-1. Measured moments The absolute moments of a measured BTC can be determined by substituting the measured concentration at the extraction well, C(T), for the lognormal probability distribution in equation 2-6. (4-1) 0)(dTTCTMNN The N th normalized moment (m N ) is determined by dividing the N th absolute moment by the 0 th absolute moment. 0MMmNN (4-2)

PAGE 108

94 For step injections the first and second normalized moments can be expressed with equation 4-3 and 4-4 (Sardin et al., 1991). An advantage of the moment equations for step injections and step decreases is that the second moment lacks the higher order term (T 2 ) and leads to decreased truncation error (Young and Ball, 2000). (4-3) dTTCm01)(1 (4-4) dTTCtm02)(12 First and second moments of nonreactive tracer The mean and variance of the travel time distribution can be determined from the first and second measured moments of the nonreactive/nonpartitioning tracer (denoted by the superscript np): 211onptTmm (4-5) 42122onponptTmTmm (4-6) where 4-5 and 4-6 are corrected for the pulse duration T o The first and second moments of the travel time distribution can then be inserted into equations 2-8 and 2-9 to determine the mean and variance of the travel time distribution. Domain average trajectory integrated NAPL content For an individual streamtube, the arrival time of a partitioning tracer is given by: Rttnpp (4-7) where the subscript p denotes partitioning tracer and the subscript np denotes nonpartitioning tracer. After Jin (1995), R is given by the following equation.

PAGE 109

95 NNnSSKR11 (4-8) with: wnnCCK (4-9) Where C n is the concentration of the tracer in the NAPL phase and C w is the concentration of the tracer in water. Note that equation 4-8 differs from Jin (1995) in that the trajectory integrated NAPL saturation along a streamtube,, has been substituted for the domain average saturation S NS N Equation 4-9 can be expressed in terms of the trajectory integrated NAPL content,, in order to be consistent with the streamtube models presented in chapter 2 by making the following substitution. S SKRn1 with: wNS1 (4-10) Assuming no correlation between R and the nonreactive tracer the expected value of 4-7 can be determined by: 22111onpopRTmTmm (4-11) The expected value ofcan then be expressed through the measured moments as: S nRSKmm111 (4-12)

PAGE 110

96 Excess spreading The traditional evaluation of partitioning tracers has been limited to analysis of the first moment (domain average saturation). Information related to the spatial distribution of the NAPL has historically been limited to soil coring (Jawitz et al., 2000; Meinardus et al., 2002; Rao et al., 1997) or partitioning tracer tests conducted at multiple sampling locations (Jawitz et al., 2000; Sillan et al., 1998). Jawitz et al. (2003a) presented equations utilizing partitioning tracer higher moments to obtain information regarding NAPL spatial distribution between pairs of injection and extraction wells. The higher moment technique is based on the concept of excess spreading. This concept is perhaps best illustrated through use of the advection dispersion equation. If we were to fit a nonreactive tracer BTC with the analytical solution to the ADE (equation 2-21injection and detection in flux) by adjusting the Peclet number, a partitioning tracer response for a homogenously distributed NAPL should then theoretically be described by the same Peclet number determined from the nonreactive tracer test and an R value equal to equation 2-5 (e.g., Valochi, 1985). Any additional spreading beyond that predicted with the Peclet number from the nonreactive tracer test is considered excess spreading potentially resulting from heterogeneity in the NAPL distribution (or any other reactive parameter). Theoretically, this excess spreading should manifest itself in the higher moments of the partitioning tracer, allowing for quantification of variability in the Lagrangian-based NAPL architecture. The excess spreading concept is illustrated below in figure 4-3.

PAGE 111

97 Pore VolumesConcentrationtheoretical value assuming uniform NAPL distribution excess spreading caused by heterogeneity in the NAPL distribution nonreactive tracer Figure 4-3. BTCs representing transport of a nonreactive tracer, a reactive tracer with homogenously distributed NAPL, and a reactive tracer with heterogeneously distributed NAPL. Moment equations for quantifying Lagrangian NAPL architecture The second moment of the retardation factor R can be expressed as: 442122122oonpnpooppRTTmmTTmmm (4-13) The second moment of the trajectory integrated NAPL content can then be expressed as: 212221nSnRSKmKmm (4-14) The first and second moments of the trajectory integrated NAPL content can then be inserted into equations 2-8 and 2-9 to determine the mean and variance of the distribution. As outlined previously in chapter 2, the mean and variance of the S S

PAGE 112

98 distribution parameters can then be combined with the mean and variance of the travel time distribution parameters to determine the mean and variance of the contribution time () distribution, which in turn is used to predict source depletion under natural gradient conditions. A summary of the steps used to determine the mean and variance of the t and distributions is shown below in figure 4-4. S Nonreactive Tracer Data Determine first and second moments of twith equations 5-5 and 5-6 Reactive Tracer Data Determine first and second moments of Rwith equations 5-11 and 5-13 Determine first and second moments of the trajectory integrated NAPL content with equations 5-12 and 5-14 Determine mean and variance of contribution time using equation 3-13 and equations 3 -15 through 3-18. Determine mean and variance of the trajectory integrated NAPL content with equations 3-8 and 3-9 Determine mean and variance of twith equations 3-8 and 3-9 Forward predict dissolution profile with equation 3-14 and 3-1 Figure 4-4. Summary of process used to predict dissolution using partitioning tracers and moment analysis. Moments for the partitioning and nonpartitioning BTCs were fit to a cubic spline function in Mathcad (figure 4-5). Cubic spline interpolation passes a curve through a set of data points such that the first and second derivatives are continuous across each point. There are two primary advantages of replacing the discrete tracer data with a continuous cubic spline function. The first advantage is that a continuous function allows for the use

PAGE 113

99 of Mathcads more accurate built in numerical integration schemes such as adaptive and Romberg integration as opposed to using a trapezoidal rule calculation in a spreadsheet type application. The second advantage is the spline function extrapolates the tracer data beyond the last sampling point in cases of incomplete data collection of the BTC tail due to substantial tailing. (Note that the level of extrapolation accomplished with the spline function may be insufficient for typical field applications where the truncation error is significantly larger than the idealized numerical simulations presented here.) Data Cubic splineinterpolation Figure 4-5. Cubic spline interpolation of arbitrary UTCHEM generated tracer data. Inverse Modeling Based Parameterization (Lognormal Distribution) An alternative parameter estimation technique to the method of moments is inverse modeling where the extraction well BTC is fit to a process model. Considering the arrival time of the nonreactive tracer as a lognormally distributed random variable allows for expression of the extraction well BTC (C f (T)) using a lognormal probability distribution

PAGE 114

100 function (pdf) for Dirac inputs (equation 4-15) and a cumulative distribution function (cdf) for step inputs (equation 4-16). )()(TfTCf (4-15) )()(TpTCf (4-16) Where f and p are used to denote a pdf and a cdf respectively. For pulse inputs (or any other type of input), 4-15 can be utilized in conjunction with a convolution integral to obtain the mean and variance of t. (4-17) TdttftTCTC0)()()( The travel of a partitioning tracer (denoted by the subscript p in equation 4-18 and 4-19) is equal to the product of t and R, meaning that equations 4-15 through 4-17 can still be utilized by making use of equation 2-10 and 2-11 to define the mean and variance of the partitioning tracer arrival time (Jury and Roth; 1990). tRp (4-18) 22ttRRp (4-19) Again note that the assumption is made that R is uncorrelated with travel time. The correlation coefficient was preserved in equation 4-19 to again emphasize that there are and infinite number of combinations of variance and correlation of individual lognormals that yield identical variances for the product of two lognormals. The exception to the previous statement in noted again as well. From inspection of 4-19, if the product of a negative correlation coefficient and the standard deviations of R and t is larger than the variance of R, the reactive tracer would exhibit less spreading than the nonreactive tracer, meaning that the best fit could not be obtained without use of the correlation coefficient.

PAGE 115

101 With R determined, the mean and variance of for uncorrelated cases can be determined using the moment generating function 2-8 and equation 4-11 and 4-14. S The inverse modeling process utilized here was to first fit the nonreactive tracer using equation 4-16 in order to determine t and t 2 Values for R and R 2 were then determined by fitting the partitioning tracer with equation 4-16 using values of t and t 2 determined from the nonreactive tracer and a and equal to equation 4-18 and 4-19. The values for R and R 2 were then inserted into the moment generating function 2-7 to determine the first and second moments of R. The first and second moments of R were then used in conjunction with equations 4-12 and 4-13 to obtain the first and second moments of. The first and second moment of, together with the first and second moments of the travel time distribution were then used to determine the first and second moments of the contribution time distribution using the equations outlined in chapter 2. S S Additional Contribution Time Parameterization Methods In simulation sets one and two the above techniques tended to over predict spreading of the dissolution profile. Because of this, two additional techniques were investigated that deemphasized velocity variability (this will be discussed in further detail in the results section). ADE based paramterization The ADE approach involved using traditional moment analysis to determine the pore volume (such that R for the nonreactive tracer was 1) and the retardation factor for the reactive tracer. The retardation factor was then used to calculate the domain average NAPL saturation, which was used in equation 2-23 in conjunction with f c to determine R for the ADE source depletion model. After retardation factors were calculated, traditional

PAGE 116

102 curve fitting techniques were used to determine Peclet numbers for the nonreactive (P nr ) and reactive (P r ) tracers. The simple estimator displayed below in equation 4-20 was used to estimate the Peclet number for the source depletion model (P SD ) based on the nonreactive and reactive tracers. Equation 4-20 is semi-empirical measure of excess spreading as defined by the nonreactive and reactive tracers: rnrSDrrnrSDPPPRPPP for P nr P r > 2 (4-20) for P nr P r < 2 where R r is the retardation factor of the reactive/partitioning tracer. Inverse modeling with no travel time variability For simulation sets 1 and 2 the inverse modeling approach over predicted dissolution profile spreading. Because of this, an approach identical to the inverse lognormal procedure outlined above, except without the second moment of travel time (travel time variability) incorporated into the second moment of the contribution time, was investigated. Using this approach, equation 2-18 can be expressed as: StfmmKm22122 (4-21) UTCHEM Simulations Prior to the NAPL dissolution simulations detailed in chapter 2, a partitioning tracer test was simulated to evaluate the ability of the techniques presented above to parameterize the source depletion models of chapter 2.

PAGE 117

103 Simulation Setup A nonpartitioning and a partitioning tracer were introduced into the domain via the injection well under steady-state flow at the same flow rate as the subsequent water flood/aqueous dissolution simulation at a concentration of 500 mg/L. A partitioning coefficient of 50 was selected after preliminary simulations indicated negligible differences in predicted dissolution behavior when utilizing alternative partitioning coefficients. Although UTCHEM has the capability to simulate tracer sorption and degradation, these effects were not evaluated here. Results Predicted domain averaged NAPL saturation Predicted domain averaged NAPL saturations using the moment, inverse lognormal, and ADE approach are displayed below in table 4-1. Recall that the ADE approach uses the traditional moment technique (Jin, 1995) to estimate the domain averaged saturation (S N ). The inverse lognormal and moment approaches outlined above yield an estimate of the trajectory integrated NAPL content (), which was converted to S S N using equation 4-10. The moment and ADE approaches yielded significantly more accurate results than the inverse lognormal approach. The moment and ADE approaches exactly predicted the NAPL saturations for each simulation in both sets and 1 and 2. The inverse lognormal approach over predicted the NAPL saturations by an average of approximately 5% for both simulation sets 1 and 2. For simulations sets 3 and 4 the moment and ADE approaches under predicted the NAPL saturations by an average of approximately 1% for simulation set 3 and 3% for simulation set 4. The under prediction of NAPL saturations as heterogeneity increases is likely due to a combination of two factors. The first factor

PAGE 118

104 that could lead to an underestimation of NAPL saturation is a decrease in hydrodynamic accessibility due to relative permeability effects as the architecture becomes more skewed towards pools. The second factor is data truncat ion due to enhanced tailing of the tracer BTCs as heterogeneity increases. This effect is likely small however as the simulations were run well into the tail portion of the BTCs, thereby minimizing truncation error. Despite the enhanced tailing of tracer BTCs and possible issues associated with hydrodynamic accessibility, the moment and ADE predictions for simulation set 3 and 4 are extremely accurate. These findings are of course influenced by the fact that the simulator produces tracer data th at is not subject to analytic al/detection limit issues and exactly conforms to the underlying assump tions (e.g., perfect linear partitioning). The inverse lognormal approach performs poorly compared to the moment based approaches, with errors of around 7% for simulations sets 3 and 4. As will be shown in subsequent sections, despite the inaccuracy of the saturation estimates, the parameters determined from the inverse lognormal approach yield excellent fits of the tracer BTCs. This suggests that the RMS object ive function yields parameter se ts that accurately fit the spreading of the tracer data yet do no t preserve the mean arrival time. Method of moments dissolution predictions Moment derived parameters describing tr acer transport for all simulations are displayed below in table 4-2. The moment determined parameters for the mean and variance of the t and R distributions were inserted into equation 4-18 and 4-19 and plotted alongside the UTCHEM simulated tracer data (figure 4-6 through 4-9) The RMS values for the moment predicted tracer BT Cs are displayed in table 4-2.

PAGE 119

105 Table 4-1. Domain averaged NAPL saturations estimated with the different techniques. Simulation Actual ADE Inverse Lognormal Moments S S Error % S Error % S Error % Set 1 1 0.0043 0.0043 0.0 0.0045 0.0045 +4.4 0.0043 0.0043 0.0 2 0.0033 0.0033 0.0 0.0036 0.0036 +8.6 0.0033 0.0033 0.0 3 0.0043 0.0043 0.0 0.0046 0.0046 +6.3 0.0043 0.0043 0.0 4 0.0030 0.0030 0.0 0.0030 0.0030 0.0 0.0030 0.0030 0.0 0.0 5.1 0.0 average Set 2 1 0.0042 0.0042 0.0 0.0043 0.0043 +2.8 0.0042 0.0042 0.0 2 0.0033 0.0033 0.0 0.0035 0.0035 +6.4 0.0033 0.0033 0.0 3 0.0033 0.0033 0.0 0.0035 0.0035 +6.4 0.0033 0.0033 0.0 4 0.0045 0.0045 0.0 0.0047 0.0047 +4.9 0.0045 0.0045 0.0 0.0 5.2 0.0 average Set 3 1 0.0044 0.0043 -2.3 0.0041 0.0041 -6.6 0.0043 0.0044 0.0 2 0.0045 0.0045 0.0 0.0047 0.0047 +4.9 0.0045 0.0045 0.0 3 0.0045 0.0044 -2.2 0.0042 0.0042 -6.3 0.0045 0.0045 0.0 4 0.0045 0.0045 0.0 0.0034 0.0034 -24.2 0.0045 0.0045 0.0 1.1 10.5 0.0 average Set 4 1 0.0065 0.0063 -3.1 0.0056 0.0056 -13.4 0.0063 0.0064 -2.2 2 0.0067 0.0067 0.0 0.0071 0.0072 +6.7 0.0068 0.0069 +1.7 3 0.0065 0.0061 -6.2 0.0064 0.0064 -1.5 0.0061 0.0061 -5.6 4 0.0062 0.0060 -3.2 0.0057 0.0057 -7.4 0.0060 0.0061 -2.3 3.1 7.2 2.9 average S S

PAGE 120

106 SS One of the problems of assuming no correlation in the moment equations is that for high travel time variances, the calculated second moment of R is either very small or negative which leads to a calculated negative second moment of (equation 4-14). This of course has no physical meaning. A negative second moment of was calculated in simulations 2 and 3 of set 4 (table 4-2). For the more homogenous simulation sets 1 and 2, the moment derived parameters significantly under predicted tracer spreading (figure 4-6 and 4-7) when inserted into the lognormal cdf. The method of moments performed better, in terms of yielding parameters that more accurately reproduced tracer spreading, as the heterogeneity of the porous media increased (figure 4-8 and 4-9). The moment derived parameters describing tracer transport were used to determine the contribution time mean and variance using the process outlined in figure 4-4. The contribution time mean and variance for all simulations are displayed below in table 4-3. The predicted dissolution profiles using the parameters in table 4-3 are displayed below in figures 4-10 through 4-13. As shown in figure 4-10 and 4-11, the moment equations provide a fairly good prediction of dissolution dynamics for simulation sets 1 and 2, predicting slightly less spreading/variability of the dissolution profile than that of UTCHEM. An under-prediction of dissolution profile spreading was expected based on the performance of the moment equations in relation to fitting the partitioning tracer BTCs, yet the under prediction of spreading for the dissolution profiles is less than that of the tracer (also see RMS values in tables 4-2 and 4-3). This suggests that there are mechanisms that cause

PAGE 121

107 Table 4-2. Tracer derived parameters using method of moments. nonreactive reactive Simulation t t 2 RMS R R 2 RMS S S 2 1-1 -0.03 0.06 0.144 0.47 0.06 0.232 -5.63 0.35 1-2 -0.03 0.06 0.236 0.21 0.01 0.380 -5.86 0.29 1-3 -0.04 0.07 0.101 0.19 0.01 0.147 -5.63 0.36 1-4 -0.06 0.12 0.039 0.19 0.02 0.072 -6.10 0.53 average -0.04 0.08 0.130 0.26 0.03 0.208 -5.81 0.38 S.D. 0.02 0.03 0.083 0.14 0.02 0.132 0.22 0.10 2-1 -0.01 0.03 0.128 0.19 0.01 0.156 -5.59 0.24 2-2 -0.01 0.02 0.131 0.15 0.01 0.115 -5.91 0.39 2-3 -0.01 0.03 0.136 1.50 0.01 0.144 -5.94 0.46 2-4 -0.01 0.02 0.132 0.20 0.01 0.147 -5.57 0.32 average -0.01 0.02 0.132 0.51 0.01 0.141 -5.75 0.35 S.D. 0.00 0.00 0.003 0.66 0.00 0.018 0.20 0.09 3-1 -0.15 0.29 0.043 0.15 0.10 0.037 -6.17 1.47 3-2 -0.16 0.32 0.032 0.19 0.02 0.063 -5.60 0.38 3-3 -0.21 0.42 0.015 0.16 0.08 0.016 -6.13 1.27 3-4 -0.23 0.46 0.007 0.13 0.14 0.023 -6.25 1.69 average -0.19 0.37 0.024 0.16 0.08 0.035 -6.04 1.20 S.D. 0.04 0.08 0.016 0.03 0.05 0.021 0.30 0.57 4-1 -0.34 0.67 0.022 0.32 0.14 0.036 -5.51 0.89 4-2 -0.42 0.86 0.047 0.33 0.01 0.022 x x 4-3 -0.43 0.86 0.021 0.27 0.05 0.014 x x 4-4 -0.37 0.74 0.015 0.23 0.07 0.008 -5.55 0.87 average -0.39 0.78 0.026 0.29 0.07 0.020 -5.53 0.88 S.D. 0.04 0.09 0.014 0.05 0.05 0.012 0.03 0.02 x = negative second moment of S

PAGE 122

108 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-6. Set 1 tracer fit (method of moments). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-7. Set 2 tracer fit (method of moments).

PAGE 123

109 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-8. Set 3 tracer fits (method of moments). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-9. Set 4 tracer fits (method of moments).

PAGE 124

110 Table 4-3. Moment derived contribution time mean and variance. SimulationRMS1-11.740.410.1231-21.50.340.1031-31.730.430.0761-41.240.660.059average1.550.460.090S.D.0.240.140.0282-11.390.260.1812-21.070.410.0362-31.050.490.0552-41.420.340.079average1.230.380.088S.D.0.200.100.0653-11.211.760.0753-21.310.70.0633-31.491.690.0633-41.122.140.174average1.281.570.094S.D.0.160.610.0544-11.651.570.1054-2xxx4-3xxx4-42.241.610.069average1.951.590.09S.D.0.420.030.025 2 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-10. Set 1 predicted source depletion (method of moments

PAGE 125

111 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-11. Set 2 predicted source depletion (method of moments). spreading of the tracers in the simulator which do not, in turn cause spreading of the dissolution profile. This is most likely due to the dispersivity parameter. In a finite difference numerical model of solute transport there is always some level of spatial averaging that must occur in order to define an average parameter value for each grid cell. The dispersivity parameter defines spatially averaged spreading that occurs as a solute move through an individual grid cell. According to the assumptions of the streamtube model, spreading of a nonreactive solute is due to differences in streamtube velocities, which in turn leads to spreading of the dissolution profile according to the equations outlined in chapter 2. In UTCHEM however, dispersivity at the scale of the individual grid cell does not lead to spreading of the dissolution profile at the grid cell scale. An alternative explanation is that dispersive mixing between streamtubes dampens spreading of the dissolution profiles This was evaluated recently by Cirpka and Kitanidis

PAGE 126

112 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-12. Set 3 predicted source depletion (method of moments). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulationnegative second moment of Snegative second moment of S Figure 4-13. Set 4 predicted source depletion (method of moments).

PAGE 127

113 (2000) in relation to bio-reactive transport. Their approach was to incorporate longitudinal mixing in each streamtube as a surrogate to transverse mixing processes. Incorporation of longitudinal dispersion into the streamtubes dampens the spreading of the BTC when compared to the advection only case. For simulations sets 3 and 4 the method of moments now over predicts spreading of the tracer BTC. This is in agreement with the dispersivity hypothesis discussed above. If the spreading evidenced in the nonreactive tracer is not entirely due to velocity variability, a more accurate fit of nonreactive tracer spreading (see figures 4-6 to 4-9) suggests an over prediction of dissolution profile spreading. As noted previously, one of the problems of assuming no correlation in the moment equations is that for cases where the reactive tracer exhibits little to no excess spreading, the calculated second moment of R is either very small or negative which leads to a calculated negative second moment of (equation 4-14). This of course has no physical meaning. A negative second moment of was calculated in simulations 2 and 3, precluding dissolution profile prediction. S S Inverse lognormal dissolution predictions The mean and variance of the travel time, retardation coefficient, and distribution determined from the inverse lognormal approach for simulation set 1 are displayed below in table 4-4. The fit tracer data is shown below in figures 4-14 through 4-17. From the RMS values in table 4-4 and figures 4-14 through 4-17, the inverse lognormal approach provides a substantially better fit of the tracer data than the method of moments. S In the more heterogeneous systems (simulation sets 3 and 4), it became increasingly difficult for the lognormal distribution to capture the later term tailing

PAGE 128

114 S For simulations sets 1 and 2, despite a highly accurate fit of the tracer data, the model predicts more spreading of the dissolution profile predicted by than the UTCHEM output. This again suggests either there is dampening of dissolution profile spreading by dispersive mixing processes not accounted for in the advection only streamtube models, or, alternatively, that dispersive mixing leads to enhanced spreading of the tracers which is erroneously interpreted as resulting entirely from velocity variability. For simulations The predicted source depletion parameters for all simulation sets along with the RMS values are displayed below in table 4-5. As discussed previously in relation to the method of moments, an accurate fit of the partitioning tracer data leads to an over prediction of dissolution profile spreading. Because of this, the second moment of the contribution time (), was calculated without the incorporation of travel time variability using equation 4-21. The source depletion parameters for no travel time variability and associated RMS values are also displayed in table 4-5. The predicted dissolution profiles are shown in figures 4-17 through 4-24. portion of the tracer BTCs. This difficulty was most evident in simulations 3-4, 4-1, and 4-2 for the partitioning tracer. Enhanced tailing of the partitioning tracer is a result of both increased velocity variability and increased variability in the trajectory integrated NAPL content. Recall that as the heterogeneity of the permeability field increases the NAPL architecture becomes more skewed to the existence of NAPL in elongated pools, enhancing tailing of the partitioning tracer. The poor fit of the tail portion of the partitioning tracer in simulation set 3-4 led to a case where there was no excess spreading of the partitioning tracer in relation to the nonpartitioning tracer, such that the variance of was zero.

PAGE 129

115 Table 4-3. Tracer derived parameters using inverse lognormal. nonreactive reactive Simulation t t 2 RMS R R 2 RMS S S 2 1-1 -0.04 0.20 0.001 0.46 0.10 0.002 -5.660 0.50 1-2 -0.04 0.19 0.002 0.21 0.05 0.003 -6.030 0.75 1-3 -0.04 0.20 0.001 0.19 0.04 0.007 -5.770 0.77 1-4 -0.06 0.21 0.003 0.18 0.05 0.009 -6.250 0.91 average -0.04 0.20 0.002 0.26 0.06 0.005 -5.928 0.73 S.D. 0.01 0.01 0.001 0.13 0.03 0.003 0.265 0.17 2-1 -0.02 0.15 0.000 0.19 0.03 0.000 -5.720 0.58 2-2 -0.02 0.13 0.000 0.15 0.03 0.000 -6.070 0.85 2-3 -0.01 0.17 0.000 0.15 0.03 0.000 -6.100 0.89 2-4 -0.02 0.15 0.000 0.20 0.03 0.000 -5.660 0.60 average -0.02 0.15 0.000 0.17 0.03 0.000 -5.888 0.73 S.D. 0.00 0.02 0.000 0.03 0.00 0.000 0.230 0.16 3-1 -0.15 0.45 0.004 0.14 0.09 0.007 -6.200 1.40 3-2 -0.15 0.45 0.004 0.18 0.06 0.004 -5.850 1.02 3-3 -0.18 0.51 0.001 0.15 0.09 0.002 -6.170 1.39 3-4 -0.19 0.50 0.003 0.16 0.00 0.008 -5.700 0.00 average -0.17 0.48 0.003 0.16 0.06 0.005 -5.980 0.95 S.D. 0.02 0.03 0.001 0.02 0.04 0.003 0.245 0.66 4-1 -0.27 0.53 0.009 0.30 0.09 0.009 -5.570 0.74 4-2 -0.31 0.56 0.002 0.27 0.07 0.002 -5.330 0.75 4-3 -0.39 0.66 0.005 0.25 0.06 0.005 -5.410 0.70 4-4 -0.31 0.60 0.004 0.18 0.14 0.004 -5.880 1.42 average -0.32 0.59 0.005 0.25 0.09 0.005 -5.548 0.90 S.D. 0.05 0.06 0.003 0.05 0.04 0.003 0.243 0.35

PAGE 130

116 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-14. Set 1 tracer fit (inverse lognormal). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-15. Set 2 tracer fits (inverse lognormal).

PAGE 131

117 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-16. Set 3 tracer fit (inverse lognormal). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-17. Set 4 tracer fit (inverse lognormal).

PAGE 132

118 sets 1 and 2 there is therefore a substantial increase in predicted dissolution accuracy when not incorporating travel time variability into the model. For simulations 1-2, and 1-3 and all of the simulation in set 2, even when not incorporating travel time variability into the second moment of the model still predicts more spreading of the dissolution profile than the UTCHEM output, suggesting that dispersive mixing acts to dampen the dissolution process. Table 4-4. Inverse lognormal derived mean and variance of contribution time. travel time variability no travel time variability Simulation RMS RMS 1-1 1.69 0.71 0.283 1.82 0.44 0.118 1-2 1.35 0.94 0.288 1.47 0.7 0.132 1-3 1.58 0.97 0.063 1.7 0.71 0.055 1-4 1.08 1.11 0.161 1.18 0.91 0.086 average 1.43 0.93 0.199 1.54 0.69 0.098 S.D. 0.27 0.17 0.108 0.28 0.19 0.034 2-1 1.25 0.73 0.322 1.33 0.58 0.178 2-2 0.9 0.98 0.538 0.97 0.85 0.396 2-3 0.89 1.05 0.454 0.97 0.89 0.302 2-4 1.46 0.76 0.219 1.53 0.6 0.070 average 1.13 0.88 0.383 1.20 0.73 0.237 S.D. 0.28 0.16 0.141 0.28 0.16 0.142 3-1 1.19 1.86 0.091 1.4 1.41 0.030 3-2 1.07 1.47 0.065 1.3 1.01 0.007 3-3 1.43 1.89 0.089 1.68 1.39 0.034 3-4 1.7 0.5 0.049 1.95 0 0.336 average 1.35 1.43 0.074 1.58 0.95 0.102 S.D. 0.28 0.65 0.020 0.29 0.66 0.157 4-1 1.66 1.27 0.050 1.92 0.74 0.048 4-2 1.48 1.3 0.044 1.77 0.75 0.190 4-3 1.84 1.37 0.532 2.17 0.7 0.175 4-4 1.97 2.02 0.099 2.27 1.42 0.070 average 1.82 1.65 0.075 2.10 1.08 0.059 S.D. 0.22 0.53 0.235 0.25 0.48 0.072 The increased accuracy of not including travel time variability for simulation sets 3 and 4 is less than simulation sets 1 and 2. This suggests that as the heterogeneity of the hydraulic conductivity field increases, the spreading of the nonreactive tracer becomes increasingly a result of macrodispersiondispersive mixing brought about by variability

PAGE 133

119 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-18. Set 1 predicted source depletion (inverse lognormal). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-19. Set 1 predicted source depletion (inverse lognormal-no travel time variability).

PAGE 134

120 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-20. Set 2 predicted source depletion (inverse lognormal). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-21. Set 2 predicted source depletion (inverse lognormal w/ no travel time variability).

PAGE 135

121 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-22. Set 3 predicted source depletion (inverse lognormal). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-23. Set 3 predicted source depletion (inverse lognormal w/ no travel time variability).

PAGE 136

122 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-24. Set 4 predicted source depletion (inverse lognormal). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-25. Set 4 predicted source depletion (inverse lognormal w/ no travel time variability).

PAGE 137

123 in the hydraulic conductivity fieldand less a result of spatially averaged mixing occurring at the individual grid cell. For simulation 3-4, the inverse lognormal approach predicted a variance of zero for the retardation coefficient, meaning that the predicted variance was zero. This prediction of no variability can be traced to the difficulty a lognormal distribution has at fitting data characterized by extensive tailing. For the case of a variance of zero, not incorporating travel time variability of course leads to a step function dissolution profile when travel time variability is subtracted (figure 4-22). S S S ADE dissolution predictions The ADE parameter estimation method was developed based on the performance of the method of moments and the inverse lognormal approach. Because the method of moments produced better estimates of the domain averaged NAPL saturation when compared to the inverse lognormal technique, the method of moments was used to determine parameters related to the first momentpore volume; domain average NAPL saturation. Because the inverse lognormal approach yielded more accurate fits of tracer spreading when compared to the method of moments, traditional curve fitting techniques were used to determine parameters related to tracer spreadingPeclet numbers. A conditional statement was also employed (see equation 4-20) to deemphasize travel time variability at low permeability field heterogeneity. The ADE tracer parameters are shown below in table 4-6 for all simulations. Tracer fits are shown in figures 4-25 to 4-28. Like the inverse lognormal approach that also employs curve fitting to obtain parameters related to tracer spreading, the ADE approach provides a better fit than the method of moments, especially at lower heterogeneity (sets

PAGE 138

124 1 and 2). Similar to the inverse lognormal approach, simulations 3-4 and 4-1 proved somewhat more difficult to fit due to substantial tailing in the BTCs. Unlike the inverse lognormal approach however, the ADE method fits the first part and tail portion of the BTC well while not fitting the middle portion well. The quality of the tracer fit is manifested in the predicted source depletion profiles where the model predicts the early time and tailing portions of the dissolution profile well, but over predicts the middle portion. Table 4-5. Tracer derived parameters for ADE method. nonreactive reactive Simulation P nr RMS P r R r RMS 1-1 43.5 0.003 18.9 1.7 0.010 1-2 50.7 0.006 32.7 1.3 0.005 1-3 45.1 0.004 31.6 1.2 0.005 1-4 36.2 0.010 23.0 1.2 0.019 average 43.9 0.006 26.6 1.3 0.010 S.D. 6.0 0.003 6.7 0.2 0.007 2-1 85.1 0.000 62.2 1.2 0.000 2-2 108.4 0.000 72.2 1.2 0.001 2-3 73.2 0.000 51.2 1.2 0.000 2-4 85.9 0.000 59.7 1.2 0.000 average 88.2 0.000 61.3 1.2 0.000 S.D. 14.7 0.000 8.6 0.0 0.001 3-1 8.1 0.006 5.4 1.2 0.010 3-2 8.0 0.006 6.3 1.2 0.005 3-3 6.3 0.004 4.5 1.2 0.004 3-4 6.2 0.006 5.3 1.2 0.018 average 7.2 0.006 5.4 1.2 0.009 S.D. 1.0 0.001 0.7 0.0 0.006 4-1 4.8 0.010 3.1 1.5 0.018 4-2 4.2 0.008 3.4 1.3 0.008 4-3 2.9 0.008 2.5 1.3 0.009 4-4 3.8 0.006 2.5 1.3 0.006 average 3.9 0.008 2.9 1.4 0.010 S.D. 0.8 0.002 0.5 0.1 0.005 The tracer determined source depletion parameters for the ADE model are displayed below in table 4-6 for all simulations along with the RMS value for the model fits. Model fits are shown in figures 4-31 to 4-34. The low RMS values for simulation sets 1 and 2

PAGE 139

125 suggests that the employed conditional statement (equation 4-20) does an effective job of deemphasizing pore scale dispersion at lower heterogeneities. The ADE approach also provides a reasonably good prediction of dissolution behavior for set 3 (see RMS values in table 4-6 and figure 4-33). As mentioned previously, the tracer fits for simulation 3-4 and 4-1, where the first and tailing portion of the tracer BTCs were fit well, manifested. itself in the predicted dissolution, where the first and tailing portions of the dissolution profile were well predicted. Simulation 4-3 proved to be somewhat of an anomaly, such that the parameterization methods entirely broke down. Despite a good fit of the tracer data, the predicted dissolution profile was significantly different than the actual dissolution profile. Similar to the other simulations in set 4, the tracers exhibited a large degree of spreading, suggesting spreading of the dissolution profile. The actual dissolution profile however was more consistent with the less heterogeneous simulations of sets 1 and 2. The initial architecture for simulation 4-3 was characterized by three large pools of similar lengths. It is possible that the close vertical proximity of these three pools, coupled with transverse vertical dispersion acted to dampen dissolution profile spreading. Aside from simulation 4-3, the ADE approach performed well, especially when considering a sophisticated numerical model has been replaced with an analytical model that can be parameterized by a combination of nonreactive and reactive tracers. Comparison of techniques The average predicted dissolution RMS value for each simulation set is shown below in figure 4-35 along with the average RMS for all simulations. Overall, each of the techniques performed well when considering that a sophisticated multiphase flow and transport simulator has been replaced with an analytical model that can be parameterized

PAGE 140

126 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-26. Set 1 tracer fit (ADE). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-27. Set 2 tracer fit (ADE).

PAGE 141

127 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-28. Set 3 tracer fit (ADE). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Model Fit (nonreactive tracer)UTCHEM SimulationModel Fit (reactive tracer)UTCHEM simulation Figure 4-29. Set 4 tracer fits. (method of moments).

PAGE 142

128 Table 4-6. Tracer determined source depletion parameters for ADE model. Simulation P R RMS 1-1 12.5 6.9 0.129 1-2 9.9 5.3 0.045 1-3 6.8 7 0.025 1-4 8 4.8 0.135 average 9.3 6.0 0.084 S.D. 2.5 1.1 0.057 2-1 9.8 4.5 0.016 2-2 24.3 3.58 0.070 2-3 13.6 4.2 0.076 2-4 12.9 4.88 0.056 average 15.2 4.3 0.055 S.D. 6.3 0.5 0.027 3-1 1.5 8.1 0.002 3-2 1.7 5.5 0.003 3-3 1.9 10.3 0.014 3-4 0.9 8.9 0.038 average 1.5 8.2 0.014 S.D. 0.4 2.0 0.017 4-1 1.7 11.3 0.052 4-2 0.77 8.3 0.785 4-3 0.33 12.6 0.031 4-4 1.3 20.9 0.079 average 1.0 13.3 0.237 S.D. 0.6 5.4 0.366 with partitioning tracers. The ADE approach, because it combined traditional moment analysis (which more accurately predicted the total NAPL mass) with curve fitting techniques (which more accurately predicted tracer spreading) and also included a conditional statement to deemphasize grid scale dispersion, had the lowest average RMS error. The inverse lognormal with no travel time variability and the method of moments had similar RMS errors for all simulations sets. The least accurate approach was the inverse lognormal approach. This is predominantly because of the first two simulation sets. In simulation sets 1 and 2 the inverse lognormal approach accurately predicted spreading of the tracers. Much of the tracer spreading however was attributable to grid scale dispersion, which does not impact the dissolution process, as opposed to

PAGE 143

129 macrodispersionspreading resulting from variability in the permeability fieldwhich does impact the dissolution process. To illustrate this, the longitudinal dispersivity was increased from the base case (0.1 m) to 1 m and to 10 m. As shown in figure 4-36, increasing longitudinal dispersivity leads to enhanced spreading of the tracer. Increasing longitudinal dispersivity does virtually nothing to the dissolution profile however, illustrating the problem that arises when interpreting tracer spreading entirely as velocity variability, as opposed to the combination of grid scale dispersion and macrodispersion. Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-30. Set 1 predicted source depletion (ADE). Conclusions and Suggestions for Further Research All of the techniques outlined above for predicting natural gradient source depletion perform well considering that a sophisticated multiphase flow and transport simulator with runtimes upwards to 30hrs has been substituted for with simple, compact

PAGE 144

130 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Source Depletion ModelUTCHEM Simulation Figure 4-31. Set 2 predicted source depletion (ADE). Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Cs Source Depletion ModelUTCHEM Simulation Figure 4-32. Set 3 predicted source depletion (ADE).

PAGE 145

131 Simulation 3Simulation 4Simulation 1Simulation 2Cumulative Pore VolumesC/Co Source Depletion ModelUTCHEM Simulation Figure 4-33. Set 4 predicted source depletion (ADE). 00.050.10.150.20.25Set 1Set 2Set 3Set 4 ADE Inv Log Inv Log MOM ADE Inv Log Inv Log (no t) MOM Figure 4-34. Root mean square errors for parameterization techniques. Points indicate average RMS for simulation set. Lines indicate average RMS for all simulations.

PAGE 146

132 00.00010.00020.00030.00040.00050.00060.00070.00080.00090.0010510152025 DL = 0.1 DL = 1 DL= 10 0100200300400500600012345 DL = 0.1 DL = 1 DL = 10 Pore VolumesCf(mg/L)Volume Fraction Dissolution Tracer Figure 4-35. Impact of increased longitudinal dispersivity on tracer response and dissolution profile. analytical models which can be parameterized with a combination of reactive and nonreactive tracers. The inverse lognormal technique does not perform as well as the other techniques at medium to low heterogeneity of the conductivity field, where the magnitude of pore scale dispersion (which does not impact the dissolution process in the simulator) is comparable to macrodispersion (which does impact the dissolution process). It is

PAGE 147

133 important to note however that the implicit assumption has been made that UTCHEM, because it has been validated with field and experimental data, represents some form of an underlying reality which we can compare with the streamtube approaches. In reality we are really just comparing two mathematical formulations, each with embedded assumptions. The issue of dissolution mechanics at the scale of the individual grid cell, that is whether or not the physical mechanisms that lead to dispersion at the grid scale impact the dissolution process, is an interesting area for further research. Soerens et al. (1998) illustrated that column scale experiments of NAPL dissolution interpreted as rate-limited (e.g., Imhoff et al., 1994) can also be interpreted as resulting from velocity variability. Based on the success of the tracers at parameterizing the simplified source depletion models, it appears that the fundamental assumption of the streamtube models, that dissolution dynamics are largely governed by the combined effects of the NAPL architecture and the velocity field are fundamentally correct, although a more systematic evaluation of the impact of neglecting mechanisms such as transverse dispersion and relative permeability is warranted. Likewise, the predicted trend of increased variability in either velocity and/or Lagrangian NAPL architecture leading to enhanced spreading of the dissolution profile appears to be fundamentally correct as well. Figure 4-37 below plots the variance of ln(K) vs. the variance of the travel time, variance of trajectory integrated NAPL content, and the variance of the contribution time. Based on the plots in figure 4-37, both and t increase initially with increased variability in ln(K). The travel time variance is highly correlated with the variance of the permeability field where the relationship between and ln(K) variance is more complex, increasing S S

PAGE 148

134 initially before appearing to asymptotically approach a maximum variance with increased heterogeneity. Previous numerical simulations have indicated a decrease in a penetration depth and enhanced spreading, quantified by spatial moments in the horizontal direction, with increased variance of the permeability field (Christ et al., 2005; Dekker and Abriola, 2000a; Lemke et al., 2004b; Mayer and Miller, 1996). There is likely some opposing interplay between penetration depth and lateral spreading of NAPL in terms of their effect on the trajectory integrated NAPL content. More spreading of the NAPL leads to an increase in variability provided there are vertical fingers present. A decrease in penetration depth compresses the architecture into relatively few streamtubes, minimizing the contrasts between vertical fingers and lateral pools which lead to a high variance. Figure 4-38 plots the variance of travel time and vs. the contribution time variance. As evident in figure 4-38, the contribution time is more highly correlated with the distribution. S S S S Additional dissolution experiments with tracers would be useful for systematically evaluating some of the assumptions of the streamtube models, especially the issue of grid scale dispersion outlined above. One such experiment was conducted here and is outlined briefly below. The experiment was conducted with the in-line flow cell system discussed in chapter 3. 10 ml of TCE were introduced into the domain, yielding the initial distribution shown in figure 4-39 below. After the NAPL was spilled, two nonreactive tracers (methanol and TBA) were introduced into the domain along with two partitioning tracers (DMB and hexanol) at a flow rate of 5 ml/min for pulse duration of 27 min. The tracer BTCs are shown in below in figure 4-40. Based on the quality of the BTC data, TBA and hexanol were chosen for

PAGE 149

135 moment analysis (figure 4-40). Due to the small separation between the nonreactive and partitioning tracers, the other techniques (ADE, inverse lognormal) were not effective. It is suggested that these techniques are applicable only for step injections where the partitioning coefficient is large enough where there is separation in the rising limbs of the BTCs. The moment derived parameters for this experiment are listed in table 4-7. 00.10.20.30.40.50.60.700.10.20.30.40.50.60.70.8 00.20.40.60.811.21.41.600.10.20.30.40.50.60.70.8 00.511.522.500.10.20.30.40.50.60.70.8 Variance ln(K)Variance contribution timeVariance NAPL traj. int. NAPL contentVariance travel time Figure 4-36. Relationship between ln(K) and contribution time parameters.

PAGE 150

136 00.511.522.500.10.20.30.40.50.60.7 00.511.522.500.20.40.60.811.21.41.6 variance traj. int. NAPL contentvariance travel timevariance contribution time Figure 4-37. Relationship between travel time,, and contribution time variability. S Figure 4-38. Initial TCE distribution for dissolution experiment with tracers.

PAGE 151

137 The experimental dissolution data along with the moment predicted dissolution profile are shown below in figure 4-41. The method of moments provides a good fit of the data. There is a slight over of total mass, likely due to experimental mass losses. Future experiments would be beneficial in order to further evaluate the streamtube models and the parameterization techniques. 0.111010010000.1110 TBA Hex 0.11101001000100000.1110 Meth TBA 23DMB Hex 0.111010010000.1110 TBA Hex 0.11101001000100000.1110 Meth TBA 23DMB Hex Figure 4-39. Tracer BTCs (top graph) and tracer BTCs selected for moment analysis (bottom graph).

PAGE 152

138 Table 4-7. Moment derived parameters for dissolution experiment with tracers. paramtervaluet-0.03t20.06S-4.67S20.783.490.84 2 Pore VolumesC/Cs Figure 4-40. Experimental and predicted dissolution data.

PAGE 153

CHAPTER 5 EXTENDING SOURCE DEPLETION MODELS TO FLUSHING BASED REMEDIATION In this chapter, the source depletion models of chapter 2 are extended to flushing based remediation. Both surfactant and cosolvent based remediation are considered. Like the simplified source depletion models, these models are perhaps of more practical utility to contaminated site managers based on the difficulty of parameterizing more sophisticated numerical multiphase flow models. The equations below are modified versions of the equations presented in Jawitz et al. (2005). The purpose of this chapter is similar to chapter 2 where the models are evaluated based on their ability to reproduce typical behavior from experimental and numerical studies. Theory Streamtube Models The streamtube source depletion models of chapter 2 can be extended to flushing based remediation by incorporating the hydrodynamics of the flushing solution which involves lagging the initialization of the flushing-based source depletion process until the flushing solution has traveled from injection well to extraction well, such that the dissolution profile for an individual streamtube is as shown in figure 5-1. Incorporation of flushing solution hydrodynamics involves two relatively straightforward steps. The first step is adding the travel time of the flushing solution (t fl ) to the contribution time calculated using equation 2-12. Recall that equation 2-12 is a simple mass balance equation on an individual streamtube, dividing the total mass in each 139

PAGE 154

140 CflCaqarrival of flushing solutionAqueous Dissolution Flushing Cflstreamtube depleted of massCumulative Pore Volumes' Cf Figure 5-1. Dissolution profiles for aqueous dissolution and flushing based remediation. streamtube by the mass discharge to arrive at the time an individual streamtube will contribute contaminant flux under a set of specific flushing solutions. The only difference between flushing and aqueous dissolution, as it relates to equation 2-12, would be a large decrease in K f due to the large increase in contaminant solubility. After adding the arrival time of the flushing solution (second term in equation 5-1), the second step is to subtract the mass that is discharged from the streamtube at aqueous solubility prior to the arrival of the flushing solution (third term of equation 5-1). The mean contribution time, referred to hear as is then given by: flfltflaqtmCCmmm1111 (5-1) Where C aq is the solubility of the contaminant in water and C fl is the solubility of the contaminant in the flushing solution. The main difference between Jawitz et al.

PAGE 155

141 (2005) and the equations presented here is that the second moment of the remedial fluid travel time is not incorporated into the second moment of This modification was introduced based on the experiments outlined below where the first pore volume experienced significantly more spreading than that of a nonreactive tracer due to density contrasts between the remedial fluid (50% ethanol and water mixture) and the resident water. The second moment of still includes the second moment of what would be a nonreactive tracer however. After Jawitz et al. (2005), the BTC for flushing based remediation can be represented using cumulative distribution functions: ))()(())(1()( ptpCtpCfTCrfflrfaqcf (5-2) where the subscript rf is used to denote remedial fluid (note that the subscript fl is used to refer to the contaminant concentration in the flushing solution where the subscript rf is used to denote the concentration of the flushing solution/remedial fluid). The term rftp1 of equation 5-2 is the fraction of streamtubes where the flushing solution has not yet reached the extraction well. This term is multiplied by C aq as it is assumed that these streamtubes are discharging a concentration equal to the aqueous solubility. The term ptprf is the fraction of streamtubes where the flushing solution has reached the extraction well yet the duration of flushing has not been long enough to entirely deplete the streamtube of mass. This term is multiplied by C fl as it is assumed that these streamtubes are discharging a concentration equal to the solubility of the contaminant in the flushing solution. For cases where C fl is much larger than C aq the third term of 5-1 can be dropped with negligible error and 5-2 can be approximated by:

PAGE 156

142 )()()(ptpCfTCflflcf (5-3) Flushing BTCs for a fixed NAPL mass, fixed flushing solution breakthrough, and increasing variability in contribution time are displayed below in figure 5-2. Similar to the aqueous dissolution model, an increase in the second moment of resulting from either increased variability in either the t or distributions, leads to more spreading/ tailing of the enhanced dissolution BTC. S 01234500.20.40.60.81 Flushing Solution m2 = 4.1 m2 = 5 m2 = 8Flushing DurationC/Co Figure 5-2. Flushing BTCs for increasing second moment of Advection Dispersion Model As discussed in chapter 2, the cdfs in equation 5-2 can be replaced with the analytical solution for the advection dispersion equation (equation 2-21) such that 5-2 can now be expressed as:

PAGE 157

143 )()()(TCTCfCTCNAPLrfcflf (5-3) where C rf (T) is the ADE describing the flushing breakthrough (with a retardation factor and Peclet number denoted by R rf and a P rf ) and C NAPL (T) is the ADE (with a retardation factor and Peclet number denoted by R NAPL and a P NAPL ) used in a manner analogous to the distribution to describe variability in the combined effects of the Lagrangian NAPL architecture and the flow field. Similar to equation 2-23, integrating 5-3 from zero to infinity yields: (5-4) NAPLrfNAPLrfRRdTTCTC0)()( which allows for expression of total mass as in equation 5-5 below. NfrfNAPLSKRR (5-5) Flushing BTCs for a fixed NAPL mass, fixed flushing solution breakthrough, and increasing variability in P NAPL are shown below in figure 5-3. Similar to the aqueous dissolution model, an increase in P NAPL leads to an increase in BTC spreading/tailing. Extensions for Dispersive Mixing Under the advection-only assumption, the concentration of the remedial agent discharging from each streamtube is zero until the arrival time of the remedial fluid, upon which there is step increase from zero to the injected concentration (figure 5-1). For example, if a 5% by volume surfactant solution was injected into the source zone, the volume fraction of surfactant discharging from an individual streamtube is assumed to be zero until the remedial fluid travel time, upon which there is a step increase in volume fraction of surfactant from 0% to 5%. The spreading that is observed in the integrated

PAGE 158

144 extraction well BTC under the advection only assumption is attributed solely to variability in the remedial fluid travel times associated with each streamtube. The major 01234567891000.20.40.60.8 Flushing Solution P = 100 P = 10 P = 2 P = 0.5Flushing DurationDimensionless Concentration Figure 5-3. Flushing BTCs for increasing second moment of computational advantage of the advection only approach is that there is no need to incorporate equations relating the solubility of the contaminant to the concentration of the remedial agent; the solubility is either at aqueous solubility (prior to flushing solution arrival time) or at flushing base solubility (after flushing solution arrival time). Any type of dispersive mixing at the leading edge of the remedial fluid front, such that transport at the streamtube scale is no longer characterized by a step increase, causes a separation in the flushing solution and contaminant BTCs relative to the advection only case. This

PAGE 159

145 enhanced separation was observed in the experiments below, requiring modifications of the above equations to incorporate mathematical relationships describing contaminant solubility as a function of remedial fluid concentration. To incorporate the effects of transverse mixing, equation 5-3 was modified as follows. For cosolvent flushing, the relationship between contaminant solubility and the ethanol fraction (f) was assumed to be log-linear(Augustijn and Rao, 1995): (5-5) faqflCfC10)( where is the cosolvency power. The ethanol fraction as a function of flushing duration (f(T)) can be defined using equation 5-6 below. )()(TCfTfrfo (5-6) Where f o is the ethanol fraction injected. A solubility enhancement term describing the normalized increase in contaminant solubility as a function of the ethanol fraction is then defined as: ofTfflcCfTfC1010))(()( (5-7) The flux averaged concentration exiting the source can then be expressed as: )(11010)()(TCCfTCNAPLfTfflco (5-8) Because of the nonlinear nature of the solubility enhancement term, a closed form solution for R NAPL is not available and R NAPL must be determined numerically after defining f(T) to ensure the appropriate mass balance. Similar to 5-5, an equation relating contaminant solubility to surfactant concentration can be defined as in 5-9 below (Pennell et al., 1997):

PAGE 160

146 surfaqsurfflmCCCC )( (5-9) where the subscript surf is used to denote surfactant and m is the slop of a solubilization line. Similar to equation 5-6, the surfactant concentration can be defined as: )()(_TCCTCrfosurfsurf (5-10) where C surf_o is the surfactant concentration in the injected remedial fluid. Similar to 5-7, a solubilization enhancement term can be defined as: osurfaqsurfaqflcsurfmCCTmCCCfTCC_)())(( (5-11) Equation 5-8 can now be defined as: )(1_)()(TComCCTmCCCfTCNAPLsurfaqsurfaqflc (5-12) Experimental Data from two cosolvent flooding experiments were modeled with the source depletion models discussed above. The experiments are reviewed only briefly below as the focus here is on the modeling effluent BTCs. Experimental Setup The same flow cell that was used for the DCA-3 aqueous dissolution experiment (see chapter 3) was used for the cosolvent flooding experiments. The porous media packing was similar to experiments TCE-1 and DCE-1 with discrete layers of 40/50 sand in background 20/30 sand. 15 milliliters of PCE for experiment one and 10 ml of PCE for experiment 2 were injected into the upper portion of the domain using a syringe pump. After PCE injection, several pore volumes of water were flushed through the system to obtain a pre-remediation measurement of contaminant flux. After the pre-remediation

PAGE 161

147 water flood, the flushing solution (50% ethanol/50% water) was injected into the flow cell at a rate of 5ml/min until the effluent PCE concentrations declined to 10% of PCE solubility in the flushing solution (C fl ), at which point water was flushed through the box until all of the ethanol had been removed in order to obtain a post-remediation measurement of mass flux. Experiment 1Experiment 2 Figure 5-4. Initial PCE distribution for experiments 1 and 2.

PAGE 162

148 Initial Architectures The initial distribution of PCE following cessation of the migration phase is shown in figure 5-4. The initial architectures were consistent with the experiments detailed in chapter 3. The darker areas in figure 5-4 are the lower permeability lenses. Prf= 20.5Rrf= 1.43Prf= 5.71Rrf= 1.02 Figure 5-5. Model fit of ethanol BTCs for experiments 1 (top graph) and 2.

PAGE 163

149 Modeling of Effluent BTCs The first step in modeling of the cosolvent flood was to fit the ethanol breakthrough with the C rf function. Model fits, along with the optimized values of P rf and R rf are displayed in figure 5-5. The fraction of streamtubes contaminated (f c ) was estimated based on the pre-remediation flux-averaged concentration, which was approximately 100 mg/L. A PCE solubility of 150 mg/L was assumed, such that f c was equal to 0.66. A cosolvency power of 3.5 for ethanol and PCE was determined based on data from Cho (2001). The model fit for experiment 1 is shown below in figure 5-6 along with the optimized values for P NAPL and R NAPL As shown in figure 5-6, the model provides an excellent fit of the observed data with an R 2 value greater than 99%. PNAPL = 4.0RNAPL = 4.3 Figure 5-6. Model fit of experiment 1 data.

PAGE 164

150 The equivalent natural gradient (or pump and treat) source depletion for experiment 1, calculated with equations 2-1 and 2-21 using the optimized value of P NAPL and a retardation factor calculated with 2-23, is shown below in 5-7. Figure 5-7. Equivalent aqueous dissolution source depletion for experiment 1 The same process outlined above was used to fit the data from experiment 2. The initial flux-averaged concentration of experiment 2 was approximately 130 mg/L, such that f c was 0.87. The model fit and optimized parameters are shown below in figure 5-8. The model again provides an excellent fit of the observed data with an R 2 value greater than 99%. The equivalent natural gradient (or pump and treat) source depletion for experiment 2 is shown below in figure 5-9.

PAGE 165

151 PNAPL = 2.7RNAPL = 3.0 Figure 5-8. Model fit of experiment 2 data. Figure 5-9. Equivalent aqueous dissolution source depletion for box 2.

PAGE 166

152 Remedial Endpoints One of the primary advantages of utilizing solubilization-based remedial technologies and orienting arrays of injection and extraction wells normal to the mean direction of groundwater flow is that after the flushing solution has completely broken through, the mathematics describing the remedial process and that of natural flowing groundwater are essentially identical. The implication of this is that the equivalent natural gradient source discharge can be predicted by scaling the flux-averaged contaminant concentration in the extraction well. For experiments 1 and 2, cosolvent injection was terminated once the effluent concentration reached 10% of its maximum solubility in the cosolvent solution (10,000 mg/L). As shown in figure 5-10 below, the natural gradient flux-averaged concentration once all of the cosolvent had been flushed out of the system was below 10% of aqueous solubility (assumed 150 mg/L). This suggests that orienting injection and extraction well arrays normal to the mean direction of groundwater flow provides an effective way to decide when to terminate remediation. Terminating the enhanced solubilization portion of the flood (cosolvent or surfactant injection) once the scaled concentration in the extraction well has declined to the target source discharge still leaves an additional pore volume of flushing solution in the source zone to account for discrepancies in mass transfer rate coefficients and residence time. Conclusions and Suggestions for Further Research The work above suggests that similar to aqueous dissolution, the simplified enhanced dissolution models can serve as computationally efficient analogs to more sophisticated numerical simulators. The models are ideal for exploring uncertainty in key input parameters, such as NAPL architecture and velocity variability. The models of Jawitz et al. (2005) were modified to incorporate experimentally observed loglinear

PAGE 167

153 01020304050607080901001101201300510152025 01020304050607080901001101201301401501601700510152025 Experiment 1 Pore VolumesConcentration (mg/L) Experiment 2 Scaled concentration (C/Cs) at transition to water flood Figure 5-10. Pre and post flux-averaged concentrations in experiments 1 and 2. effects thought to result from mixing processes at the leading edge of the injected remedial fluid front. The above work would benefit from simulations of cosolvent or surfactant based remediation and additional experiments with which to compare the simplified models. Similar to the aqueous dissolution models, a more in depth analysis of the implications of neglecting dispersion, remedial fluid properties, relative permeability, and mass transfer rate coefficients is an area for further research. The difficulty in comparing UTCHEM surfactant simulations with equation 5-12 is that UTCHEM uses a ternary phase diagram to determine contaminant solubility (Delshad et al., 1996), such that the solubility of the contaminant is dependent on the local NAPL saturation in addition to the surfactant concentration. The MALVOR model

PAGE 168

154 (Dekker and Abriola, 2000b; Lemke and Abriola, 2003; Lemke et al., 2004b) and a modified MT3D/MODFLOW code (Saenton et al., 2002) contain formulations similar to 5-9 and are perhaps more appropriate for evaluating the above models. Table 5-1. Fitted model parameters from UTCHEM surfactant simulations. SetSimulationmPNAPLRNAPL210.72108.81.7220.6681.41.6310.562.9330.413.23.2410.661.61.9420.53.62.7 Despite the fact that UTCHEM does not employ a relationship similar to 5-9, equation 5-12 was fit to surfactant flood simulations conducted on two spills from each of simulation sets 2, 3, and 4 for comparison purposes. Parameters describing the ternary phase diagram were adapted from a UTCHEM example file previously developed to model the column data of Mayer et al. (1999). Because the ternary phase diagram describes the relationships between NAPL solubility and the local surfactant concentration and not equation 5-9, the parameter m was treated as a fitting parameter. Table 5-1 displays the results from the curve fitting process. Figure 5-11 displays the model fits. As shown in figure 5-11, the curve fits are reasonable despite the fact that UTCHEM does not include a relationship similar to equation 5-9 to model surfactant induced solubility enhancement. Similar to the aqueous dissolution simulation, simulation 1 of set 4 exhibits a bimodal behavior with two distinct portions of the contaminant BTC, one associated with the residual component and another associated with the late term pooling portion.

PAGE 169

155 Model Fit (surfactant breakthrough)UTCHEM SimulationModel Fit (TCE breakthrough)UTCHEM simulation Set 2; Sim1 Set 2; Sim1Set 3; Sim1Set 3; Sim3Set 4; Sim1Set 4; Sim2Pore VolumesC/Cs Figure 5-11. Model fits of UTCHEM surfactant floods.

PAGE 170

CHAPTER 6 SYNTHESIS: A FLUX BASED FRAMEWORK FOR MANAGING DNAPL CONTAMINATED SITES Defining Efficiency and Favorability Three Definitions of Efficiency Prior to conducting any type of contaminated site remediation, it is prudent to define the desired outcome, or the objectives of the remedial effort. Ideally, these objectives should be quantifiable, such that upon completion of remediation the employed technology can be evaluated for its effectiveness in terms of meeting the defined performance objectives. A recent EPA report (EPA, 2003) grouped source zone remediation performance metrics into three categories. Type I metrics are measurements that can be reliably acquired and are commonly used such as the total mass of DNAPL removed from the subsurface and changes in resident groundwater concentrations. Type II metrics are metrics that can sometimes be measured, but are not in wide use such as remaining DNAPL mass and DNAPL architecture. Type III metrics are metrics that are theoretically possible and under development such as mass flux and mass discharge. The above performance metrics can be generalized as metrics relating to resident groundwater concentrations in the source zone, metrics relating to contaminant mass, and metrics relating to contaminant flux or source discharge. As discussed in the opening chapter, a remedial objective of reducing local resident groundwater concentrations below the regulatory MCL is almost always, without question, technologically infeasible. 156

PAGE 171

157 When utilizing mass removal based performance metrics, an efficiency assessment for a particular technology would quantify the mass removed per unit work input (e.g., financial expenditures, number of pore volumes, etc. etc.). As an example consider three hypothetical source zones, source zone A, B, and C, which are characterized by sequentially increasing degrees of spatial variability in NAPL architecture (figure 6-1) (Note that the curves below were generated using the ADE model outlined in chapter 2). Implementing a pump and treat system and considering only mass reduction efficiency, source zone A would be considered the most favorable, yielding the largest reduction in mass for a given work input (figure 6-1). Cumulative Work InputCumulative Mass Removed Source Zone ASource Zone BSource Zone C Figure 6-1. Comparison of source zones A, B, & C using mass reduction efficiency.

PAGE 172

158 If the ultimate goal is a reduction in risk however, the cumulative mass removed tells us little about risk reduction unless it can be correlated in someway to the expected reduction in mass discharge. An alternative definition of efficiency to consider then is the reduction in source discharge for a given reduction in mass. Using this definition of efficiency, source zone C would now be considered the most favorable where source zone A would now be considered least favorable (figure 6-2); which is the opposite conclusion one would make when considering mass reduction efficiency. Because cumulative mass removed is the time integral of source discharge, the most favorable systems in terms of mass removed per unit input of work are, by definition, the least favorable in terms of source discharge reduction per unit reduction in mass. Cumulative Reduction in MassCumulative Reduction in Source DischargeSource Zone ASource Zone BSource Zone C Figure 6-2. Comparison of source zones A, B, & C using reduction in source discharge per reduction in mass.

PAGE 173

159 From the perspective of a site manager who is primarily concerned with risk reduction, the definition of efficiency which is most applicable is the reduction is source discharge per unit input of work. Using this definition of efficiency, the question of which source zone is more favorable becomes time dependant (figure 6-3). At early stages, the more heterogeneous source zone C is more favorable, realizing a larger reduction in source discharge per unit work input. At later stages, as the high flux producing mass is preferentially removed and source zone C enters into its later term tailing stage, source zone A becomes more favorable, realizing a larger reduction in source discharge per unit work input. Cumulative Work InputCumulative Reduction in Source DischargeSource Zone ASource Zone BSource Zone C Figure 6-3. Comparison of source zones A, B, & C using reduction in source discharge efficiency.

PAGE 174

160 The implication of the time dependency of the source discharge efficiency metric is that for increased aging (i.e. natural gradient dissolution occurring after the initial contamination event and before remediation), the relationship between mass and source discharge becomes less favorable for highly heterogeneous systems and more favorable for relatively homogenous systems. As source zone heterogeneity increases, more mass is shifted to the late term tailing portion of the natural gradient (or pump and treat) dissolution profile, such that there is a large initial decrease in source discharge for a relatively small reduction in mass (figure 6-2; figure 6-3). Relatively early in the aging process of highly heterogeneous systems, the high flux producing mass is consumed by the natural gradient dissolution process, such that the relationship between mass and source discharge becomes less favorable with aging (figure 6-4). The opposite is true for the more homogenous source zone, source zone A; with increased aging the initial stage of constant source discharge is consumed by the natural gradient process, such that the relationship between mass and flux becomes more favorable with aging (figure 6-4). For source zone B, which for no aging is defined by an approximate 1:1 relationship between mass and source discharge, aging has minimal impact on the relationship between mass and flux. Not considering time scale dependency and expressing mass and source discharge reduction in a dimensionless form, as in figure 6-4, becomes extremely limiting when considering the implications of aging on source zone remediation. From figure 6-4, the relationships between reductions in mass and source discharge moves closer to the 1:1 line with increased aging. Assuming that source zone C remains a good candidate for remediation after it has aged appreciably based solely on the fact that the relationship

PAGE 175

161 Cumulative Reduction in MassSource Zone CSource Zone BSource Zone ACumulative Reduction in Source Discharge Increased AgingIncreased Aging Figure 6-4. Aging effects on relationships between reduction in mass and reductions in source discharge for source zone A, B, & C.

PAGE 176

162 between mass and flux is still favorable would be a mistake. Considering aging effects on source discharge efficiency (source discharge reduction per unit work input) one would come to a completely different conclusion regarding the appropriateness of source zone remediation. The scenarios considered in figure 6-4 are plotted in terms of source discharge efficiency below in figure 6-5. Considering source discharge efficiency, it is apparent that source zone C becomes highly unfavorable relatively early in the aging process as all of the high flux producing mass is removed and the system enters into its late term tailing stage where there is little reduction in source discharge per unit input of work. The trending of the mass reduction/source discharge reduction relationships towards the 1:1 line is primarily a property of the ADE and the lognormal functions, lacking a physical meaning that is readily deducible from currently available literature. Mathematically, this suggests that the latter portion of the ADE and lognormal BTCs are defined by exponential decay. Recall from chapter 2 that a 1:1 relationship between mass and source discharge is defined by an exponential source discharge profile where the source discharge (S D ) is defined by: oDoMTSoDDeSTS)( (6-1) where S Do is the initial source discharge [M/T] and M o is the initial mass [M]. Equation 6-1 is useful for demonstrating that a 1:1 relationship between mass and source discharge can still be characterized by highly unfavorable source discharge reduction efficiency if the ratio of the initial source discharge to the initial mass is very high. A suite of source discharge profiles are shown in figure 6-6 for a fixed initial mass and alternate initial source discharges. In cases where the ratio of initial source discharge to the initial mass is

PAGE 177

163 Cumulative Pore VolumesCumulative Reduction in Source DischargeSource Zone ASource Zone BSource Zone C Increased Aging Increased Aging Increased Aging Figure 6-5. Aging effects on source discharge efficiency for source zone A, B, & C.

PAGE 178

164 Source DischargeCumulative Pore Volumes Figure 6-6. Source discharge profiles for a fixed initial mass and alternative initial source discharge. very low, the source discharge profile closely resembles that of a step functionthat is the least favorable mass reduction/source discharge reduction relationship. This is essentially true for all relationships between mass and source dischargewhen the ratio of the initial source discharge to the initial mass is very low, there is minimal discrepancy in the source discharge efficiencies regardless of what the mass reduction/source discharge reduction relationship is. Combining Three Efficiency Definitions The dimensionless plots of reduction in source discharge vs. reduction in mass are limiting in that they do not consider time scale dependency and secondly, they are dimensionlessthat is they are divorced from the site specific target for source discharge reduction and absolute values of source discharge reduction. Clearly, if the source discharge is already below the site specific target, an expensive source zone remediation effort is unnecessary. By incorporating tie lines of equivalent work input (e.g. pore volumes; financial expenditures) and the target reduction in source discharge, the plot

PAGE 179

165 displayed in figure 6-2 now encompasses the three definitions of efficiency (mass reduction/work input; source discharge reduction/mass reduction; source discharge reduction/work input) and is connected tangentially to information relating to the site specific target reduction in source discharge (contaminant fate and transport, exposure pathways) (figure 6-7). Instead of incorporating the site specific target for source discharge reduction as in figure 6-7, figure 6-2 could be expressed in a dimensional form (i.e. source discharge and mass reduction take on actual values) with the target source discharge denoted on the y-axis in some fashion. Figure 6-7. Figure 6-2 modified to include net work input and site specific flux-based performance metric. Flux Based Characterization Source Identification (Pre-Reactor Based Characterization) The utilization of tracers for the parameterization of the source depletion models is expected to be implemented in the later stages of site characterization where a source has

PAGE 180

166 been identified and preliminary contaminant flux measurements have indicated that the source discharge is large enough to warrant potential remediation. The source zone during tracer analysis is idealized as functioning in an analogous manner to a chemical reactor, with inflows and outflows to the reactor controlled and analyzed for relevant chemical constituents on site. There are of course several preliminary, less cost intensive steps prior to the establishment of reactor-like conditions. For discussion purposed I will group all of the traditional phase I & II type characterization (document (prior use) review, soil coring, groundwater samples, etc. etc.) into pre-reactor characterization. Transitional Reactor-Based Characterization: Mass Flux as a Remedial Screening Tool and Qualitative Estimator of Aging The decision to move from pre-reactor characterization to more cost intensive reactor-based characterization should be determined based on the perceived level of risk posed by a given site, which, as discussed previously, involves an estimate of source discharge coupled with predictive contaminant fate and transport analysis to determine exposure pathways and potential contaminant receptors. Mass flux measurements should function not just as a remedial performance assessment tool, but as an integral component in the decision to conduct source zone remediation and in the remedial design process. Technologies such as the passive flux meter (PFM) (Hatfield et al., 2004) which are capable of providing flux-based measurements without reactor-based infrastructure and without generating large volumes of sampling derived waste, are ideal for preliminary flux measurements and should function as a transitional characterization technology between pre-reactor characterization and reactor-based characterization. Aside from an estimate of source discharge, the PFM also provides a qualitative estimate of aging, which is a critical component to consider when deciding whether to

PAGE 181

167 implement a source zone remediation technology. Figure 6-8 below displays the source discharge profile and down gradient flux planes for an arbitrary source zone (simulation conducted using UTCHEMsee chapter 2). After the initial spill, the flux plane is well connected along the initial NAPL migration pathway. As the source zone ages, the residual component is consumed by natural gradient dissolution, such that the flux plane becomes discontinuous, characterized by isolated hot spots associated with up-gradient pools. Utilized in this way, spatially resolved flux information functions as a valuable screening tool, especially for flushing based remediation where successful remediation is dependent upon the same flux generation mechanisms as natural gradient dissolution (e.g. hydrodynamic accessibility to the NAPL interface). Put in another way, high flux measurements at the down-gradient flux measurement plane connote up-gradient flux producing mass that can be effectively targeted and removed with flushing based technologies. Discontinuous, lower flux measurements connote isolated pools which are likely to be significantly more difficult to remediate, requiring several additional pore volumes to achieve appreciable reductions in flux, which perhaps suggests implementation of a more passive technology with a lower flux assimilation capacity, such as enhanced bioremediation. In relation to predicting remedial performance it is important to emphasize the difference between spatially resolved mass flux measurements and spatially integrated source discharge measurements. A highly aged source zone that spans a large cross sectional area could have a high source discharge yet be characterized by low point flux measurements. High flux and not high source discharge indicates up-gradient

PAGE 182

168 hydrodynamically accessible mass that can be effectively targeted with flushing based remedial technologies. 010203040506070809010001020304050Dissolution Time (Yrs) Source Discharge (g/day) Mixed length residual Laterally Ext. pooling pools 40% initial mass20% initial flux Figure 6-8. Source discharge profile and down-gradient flux planes at selected time steps. Reactor-Based Characterization If the transitional flux measurement is high enough to suggest potential source zone remediation, the next step should be the establishment of reactor-like conditions. The establishment of reactor-like conditions, such that the inflow and outflow to the source zone are controlled, provides immediate source stabilization where continued contaminant mass discharge from the source zone to the dissolved plume is arrested. It also provides an integrated measurement of source discharge (Einerson and Mackay, 2001) in addition to the PFM. The question of whether the integrated measurement provides a more accurate estimate of source discharge, such that it should serve as the

PAGE 183

169 final step in the decision process of whether source zone remediation is warranted, is dependant upon the spatial variability of the flux plane. If the spatial variability is high enough where the PFM misses hot spots that are between sampling locations, the integrated measurement is likely to provide a more accurate measurement of the total source discharge. I have not evaluated this question in any systematic fashion here. If it is eventually established that the PFM and integrated steady state measurement are comparable, then moving immediately to the reactor-based characterization discussed below, as opposed to first obtaining an integrated measurement of source discharge, may be appropriate. Of critical importance in remedial design is the successful delivery of remedial constituents to the NAPL interface. In the flux-based paradigm, the delivery of remedial fluids to the high flux producing components of the source zone is especially important. The mapping of flux producing mass, a priori, in terms of its location in relation to injection and extraction wells can be accomplished through utilizing the spatially resolved flux information from the PFM and assuming the high flux producing mass is immediately up gradient, or, perhaps more accurately, injecting a unique nonreactive tracer into each injection well and using the tracer response at the extraction wells, together with the contaminant concentration of interest, to map the location of high flux producing mass. The utilization of spatially resolved flux information for flux mapping requires two key assumptions (figure 6-9). The first assumption is that the flux producing mass is immediately up-gradient. The second assumption is that remedial fluids can be focused towards the high flux producing mass by injecting fluids immediately up-gradient and extracting immediately down-gradient (figure 6-9). The above assumptions

PAGE 184

170 are likely less tenable as heterogeneity in the porous media increases, increasing the likelihood for more tortuous pathways between injection and extraction wells. Trace flux back to source directly up-gradient?Assumption 1: Flux producing mass is immediately up-gradientAssumption 2: Pairs of injection and extraction wells are in direct communication with each other High flux producing NAPL mass injectorsextractors Figure 6-9. Key assumptions required when using spatially resolved flux information for flux mapping. A more rigorous flux mapping approach involves the injection of a unique nonreactive tracer into each injection well. The injection of unique tracers allows for a measurement of communication between injection and extraction wells and a pseudo flux measurement based on the contaminant concentration measured in the extraction well (figure 6-10). Based on the extraction well data depicted in figure 6-10, the majority of the flux is discharging out of extraction well 2 (EW2)which is in good communication with injection well 2 (IW2), such that the preferential injection and withdrawal into and out of IW2 and EW2 is likely to lead to an increase in flux reduction efficiency. Likewise, there appears to be virtually no flux producing mass in communication with

PAGE 185

171 EW3, suggesting that minimal, if any remedial fluids should be injected and withdrawn into and out of IW3 and EW3. A key assumption to the above statement is that the remedial fluid travels on similar paths as the nonreactive tracer. The influence of remedial fluid properties on the applicability of using nonreactive tracers to map communication between injection and extraction wells is an area for further research. 0510152025300246 05101520250246 0510152025300246EW1 Unique tracers used for mapping communication between injector and extractor pairs (after initial baseline flux measurement) PCEIW1 tracerIW2 tracerIW3 tracerUse NAPL EW concentration for pseudo flux measurementC(mg/L)Pore VolumesEW2EW3 injection wells extraction wells Figure 6-10. Use of unique nonreactive tracers for flux mapping. In addition to the injection of unique nonreactive tracers, partitioning tracers can also be injected in order to provide an estimate of the total NAPL mass (Jin et al. 1995) and to parameterize the source depletion models outlined in chapter 2. The question of whether or not to inject partitioning tracers once reactor like conditions have been established is primarily an economic one. Based on the arguments outlined above, flushing based remediation should be conducted based on high flux measurements, which

PAGE 186

172 connote up-gradient high flux producing mass which can be effectively targeted. Information derived from partitioning tracers would help to more accurately estimate the number of pore volumes of flushing solution required to meet the source discharge reduction objective. Instead of utilizing partitioning tracers, uncertainty relating to the NAPL mass and spatial distribution could be incorporated into the design process by allowing for additional flushing solution. Again, this is primarily an economic question where the cost of conducting a partitioning tracer test must be evaluated in terms of whether the enhanced prediction capabilities of the remedial process provided by the tracers allow for savings in the volume of flushing solution purchased which would offset the cost of the tracer test. Remedial Endpoints and Site Closure An integrated measurement of source discharge, as outlined above, would of course require that the extraction wells be oriented in the mean direction of groundwater flow. It is suggested that orienting the well configuration in the mean direction of groundwater flow and implementing flushing based remediation significantly simplifies the decision of when to terminate remediation. After complete remedial fluid breakthrough, such that the concentration of the remedial agent (e.g. surfactant or cosolvent) extracted equals that injected, the mathematics describing the remedial process are identical to those describing natural gradient dissolution. The implication of this is that the concentration of the contaminant in the extraction well during remediation can be used to predict the natural gradient source discharge by accounting for discrepancies in mass transfer rate coefficients, solubility, and velocity between remediation and natural gradient conditions. It is suggested that terminating remedial fluid injection and transitioning to water injection once the extracted contaminant concentration is equal to the scaled site specific

PAGE 187

173 source discharge target is a conservative definition of a remedial endpoint as there is still an entire pore volume of flushing solution in the source zone that must be extracted. It is also suggested that in a mass discharge type framework, it is appropriate to maintain reactor like conditions (source stabilization) until the system returns to natural gradient conditions and it is demonstrated that the source discharge has been reduced below the target level. Figure 6-7 as a Conceptual Site Model and Design Tool The characterization outlined above is sufficient to parameterize figure 6-7 using the techniques outlined in subsequent chapters. Figure 6-7 can function as a useful conceptual site model with which to evaluate remedial alternatives. The different mass reduction/source discharge lines previously defined as representing different source zones with different degrees of NAPL architecture variability can be used instead to represent uncertainty in the mass reduction/source discharge relationship at a single site. As an example, consider the following scenarios displayed in figure 6-11. In the first scenario, a large increase in plume attenuation is possible with the implementation of a relatively passive technology. With the large increase in attenuation capacity, the site specific target for source discharge reduction is substantially lowered such that the time scale of plume management is small enough where it becomes the more economical remedial alternative. In the second case, the capacity for enhancing attenuation is limited, such that the current source discharge is well above the target level. Plume management, if implemented, would need to be in place for a substantial time frame. In addition, high flux measurements indicate accessible mass that can be effectively targeted and removed with flushing based remedial technologies. In this case, source zone remediation is preferable.

PAGE 188

174 Possible to significantly increase Plume Attenuation Capacity Such that Plume Management is more cost effective Limited Potential for Increasing Plume Attenuation Capacity Such that Source Zone Remediation is more cost effective site specific target reduction in source discharge Figure 6-11. Conceptual example of using figure 6-5 in the remedial design process. Conclusions In this chapter I have attempted to outline a general flux-based framework for addressing DNAPL contaminated sites. Critical to this framework is that flux based measurements function not just as remedial performance assessment tools, but as integral components of the decision to conduct source zone remediation and in the actual design of source zone remediation. The models and parameterization techniques outlined in the previous chapters are designed to fit into this framework.

PAGE 189

LIST OF REFERENCES Annable, M. D., P. S. C. Rao, K. Hatfield, W. D. Graham, A. L. Wood, and C. G. Enfield, Partitioning tracers for measuring residual NAPL: Field-scale test results, J Environ Eng-Asce, 124(6), 498-503, 1998. Augustijn, D. C. M., and P. S. C. Rao, Enhanced removal of organic contaminants by solvent flushing, Acs Sym Ser, 607, 224-236, 1995. Berglund, S., Aquifer remediation by pumping: A model for stochastic-advective transport with nonaqueous phase liquid dissolution, Water Resour Res, 33(4), 649-661, 1997. Berglund, S., and A. Fiori, Influence of transverse mixing on the breakthrough of sorbing solute in a heterogeneous aquifer, Water Resources Research, 33(3), 399-405, 1997. Bradford, S. A., K. M. Rathfelder, J. Lang, and L. M. Abriola, Entrapment and dissolution of DNAPLs in heterogeneous porous media, J Contam Hydrol, 67(1-4), 133-157, 2003. Bradford, S. A., R. A. Vendlinski, and L. M. Abriola, The entrapment and long-term dissolution of tetrachloroethylene in fractional wettability porous media, Water Resour Res, 35(10), 2955-2964, 1999. Brooks, M. C., M. D. Annable, P. S. C. Rao, K. Hatfield, J. W. Jawitz, W. R. Wise, A. L. Wood, and C. G. Enfield, Controlled release, blind tests of DNAPL characterization using partitioning tracers, J Contam Hydrol, 59(3-4), 187-210, 2002. Brusseau, M. L., N. T. Nelson, M. Oostrom, Z. H. Zhang, G. R. Johnson, and T. W. Wietsma, Influence of heterogeneity and sampling method on aqueous concentrations associated with NAPL dissolution, Environ Sci Technol, 34(17), 3657-3664, 2000. Brusseau, M. L., Z. H. Zhang, N. T. Nelson, R. B. Cain, G. R. Tick, and M. Oostrom, Dissolution of nonuniformly distributed immiscible liquid: Intermediate-scale experiments and mathematical modeling, Environ Sci Technol, 36(5), 1033-1041, 2002. 175

PAGE 190

176 Cain, R. B., G. R. Johnson, J. E. McCray, W. J. Blanford, and M. L. Brusseau, Partitioning tracer tests for evaluating remediation performance, Ground Water, 38(5), 752-761, 2000. Cho, J., Characterization of spatial NAPL distribution, mass transfer and the effect of cosolvent and surfactant residuals on estimating NAPL saturation using tracer techniques, 2001. Christ, J. A., L. D. Lemke, and L. M. Abriola, Comparison of two-dimensional and three-dimensional simulations of dense nonaqueous phase liquids (DNAPLs): Migration and entrapment in a nonuniform permeability field, Water Resour Res, 41(1), 102-117, 2005. Cirpka, O. A., E. O. Frind, and R. Helmig, Numerical simulation of biodegradation controlled by transverse mixing, Journal of Contaminant Hydrology, 40(2), 159-182, 1999. Cirpka, O. A., and P. K. Kitanidis, An advective-dispersive stream tube approach for the transfer of conservative-tracer data to reactive transport, Water Resources Research, 36(5), 1209-1220, 2000. Cvetkovic, V., G. Dagan, and H. Cheng, Contaminant transport in aquifers with spatially variable hydraulic and sorption properties, P Roy Soc Lond a Mat, 454(1976), 2173-2207, 1998. Dagan, G., Flow and transport in porous formations, Springer-Verlag, Berlin ; New York, 1989. Dekker, T. J., and L. M. Abriola, The influence of field-scale heterogeneity on the infiltration and entrapment of dense nonaqueous phase liquids in saturated formations, J Contam Hydrol, 42(2-4), 187-218, 2000a. Dekker, T. J., and L. M. Abriola, The influence of field-scale heterogeneity on the surfactant-enhanced remediation of entrapped nonaqueous phase liquids, J Contam Hydrol, 42(2-4), 219-251, 2000b. Delshad, M., G. A. Pope, and K. Sepehrnoori, A compositional simulator for modeling surfactant enhanced aquifer remediation .1. Formulation, J Contam Hydrol, 23(4), 303-327, 1996. Eberhardt, C., and P. Grathwohl, Time scales of organic contaminant dissolution from complex source zones: coal tar pools vs. blobs, J Contam Hydrol, 59(1-2), 45-66, 2002. Einarson, M. D., and D. M. Mackay, Predicting impacts of groundwater contamination, Environ Sci Technol, 35(3), 66a-73a, 2001.

PAGE 191

177 Environmental Protection Agency (EPA), The DNAPL Remediation Challenge: Is There a Case for Source Depletion?, EPA/600/R-03/143, E.P.A., 2003. Falta, R. W., Modeling sub-grid-block-scale dense nonaqueous phase liquid (DNAPL) pool dissolution using a dual-domain approach, Water Resour Res, 39(12), -, 2003. Falta, R. W., C. M. Lee, S. E. Brame, E. Roeder, J. T. Coates, C. Wright, A. L. Wood, and C. G. Enfield, Field test of high molecular weight alcohol flushing for subsurface nonaqueous phase liquid remediation, Water Resour Res, 35(7), 2095-2108, 1999. Fiori, A., Finite Peclet extensions of Dagan's solutions to transport in anisotropic heterogeneous formations, Water Resources Research, 32(1), 193-198, 1996. Gelhar, L. W., Stochastic subsurface hydrology, Prentice-Hall, Englewood Cliffs, N.J., 1993. Geller, J. T., and J. R. Hunt, Mass-Transfer from Nonaqueous Phase Organic Liquids in Water-Saturated Porous-Media, Water Resour Res, 29(4), 833-845, 1993. Gerhard, J. I., and B. H. Kueper, Capillary pressure characteristics necessary for simulating DNAPL infiltration, redistribution, and immobilization in saturated porous media, Water Resour Res, 39(8), 210-220, 2003a. Gerhard, J. I., and B. H. Kueper, Influence of constitutive model parameters on the predicted migration of DNAPL in heterogeneous porous media, Water Resour Res, 39(10), 220-230, 2003b. Gerhard, J. I., and B. H. Kueper, Relative permeability characteristics necessary for simulating DNAPL infiltration, redistribution, and immobilization in saturated porous media, Water Resour Res, 39(8), 230-245, 2003c. Ginn, T. R., Stochastic-convective transport with nonlinear reactions and mixing: finite streamtube ensemble formulation for multicomponent reaction systems with intra-streamtube dispersion, Journal of Contaminant Hydrology, 47(1), 1-28, 2001. Hatfield, K., M. Annable, J. H. Cho, P. S. C. Rao, and H. Klammler, A direct passive method for measuring water and contaminant fluxes in porous media, J Contam Hydrol, 75(3-4), 155-181, 2004. Illangasekare, T. H., J. L. Ramsey, K. H. Jensen, and M. B. Butts, Experimental-Study of Movement and Distribution of Dense Organic Contaminants in Heterogeneous Aquifers, J Contam Hydrol, 20(1-2), 1-25, 1995.

PAGE 192

178 Imhoff, P. T., S. N. Gleyzer, J. F. Mcbride, L. A. Vancho, I. Okuda, and C. T. Miller, Cosolvent Enhanced Remediation of Residual Reuse Nonaqueous Phase Liquids-Experimental Investigation, Environ Sci Technol, 29(8), 1966-1976, 1995. Imhoff, P. T., P. R. Jaffe, and G. F. Pinder, An Experimental-Study of Complete Dissolution of a Nonaqueous Phase Liquid in Saturated Porous-Media, Water Resour Res, 30(2), 307-320, 1994. Jawitz, J. W., M. D. Annable, C. J. Clark, and S. Puranik, Inline gas chromatographic tracer analysis: An alternative to conventional sampling and laboratory analysis for partitioning tracer tests, Instrum Sci Technol, 30(4), 415-426, 2002. Jawitz, J. W., M. D. Annable, G. G. Demmy, and P. S. C. Rao, Estimating nonaqueous phase liquid spatial variability using partitioning tracer higher temporal moments, Water Resour Res, 39(7), 112-115, 2003a. Jawitz, J. W., M. D. Annable, and P. S. C. Rao, Miscible fluid displacement stability in unconfined porous media: Two-dimensional flow experiments and simulations, J Contam Hydrol, 31(3-4), 211-230, 1998a. Jawitz, J. W., M. D. Annable, P. S. C. Rao, and R. D. Rhue, Field implementation of a Winsor type I surfactant/alcohol mixture for in situ solubilization of a complex LNAPL as a single phase microemulsion, Environ Sci Technol, 32(4), 523-530, 1998b. Jawitz, J. W., D. P. Dai, P. S. C. Rao, M. D. Annable, and R. D. Rhue, Rate-limited solubilization of multicomponent nonaqueous-phase liquids by flushing with cosolvents and surfactants: Modeling data from laboratory and field experiments, Environ Sci Technol, 37(9), 1983-1991, 2003b. Jawitz, J. W., R. K. Sillan, M. D. Annable, P. S. C. Rao, and K. Warner, In situ alcohol flushing of a DNAPL source zone at a dry cleaner site, Environ Sci Technol, 34(17), 3722-3729, 2000. Jensen, K. H., K. Bitsch, and P. L. Bjerg, Large-Scale Dispersion Experiments in a Sandy Aquifer in Denmark Observed Tracer Movements and Numerical-Analyses, Water Resour Res, 29(3), 673-696, 1993. Jin, M. Q., M. Delshad, V. Dwarakanath, D. C. Mckinney, G. A. Pope, K. Sepehrnoori, C. E. Tilburg, and R. E. Jackson, Partitioning Tracer Test for Detection, Estimation, and Remediation Performance Assessment of Subsurface Nonaqueous Phase Liquids, Water Resour Res, 31(5), 1201-1211, 1995. Jury, W. A., and K. Roth, Transfer Functions and Solute Movement through Soil, Berkhauser Verlag, 1990.

PAGE 193

179 Kaluarachichi, J. J., and J. C. Parker, Multiphase Flow with a Simplified Model for Oil Entrapment, Journal of Contaminant Transport, 5, 349, 1990. Kapoor, V., L. W. Gelhar, and F. MirallesWilhelm, Bimolecular second-order reactions in spatially varying flows: Segregation induced scale-dependent transformation rates, Water Resources Research, 33(4), 527-536, 1997. Kapoor, V., C. T. Jafvert, and D. A. Lyn, Experimental study of a bimolecular reaction in Poiseuille flow, Water Resources Research, 34(8), 1997-2004, 1998. Klenk, I. D., and P. Grathwohl, Transverse vertical dispersion in groundwater and the capillary fringe, J Contam Hydrol, 58(1-2), 111-128, 2002. Kueper, B. H., and J. I. Gerhard, Variability of point source infiltration rates for two-phase flow in heterogeneous porous media, Water Resour Res, 31(12), 2971-2980, 1995. Kueper, B. H., D. Redman, R. C. Starr, S. Reitsma, and M. Mah, A Field Experiment to Study the Behavior of Tetrachloroethylene Below the Water-Table Spatial-Distribution of Residual and Pooled Dnapl, Ground Water, 31(5), 756-766, 1993. Kuhn, T. S., The structure of scientific revolutions, University of Chicago Press, Chicago, 1970. Lemke, L. D., and L. M. Abriola, Predicting DNAPL entrapment and recovery: the influence of hydraulic property correlation, Stoch Env Res Risk A, 17(6), 408-418, 2003. Lemke, L. D., L. M. Abriola, and P. Goovaerts, Dense nonaqueous phase liquid (DNAPL) source zone characterization: Influence of hydraulic property correlation on predictions of DNAPL infiltration and entrapment, Water Resour Res, 40(1), 78-94, 2004a. Lemke, L. D., L. M. Abriola, and J. R. Lang, Influence of hydraulic property correlation on predicted dense nonaqueous phase liquid source zone architecture, mass recovery and contaminant flux, Water Resour Res, 40(12), 119-127, 2004b. Mackay, D. M., P. V. Roberts, and J. A. Cherry, Transport of Organic Contaminants in Groundwater, Environ Sci Technol, 19(5), 384-392, 1985. Mayer, A. S., and C. T. Miller, The influence of mass transfer characteristics and porous media heterogeneity on nonaqueous phase dissolution, Water Resour Res, 32(6), 1551-1567, 1996.

PAGE 194

180 Mayer, A. S., L. R. Zhong, and G. A. Pope, Measurement of mass-transfer rates for surfactant-enhanced solubilization of nonaqueous phase liquids, Environ Sci Technol, 33(17), 2965-2972, 1999. Meinardus, H. W., V. Dwarakanath, J. Ewing, G. J. Hirasaki, R. E. Jackson, M. Jin, J. S. Ginn, J. T. Londergan, C. A. Miller, and G. A. Pope, Performance assessment of NAPL remediation in heterogeneous alluvium, J Contam Hydrol, 54(3-4), 173-193, 2002. Miller, C. T., M. M. Poiriermcneill, and A. S. Mayer, Dissolution of Trapped Nonaqueous Phase Liquids Mass-Transfer Characteristics, Water Resour Res, 26(11), 2783-2796, 1990. Nambi, I. M., and S. E. Powers, NAPL dissolution in heterogeneous systems: an experimental investigation in a simple heterogeneous system, J Contam Hydrol, 44(2), 161-184, 2000. Nelson, N. T., and M. L. Brusseau, Field study of the partitioning tracer method for detection of dense nonaqueous phase liquid in a trichloroethene-contaminated aquifer, Environ Sci Technol, 30(9), 2859-2863, 1996. NRC, Alternatives for ground water cleanup, National Academy Press, Washington, D.C., 1994. Oostrom, M., C. Hofstee, R. C. Walker, and J. H. Dane, Movement and remediation of trichloroethylene in a saturated heterogeneous porous medium 1. Spill behavior and initial dissolution, J Contam Hydrol, 37(1-2), 159-178, 1999. Parker, J. C., R. J. Lenhard, and T. Kuppusamy, A Parametric Model for Constitutive Properties Governing Multiphase Flow in Porous-Media, Water Resour Res, 23(4), 618-624, 1987. Parker, J. C., and E. Park, Modeling field-scale dense nonaqueous phase liquid dissolution kinetics in heterogeneous aquifers, Water Resour Res, 40(5), 247-264, 2004. Pennell, K. D., A. M. Adinolfi, L. M. Abriola, and M. S. Diallo, Solubilization of dodecane, tetrachloroethylene, and 1,2-dichlorobenzene in micellar solutions of ethoxylated nonionic surfactants, Environ Sci Technol, 31(5), 1382-1389, 1997. Powers, S. E., L. M. Abriola, and W. J. Weber, An Experimental Investigation of Nonaqueous Phase Liquid Dissolution in Saturated Subsurface Systems Steady-State Mass-Transfer Rates, Water Resour Res, 28(10), 2691-2705, 1992.

PAGE 195

181 Powers, S. E., L. M. Abriola, and W. J. Weber, An Experimental Investigation of Nonaqueous Phase Liquid Dissolution in Saturated Subsurface Systems Transient Mass-Transfer Rates, Water Resour Res, 30(2), 321-332, 1994. Powers, S. E., I. M. Nambi, and G. W. Curry, Non-aqueous phase liquid dissolution in heterogeneous systems: Mechanisms and a local equilibrium modeling approach, Water Resour Res, 34(12), 3293-3302, 1998. Rao, P. S. C., M. D. Annable, R. K. Sillan, D. P. Dai, K. Hatfield, W. D. Graham, A. L. Wood, and C. G. Enfield, Field-scale evaluation of in situ cosolvent flushing for enhanced aquifer remediation, Water Resour Res, 33(12), 2673-2686, 1997. Rao, P. S. C., and J. W. Jawitz, Comment on "Steady state mass transfer from single-component dense nonaqueous phase liquids in uniform flow fields" by T.C. Sale and D.B. McWhorter, Water Resour Res, 39(3), 142-157, 2003. Rao, P. S. C., J. W. Jawitz, C. G. Enfield, R. W. Falta, M. D. Annable, and A. L. Wood, Technology integration for contaminated site remediation: Cleanup goals and performance metrics., in Ground Water Quality 2001, vol., pp. 410-412, Sheffield, U.K., 2001. Rathfelder, K., and L. M. Abriola, The influence of capillarity in numerical modeling of organic liquid redistribution in two-phase systems, Adv Water Resour, 21(2), 159-170, 1998. Russo, D., Stochastic-Analysis of Simulated Vadose Zone Solute Transport in a Vertical Cross-Section of Heterogeneous Soil During Nonsteady Water-Flow, Water Resour Res, 27(3), 267-283, 1991. Saenton, S., T. H. Illangasekare, K. Soga, and T. A. Saba, Effects of source zone heterogeneity on surfactant-enhanced NAPL dissolution and resulting remediation end-points, J Contam Hydrol, 59(1-2), 27-44, 2002. Sale, T. C., and D. B. McWhorter, Steady state mass transfer from single-component dense nonaqueous phase liquids in uniform flow fields, Water Resour Res, 37(2), 393-404, 2001. Sardin, M., D. Schweich, F. J. Leij, and M. T. Vangenuchten, Modeling the Nonequilibrium Transport of Linearly Interacting Solutes in Porous-Media a Review, Water Resour Res, 27(9), 2287-2307, 1991. Schroth, M. H., S. J. Ahearn, J. S. Selker, and J. D. Istok, Characterization of miller-similar silica sands for laboratory hydrologic studies, Soil Sci Soc Am J, 60(5), 1331-1339, 1996. Schwille, F., Dense chlorinated solvents in porous and fractured media : model experiments, Lewis Publishers, Chelsea, MI, 1988.

PAGE 196

182 Sillan, R. K., M. D. Annable, P. S. C. Rao, D. P. Dai, K. Hatfield, and W. D. Graham, Evaluation of in situ cosolvent flushing dynamics using a network of spatially distributed multilevel samplers, Water Resour Res, 34(9), 2191-2202, 1998. Soerens, T. S., D. A. Sabatini, and J. H. Harwell, Effects of flow bypassing and nonuniform NAPL distribution on the mass transfer characteristics of NAPL dissolution, Water Resour Res, 34(7), 1657-1673, 1998. Soga, K., J. W. E. Page, and T. H. Illangasekare, A review of NAPL source zone remediation efficiency and the mass flux approach, J Hazard Mater, 110(1-3), 13-27, 2004. Taylor, T. P., K. D. Pennell, L. M. Abriola, and J. H. Dane, Surfactant enhanced recovery of tetrachloroethylene from a porous medium containing low permeability lenses 1. Experimental studies, J Contam Hydrol, 48(3-4), 325-350, 2001. Tompson, A. F. B., R. Ababou, and L. W. Gelhar, Implementation of the 3-Dimensional Turning Bands Random Field Generator, Water Resour Res, 25(10), 2227-2243, 1989. Travis, C. C., and C. B. Doty, Can Contaminated Aquifers at Superfund Sites Be Remediated, Environ Sci Technol, 24(10), 1464-1466, 1990. Walker, R. C., C. Hofstee, J. H. Dane, and W. E. Hill, Surfactant enhanced removal of PCE in a nominally two-dimensional, saturated, stratified porous medium, J Contam Hydrol, 34(1-2), 17-30, 1998. Zhu, J. T., and J. F. Sykes, Simple screening models of NAPL dissolution in the subsurface, J Contam Hydrol, 72(1-4), 245-258, 2004.

PAGE 197

BIOGRAPHICAL SKETCH Adrian Fure grew up on the south shore of Lake Superior in Marquette, Michigan. After accepting a hockey scholarship to Michigan Technological University he somewhat arbitrarily decided to study environmental engineering based on his enjoyment of the outdoors. In his third year of study, the first year of actual environmental engineering coursework, he had the good fortune of taking Marty Auers Surface Water Quality Modeling class. It was in this class that he found his passion for environmental engineering, deciding he would attend graduate school. Upon graduating from Michigan Tech he accepted an Alumni Fellowship to the University of Florida to pursue his doctoral degree under the auspices of fellow Michigan Tech alumni Dr. Mike Annable and renowned flag football offensive guru Dr. Jim Jawitz. 183