. .. 25:93.}.3ia .a 5...... . a .93. a 3M... .u .r: .fi . . 2.1. 7. to... .. r. .. . i araanwmnmua . 1%.-..gmfi 2343.... d... .. an? mu} a... . 1% .at.....afiu..fifiu...§afisi . 3...... )9... m3 .3...“ ,...v.........:fi . .mnasvifl! In . . 1 . $5.133“. .» $3.... 5*. 1.... I: .52....“ .gfimrr s: n 321.13.; .5». ,. H.523... . agitdv. ghivndgxlfifi .3 mt unfit}. I. .51 2m: n. it hfhrihulzv EL ”9.5.1.3... .21.»... gapiggfl.fl I. 2.“...Q aéhfihufiavm N ‘flfisagi irafifikn .3 O... vs... at... . n| .t fell. III...” .3 5 .. 9339411.:535 .Io 5.3%....Stlal itch-I Vt. E‘s-aunt: ,5?‘ '3‘}: i: 1. saw}... 3 .fih a: . 2.59%..» 3! .1 \ 12.3.1533}... ll.‘.l..i.ih.d Eli-4“! 3:35.213: a}... 1. 1.1:... sci-35:51.... at)... if 111.25.: .2711! 13.2531... the 3:33...‘ gird. {Bui‘xil 3... I)‘:..3¢) a... :3: i’l‘!’ . it... .2. . .: a...» ..mu...... “1.3.5.“... 1...... 37:1}... . . it... 51...... a .. .. . . .1... . villaiz. 5 3.2... .388: :3: :01)! .thn C521}! 1 E3. 21.5.! a: tag-53.: a... 5 .11 151.315 .9. 53.. 51,135. xiii...» .15}... xx 1 11.5.1.3. kit-h? .1. 1.! 11; |\Y.\.lv.l : 11 . ..x . . . 3......1331Ii . . , . :. .. £3... THESIS \ 2 007/ This is to certify that the thesis entitled The Influence of Detailed Aquifer Characterization on Groundwater Flow and Transport Models at Schoolcraft,Michigan presented by Christopher John Hoard has been accepted towards fulfillment of the requirements for Master ' 5 degree in Environmental Geosciences 7).. flwflmlk. / Major professor Date /"3/'O; 0.7639 MS U is an Affirmative Action/Equal Opportunity Institution LIBRARY Michigan State University PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 6/01 c:/CIRC/DateDue.p65-p.15 THE INFLUENCE OF DETAILED AQUIFER CHARACTERIZATION ON GROUNDWATER FLOW AND TRANSPORT MODELS AT SCHOOLCRAFT, MICHIGAN By Christopher John Hoard A THESIS Submitted to Michigan State University in partial fialfillment of the requirements for the degree of MASTER OF SCIENCE Department of Geological Sciences 2002 ABSTRACT THE INFLUENCE OF DETAILED AQUIFER CHARACTERIZATION ON GROUNDWATER FLOW AND TRANSPORT MODELS AT SCHOOLCRAFT, MICHIGAN By Christopher John Hoard Aquifer characterization is a key aspect of any groundwater study, however it is difficult to determine when an aquifer has been adequately characterized. This study explored several methods of describing aquifer properties and the value of knowing those properties for developing accurate groundwater flow and transport simulations. Data collected from an aquifer bioremediation study in Schoolcrafl Michigan was used for this research. The value of incorporating various levels of aquifer characterization was assessed using groundwater flow and transport simulations. The relationship between several aquifer properties including grain-size, ground penetrating radar velocity, and hydraulic conductivity were explored using rock physics approaches. Regional aquifer porosity estimates were estimated using ground penetrating radar tomography, and the impact of these porosity values were evaluated in solute transport models. The approaches developed through this research can be applied to future studies to help design cost-effective sampling strategies for detailed aquifer characterization. ACKNOWLEDGEMENTS There are many people that I need to thank for helping me complete this thesis. I should start off by thanking the guy who made it all possible and that stuck by me through the long process. Through his support, Dave Hyndman has helped me to develop into a better scientist and a better writer. Ran Bachrach was another individual whose assistance was invaluable. He helped me learn the language of MatLab as well as aid me with aspects of the geophysical and rock physics work. Dave Wiggert was another member of my committee that helped by making significant contributions to the manuscript, and helping me see the bigger picture of the research. The research was funded by the Michigan Department of Environmental Quality (contract #Y403 86). There is a long list of others that have helped with me through graduate school. I would like to thank my family for constantly reminding me that I wasn’t quite finished. You don’t have to remind me anymore cause I am finally done. My family was very supportive through the whole process and I wouldn’t have gotten through this without them in my corner. I also want to thank Big A for being a good fi'iend through all this. Now we got more time to watch NASCAR. I also gotta thank the Mosser, Joel, Robin, Linker and Brian for the misadventures we had that kept us sane through this crazy thing called grad school. What a long strange trip its been... iii TABLE OF CONTENTS LIST OF TABLES LIST OF FIGURES 1. EVALUATION OF GROUNDWATER FLOW AND TRANSPORT MODELS INTRODUCTION SITE BACKGROUND REMEDIATION SYSTEM OPTIMIZATION FIELD TRACER TEST METHODS GROUNDWATER FLOW AND TRANSPORT MODELING RESULTS II. PHYSICAL CHARACTERIZATION INTRODUCTION GRAIN-SIZE ANALYSIS RADAR TOMOGRAPHY III. CONCLUSIONS BIBLIOGRAPHY iv V vi 1 12 22 36 36 36 48 61 LIST OF TABLES TABLE 1 Description of data used for each flow and transport simulation. . . 1 5 TABLE 2 Summary of variogram analysis for various well scenarios ................. 19 FIGURE 1 FIGURE 2 FIGURE 3 FIGURE 4 FIGURE 5 FIGURE 6 FIGURE 7 FIGURE 8 FIGURE 9 FIGURE 10 FIGURE 11 FIGURE 12 FIGURE 13 FIGURE 14 FIGURE 15 LIST OF FIGURES Images in this thesis are presented in color Location of plume and study area. (Taken from Hyndman et al., 2000). . ..4 Well layout for study area ......................................................... 7 Example of 50% quantile arrival times for well 9-75 ........................ 13 Dimensions of model grid a) plan view of model grid b) cross-sectional view of model grid ................................................................ 14 Comparison of histograms for normal and log-transformed hydraulic conductivity data ................................................................... 1 7 Example variogram for the 7 well hydraulic conductivity case ............. l 8 Relationship between hydraulic conductivity and depth in the aquifer illustrating separate zones used for hydraulic conductivity interpolations ........................................................................................ 20 Spatial distribution of log hydraulic conductivity field, including the location of kriging zone boundaries ............................................ 21 Concentration history of observation wells 9 - 13 for 1 , 3 , 7 , and 7 well deviated and scenarios ............................................................. 23 Comparison of transport simulation residual plots a) 10% quantile comparison b) sum of quantile comparison c) total mass comparison. . ...26 Estimation variance and kriged log hydraulic conductivity outcome for two different 3 well scenarios .................................................... 29 Estimation variance and kriged log hydraulic conductivity outcome for the 7 well case .......................................................................... 30 Example of 50% quantile evaluation for monitoring wells 9 and 11 ...... 32 Comparison of 50% quantile to total concentration deviation for monitoring well 11 ................................................................ 35 Location of grain size samples along delivery well gallery .................. 37 vi FIGURE 16 FIGURE 17 FIGURE 18 FIGURE 19 FIGURE 20 FIGURE 21 FIGURE 22 FIGURE 23 FIGURE 24 FIGURE 25 FIGURE 26 FIGURE 27 Plot of relationship between measured hydraulic conductivity, depth and efl‘ective grain size ................................................................. 39 Plot of relationship between measured porosity, depth and effective grain size ................................................................................... 40 Plot of relationship between measured hydraulic conductivity, depth and effective grain size ................................................................. 41 Kozeny-Carman estimated permeability vs. measured permeability. . .....43 Kozeny-Carman estimated permeability vs. measured permeability using the d10 efl‘ective grain size ........................................................ 45 Kozeny-Carman estimated permeability vs. measured permeability using the d20 effective grain size ....................................................... 46 Kozeny-Carman estimated permeability vs. measured permeability using the d30 effective grain size ...................................................... 47 Map of GPR tomography planes using a normalized velocity. Labels indicate well locations ............................................................. 51 Relationship between GPR velocity, porosity and hydraulic conductivity ......................................................................... 52 Plot of Hashin-Shtrikman bounds and results of various mixing theory methods to calculate water content ............................................... 56 Radar derived porosity for delivery well gallery (DI-D15) and the comparison of those values to measured values ............................... 57 Simulated and observed tracer concentrations for radar derived porosity, measured porosity and constant porosity ........................................ 59 vii CHAPTER I EVALUATION OF GROUNDWATER FLOW AND TRANSPORT MODELS INTRODUCTION Aquifer characterization is an essential component of groundwater modeling and remediation studies. Poor aquifer characterization can lead to inadequate models and poor remediation design decisions. The remediation design may be either over- engineered, which wastes money on excessive treatment costs, or the design may not effectively remediate the aquifer due to inadequate characterization. A key question becomes: how do we know if a site has been adequately characterized? The answer is not simple because it can be argued that large amounts of data are needed to adequately represent the heterogeneity of an aquifer. The purpose of this paper is to examine different methods of evaluating the value of hydraulic conductivity data on groundwater flow and transport simulations. Including more hydraulic conductivity data in simulations should improve the predictive power of the simulations, but at some point the improvement in model performance should plateau with the addition of further data. Recent tracer studies have shown the value of identifying heterogeneity in aquifer systems. Vereecken and others (2000) performed a large-scale tracer test on heterogeneous fluvial deposits at the Krauthausen site, in Germany. Investigations revealed that estimated dispersivities of the bromide tracer test, based on the theory by Gelhar and Axness (1983), did not match the observed dispersivities. This was most likely from heterogeneities on the site that caused loss of tracer mass and uncertainty in hydraulic conductivity estimates. The few studies that have had detailed aquifer characterization have shown that aquifer heterogeneity plays a major role in solute transport. For example, studies carried out at Base Borden, Ontario (Mackay et a1. 1986; Freyberg, 1986; Sudicky, 1986), Cape Cod (LeBlanc et al. 1991; Garabedian et al.1991; Hess et al. 1992) and the MacroDispersion Experiment Site, near Columbus, Mississippi, (Boggs et a1. 1992; Rehfeldt at al.1992) all illustrate the effect of heterogeneity on flow and transport. Most remediation design studies do not have the luxury of working with highly characterized aquifers. Yet few, if any, studies have explored the value of hydraulic conductivity data for transport simulations using field data. Several papers have addressed the issue of data worth and sampling design. Many researchers have used Bayesian statistics to evaluate the worth of hydrogeologic data. Most of these studies are cost-risk evaluations that tie economic decision making with hydrogeologic uncertainty, thus, addressing the value of data in a monetary sense. The basic premise in most Bayesian schemes is that if the cost of a sample does not significantly lower the cost of treatment or remediation design, then the sample is not worth taking. This methodology has been used to examine hydrogeologic data worth in design of waste management facilities and pump-and-treat remediation plans (Freeze et a1. 1992; James and Freeze 1993; James and Gorelick 1994). Others have examined the data worth issue in designing optimal data sampling strategies to improve model predictions. Loaiciga (1989) used a mixed integer programming method to determine an optimal sampling network for a landfill in Ohio. In that case, model predictions were improved using parameter estimation as well as the increasing sample density following the analysis. Wagner (1995) performed a similar study using two types of optimization algorithms to evaluate the cost of a sample versus the information obtained from a sample. In addition, he examined the cost of different hydrogeologic data types and the effect of each those data types had at reducing model prediction uncertainty. Wagner (1999) built on his previous work by developing a four- part methodology that identifies the most cost effective sampling strategy for groundwater management studies. The majority of the previous data worth analysis have been performed on synthetic data sets. This paper presents a case in which we evaluate the impact of varying amounts of hydrogeologic data on groundwater flow and transport model performance. The influence of adding hydrogeologic data is demonstrated through groundwater flow and transport simulations. A tracer test is simulated using a subset of the available hydraulic conductivity data. The same simulation is performed again, only using a larger subset of data. Data from well pairs are added incrementally until all the available data are incorporated in the estimated hydraulic conductivity field. The results of the solute transport models are then compared against each other, as well as the measured data from the tracer test. Several different methods were used to evaluate solute transport performance. Simulated transport concentrations were evaluated against the measured values using median quantile arrivals, sum of 10-90% quantiles, and the difference between total mass of tracer simulated and measured (concentration deviation). These analyses can be used to evaluate the value of detailed aquifer characterization as well as explore the best method to compare simulation results. SITE BACKGROUND The study area is in the East side of the village of Schoolcrafi, in Kalamazoo County, Michigan (Figure 1). Schoolcraft Plume A is a carbon tetrachloride (CT) plume affecting a region of the unconfined aquifer in Schoolcraft. The source of CT .5. Figure 1: Location of plume and study area. (Taken from Hyndman et al., 2000) contamination is speculated to be an abandoned grain elevator near the center of the village where the contaminant was used as a pest filmigant (Mayotte, 1996). Carbon tetrachloride is a suspected human carcinogen and has been shown to cause acute liver toxicity in laboratory animals (Sittig, 1985). The extent of the plume is approximately 1.2 km long, 90 m wide and lies between 8 and 26 meters below ground surface (bgs) (Hyndman et al., 2000). The aquifer below the village of Schoolcraft is composed of glaciofluvial sediments left behind by the retreating Laurentide ice sheet. The Lake Michigan lobe of the Laurentide ice sheet retreated across Schoolcrafi to the location of the Kalamazoo moraine (Monaghan et al., 1986). At this point, it is theorized that melt water ponded against the combination of ice and moraine until it breached the moraine and spilled out across the area (Steinman, 1994). The material deposited by this massive flood was termed the Prairie Ronde fan by Steinman (1994). The contaminated aquifer in the study area is composed of the outwash material from this flood. Directly underlying the coarse aquifer material is a layer of dense clay. This clay layer kept the contaminant from migrating into deeper parts of the aquifer. The origin and extent of the clay layer is uncertain, but at this remediation site all wells drilled encountered clay at approximately 90 ft bgs. REMEDIATION STUDY This study is part of a larger study aimed at developing new, widely applicable strategies to bioremediate contaminated groundwater. The Schoolcrafi site provided an ideal place to test the new bioremediation design because a microbe was identified that could degrade CT into nonharmful products, there were adequate levels of nitrate available as an electron acceptor, and the aquifer material was relatively homogenous with only two orders of magnitude range in hydraulic conductivity. The design approach involved creating an active bioremediation zone in the aquifer by injecting microbes, substrate and other necessary constituents to enhance microbial activity. This zone of active degradation was termed a biocurtain (Hyndman et al., 2000). The biocurtain was established by injecting microbes, substrates and other constituents across a transect of wells that are perpendicular to natural groundwater flow gradient (Figure 2). Natural groundwater gradient then carries contaminated water through the bioreactive zone where contaminants are degraded. Successful degradation of the contaminant relied on efficient delivery of microbes, substrates and other constituents across the transect. Failure to completely deliver material across the transect would result in gaps in the biocurtain, and contaminated water could pass through gaps. Therefore, it was important to identify heterogeneity that might effect delivery across the transect. SYSTEM OPTIMIZATION Groundwater flow and transport models linked to an optimization function were used to design a well layout and pumping scheme for the remediation system (Hyndman et al., 2000). The goal of the optimization was to determine flow and well field properties (e.g. number of wells, well spacing, duration and rates of pumping, and duration of flow reversal) that achieved complete delivery of microbes and necessary nutrients while minimizing design cost. On the basis of these results, a one meter well spacing for the delivery wells (Dl-DlS Figure 2) and a weekly 6 hr pumping strategy was chosen for the projected 25 yrs of system operation. Direction of N groundwater 21 “0W 4 8 0:5 a W E O D 014 .P11 013 1c} 18 S I 012 q- 7 011- 0 P9 31 D 010 ‘30 I . 09. 13 20 24 2 I . I: . 3 D: 10:: 15 P7 '3 DB 0 P6 05 - 102 19 . ' 28 D4 29 D 2:? 25 g. 5 [)3 ' P10 '3 " 9 14 D‘- D (,1. 27 a P12 0 18 0 Well Symbols I Injection/Extraction delivery wells 9 Nested 2.54mi multilevel wells with 0.304m screens at depths of 17 9.1m.12.2m.15.2m.1B.3m.21.3m,24.4m.and 27.4m :1 A Nested 2.54m multilevel wells with 0 .809m screens centered 22 at depths of 24.4m.25 .9m.and 27.4m 10.16cm wells with multilevel diffusion cell sampling Nested 2.54mi multilevel wells with 0.809m screens at depths of 10.7m.13.7m.1B.8m,19.8m.and 22.9m 2" Observation Wells 0 4 8 E Figure 2: Well layout for study area. 12 Meters The pumping strategy developed included a recirculation loop in which the even numbered delivery wells would inject the solutes, while the odd numbered wells would extract formation water, in effect, pulling the injected material across the delivery well gallery. The initial pumping period was for 5 hrs at a flow rate of 9.1 m3/hr. This was followed by 1 hr of reversal flow in which the odd numbered wells would begin injecting and the even numbered wells would extract water. The following week this schedule would be reversed, so the odd number wells would inject for the first 5 hrs and the even would inject for the 1 hr reversal. This alternating scheme was designed to evenly distribute injected material across the delivery wells. FIELD TRACER TEST Following installation of the field system, a tracer test was performed to evaluate the designed pumping strategy and natural gradient flow through the observation gallery. For this study two tracers, fluoroscein and bromide, were chosen because they are easily detected, do not affect microbial growth, do not affect contaminant degradation and are relatively conservative (Hyndman et al., 2000). The proposed strategy was to inject tracer into the even-numbered delivery wells at a combined pumping rate of 9. l m3/hr (40 gpm). Bromide was injected at an average concentration of 16 ppm while fluoroscein was injected at approximately 130 ppb. At the same time, water was extracted through the odd-numbered wells at a total rate of 9. lm3/hr. After 5 hrs, the flow was reversed for 1 hr and then natural gradient was allowed to carry the tracer through the sampling grid for the remainder of the week. During the flow reversal stage, extracted water contained tracer fiom the previous 5 hr injection. Additional bromide and fluoroscein was added in the recirculation loop, bringing the average concentration of the bromide to 22 ppm and fluoroscein to 215 ppb for the recirculated water for the final hour of injection. A period of 21 days of natural gradient flow followed the tracer injection to allow the tracer slug to move through the monitoring wells without added disturbance due to pumping. Additional bromide was injected into the aquifer for all remaining pumping events to tag all injected water used for pH adjustment. Results of the tracer test revealed the pumping strategy was eflective for distributing material across the aquifer and the average of the measurements was quite similar to the predicted concentrations (Hyndman et al., 2000). A complete description of the tracer test is found in Hyndman et a1. (2000). The initial attempts at modeling the tracer test accurately predicted averaged observed tracer concentrations, but had significant difficulties in matching individual measured concentration histories (Hyndman et aL, 2000). At the time of the tracer test, only a small subset of the hydraulic conductivity data were available. The discrepancy between individual simulated and observed concentration histories was most likely caused by uncertainties in well location or significant heterogeneity along the transect that were not identified with the available data. As more data became available, transport simulations improved providing a better match between the simulated and observed concentration histories. The effect of incorporating additional hydraulic conductivity data in transport models had not been quantified, which is the motivation for this research project. METHODS Seven sediment cores were collected in the region of the delivery well gallery (Figure 2) using the Waterloo cohesionless continuous sand sampler (Dybas et al. 1998). These borings were taken in the locations of the even numbered delivery wells between the depths of 9.1 and 24.4 m. The sampler operates by advancing 1.5 m long, 5 cm diameter plastic core tube ahead of the auger bit. A plunger-type device was inserted in each tube to maintain vacuum as the sampler advanced into the aquifer, thus improving core sample recovery. Incomplete core recovery often resulted in sediment shifting in the sample tube, resulting in uncertainty with the exact depth of the sediment sample. The cores collected from the aquifer were used to estimate the hydraulic conductivity of the aquifer. Core samples were cut into 20 to 25 cm sections and 220 of these sections were tested using constant-head permeameters. Sediment was removed from selected core sections and placed into sample bags. Following oven drying at 60° C for 24 hours, the samples were allowed to air dry at room temperature for 48 hours. Approximately 300 grams of sediment were weighed and placed into a permeameter experiment cell on top of a porous conductivity stone. The sample was compacted by placing a plug on top of the sample and lmmmering that with a drop weight. The drop hammer was 1.74 kg and dropped fiom a height of 7.54 cm 10 times. A second porous conductivity stone was placed on top of the packed sediment to seal the experimental cell. To prevent water leaks fi'om the permeameter, PT PE pipe thread and seal tape was wrapped around the edges of the stone. This created a tight seal between the experiment cell and the stone. The sample was then placed into the saturation apparatus. Approximately five pore volumes of carbon dioxide were flushed through the sample to remove air. Then the sample was saturated from the bottom with deaired water from the Schoolcraft site. The hydraulic conductivity was measured by passing the deaired Schoolcraft water through the sample at a constant head gradient. Results of the 10 permeameter analysis yielded hydraulic conductivity values ranging fi'om 0.0011 cm/s to 0.1056 cm/s with a mean value of 0.0028 crn/s. Given the close well spacing (1 m separation) in the delivery well gallery, even slight deviation of the well bores from vertical could significantly affect delivery and transport. For example, two wells deviated away from each other would not act the same as two perfectly straight wells. The deviated wells require an increased volume of pumping to achieve the same level of breakthrough as the straight wells, because the distance between the wells is larger than if they were perfectly straight. Deviation of the delivery well bores from vertical was measured using a borehole colloidoscope, which used a compass to determine the orientation of the well and a visual dipmeter to measure the angle from vertical. The tool was lowered in 1.5 m increments, at each location the measured direction and angle were recorded. Between colloidoscope measurements, the amount of well bore deviation was interpolated using a linear interpolation. Some of the wells were estimated to deviate as much as 32 cm fiom a vertical line. However, due to the inexact nature of these measurements, there is still uncertainty in the accuracy of the deviation measurement. There was an additional issue with the flame of reference used when applying the deviation measurements to model coordinates. The possibility of the deviation measurements being 180° off was tested by comparing transport simulations for the two different cases. The deviations that best matched the measured data were used for the fixture modeling. The performance of solute transport simulations were evaluated using a variety of methods including concentration quantile arrival times. The quantile arrival times were calculated using a Matlab script that estimated the mass under the tracer concentration 11 curve. Based on the integrated mss under the concentration curve, the time required for a certain percentage (e.g. 50%) of tracer mass to reach an observation point was calculated. For example, the 50% quantile is the time it takes for 50% of the tracer mass to be observed at a location. Figure 3 shows the calculated 50% quantile arrival times (red and blue X) for well 9-75. MODELING In order to understand the effects of heterogeneity on the groundwater flow system, a series of groundwater flow and transport models were developed. Three- dirnensional flow and transport models were constructed using the Groundwater Modeling System (GMS), a pre and post processing interface for MODF LOW (McDonald and Harbaugh, 1996) and RT3D (Clement, 1997). A control volume with the dimensions 102 m by 57 m by 27.44 m was used as the model domain (Figure 4). Specified head boundaries were used at two ends of the model to simulate the natural head gradient (0.0011) measured at the site (Hyndman et al., 2000). The bottom and other two sides of the model domain were considered no-flow boundaries. The model domain was discretized into a network of 132 columns 86 rows and 23 layers. Fine cell spacing, about 20 cm by 20 cm, was used in the 8 m by 20 m region around the delivery well gallery to accurately model the flow dynamics in that zone. The cell spacing increases away fi'om the delivery well gallery at just less than 1.5 times the previous cell’s length in order to maintain stability of the model solution. The model was vertically discretized into 22 equally spaced 1 m layers with a 5.44 m thick top layer, which was mostly composed of the unsaturated zone. The 1 m 12 .23 :03 com 85.: FEB 035.5 $3 mo 29:an ”m 253m 95h 9N up 2. fir Nr 2. w w v N a ‘ ‘ I i- . i ’ I I ‘ ‘ ¢ I l m t l GP . . 3 9 m .. ON o L m o oi u _. I ma I l on . $.58 :8 8.225 x . an 05550 $9” U05m~0§ X bun-NC. gun-WNQ: . .88... Sign.» .1 n F I n n n n n n D! nhua l3 2a a Bee me 263 fiuowoormmob 3 Saw .ovofimo 303 Sam 3 Saw E88 mo mac—€085 .v Raw—m E VKN e E N3 firm As 14 spacing was chosen to provide a reasonable description of the hydraulic conductivity distribution for the flow and transport simulations. The goal of the flow and transport modeling was to examine the effect of different levels of aquifer characterization on the accuracy of flow and transport simulations. To explore this, a series of simulations were developed that incrementally add more hydraulic conductivity data to the simulations. Adding more data to the model should improve the model predictions as well as reduce uncertainty of estimated hydraulic conductivity values. Six separate simulations were formulated to represent varying stages of aquifer characterization. The first modeled scenario is based on the hydraulic conductivity data obtained from one well core in the center of the delivery well gallery (Figure 2). The second scenario uses the hydraulic conductivity data of two well cores, at the ends of the delivery well gallery. Scenario 3 utilizes the hydraulic conductivity data from three different well borings, the two end wells and the center. In the fourth scenario, hydraulic conductivity from 5 well cores is used in the flow model, using the 2 end wells and 3 wells in the center. The fifth scenario modeled incorporates the hydraulic conductivity Table 1: Description of data used for each flow and transport simulation. Scenario Well Data Used Interpolation Scheme 1 D8 Layered Average 2 D2, D14 Layered Average 3 D2, D8, D14 Zonal Kriging 4 D2, D4, D8, D12, D14 Zonal Kriging 5 D2, D4, D6, D8, D10, D12, D14 Zonal Kriging 6 D2, D4, D6, D8, D10, D12, D14 and deviation Zonal Kriging data of all 7 well cores. The sixth scenario uses hydraulic conductivity data fi'om the 7 well cores, but it also includes the influence of borehole deviation on the pumping wells. 15 Hydraulic conductivity values were assigned to model grids through a layered averaging of values or by interpolation values across the grid by kriging depending on the amount of assumed data available. The layered average technique was used for the first and second model scenarios, because inadequate data were available to fit a variogram to those sparse data sets. Hydraulic conductivity estimates within each layer were averaged and assigned to the appropriate layer. The hydraulic conductivity values for scenarios 3 through 6 were interpolated across the grid using a zonal kriging approach. Prior to the variogram analysis, histograms of the data sets were generated to examine the distribution of the data set. For all cases, the histograms were significantly skewed (Figure 5) so the analysis was done on log-transformed hydraulic conductivity data, which resulted in a more normal distribution. Since kriging, assumes a normal distribution for the data the log transformation was chosen for all kriging cases. A normal scores transform (Deutsch and Journel, 1998) was also performed on the data for one case, to convert the skewed distribution to a normal distribution. Kriging the normal scores transform data did not produce significantly different results from the log transformed data. The log transformation was not as computationally intensive, therefore it was used instead of the normal scores transform for all remaining results. Prior to kriging, experimental variograms were produced for the data set of interest. An omnidirectional experimental variogram, as well as variograms for the horizontal and vertical directions were generated. An exponential model was used to fit the experimental variograms for all cases. Anisotropy was applied to the model variogram to fit the vertical component of the data. Figure 6 provides an example 16 .83 £303.80 0:363 3&0805502 Ea 3E5: 000 mEEwSmE mo acetamfioo um Emmi 8.0. 8.7 02.. one. 00.? 01.? 36. 8.7 054.. 00.7 00.? SN. 94.. and. as. =a as Shoes 3300:0000 0:080»: vocaommg ms 3.0 9.0 8.0 8.0 8.0 5.0 8.0 8.0 8.0 3.0 8.0 8.0 «0.0 5.0 3% =a 80 0805203888 2323: 17 .030 33,503.80 0:353 :03 A. 05 .8.“ 05.00:? can—«£6 8:03 18 variogram fi'om the 7 well case. A summary of the model variogram parameters used for each scenario is found in Table 2. It should be noted that the same hydraulic conductivity field is used for both the 7 well case that incorporates borehole deviation and the 7 well case that does not. Table 2: Summary of Variogram Analysis for Various Well Scenarios 3 well 5 well 7 well Nugget 0.0065 0.0063 0.01 1 Horizontal Correlation Length (m) 7.03 10.12 13.8 Vertical Correlation Length (m) 1.69 4.76 5.77 Variance 0.044 0.057 0.1 1 Three distinct hydraulic conductivity zones are apparent in Figures 7 and 8; the top zone (0-15.5 m bgs) has relatively low hydraulic conductivity, the middle zone (15.5 — 22.3 m bgs) has slightly higher hydraulic conductivity, and the deep zone (22.3 — 27.4 m bgs) has high hydraulic conductivity. Because of this, the zones were separately kriged in an effort to maintain stationarity. The data collected between 22.3 and 27.4 m bgs were kriged to the model grid using the variogram model obtained from the aforementioned procedure. The same was done with the middle and upper hydraulic conductivity zones as well. The three, kriged hydraulic conductivity zones were transformed back to the correct values by taking the inverse log of the estimates. Next, the three zones were combined to fit the entire model grid using a Matlab script and imported into the groundwater flow model. 19 down—0935 33:03:00 BEE—0.3 do.“ 09.: 800.». 00838 @0882: 00:33 05 E 5000 was 33320000 0:580»: 8250: 0280023“ ”n Semi o a o o 0 o coco .0 00 .mi FN o o.» o o o o. H: MN 0 0 O” 0 O M 0 .. 0&5 emoanSBE ouo‘. o p o o o». 0 mi MN 1 0 > #9 N r o—. w w v N ‘ 3 .0 T C «m L m m a a I w I a g a a 3 a. M,“ 5% 7 c bu. A.” ‘ C 0 s s x .0 ( .1 Q ,m.‘ .. E x l ,. fl 3 e I i I N O .0,” C 3 :9 a .. m p a a 0 0 fi .2. Q t. G r, a .00 a a 0 a a c... w; a . 7 v . ,, , j 3 a . W M m e c a 333:8 0.892 no. .6 8.5250 3:30 am. :2 “16°C! 21 The pumping stress periods of the flow model mimic those of the tracer test mentioned previously. A flow of 9. l m3/min was extracted fi'om odd numbered wells and distributed across the even numbered injection wells. To correctly model the flow to wells, a conductivity-weighted average was used to calculate the flux from each layer entering or exiting the well. The conductivity-weighted average was performed independently for each model scenario. The flux from the aquifer is naturally distributed according to hydraulic conductivity so the flux-weighted average accounts for this behavior. The flux-weighted average was needed because of the large difference in hydraulic conductivity between the region at the top of the well screen (2 9.1 m) and the region near the bottom of the well screen ("~" 24.4 m). Concentrations of the injected bromide tracer were measured throughout the duration of the 6-hour pump event. Since the concentration was not constant throughout the tracer test, the transport model stress periods needed to be adjusted to match the measured concentration history for the injection wells. The 21 day period of no tracer injection was followed by injection of bromide at concentrations of approximately 130 ppm for each weekly feeding event. This mass needed to be accounted for in the simulations since it mixed with the original tracer injection in the lower conductivity shallow observation points. As a result 52 stress periods were modeled starting with the tracer test and including the initial inoculation of microbes and the subsequent weekly feeding events. This concentration data was used as the inputs for the RT3D simulations. RESULTS Examining the concentration history (Figure 9) provides a qualitative way to interpret the results of transport simulations. All concentrations presented in this figure 22 Well Number 2 S o— dor—«:08 can got :03 5 v5 . b . m . _ no.“ 2 - a m=03 8.333.903 @803 eons—0:00:00 5 arm AgaDVoEF oneroonoweoooovonoomovowo l I l J o x x x 0. m6 In x 8 O 1- O x x» lv 01- x a“ md x O!- Q Q Q Q QQ Q a. x . x x x x I ll llF x x o k e». 00 x x x x x . I In: in? 3538?: II . o =an II in II xx 90 005305. x x m. 3.: I x w on mm mm 80 589 23 are normalized as C/Co, The difference in hydraulic conductivity with depth is reflected in the arrival times of the concentration peaks, with tracer arriving much faster at the deep observation points (75 it) than at the shallow observation points (35 and 45 it). This is further evidenced by the fact that one distinct peak is exhibited in the deeper zones while the shallower observed depths show an initial peak merged with a second larger peak of tracer. The second rise seen, mainly in depths 35 and 45, is attributed to the injection of additional bromide used to tag microbial feeding events that followed inoculation. Since, the shallow zones had a higher residence time, associated with low hydraulic conductivity, the initial peak did not have enough time to fully pass the observation point before the additional tracer was added to the system. This blending of tracer peaks is not seen in the deeper observation points ofthe aquifer because the residence times in that zone was much lower so the initial tracer peak traveled past the observation point before the second injection of tracer reached that observation point. Based on the concentration histories, it is clear that there is a significant difference between the cases using one and two well cores of hydraulic conductivity data and the other cases shown. The 3, 5 and 7 well scenarios do a significantly better job at representing the measured concentration histories than the one well case. For all cases except depth 35 ft, the tracer simulated by the one well case arrives later and is spread across a wider time range than all the other cases. Well 10 matched the measured values closer tlmn the other wells for the one well case, because the hydraulic conductivity for this case came from D8, which lies directly up gradient of well 10 (Figure 2). The 3 and 7 well case transport simulations are nearly the same for all observation locations, although incorporating the effect of borehole deviation does significantly affect 24 the majority of the observation locations. The exception to this is well 12, which indicates almost no influence of borehole deviation on the transport simulations. This is because the delivery wells directly upgradient of monitoring well 12 are not significantly deviated. For the shallow aquifer zone, well 11 shows the most significant difference between deviated and non-deviated simulations. The deviated simulation matches the measured tracer much better because the upgradient wells are deviated towards well 11. Simulated tracer history including deviation for well 9 does a much better job of matching the measured history. In that case the delivery wells upgradient of well 9 are deviated away fiom well 9 which slows the tracer peak and better matches the measured concentration history. Although it is apparent fiom this series of plots that adding more information to simulations improves model predictions, a quantitative approach was needed to compare the results, which is why quantile and concentration deviations were used to evaluate the model simulations. Three different methods were used to assess the value of adding information to groundwater flow and transport models: 50% concentration quantile arrival times, sum of 10 - 90% concentration quantile arrival times, and total concentration deviations. Each method produced similar results as evidenced by Figure 10, which illustrates the change in residual as a fimction of the number of well cores used in the model. The results of the quantitative methods chosen to explore the effects of heterogeneity are found in Figure 10. The most obvious feature in Figure 10 is the sharp drop in residual where the model uses 3 well cores as opposed to two well cores of hydraulic conductivity data. This is probably the effect of being able to interpolate the hydraulic conductivity field using kriging instead of using a layered averaging approach 25 demteqaoo 038 .88 Au eaten—Eco 0553: mo 83 3 :85qu0 26:36 $3 Q 303 3:28.. cote—58mm ton—05b mo :BEQEoU 5— 23E 1.— sois .352 a o u 0 n c n u . 1 4 4 l4 8" 7 1 /// 82 I/l. // l l - L ./ gr / m 1 saw 0 / m /_ 1 8am ' / m / s / 1 8:“ T 1... 1 82 a //i l/J » 80a 8230 8952 82.3 .852 a o a o n q n n p o m 4 n w — J . 8 1 08 n .1, I’ll. ..l ll 0 1|; lllllllth'lllllvllllll Ill-['0 L 8 Ill! llllllll J’ 180 2,, _ 2., 1 8w , M08. .., , . .. ..1 .— .r 1 ONF / . ..M 1,. ..oom— .... I. M ,, /, M a 1... Ova m 1. 8: m .... w H M a w 2., 8w _ /. f 1. 89 , t M If 1M 8F 1 _. ,2 /1 1. g— y MOON _ 1, M . 11.. . 8% L ova 26 to assigning hydraulic conductivity. In this case, 3 wells were within one correlation length which may adequately represent the aquifer due to the significant lateral correlation lengths at the site. The performance of the one and two well scenarios could vary significantly depending on which well locations were chosen to represent the hydraulic conductivity of the aquifer. However, it is unlikely that any averaged result, from a small subset of available data would adequately describe the aquifer. For a layered average case to be appropriate, each of the concentration histories for sample wells 9 10 and 11 (Figure 2) 1 m down gradient from delivery well should be identical. This is not reflected in the measured data at this site (Figure 9). Differences in the concentration histories for those wells are likely due to borehole deviation and lateral variability in aquifer properties. Aquifer heterogeneity is likely the most significant factor causing the differences in measured data for the three down gradient monitoring wells. For this site, 3 well cores of hydraulic conductivity data were enough to capture the larger scale heterogeneities and spatial distribution of hydraulic conductivity. Even with the addition of further hydraulic conductivity data (5 and 7 well cases), there was not much influence on the transport simulations. This is evidenced by the plateau of residuals from the 3 to the 7 well cases (Figure 10). However, variogram analysis of the 5 and 7 well cases produced significantly larger correlation length estimates in both the horizontal and vertical directions (Table 2). In the horizontal direction the estimated correlation length is doubled from the 3 well to the 7 well case. The estimation variance of the kriged conductivity field clearly illustrates the i value of incorporating additional data as well as the efl‘ect of spatial location of the data 27 points. Figure 11 shows the kriged result and the estimation variance for two different 3 well scenarios. Figure 11a demonstrates the effect of incorporating data from wells D2, D8 and D14 (Figure 2), which is the same 3 well scenario used for the analysis of transport simulations in the rest of the paper. Figure 11b uses data from wells D4, D8, D12. Based on the estimation error of the two cases, Figure 11b has a lower estimation error between the wells because well spacing is closer (4 m) for 11b as opposed to 6 m for 11a. The uncertainty in the kriging estimates is reduced by having closely spaced samples. Similarly, increasing the amount of data used in the hydraulic conductivity interpolation significantly reduces the uncertainty of the interpolation outcome. Comparing Figure 11a to Figure 12 reveals that even though the hydraulic conductivity outcomes are similar, the estimation variance is much different for the two cases. The 7 well case (Figure 12) lms a much lower estimation error associated with the hydraulic conductivity estimates than the 3 well case (Figure 11a). This is expected because the data points are only 2 m apart for the 7 well case instead of6 m apart in the 3 well case (Figure 11a). Despite the reduction in estimation error using 7 wells, the transport simulation results for the 7 well case did not considerably improve over the 3 well case (Figure 10). This is due to the fact that for the 3 well case, the horizontal correlation length was long (7 m) and because wells selected for the 3 well case provided a reasonable representation of the aquifer heterogeneity. Therefore, incorporating additioml data did not significantly affect the results of kriging, as evidenced by similar hydraulic conductivity fields produced by the 3 and 7 well case in Figure 11b and 12. 28 .8588 :03 m ~85&% 93 .50 08038 33330.80 0:35.»: me— 0005 05 35:3, eeuefiumm M: 2:05 3:53. we; 0.0 _.0. 0 0 _.0. 0. «N00 0N0.0 v00. 0 3.0 . 90.0 35.53 35.53 coqudam 028.5; cogs—Em 29 Estimation variance k Variance 0.06 0.053571 0.047142 . 0.04071 7 0.03428 '. ‘. 0.02785 ‘1 0.02142 0.015 A L V 7 12m Log K(cm/s) Log hydraulic conductivity field Figure 12: Estimation variance and kriged log hydraulic conductivity outcome for the 7 well case. 30 The addition of well borehole deviation also causes a reduction in both quantile and concentration deviation residuals (Figure 10). The composite residual plots make it appear that the incorporation of borehole deviation greatly improves the simulation results. However, by adding efl‘ects of borehole deviation, residuals at some sampling points increase, indicating poor model prediction. At the same time, some residuals at other sampling points decrease, indicating improved model prediction. For example, the 50% quantile arrivals for sampling points 9-35, 9-55 and 9-65 (Figure 13) all get worse with the addition of borehole deviation; however, for well 13 all sampling points improved (Figure 11). In general the, the sampling points in the shallower depths performed worse or stayed the same, while the deeper sampling points improved by including borehole deviation data. The effect of borehole deviation was expected to be greater at depth, which is consistent with the model results. The most likely reason for the poor performance of the shallow sampling points is the lack of measured hydraulic conductivity data available to compare to the simulated values. The shallow depths (35 and 45 it) were sparsely sampled (Figure 9). For calculation of quantiles and concentration mass, the measured concentration history is interpolated using a trapezoidal area calculation. As a result, accuracy of the mass under the curve heavily relies on a high tracer sampling frequency for accurate results. The sample frequency was very low for the shallow depths, which introduces significant uncertainties in both the quantile analysis and the concentration deviation analysis. In addition, uncertainty associated with the method of repacked testing of hydraulic conductivity estimates may also be causing the poor behavior of the slmllow observation points. The repacked method results in a homogenized sample that does not 31 .2 05 a £03 mew—SEQ: com 5.53.56 055:: .x.% .«e 20:83 22 2:03 3.9338952 m o h o n 4 n N _. . . .2 . cw . . i o MW . 57/, , ..c o o: in r G. uuuuuuuuuuuuuuu m nnnnnnnnnnnnnnn fix. it w, u .o . ./m. M0 4-. 1 u a .Q. Av I . ..o m . m g p - $5... 5. o m. . -upm . is: x . -2 3.9 o .9 8.2 -e. . L 3.2 .. 3,1,2 90 3.3 H p 1» w _ h — ON .sildionscomtficoo 8355.355: 0 0 n 0 n v N p q 1 11 A] 4 14 0 Av \RxW ............. {ax uuuuuuuuuuuuuu .a “v-11 , - 1 N . o , - - - - b - - - 7 ...... é .o.. ,. . n... v M. - m I 1 0 J DNIQ ¢. V 1 OF 3.0 -m. no.0 t 01.0 C O. 00.0 A. F n h P _ Lr NF 00.2.5 *3 .0 80.39.30 SIMPFOH P “"3 32 necessarily reflect the in situ properties of that sample. For poorly sorted grain-size distributions, the smaller grains will fill void space between larger grains. In effect, reducing the connectedness of the pores pace and reducing the effective hydraulic conductivity. The discrepancy between the deep and shallow model behavior may also indicate a problem with the accuracy of the deviation measurement. There is significant uncertainty associated with the borehole deviation data collected with the borehole colloidoscope. Until a more accurate method of measuring well deviation can be performed this issue will remain uncertain. Regardless, the dramatic effect deviation has on some simulation results warrants consideration of well deviation especially in highly detailed, small-scale modeling exercises like this one. Borehole deviation might also be a factor influencing the observation points in the aquifer. Borehole deviation information is not available for observation wells because the diameter of the wells was too small (2.54 cm) for the colloidoscope to measure. Because there is no way to measure the deviation of the observation points we have to assume they are straight. Since this is unlikely, observation deviation becomes another factor contributing uncertainty to simulation results. It appears that the concentration deviation method of calculating model residuals does the best job of comparing model results to measured values. The quantile calculation tends to be less accurate because it searches for the time when 50% of the mass has been sampled. For example, in well 11-45 (Figure 9) the 7 well case with deviation clearly does a better job at representing the concentration history than the other simulated scenarios. However, the 50% quantiles are nearly the same (Figure 14) 33 because 50% of the simulated mass for the other scenarios passes at the same time as the 7 well case with deviation. In contrast, the concentration deviation method examines the total mass under the curve, for which the 7 well case with deviation reduces the residual much more than the other scenarios, illustrated in Figure 14. Likewise, the total concentration deviation, when compared to the sum of the quantiles method, does a better job of representing simulation performance. For well 11- 45 (Figure 9) the all of the quantiles calculated for this case (10 through 90%) for the 3, 5 and 7 well cases will be similar to the quantiles for 7 well case with deviation. Even though the 3, 5 and 7 well case simulation results are much worse than the 7 well case with deviation, the quantiles do not reflect this as well as the total mass method. The only potential problem with using the total concentration deviation method is the same problem associated with all of the methods used, low sampling fi'equency of the measured tracer arrival. The concentration deviation method does a better job of comparing model behavior to measured results than the other two methods evaluated for this study. 34 .: :03 3% 8:006: 8955080 .88 8 055:0 $8 me "—85:88 ”3 950$ 83533252 8953.352 o a a. m .1. v n N a a a a .1. n u p q T 1 q . 4 4 ON q a _ 1 a 5 O 1.. . o .. o .. o .9 ..... o 1 1 3 1 1 1 N o ttttttttttt m uuuuuuuuuuuuuuu 1.0. 111111111111111 m. 1 8 £11 1.. ... 1 4 1 1.0 1 1 o 111111 1 0: , AV 1 a - .1 . - 4 1 . .1. o #00 1 w 1 0 1111111111 111w111111111111111m11111111111111 ,, W 2111 11m c a m ., W . . 1 , 1 m 1 . . 3; 89 p x , 1. 0 p , a ,. a 1 1 N. ,. W .1 4 1 1 1 a 1 ON? W 1 ., . ,1...er U x .lo. . m. -.l. ..1 z r 1 i . a . 1 .. ,. Mo: 1 ,1 . 1 a: x .p b ,a :4 \ 1 ., 18: 1 .. 1 : as: ,3. 8.111111: 9.: o ... 8.: -4 8.: -m1 . . 1 on: 1. 1.8. 3.: 1 a 12 9.: .0 9.: o 001—... . 061:. H p h p p b L? g A p _ _ _ — mSEspo 8.92.8080 .0 55980 8.230 88 B 805980 or 35 CHAPTER II PHYSICAL CHARACTERIZATION INTRODUCTION Several physical property relationships were explored during the characterization of the aquifer. Measured property values of porosity, permeability, grain size and ground-penetrating-radar (GPR) velocity were compared and used to further characterize the aquifer. Empirical relationships, like Kozeny—Carman, were explored to obtain hydraulic conductivity estimates from grain-size information. Since grain-size distributions are typically easier to measure than permeability estimates, any grain size relation that yields reliable permeability estimates would be valuable. In addition, several planes of crosswell radar tomography were collected and it was hoped that a relationship could be found between hydrogeologic properties and GPR velocity. Such a relationship, would allow the GPR to be used to improve estimates of hydro geo lo gic properties where no direct measurements were available. GRAIN-SIZE ANALYSIS During the aquifer characterization process, grain-size analysis was performed on subsections of the core samples collected from the delivery well gallery. Grain-size sample locations are shown in Figure 15. Grain-size samples were sieved using a series of 11 mesh sizes (2000, 840, 600, 500, 425, 300, 250, 212, 180, 150, 75 microns). This produced a total of twelve different grain-size categories, because everything less than 75 microns was collected in a pan below the series of sieves. It was hoped that a method to estimate hydraulic conductivity fiom grain-size information would provide useful results, since it was easier and faster to perform the grain-size measurements than the permeameter measurements. 36 De pth Location ofGtain Size Samples -8 r 1 O -10 . . o o ' -12 ~ . . I -14 L . . O O '16 * ' o o 1 ' ' i -18 - o . . ' o 3 ’ -20 . ’ o o ' O -22 e . : . -24 ~ g ' 0 ~26 _ : 2 0 I _28 g 1 1 1 1 1 1 D2 D4 D6 D8 D10 D12 D14 Well Number Figure 15: Location of grain size samples along delivery well gallery. 37 Exploratory data analysis of the physical properties of the aquifer examined the relationships between depth, porosity, effective grain size and hydraulic conductivity. The effective-grain size (Defl) is defined as: _ 21D13n1’ Defl— 2,193". (1.1) where D,- is the discrete grain size and n1- is the fraction of grains from a sample for a particular discrete grain size (Mavko et al., 1998). The discrete grain size categories used were based on the sieve sizes used in the sieve analysis. The effective particle size approach is used primarily for systems with poorly sorted grain size distributions (Mavko et al., 1998), which is the case at Schoolcraft where particle sizes range from gravel to silt size. Figure 16, 17 and 18 illustrate the relationships for the physical properties of the aquifer. For this site, high hydraulic conductivity sediments found deep in the aquifer, correspond to low porosity, and large effective grain sizes. A Kozeny-Carman type equation was used in an attempt to estimate permeability from grain-size information 11 = 19¢13111,,,.2 (1.2) where x is the permeability, B is the geometric fitting factor, ¢ is the porosity and def is the effective grain size (Mavko et al., 1998). This form of the Kozeny-Carman equation was used because it was robust and had few assumptions associated with the geometric configuration of the sample (e.g. tortuosity, packing, and sorting). Any problems associated with the geometric configuration of grains or pores were taken care of through the B term of the equation 1.2. The geometric fitting factor B was determined through a least squares regression between the measured permeability data fiom the permeameter 38 Comparison of Depth to Hydraulic Conductivity of samples .8- _ 03” 10 - . .3 -12- c, a... -14- .0 O -16— ’ ea E o ., , g 2-18— 0 . ‘d O 8 “O. C; -20- 011 ' ’ 0 ‘2T .. O O -24~ 3 o a ,1. -23- ;. . . . O . . O: '28 l l l I l 0.02 0.04 0.06 0.08 0.1 Hydraulic Conductivity (cm/s) 0.12 Deff (m) x 10‘3 0.6 0.4 H 0.2 Figure 16: Plot of relationship between measured hydraulic conductivity, depth and effective grain size. 39 Comparison of Depth to Porosity of samples Deff (m) x 103 .8.— .. 1 -10- . . . O . 4.8 m ‘ '12‘ 1 1.6 e o3 ’ '14— . O — -14 a _ a A46 3 : - «1.2 g C . a . . . 3"” ' ' . 1 “o a ‘ . -201— . 9’ . - ~0.8 . . -22— . . g ‘ . 0.6 -24— a. . ‘ ’ 0.4 -26- . . O . o 'g ’0 . 1 02 -28 I I 1 1 0.3 0.35 0.4 0.45 Porosity Figure 17: Plot of relationship between measured porosity, depth and effective grain size. 40 Comparison of Porosity to Hyd'aulic Conductivity of samples 0.46 — 0.44 L . . O O 0.42 — ¢ 0 of . o 0 4 — . ‘Wéi . .. ‘ J o ,1 1., :1. g 0.38 - . . 0 O (113 ‘“ 10 e» g 0.36 9 " . g Q A; o 0 34 3 o g Q o P ' o 0.32 — . o o 0.3 - o . 0.28 I I I I l 0 0.02 0.1 .04 0.06 0.08 Hydraulic Conductivity (cm/s) Figure 18: Plot of relationship between measured hydraulic conductivity, porosity and effective grain size. 41 Dei‘f(m)x10'a 1.0 ’1.6 - 0.8 0.6 ' 0.4 0.2 analysis and the estimated permeability fi'om the Kozeny—Carman equation. The B term takes into account the effect of grain packing and other effects associated with geometric configuration of grains and pores of the sample. The porosity values were estimated using the relation ¢=1—PA (1.3) ps where 4) is the porosity, p11, is the bulk density of the sample and p, is the sediment particle density of the sample (Freeze and Cherry 1979). This was done because no direct measure of the porosity was ever taken and because bulk density and particle density values were available for each sample. The bulk density values used for this analysis were measured during the repacked permeameter tests. The mass of the saturated sample and the volume the bulk sample occupied in the test chamber were known, so the bulk density was calculated as the mass of the saturated sample divided by the bulk volume. The particle density was measured using water displacement of a sediment sample. A known volume of water was placed into a graduated cylinder. Then a known mass of oven-dried sediment was poured into the graduated cylinder and the volume change was recorded. The particle density was then calculated as the mass of the sediment divided by the volume of water the sediment displaced. Given that the bulk density was calculated fiom a repacked sample, it is unlikely that the porosity measured in the lab is the same as the in situ porosity. Initially, the Kozeny-Carman equation did not match the measured permeabilities well (Figure 19). Comparison of Kozeny-Carman permeability estimates to permeameter derived estimates yielded a correlation coefficient of 0.5794. However, closer examination of Figure 19 revealed that the effective grain size (note color scheme on 42 k estimated Effective Grain Size x 1040 Estimated vs. Measwed Permeability x 10° .2 - correlation COOIIOCIOM 0.5794 1.6 1.4 ~ < 1.2 O l. . 1 . L 0.8 0 0.2 o 0.6 0.8 1 71.2 k measured x 10 to Figure 19: Kozeny-Carman estimated permeability vs. measured permeability. 43 Figure 19) was the dominating variable in the Kozeny-Carman equation. Three distinct color bands can be seen linking the large effective grain sizes to high permeability values and small effective grain sizes to low permeability values. This demonstrated the need for a different method of describing the effect of the grain-size distribution on permeability. Instead of performing the effective grain size calculation on the whole grain-size distribution, this calculation was performed on the grains that were 10% of the distribution and finer (d10). Figure 20 shows the result of using the effective d10 approach. Using the d10 effective grain size the correlation coefficient increased to 0.84. This was done with the d20 and d30 effective grain size as well, and similar results were achieved (Figures 21 and 22 respectively). The d20 effective analysis yielded a correlation coefficient of 0.81 and the d30 effective yielded a correlation coefficient of 0.83. Therefore, when applying the Kozeny-Carman permeability equation to this data, a d30 or smaller effective grain size should be used for reasonably close results to the values estimated using the permeameter measurements. There is also uncertainty associated with the permeability values estimated through the repacked permeameter process. During the repacked permeameter measurement, the grains in the sample are essentially mixed and homogenized. Therefore, any sorting as a result of deposition of the aquifer material is no longer represented with this method. Homogenizing the sample causes smaller grain sizes to fill in the pore space between larger grains. So the smaller grains become the dominant factor controlling the estimate of hydraulic conductivity. This behavior is reflected in the fact that by using an effective grain size that takes into account both large and small grain 44 Estimated vs. Measured Permeability using Effective d10 Correlation Coefficient 0.8438 It estimated 0 0.02 0.04 0.06 0.08 0.1 0.12 11 measured Figure 20: Kozeny-Carman estimated permeability vs. measured permeability using the d10 effective grain size. 45 Estimated vs. Measured Permeability using Effective d20 Depth (m) 0.12 — m 1’. i 26 Correlation Coefficient 0.8183 " 2‘ 0.1 L '1' < 22 0.08 b - ~ 20 . 0.06 — - . 8 .1: O O o O - - 16 0.04 - . e a 14 0.02 — 12 aw i‘x‘. . 7 " '6’ ‘ 10 or / 1 1 1 1 1 1 L L 1_ 1 h 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 k measured Figure 21: Kozeny-Carman estimated permeability vs. measured permeability using the d20 effective grain size. 46 Estimated vs. Measured Permeabiity using Effective d30 09 . Correlation Coefficient 0 0.8398 / “T“. . . . 0 0.01 0.02 0.03 0.04 Figure 22: Kozeny-Carman estimated permeability vs. measured permeability using the d30 effective grain size. 0.05 k measured 47 0.06 0.07 0.08 0.09 0.1 Depth (m) . 26 ~22 ~20 sizes, the Kozeny-Carman permeability estimates did not match the permeameter estimated permeabilities well. Meanwhile, using a srmll effective grain size (< D30), the Kozeny-Carman permeability estimates matched repacked permeameter estimates relatively well. Because repacked permeameter measurements might not accurately reflect the in situ value of permeability for some samples, the Kozeny-Carman equation that incorporates the total effective grain size of a sample may better represent the in situ values of aquifer permeability. This could be tested in future work through comparison of transport model results between two simulations that used the hydraulic conductivity estimates derived from both methods. RADAR TOMOGRAPHY Several studies have explored using geophysics to estimate hydraulic properties of aquifers. Despite this, there are few examples that have been able to successfully infer hydrogeologic properties fiom geophysical measurements. This is mainly because relationships between geophysical and hydrogeologic properties are usually non-unique and uncertain. Significant progress has been made using geophysical measurements to estimate hydraulic properties and identify aquifer structure (Bourbie and Zinszner, 1985; Rubin et al., 1992; Hyndman et al., 1994; Hyndman and Gorelick, 1996; Rea and Knight, 1998; Hubbard et al. 1999; Ezzedine et al., 1999; Hyndman et al., 2000; Hubbard et al., 2001). Electromagnetic waves have also been applied to obtain values of soil water saturation (Topp et a1. 1980; Alharthi and Lange, 1987; Hubbard et al., 1997; Chan and Knight, 1999), which is equal to the porosity for saturated aquifers. The goal of geophysical research performed in Schoolcrafl was to relate GPR velocity to porosity and use those porosity estimates to improve solute transport models. 48 Although research using GPR has provided encouraging results, no research has incorporated porosity estimates obtained from GPR into solute transport models. Porosity is a key parameter in solving the advective flow equation shown in equation 1.4: V, = Iii (1.4) 11, d] where VJr is the average linear velocity, K is the hydraulic conductivity, n, is the effective porosity and dh/dl is the hydraulic gradient (Fetter 1994). Modeling advective solute transport requires accurate porosity values, although it is rarely measured at field sites, so homogeneous values between 25 % and 40 % are often chosen for transport simulations. The heterogeneous porosity estimates obtained from GPR should provide a more accurate representation of the aquifer porosity than an approximate value between 25% and 40%. Several planes of cross-well radar tomography were shot across the well grid using 100 MHz antennas. Data was collected from 6 m bgs to the bottom of the well casing which varied from 22.5 m bgs to 27.5 m bgs depending on the well used. This was done because the water table was at approximately 5 m bgs at this site and we were only interested in the saturated zone of the aquifer. Shots were collected at increments of 0.5 m to avoid spatial aliasing of the data. The first arrival of each shot was picked using the pulseEKKO software (Sensors and Software, 1997). The first arrival picks were used as input for the straight ray inversion code provided in the pulseEKKO software package (Sensors and Software, 1997). When the inversion convergence criterion was met, the tomogram was displayed for visual inspection. Of the several planes of radar that were processed only five planes were used. It is speculated that problems with possible well deviation and picks of first arrivals may have caused the other planes of radar data to fail to converge in the inversion code. The 49 five planes of GPR data that provided reasonable results were from wells Dl-D15, P9- P10, PlO—P6, P11-P6 and P12-P6 shown in Figure 23. Comparing radar velocity to hydraulic conductivity revealed an unexpected relationship. High radar velocities correlated to high hydraulic conductivities as well as low porosities (Figure 24). However, for well-sorted sandy aquifers, like the aquifer in Schoolcraft, it is expected that high hydraulic conductivities correlate to high porosity (Freeze and Cherry, 1979). At this site, the aquifer is composed of very coarse sediment (gravel and coarse sand) at the bottom of the aquifer. Because of the range of sediment sizes the smaller sand grains may be filling in the pore space between the larger gravel. Although the porosity decreases, the size of the connecting pore throats are large and can communicate water extremely well resulting in high hydraulic conductivities in this deep zone. Therefore, the relationship at Schoolcratt is not the expected direct relationship between hydraulic conductivity and porosity. GPR velocity is inversely proportional to the square root of the dielectric constant given in the equation 2 0.5 V0”: -2—[ 1(1) +1] (1.5) ,us 8a) where VGPR is the GPR velocity, p is the magnetic permeability, a is the dielectric permittivity, ais the conductivity, and (0 is the angular frequency of the wave (Reynolds 1997). This is often simplified to the following in low loss mediums VGPR = yJ; (1'6) where VGPR is the GPR velocity, c is the wave velocity in air, and s is the dielectric constant of the medium (Reynolds, 1997). 50 Figure 23: Map of GPR tomography planes using a normalized velocity. Labels indicate well locations. 51 WWWW-M ul‘ 0... "8 .'. I. .. . I ' ..r"n. ’ . I. I O Ow. I C Ow... : I-.. ..U "' i ! . .1... O . - .0 i! O h. . 1" ‘i v m M Rel 'onship between GPR velocity porosity and hydraulic conductivity. Figure 24 52 The dielectric constant of sands depends primarily on the volumetric content of water, and for saturated systems the volumetric content of water is equal to the porosity. Because of the large contrast between the dielectric constant of water and sand, 80 and 3.5 respectively (Reynolds, 1997), and if a binary system of quartz and water is assumed, then some form of mixing theory may be used to estimate the volumetric percent of sand grains and water. To estimate porosity from the GPR tomograms, the GPR velocities need to be converted to dielectric constants using equation 1.6. The dielectric constants obtained from this equation are effective dielectric constants influenced by the respective amounts of sand and water in the system. In order to estimate the porosity, the individual effect that each constituent had on the effective dielectric constant was determined using three methods: Topp’s relation, differential effective medium theory (DEM) and the complex refractive index method (CRIM). Topp’s relation is based on a regression fiom a variety of soil samples and shown as a two part equation: 19, = -5.3x10‘2 +2.92x10’2xa -5.3x10“1c,,2 +5.3x10'6xa3 (1.7) x, =3.03+9.3a, +146.0a,,2 —76.70,3 (1.8) where 6,. is the volumetric water content and It}, is the apparent dielectric constant (Topp et al., 1980). This equation does not work well with soils that have hrge amounts of clay due to bound water in clay minerals. Differential effective medium theory is a mixing theory that assumes a binary system of a host material containing spherical inclusions of another material. In this case, 53 it can be approached in two ways, using either water or quartz as the host material. These two approaches provide an upper and lower estimate for the porosity using the equation 1_y=[£2_8DEM(y)]|: 81 ]3 (19) £2 - 51 8.05M 0’) where y is the volumetric content of spherical inclusions, 6.051.102) is the apparent dielectric constant, 8; is the dielectric constant of the material 1 and 82 is the dielectric constant of material 2 (Berryman, 1995). For this exercise a value of 3.5 was used for the dielectric constant of sand and a value of 80 was used for the dielectric constant of water (Reynolds, 1997). The complex refractive index method is another mixing theory that determines the percent fiaction of materials based on a composite value of a physical property, in this case the dielectric constant. The equation is given as: JF=f,,/Z+f,,/Z (1.10) where x. is the apparent dielectric constant, k; is the dielectric constant of material 1, x2 is the dielectric constant of material 2, fl is percent fraction of material 1 and f2 is percent fraction of material 2 (Mavko et al., 1998). Again 3.5 was used as the dielectric constant of sand and 80 was used for the dielectric constant of water. None of these methods assume any prior knowledge of the geometry of water and sand system. Because of this, a set of bounds called the Hashin-Shtrikman bounds, are used to define the narrowest possible range of values without any previous knowledge of the geometry of the constituents (Mavko et al., 1998). For a two-component system this is defined as: 54 £i=sl+ f2 (1°11) (£2 —6‘1)-l+§% where a, is the dielectric constant of the material 1, 82 is the dielectric constant of material 2, f 1 is percent fraction of material 1, and f 2 is percent fi'action of material 2 (Mavko et al., 1998). Figure 25 shows the Hashin-Shtrikman bounds plotted along with the results of the three different methods used to calculate porosity. All cases fall within the bounds so the results are reasonable by this measure. Based on this figure, the DEM upper model and the CRIM model seem to do a better job of predicting the porosity than either the Topp relation or the DEM lower model. The Topp relation begins to fail as it approaches a volumetric water content of approximately 50% because it begins to go outside the Hashin-Shtrikman bounds near this point. The DEM lower bound does not provide reasonable values of porosity for this site. Values of porosity derived from the upper DEM model were calculated fi'om each of the five planes of radar data. The mean of the porosity for each tomography plane was adjusted to the mean measured porosity using the repacked samples (37.5%), which involved adjustments of 2.4 to 7.1 %. These mean-adjusted porosity estimates were then interpolated across the transport model grid using ordinary kriging. Figure 26 provides an example of the GPR derived porosity for the plane of radar between wells D1 and D15. Following the kriging, the interpolated porosity was input into the BTN package of RT3D and the tracer test discussed previously in the report was simulated. This result was compared against a transport simulation that used the mean of the measured porosity value as a constant throughout the model grid as well as a kriged estimate of the measured porosity values. The results of the three simulations are seen in Figure 27. 55 Hasrin-Shtrilrman Bounds compared to Porosity Estimates 80 L l I l r l l l r ' «3— Hashin-Shtriicnan Upper ”L." . 4}— Hashin-Shtrlianan Lower a: 70 h .31-} t." .. ,1” 60 ‘ 0'.” 7‘ 1:) _, J" ,1 u a ‘3' g DEM Upper bound :3" ,j' 8 so - , 0‘ , o e" .0 it? ‘ ‘. {if} :1 a: 4° ' ,x CRIM Method 4;. - 5 :1“ 15' Topp Relation (j 330 ‘- \? i . ‘ 'f _ a); i 20 _ 19p DEM lower bound mid L .. ' 1k . ”xv,“ {1:5 _ 101' \qu _ , ‘ ‘~ U -..~.-.11+‘H.* "V“ xii-“v.1 N““"" ‘ 0 l l l 1 l i l l l 0 0.1 0.2 0.3 0.4 0 5 0.6 0.7 0.8 0.9 1 Saturation (96) Figure 25: Plot of Hashin-Shtrikman bounds and results of various mixing theory methods to calculate water content. 56 D1 D15 Radar derived porosity from DEM upper relation. Depth (m) 0.46 [- 22 Correlation Coefficiem 0.1259 0.44 — . . 20 . i 0.42 - . . O 0 g 0.4 - 1. e .. e: ,3} .. C . . 0.381 ‘F e e 5; 1 1‘ C . g .. 0.36 . O o 34 - O . 0.65 I l I I l I "6.32 0.34 0.38 0.33 0.4 0.42 0.44 Radar Porosity Comparison of measured porosity to radar derived porosity. Figure 26: Radar derived porosity for delivery well gallery (DI-DIS) and the comparison of those values to measured values. 57 Based on Figure 27, using the DEM derived porosity as opposed to a constant mean porosity or the kriged estimate of the measured porosity does not significantly affect the simulation results. The concentration histories of the three different transport simulations do not show a significant difference. In addition, there is poor correlation between the measured and GPR derived porosity values (Figure 26). This suggests that either one or both of the porosity estimates are not accurate representations of the aquifer. The problems with GPR derived porosity could be due to several factors, including the effect of borehole deviations, poor processing of the radar tomography or problems with the assumption that the aquifer is a binary system. If the system has large amount of clay material or material that has significantly different dielectric constant than sand then the porosity values derived from binary mixing theories are no longer valid. Uncertainty with the measured porosity based on the repacked bulk density measurement, discussed in the previous section, may also be a factor. In addition to the porosity analysis using GPR, correlation lengths of GPR velocity were generated and compared to the correlation lengths of the hydraulic conductivity of the aquifer. The correlation lengths generated using GPR (25.6 m horizontal and 3.97 m vertical) were similar to those of the hydraulic conductivity (13.8 m horizontal and 5 .77 m vertical). GPR provided reasonable estimates of the aquifer correlation structure at this site, which indicates that such geophysical approaches add value to aquifer characterization studies. Knowledge of the correlation lengths may help with the design of aquifer sampling strategies that will adequately sample the heterogeneity of an aquifer. 58 O N AguovoEFo ouowoomomor ooovomooooeomo 2 N— o— >8. PH 5%an ESE—8 v5 hump—om “025308 5&an @3th have .8.“ muoug=oo=8 .593 32030 was vow—«EEG KN BamE no caisfim ameoa 3388 II F cesafim $8888.0on II o =3.§:ooaoo 598p. uptempo x doze—23m .EmEom @25332 l 0.0 _. mm av .85 59 Despite the problems encountered with using GPR to improve transport simulations at this site, there is potential for this technique to aid in aquifer characterization. GPR tomography can provide estimates of aquifer characteristics for portions of the aquifer that are not directly sampled by coring. Indirect measurements like GPR are often faster and less expensive to collect, but in order to provide meaningful information about the aquifer uncertainty in data acquisition and data processing must be minimized. Using GPR to obtain accurate porosity measurements will help improve aquifer characterization, which should lead to improved groundwater flow and transport models. 60 CHAPTER III CONCLUSIONS It is widely known that aquifer characterization is an integral part of any groundwater modeling study. However, costs are often prohibitive for performing the dense spatial sampling of hydraulic characteristics that is needed to adequately identify heterogeneity in aquifer systems. Once the spatial correlation structure is adequately represented, then interpolation of hydraulic conductivity estimates through kriging provides a good representation of aquifer heterogeneity. Based on the information obtained at the Schoolcraft Study area, three well cores of hydraulic conductivity data appear to provide enough information to adequately model the system. Even though the uncertainty of estimated hydraulic conductivity values was reduced by adding more hydraulic conductivity data, little change was seen in the performance of flow and transport values from the 3 well to the 7 well cases. This result is specific to this site and likely due to the long correlation lengths of hydraulic conductivity data found at this site. For highly detailed, small-scale modeling investigations, the effects of borehole deviation are not trivial. Borehole deviations have the potential to significantly affect local flow around wells. An accurate method of measuring borehole deviation should be used to provide more accurate modeling results, as well as to account for flow behavior in sensitive remediation systems like the one at Schoolcraft, Michigan. Methods like quantile analysis and concentration deviations can be used to quantify the ability of simulations to match measured results. It is important to note that to perform these armlyses, a relatively high sampling-fiequency is needed for the measured data This ensures accurate calculation of quantiles, as well as total mass of 61 solute under the curve. Of the three methods examined the total concentration deviation seemed to provide the best method of comparing simulations to measured results. The total concentration deviation better represents the actual concentration history. Certain cases occur when simulated concentration quantiles match the measured quantiles relatively well, while the simulated tracer mass is much different than the measured tracer mass. Therefore, care should be taken when applying these techniques to evaluate transport simulations. Thorough examination of the physical properties of the aquifer at Schoolcraft revealed many important relationships. Permeability calculated from grain-size information provided accurate results when the effective grain-size used was limited to d30 or lower. This is likely because the small grain-sizes are controlling the connecting pore throats and thus controlling the permeability. Another important relationship discovered at this site was that high hydraulic conductivity correlated to high radar velocity (low porosity), which contrary to the expected behavior for well-sorted sandy aquifers. This type of data analysis should be applied to any aquifer characterization study. Simply comparing the available data proved to be valuable in recognizing the important relationships between the aquifer properties at this site. Additionally, GPR has potential to improve flow and transport models by estimating an aquifer property that is not directly measured throughout the aquifer. The dense spatial sampling that geophysical methods provide could greatly increase the information about the aquifer being studied. Provided that the uncertainty associated with acquisition of the geophysics and the processing of the data are minimal. For this site, the porosity estimation should have made a significant difference in the transport 62 modeling. However, uncertainties associated with the radar data likely mask the effect of using the radar derived porosity. This research project provides an example of how flow and transport models based on poor aquifer characterization (layered average case) compare to models developed fi'om a well characterized aquifer (7 well case with deviation). A point exists at which model results do not significantly improve from incorporating additional data Once this point is reached, research efforts should change fi‘om characterization to design. Future aquifer characterization efforts should focus on collecting enough data to adequately determine the correlation lengths of hydraulic conductivity in the aquifer. Based on the correlation structm'e of the data and the scale of the study, the data collection strategy can be adjusted to minimize uncertainty of hydraulic conductivity estimates, while at the same time reducing the amount of hydraulic conductivity data needed to characterize the aquifer. Sites that are much more heterogeneous than the Schoolcraft site likely require more than 3 wells of hydraulic conductivity data to adequately simulate flow and transport in the aquifer. Correlation structure could also be estimated using additional information like GPR, or grain-size measurements. Incorporating different aquifer properties and using the resulting correlation structure of those properties, a sampling strategy could be developed that leads to adequate aquifer characterization. 63 BIBLIOGRAPHY Alharthi, A., and J. Lange, Soil water saturation; dielectric determination, Water Resources Research, 23 (4), 591-595, 1987. Berryman, J .G., Mixture theories for rock properties, in A Handbook of Physical Constants, in AGU reference shelf, [-3, edited by T.J. Ahrens, pp. 3 v., American Geophysical Union, Washington, DC, 1995. Boggs, J.M., S.C. Young, L.M. Beard, L.W. Gelhar, K.R. Rehfeldt, and EB. Adams, Field study of dispersion in a heterogeneous aquifer; 1, Overview and site description, Water Resources Research, 28 (12), 3281-3291, 1992. Bourbie, T., and B. Zinszner, Hydraulic and acoustic properties as a function of porosity in F ontainebleau Sandstone, JGR. Journal of Geophysical Research. B, 90 (13), 11,524-11,532, 1985. Chan, C.Y., and RJ. Knight, Determining water content and saturation fi'om dielectric measurements in layered materials, Water Resources Research, 35 (1), 85-93, 1999. Clement, T.P., RT3D - A modular computer code for simulating reactive multi-species transport in 3-dimensional groundwater aquifers, Pacific Northwest National Laboratory Report, PNNL-11720, 1997. Deutsch, C.V., and A.G. Journel, GSLIB : geostatistical software library and user '3 guide, x, 369 pp., Oxford University Press, New York ; Oxford, 1998. Dybas, M.J., M. Barcelona, S. Bezborodnikov, S. Davies, L. Fomey, H. Heuer, O. Kawka, T. Mayotte, T.L. Sepulveda, K. Smalla, M. Sneathen, J. Tiedje, T. Voice, D.C. Wiggert, M.E. Witt, and OS. Criddle, Pilot-scale evaluation of bioaugmentation for in-situ remediation of a carbon tetrachloride-contaminated aquifer, Environmental Science and Technology, ES and T, 32 (22), 3598-3611, 1998. Ezzedine, S., Y. Rubin, and J. Chen, Bayesian method for hydrogeological site characterization using borehole and geophysical survey data; theory and application to the Lawrence Livermore National Laboratory Superfund site, Water Resources Research, 35 (9), 2671—2683, 1999. Fetter, C.W., Applied hydrogeology, 691 pp., Macmillan ; New York, 1994. Freeze, R.A., and J .A. Cherry, Groundwater, xvi, 604 pp., Prentice-Hall, Englewood Cliffs, N.J., 1979. Freeze, R.A., B. James, J.W. Massmann, T. Sperling, and L. Smith, Hydrogeological decision analysis; 4, The concept of data worth and its use in the development of site investigation strategies, Ground Water, 30(4), 574-588, 1992. F reyberg, D.L., A natural gradient experiment on solute transport in a sand aquifer; 2, Spatial moments and the advection and dispersion of nonreactive tracers, Water Resources Research, 22 (13), 2031-2046, 1986. Garabedian, S.P., D.R. LeBlanc, L.W. Gelhar, and M.A. Celia, Large-scale natural gradient tracer test in sand and gravel, Cape Cod, Massachusetts; 2, Analysis of spatial moments for a nonreactive tracer, Water Resources Research, 27(5), 911- 924, 1991. Gelhar, L.W., and CL. Axness, Three-dimensional stochastic analysis of macrodispersion in aquifers, Water Resources Research, 19 (1), 161-180, 1983. Harbaugh, A.W., and M.G. McDonald, User's documentation for MODFLOW-96, an update to the US. Geological Survey modular finite-difference ground-waterflow model, 56 pp., US. Geological Survey, Reston, Va, 1996. Hess, K.M., S.H. Wolf, and M.A. Celia, Large-scale natural gradient tracer test in sand and gravel, Cape Cod, Massachusetts; 3, Hydraulic conductivity variability and calculated macrodispersivities, Water Resources Research, 28 (8), 2011-2027, 1992. Hubbard, S.S., J. Chen, J. Peterson, E.L. Majer, K.H. Williams, D.J. Swift, B. Mailloux, and Y. Rubin, Hydrogeological characterization of the South Oyster Bacterial Transport Site using geophysical data, Water Resources Research, 37 (10), 2431- 2456, 2001. Hubbard, S.S., Y. Rubin, and E. Majer, Spatial correlation structure estimation using geophysical and hydrogeological data, Water Resources Research, 35 (6), 1809- 1825, 1999. Hubbard, 8.8., Y. Rubin, and EL. Majer, Ground-penetrating-radar-assisted saturation and permeability estimation in bimodal systems, Water Resources Research, 33 (5), 971-990, 1997. Hyndman, D.W., M.J. Dybas, L. Fomey, R. Heine, T. Mayotte, M.S. Phanikumar, G. Tatara, J. Tiedje, T. Voice, R. Wallace, D. Wiggert, X. Zhao, and CS. Criddle, Hydraulic characterization and design of a full-scale biocurtain, Ground Water, 38 (3), 462-474, 2000a. Hyndman, D.W., and SM. Gorelick, Estimating lithologic and transport properties in three dimensions using seismic and tracer data; the Kesterson Aquifer, Water Resources Research, 32 (9), 2659-2670, 1996. 65 Hyndman, D.W., J.M. Harris, and SM. Gorelick, Coupled seismic and tracer test inversion for aquifer property characterization, Water Resources Research, 30 (7), 1965-1977, 1994. Hyndman, D.W., J .M. Harris, and SM. Gorelick, Inferring the relation between seismic slowness and hydraulic conductivity in heterogeneous aquifers, Water Resources Research, 36 (8), 2121-2132, 2000b. James, B.R., and R.A. Freeze, The worth of data in predicting aquitard continuity in hydrogeological design, Water Resources Research, 29(7), 2049-2065, 1993. James, BR, and SM. Gorelick, When enough is enough; the worth of monitoring data in aquifer remediation, in AGU 1994 spring meeting, edited by Anonymous, pp. 178, American Geophysical Union, Washington, DC, United States, 1994. LeBlanc, D.R., S.P. Garabedian, K.M. Hess, L.W. Gelhar, R.D. Quadri, K.G. Stollenwerk, and W. W. Wood, Large-scale natural gradient tracer test in sand and gravel, Cape Cod, Massachusetts; 1, Experimental design and observed tracer movement, Water Resources Research, 27(5), 895-910, 1991. Loaiciga, H.A., An optimization approach for groundwater quality monitoring network design, Water Resources Research, 25 (8), 1771-1782, 1989. Mackay, D.M., D.L. Freyberg, P.V. Roberts, and J.A. Cherry, A natural gradient experiment on solute transport in a sand aquifer; 1, Approach and overview of plume movement, Water Resources Research, 22 (13), 2017-2029, 1986. Mavko, G., T. Mukerji, and J. Dvorkin, The rock physics handbook: tools for seismic analysis in porous media, x, 329 pp., Cambridge University Press, Cambridge ; New York, 1998. Mayotte, T.J., M.J. Dybas, and OS. Criddle, Bench-scale evaluation of bioaugmentation to remediate carbon tetrachloride-contaminated aquifer rmterials, Ground Water, 34(2), 358-367, 1996. Monaghan, G.W., G.J. Larson, and GD. Gephart, Late Wisconsinan drift stratigraphy of the Lake Michigan Lobe in southwestern Michigan, Geological Society of America Bulletin, 97(3), 329-334, 1986. Rea, J ., and R. Knight, Geostatistical analysis of ground-penetrating radar data; a means of describing spatial variation in the subsurface, Water Resources Research, 34 (3), 329-339, 1998. 66 Rehfeldt, K.R., J.M. Boggs, and L.W. Gelhar, Field study of dispersion in a heterogeneous aquifer; 3, Geostatistical analysis of hydraulic conductivity, Water Resources Research, 28 (12), 3309-3324, 1992. Reynolds, J .M., An introduction to applied and environmental geophysics, ix, 796 pp., John Wiley, Chichester ; New York, 1997. Rubin, Y., G. Mavko, and J. Harris, Mapping permeability in heterogeneous aquifers using hydrologic and seismic data, Water Resources Research, 28 (7), 1809-1816, 1992. Sittig, M., Handbook of toxic and hazardous chemicals and carcinogens, xxvi, 950 pp., Noyes Publications, Park Ridge, N.J., U.S.A., 1985. Steinman, W.K., Geological and hydrogeological characterization of the Schoolcraft Aquifer, Schoolcraft, Michigan, Masters thesis, Western Michigan University, Kalamazoo, 1994. Sudicky, E.A., A natural gradient experiment on solute transport in a sand aquifer; spatial variability of hydraulic conductivity and its role in the dispersion process, Water Resources Research, 22 (13), 2069-2082, 1986. Topp, G.C., J .L. Davis, and AP. Arman, Electromagnetic determination of soil water content; measurements in coaxial transmission lines, Water Resources Research, 16(3), 574-582, 1980. Vereecken, H., U. Doering, H. Hardelauf, U. Jaekel, U. Hashagen, O. Neuendorf, H. Schwarze, and R. Seidemann, Analysis of solute transport in a heterogeneous aquifer; the Krauthausen field experiment, Journal of Contaminant Hydrology, 45 (3-4), 329-358, 2000. Wagner, B.J., Sampling design methods for groundwater modeling under uncertainty, Water Resources Research, 31 (10), 2581-2591, 1995. Wagner, B.J., Evaluating data worth for ground-water management under uncertainty, Journal of Water Resources Planning and Management, 125 (5), 281-288, 1999. 67