Published online 3 October 2006
Published in Vadose Zone J 5:1129-1142 (2006)
DOI: 10.2136/vzj2006.0015
© 2006 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Estimating Large-Scale Fracture Permeability of Unsaturated Rock Using Barometric Pressure Data
Yu-Shu Wu,
Keni Zhang* and
Hui-Hai Liu
Earth Sciences Division, Lawrence Berkeley National Lab., 1 Cyclotron Rd., Berkeley, CA 94720
* Corresponding author (kzhang{at}lbl.gov)
Received 20 January 2006.
 |
ABSTRACT
|
|---|
We present a three-dimensional modeling study of gas flow in the unsaturated fractured rock of Yucca Mountain. Our objective was to estimate large-scale fracture permeability, using the changes in subsurface pneumatic pressure in response to barometric pressure changes at the land surface. We incorporate the field-measured pneumatic data into a multiphase flow model for describing the coupled processes of liquid and gas flow under ambient geothermal conditions. Comparison of field-measured pneumatic data with model-predicted gas pressures is found to be a powerful technique for estimating the fracture permeability of the unsaturated fractured rock, which is otherwise extremely difficult to determine in field studies with large scales of interest. In addition, this study demonstrates that the multidimensional flow effect on estimated permeability values is significant and should be included when determining fracture permeability in heterogeneous fractured media.
Abbreviations: CFu, Crater Flat undifferentiated unit • CHn, Calico Hills nonwelded unit • ESF, Exploratory Studies Facility • PTn, Paintbrush Tuff nonwelded unit • TCw, Tiva Canyon welded unit • TSw, Topopah Spring welded unit • UZ, unsaturated zone
|
INTRODUCTION
|
|---|
DURING THE PAST TWO DECADES, as the proposed site of a geological repository for storing high-level radioactive waste, the unsaturated zone (UZ) of the highly heterogeneous, fractured tuffs at Yucca Mountain, Nevada has been extensively investigated. Driven by the need to conduct long-term performance assessment of the repository, many site characterization studies have been performed, and various types of data have been collected. Because of the large temporal and spatial scales involved, it has been recognized that the development and use of numerical models are necessary to make quantitative evaluations and long-term predictions of flow and transport processes in the UZ under future climates. This need has motivated a continual research effort to develop and apply large mountain-scale flow and transport models (e.g., Wu et al., 1999, 2002) for future repository performance analysis. Thus, most site characterization investigations have been conducted (i) to understand unsaturated flow and transport processes in the UZ and (ii) to estimate various hydrological parameters required for model input and predictions (e.g., Rousseau et al., 1999; Wu et al., 2003).
Even with the significant progress made in characterizing the Yucca Mountain UZ system since the 1980s, the complexity in both site geological conditions and physical processes has challenged quantitative characterization efforts. In particular, determining the appropriate model input parameters on the temporal and spatial scales relevant to understanding the nuclear waste disposal system and assessing the repository performance remains a very difficult task. In a continual research effort, many core samples have been taken from boreholes, tunnels, and outcrops for different hydrogeological units or layers. These samples have been very useful for making laboratory determinations of rock-matrix hydrological properties. However, resolving how to properly estimate fracture flow properties at large spatial scales is still a challenging task. This is because fracture flow properties, such as fracture permeability, are related to flow processes occurring at the large spatial scales (10
100 m) required as input for model predictions. These large-scale model parameters for fractures are generally more difficult to measure at the site than those for the rock matrix.
Estimation of large-scale flow parameters for the unsaturated rock at Yucca Mountain relies primarily on inverse modeling studies, which incorporate field- and laboratory-measured moisture data to obtain parameter sets (e.g., Bandurraga and Bodvarsson, 1999; Ghezzehei and Liu, 2004). However, past studies using inverse modeling have concluded that fracture flow properties are not very sensitive to measured moisture data from core samples of the site. For example, measured liquid-saturation and water-potential data for the rock matrix (although very useful for determining matrix properties) do not provide much information about fracture properties. This is because under ambient conditions associated with low infiltration and arid climates, fractures contain little moisture (i.e., they are "dry"). At the same time, as shown in field tests, the fracture system is well connected and highly permeable, with a saturated fracture hydraulic conductivity many orders of magnitude higher than infiltration rates. Under such conditions, slight changes in moisture conditions, caused by the ambient infiltration, have little impact on fracture flow responses. On the other hand, the dry, large pore space in fractures makes almost entire fracture apertures available for gas flow. Therefore, pneumatic signals and data from air injection tests at the Yucca Mountain site are found to be more sensitive to fracture flow properties and have been used as the main data source for determining fracture flow parameters (e.g., Ahlers et al., 1999).
Airflow or gas (a mixture of air and water vapor) movement through the unsaturated zone is driven by changes in barometric pressure, temperature-induced density differences, wind effects, and topography (Rousseau et al., 1999; Weeks, 1987). Changes in barometric pressure at the land surface between day and night or seasonally result in corresponding changes in pneumatic pressures at different depths inside the UZ. For example, substantial airflow was observed in two wells drilled at Yucca Mountain, showing air exchange between the subsurface and the atmosphere, with flow into the wells during winter and out of the wells during summer (Weeks, 1987). The process is controlled by geothermal gradients or density differences. Subsurface pneumatic responses to surface barometric pressure fluctuations reflect the occurrence of gas flow through highly permeable subsurface fractures or permeable porous media. With more flow resistance at deeper depth for downward gas flow from the ground surface, subsurface pressure signals are amplitude-attenuated and time-lagged relative to the surface pressure signals. Therefore, naturally occurring gas-pressure variations provide a good indication of how well the formation is able to transmit gas flow, a measure then used to derive fracture permeability.
A number of studies in the literature use gas flow data or barometric pressure cycles in characterizing the vadose zone. Among the early research efforts, Weeks (1978) presented a systematic study (including a model formulation and numerical code as well as field testing results) for determining vertical permeability to air in a layered unsaturated zone. In 1987, he further provides a conceptual model and site-specific analysis of gas flow in the Yucca Mountain UZ caused by topographic effects.
In the literature, there are several types of gas flow data sources used in characterizing unsaturated porous media. The first type is the measurement of naturally occurring barometric pressure fluctuations on the surface and in the subsurface (e.g., Ahlers et al., 1999; Neeper, 2002). The second approach, called the air-injection permeability test (air-k test), involves the injection of air into boreholes or wells and monitoring of gas pressure or gas flux changes at and near injectors (LeCain, 1999; Illman and Neuman, 2000, 2001, 2003; Vesselinov et al., 2001a, 2001b; Illman and Tartakovsky, 2005; Illman, 2005). Both types of data have been applied to characterizing the Yucca Mountain site, with the air-k tests extremely valuable for calculating fracture permeability within a small-scale (1 m) spatial domain (Huang et al., 1999). In addition, Unger et al. (2004) reported on an indirect approach of using radon gas concentration data, measured in an underground tunnel at Yucca Mountain, to estimate large-scale fracture properties. It should be mentioned, however, that the application of subsurface pneumatic responses to surface barometric signals is in most cases limited to one-dimensional vertical flow scenarios (Weeks, 1987; Ahlers et al., 1999; Neeper, 2002). In particular, existing studies are suitable primarily for handling homogeneous or horizontal layered unsaturated formations only.
We describe a comprehensive modeling effort to estimate large-scale fracture permeability using pneumatic data measured from boreholes of the Yucca Mountain UZ. Our modeling approach, built on the current three-dimensional (three-dimensional) mountain-scale UZ flow model (Wu et al., 2003, 2004), incorporates pneumatic data into a modeling analysis of two-phase liquid and gas flow under ambient geothermal conditions. The gas flow modeling studies are performed under present-day infiltration conditions using the site-specific geological model and characterization data. Calibration of model-predicted gas pressures against field-measured pneumatic data leads to a methodology for estimating fracture permeability in the unsaturated fractured rock, an important parameter that would otherwise be difficult to determine at large scales (Ahlers et al., 1999).
This study represents a continuation of the work by Ahlers et al. (1999) to characterize large-scale fracture permeability in the Yucca Mountain UZ. However, even though the methodology in this study is similar to that of Ahlers et al. (1999), there are significant differences between the two studies and their results. Specifically, the current study is focused on estimating large-scale fracture permeability through three-dimensional modeling studies, while previous investigations (Ahlers et al., 1999) consisted primarily of one- and two-dimensional model results, as well as some three-dimensional modeling efforts with limited spatial and temporal scale. The main objective of Ahlers' study was to develop a methodology for characterizing subsurface pneumatic responses at the Yucca Mountain. Our study differed in that it was based on (i) the updated geological framework model, (ii) the current hydrological UZ flow conceptual model, (iii) the newly designed three-dimensional numerical grid model, and (iv) the updated fracture and matrix properties. More specifically, gas flow within the Yucca Mountain UZ is shown to be a three-dimensional phenomenon in this study. In general, model dimensionality has a significant influence on model-estimated permeability values, and multidimensional flow effects should be accounted for when simulating gas flow in heterogeneous fractured media.
|
HYDROGEOLOGICAL SETTING, PNEUMATIC DATA, AND CONCEPTUAL MODEL
|
|---|
The geological setting (implemented in the three-dimensional pneumatic model of this study) for describing the Yucca Mountain hydrogeological condition is based on the current site-scale UZ flow and transport model (Wu et al., 2003), which is in turn built on a geological framework model (BSC, 2004a) for Yucca Mountain. This section briefly discusses the geological model, the measured pneumatic data, and the gas flow conceptual model used in this modeling study.
Geological Model
Figure 1 shows a typical vertical west–east cross section in the proximity of the repository, where the UZ is between 500 and 700 m thick and overlies a relatively flat water table. Geologically, Yucca Mountain is a structurally complex system of Tertiary volcanic, layered, anisotropic, and fractured volcanic rocks (Scott and Bonk, 1984). These volcanic formations consist of alternating layers of welded and nonwelded ash flow and air-fall tuffs. The primary geological formations, from the land surface downward, are the Tiva Canyon, Yucca Mountain, Pah Canyon, and Topopah Spring tuffs of the Paintbrush Group. Underlying these are the Calico Hills Formation and the Prow Pass, Bullfrog, and Tram tuffs of the Crater Flat Group (Buesch et al., 1995).
View larger version (47K):
[in this window]
[in a new window]
|
Fig. 1. Schematic showing typical vertical profiles of hydrogeological layers and units, as well as the conceptualized barometric pressure signals within a typical east–west cross section of the unsaturated zone flow model domain (see Fig. 2 for the cross-section location).
|
|
View larger version (59K):
[in this window]
[in a new window]
|
Fig. 2. Plan view of the three-dimensional model domain, showing model boundary, horizontal grid layer, major fault locations, and selected borehole locations.
|
|
For hydrological investigations, the UZ geologic formations have been categorized into several hydrogeological units based primarily on their degree of welding (Montazer and Wilson, 1984). These units are classified as the Tiva Canyon welded (TCw) hydrogeological unit; the Paintbrush Tuff nonwelded unit (PTn), consisting primarily of the Yucca Mountain and Pah Canyon bedded tuffs; the Topopah Spring welded (TSw) unit; the Calico Hills nonwelded (CHn) unit; and the Crater Flat undifferentiated (CFu) unit. Table 1 lists the geological units and layers for different hydrogeological units and the associated grid-layer information for the numerical model. These hydrogeological units and layers are generally three-dimensionally distributed and vary significantly in thickness and slope across the model domain, as shown in Fig. 2.
View this table:
[in this window]
[in a new window]
|
Table 1. Lithostratigraphy and correlations of model grid layer and hydrogeological unit used in the three-dimensional pneumatic model.
|
|
Pneumatic Data
As part of the Yucca Mountain site characterization effort, several deep boreholes instrumented in the UZ are continuously monitored to record changes in pneumatic pressure at different depths (Rousseau et al., 1999). Gas pressures are also measured at the land surface, which is referenced to atmospheric pressure at the time measurements in boreholes are taken. Measurements in subsurface boreholes are conducted through isolated pressure transducers or pressure monitoring ports. Monitoring points in these pneumatic monitoring boreholes are distributed primarily along the shallow TCw, PTn, and TSw units, and a few of the boreholes have monitoring points below the bottom of the TSw.
Pneumatic pressure data have been collected by USGS scientists at a number of boreholes, including UZ-1, UZ#4, UZ#5, NRG-6, NRG-7a, SD-7, SD-9, SD-12, and UZ-7a, with most of the measurements performed in 1995 and 1996. Locations of these boreholes are shown in Fig. 2. Available gas pressure data cover time periods ranging from slightly less than 6 mo up to 1 yr. The longer-period records include downhole pressure measurements that cover most of the annual barometric pressure cycle, a large portion of which was not or insignificantly disturbed by interference from construction of the underground tunnels. Boreholes UZ-1, NRG-6, NRG-7a, UZ#4, and UZ#5 were instrumented with downhole pressure transducers that measure absolute pneumatic pressure. Detailed discussion of the field measurements and test methods of gas-pressure pneumatic data can be found in Rousseau et al. (1999).
Figure 3 presents typical pneumatic pressure records for instrument stations in the monitored Borehole NRG-7a. Figure 3 indicates that sensors in the TCw record little to no amplitude attenuation or phase lag compared with the surface barometric signals. Results from spectral analyses of in situ pneumatic pressure responses to the synoptic pressure variations indicate that the phase lags range from a few hours to several tens of hours at different depths (Rousseau et al., 1999). As shown in Fig. 3, sensors in the PTn unit record increasing attenuation and lag with increasing depth. In most boreholes, sensors in the TSw unit monitor a similar amount of attenuation and a practically indistinguishable lag over the entire vertical interval of the TSw, because of the high density of fractures within the TSw unit. Horizontally within the TSw, the observed attenuation varies and appears to depend on the thickness of the overlying PTn (Ahlers et al., 1999).
View larger version (26K):
[in this window]
[in a new window]
|
Fig. 3. Pneumatic pressure data from Borehole NRG-7a for typical pneumatic pressure responses at the site. Data from sensors located in the middle TSw, upper TSw, PTn, and TCw geological layers, and land surface are shown.
|
|
The gas pressure signal consists of multiple components, including daily, seasonal, and annual changes as well as the interference from each other. Seasonal components of the annual barometric pressure cycle are characterized by high-frequency, large-amplitude synoptic-pressure signals during the fall and winter months and by lower-frequency, small-amplitude synoptic-pressure signals during the spring and summer months. The large amplitude can be up to five times larger than the daily (diurnal) signal. Mean atmospheric pressure values are also higher during the fall and winter months than during the spring and summer months. Note that with the construction of an underground tunnel, the Exploratory Studies Facility (ESF) of Yucca Mountain, the underground barometric records were disturbed to a certain extent by additional barometric signal sources. Comparison of the barometric signals in the tunnel and at the land surface shows that they are nearly identical. As construction of the tunnels brought signal sources close to the pneumatic monitoring boreholes, the downhole pressure in nearby borehole signals changed. These changes will be enhanced at the locations near faults or fractures that intersect the tunnels (Ahlers et al., 1999). Such interference effects from multiple sources at the surface and in the tunnels were observed in nearby boreholes.
In this study, to reduce the uncertainties associated with the impact of ESF tunnel construction on pneumatic responses within the UZ, only the in situ pressure data from four boreholes (NRG-7a, UZ-7a, SD-7, and SD-12) are used for estimating fracture permeability in model calibration. In particular, all calibrations and analyses are made for the periods before data at any of the boreholes are affected by penetration of the tunnel. Table 2 shows sensor locations (elevations and depths) and observation periods used in this study. Each borehole has four or five observation points. At least one sensor in each borehole is located in each of the TCw, PTn, and TSw units.
View this table:
[in this window]
[in a new window]
|
Table 2. Sensor locations and observation periods for pneumatic pressure measurements at the four boreholes used in this study.
|
|
Conceptual Model and Physical Processes
Variations in subsurface gas pressures are caused by bulk gas flow, mainly through relatively dry and well-connected fractures within the UZ. The gas flow is driven primarily by atmospheric barometric-pressure fluctuations on the land surface and geothermal-gradient-induced density differences, in addition to some minor effects (e.g., wind speeds over the mountain and atmospheric tide responses). As shown in Fig. 1, a barometric pressure cycle signal on the land surface is rapidly transmitted down into the mountain through the highly permeable and well-connected fractures of the TCw unit, with little flow resistance (see Fig. 3). The low-permeability matrix may have negligible impact on rapid gas transmission through this top unit. However, moving down into the nonwelded PTn, the large matrix porosity, higher matrix permeability and few fractures of this unit provide significant storage volume as well as resistance to downward gas flow, leading to large attenuation and phase lag in response to surface changes. Further down into the densely fractured TSw unit, similar patterns in gas-pressure changes to those on the land surface may be found. The significant attenuation or lag, seen in the TSw unit, results mainly from effects or energy loss while traveling through the PTn unit.
In general, the daily or seasonal gas-pressure fluctuation at the land surface causes corresponding changes in the unsaturated subsurface. The typical behavior of subsurface pneumatics consists of three different characteristic periods (Ahlers et al., 1999), as shown in Fig. 3 by measured borehole data. Short periodic cycles correspond to daily (every half day) heating and cooling events in the atmosphere, as well as tidal effects. Intermediate period variations, on the order of days to weeks, result from weather frontal systems as they move across the mountain. The longest period occurs yearly, caused by seasonal temperature variations in the local atmosphere. Note that pneumatically static pressure decreases with increasing elevation and that the surface barometric pressures have the lowest mean value. As the surface pressure signal propagates into the subsurface, the amplitude of the signal gradually decreases, and the phase of the signal will be delayed. Furthermore, faults and more permeable fractured zones, where they exist, will provide a shortcut for transmitting rapid gas flow deeper into the UZ.
In addition to pneumatic processes, the ambient UZ system is also subject to other hydrologic, geochemical, and geothermal processes. Because gas flow caused by barometric pressure variations is a short-time phenomenon relative to other processes, only moisture and heat flow are included as relevant in our modeling study. The barometrically induced transient gas flow is analyzed under steady-state water and heat-flow conditions under the present-day infiltration and ambient geothermal conditions. In particular, the following two conceptualizations and assumptions are made in this study. First, ambient water and heat flow in the UZ system is at a quasi-steady-state condition, subject to spatially varying steady-state infiltration on the ground surface. Second, hydrogeological units and layers, as defined by the geological model, are internally homogeneous, unless interrupted by faults or alterations.
|
MODELING APPROACH AND NUMERICAL MODEL
|
|---|
The three-dimensional nature and complexity of the UZ geological system, as well as the coupling of liquid and heat flow, make it necessary to use a numerical modeling approach for conducting gas flow analyses. In this section, we describe the modeling approach used for simulating transient gas flow and for handling fracture–matrix interaction, the numerical scheme and codes, numerical model grids, and input parameters. We also discuss treatment of initial and boundary conditions used in the modeling study.
Modeling Approach and Numerical Code
The dual-permeability concept, as built into the UZ flow and transport model (Wu et al., 2003), is used in this study to simulate transient gas flow as well as matrix–fracture interaction in the unsaturated fractured rock. In this approach, global fluid flow is considered to occur not only between fractures but also between matrix blocks. In addition, the fluid and heat flow between fractures and the matrix is evaluated using a quasi-steady-state approximation (Warren and Root, 1963; Pruess and Narasimhan, 1985).
Model calibration and simulation of gas flow in this study were performed using the EOS3 module of the TOUGH2 code (Pruess et al., 1999; Wu et al., 1996). In the TOUGH2 code, gas and liquid two-phase fluid flow and heat transfer are described using the general multiphase Darcy's Law associated with energy and mass transfer through porous media. An integral finite-difference scheme is used for spatial discretization, and time discretization is performed with a backward first-order, finite-difference scheme. The resulting discrete nonlinear algebraic equations are written in a residual form and solved using Newton–Raphson iterations.
Numerical Model Grid
The three-dimensional numerical model grid used in this study is shown in plan view in Fig. 2. The three-dimensional model grid was generated from an integral finite-difference scheme (Pan et al., 2000) using an irregular, unstructured, three-dimensional control-volume spatial discretization. As shown in Fig. 2, the three-dimensional model grid covers a model domain approximately 20 km2 in area. The model grid consists of 980 mesh columns of fracture and matrix continua per grid layer, for a total of 86 440 gridblocks and 350 000 connections in a dual-permeability mesh. Vertically, the model grid has an average of 45 computational grid layers. This model grid was designed initially for modeling ambient heat flow within the UZ (Wu et al., 2003). Note that in the three-dimensional model, faults are explicitly represented by vertical 30-m-wide zones.
Model Input Parameters
Most rock and fluid-flow parameters, except fracture permeability, are taken from several related research reports of UZ flow investigations (Wu et al., 2003; Ghezzehei and Liu, 2004; Pan and Liu, 2004). These parameters include matrix and fracture porosities, matrix permeability, and matrix and fracture van Genuchten parameters. Table 3 lists fracture permeability and porosity for all the model layers. Fracture permeability data of Table 3 consists of (i) uncalibrated or initially estimated values of Column 2 (Pan and Liu, 2004), (ii) results of one-dimensional model inversions of Column 3 (Ghezzehei and Liu, 2004), and (iii) the three-dimensional calibrated results of Column 4 from this study. Among the three types of fracture permeability data, the permeability values in Column 2 are estimated using air-injection tests and used as initial conditions for one-dimensional model inversions, while one-dimensional model results of Column 3 are used for initial conditions for three-dimensional model calibration. Note that ranges of air injection testing–derived fracture permeability values of Column 2 of Table 3 for geological units of Yucca Mountain UZ are similar to those obtained from cross-hole injection tests (Table 2; Illman and Neuman, 2001) at a different site.
The fracture porosity is estimated from gas tracer tests conducted in the unsaturated zone of Yucca Mountain (Ghezzehei and Liu, 2004), and its average value is at 1% (Column 4 of Table 3) for most of the upper geological units above the CHn. Note that the fracture porosity values, as used in the current study in Table 3, are on the same order of magnitude as estimated air-filled porosity in the studies by Illman and Neuman (2001) for a different site and by Vesselinov et al. (2001b) for another site. Except in the lower CHn units, Table 3 shows that the Yucca Mountain tuff has relatively lower fracture porosity values. Temperature- and pressure-dependent fluid properties, such as density, viscosity, and specific enthalpy, are calculated internally using the formulation embedded in the TOUGH2 code.
Boundary and Initial Conditions
In this study, UZ liquid and heat flow are assumed to be at a steady state under ambient conditions. This assumption establishes boundary and initial conditions for the gas flow analysis by running the three-dimensional model to steady state under the present-day infiltration rate. The model uses the ground surface of the mountain (or the tuff–alluvium contact in areas of significant alluvial cover) as the top model boundary and the water table as the bottom model boundary. Both the top and bottom boundaries of the model are treated as Dirichlet-type conditions with spatially varying temperatures and pressures that remain unchanged through time. In addition, surface water recharge, as described by the net infiltration map, is applied using a source term in the fracture gridblocks.
All lateral side boundaries, as shown in Fig. 2, are treated as no-flow (closed) boundaries, which allow flow only along the vertical plane. The water table, the bottom boundary of the three-dimensional model, is shown to be a relatively flat surface, and for most model areas, the flat portion of the water table is about 730 m above sea level. The gas pressures are specified at 92 kPa at an elevation of 730 m, while surface gas pressures are determined by running the TOUGH2 code to steady state under given temperature, bottom pressure, and surface-infiltration conditions. This is necessary to generate a steady-state, equilibrated pneumatic-pressure boundary condition to avoid artificial airflow or circulation, which may occur if nonequilibrated pressures are imposed on both the top and bottom boundaries.
The top boundary temperature condition is determined by correlating average atmospheric temperature with surface elevations. Measured mean surface temperatures are used with a linear equation that correlates surface temperature with elevation. The surface temperatures Ts (°C) at any elevation Z (m) are then computed as constants according to the following equation (Wu et al., 1999):
 | [1] |
The initial estimates of the temperature distribution at the water table were made by contouring the temperature data, measured from 25 boreholes, at an elevation of 730 m (Sass et al., 1988; Rousseau et al., 1999). Because the water table is not perfectly flat over the model main, the actual water table temperatures are determined by linearly interpolating the values at 730 m and the model surface boundary elevation.
Water recharge imposed on the UZ model surface boundary is described by the steady-state net infiltration map. The mean infiltration map, estimated by studies of site climate and infiltration (BSC, 2004a, 2004b), is used. A plan view of the spatial distribution for the infiltration, as interpolated onto the three-dimensional model grid of Fig. 2, is shown in Fig. 4. The figure shows a flux distribution for the present infiltration, with higher infiltration rates in the northern part of the model domain and along the mountain ridge east of the Solitario Canyon fault. The average annual net infiltration rate is 3.6 mm yr–1 over the model domain.
View larger version (66K):
[in this window]
[in a new window]
|
Fig. 4. Plan view of net infiltration distributed over the three-dimensional model domain of the present-day mean infiltration, used for gas flow analysis.
|
|
With the surface recharge, surface and bottom temperature, and bottom gas pressure conditions specified as described, the three-dimensional model is run to steady state. The steady-state simulation results give estimates of the initial condition and gas pressures at each surface gridblock, which corresponds to a pneumatic static condition. The surface, time-dependent, barometric gas pressure condition is then imposed by adding a time-varying pressure change to the steady-state gas pressure at each surface block. The time-varying pressure change is defined by taking the actual barometric pressures, measured at one surface location, and subtracting the time-average gas pressure value over the entire measurement period. This approach for specifying rapidly varying gas pressures at the surface is equivalent to assuming the same pressure time-varying patterns or cycles for all the surface nodes or assuming static pneumatic conditions in the atmosphere. Nevertheless, the actual gas pressure values imposed are still correlated to elevations of the location through superposing onto the steady-state gas pressure, estimated by the steady-state flow simulation.
|
MODEL RESULTS AND ANALYSES
|
|---|
As discussed above, the main objective of the three-dimensional model calibration to pneumatic data is to aid in estimating large-scale fracture permeability for the three-dimensional UZ system. Past investigations (e.g., Ghezzehei and Liu, 2004; Wu et al., 2003; Ahlers et al., 1999) have found that transient time-dependent pneumatic data are among the most important data sources for estimating large-scale fracture permeability. By contrast, many other types of data, such as liquid saturation, water potential, temperature, or chloride data, are found to be relatively insensitive to constraining fracture properties. This study uses matrix and other rock properties (except fracture permeability), as determined from core samples, field observations, and steady-state moisture data (Ghezzehei and Liu, 2004; Pan and Liu, 2004; Wu et al., 2003). In addition, this study separates fracture permeability from fracture porosity in pneumatic diffusivities in different layers, using reliable fracture porosity data determined from different studies using gas tracers and other methods (Ghezzehei and Liu, 2004). In general, fracture permeability and porosity cannot easily be estimated separately using pneumatic data alone under nearly single-phase gas flow condition (Ahlers et al., 1999).
Initial Estimates
The first step is to estimate fracture permeabilities using inverse modeling to match the observed pressure signals. Since an automatic inversion involves a great number of forward runs and demands intensive computational effort, one-dimensional models (corresponding to selected boreholes where the pressure signal data are available, including NRG#5, NRG-7a, SD-7, and SD-12) are employed for the automatic inversion, using the iTOUGH2 code (Finsterle, 1999). Note that because airflow in the unsaturated fractured rock is a three-dimensional diffusive process, the one-dimensional models cannot in general capture the process accurately. Here, one-dimensional models are used for providing initial estimates only; that is, fracture permeability values obtained from the one-dimensional models will be adjusted in the following three-dimensional modeling study considering multidimensional flow effects.
Grids for the one-dimensional models are directly extracted from the three-dimensional model grid, and consequently the one-dimensional models in the vertical direction have the same spatial scale as that in the three-dimensional model. For a given one-dimensional model representing a borehole, the top boundary condition for the airflow is described using observed time-varying pneumatic pressure over a time period of 240 d. Before inversion, the combination of a steady-state liquid water flow field and pneumatically static conditions are determined and used as initial conditions for the inversion, similar to the approach used in specifying boundary and initial conditions for the three-dimensional model. A previous study by Ahlers et al. (1999) indicated that after 30 d, the initial conditions for the pneumatic pressures have an insignificant effect on simulation results. To exclude possible initial-condition effects, we matched the data observed after 30 d from simulation start for both one-dimensional and three-dimensional modeling studies.
Automatic inversion is applied to determine fracture permeabilities for the top two hydrologic units of TCw and PTn only (i.e., fracture permeabilities were estimated directly by the ITOUGH2 code) using simultaneous inversion of the pneumatic data measured from several boreholes in these two units. Fracture permeability values were estimated for each model layer in these units (Table 1). The objective function used in the inversion is defined as a summation of the square of the difference between simulated and measure air pressure values. A different analysis approach, as discussed below, was developed and employed to estimate permeability values for the lower TSw unit (below the PTn). This is because pneumatic data not only are limited, but also show little attenuation in gas-pressure signals across the TSw unit, which are difficult to match by simultaneous and automatic inversion using gas-pressure data in association with TCw and PTn units. However, the lack of significant attenuation in the TSw unit is considered an important feature because it implies that the TSw fractures are highly permeable and well connected. The calibrated fracture permeabilities for the model layers in the TSw unit must be consistent with the observation. Therefore, fracture permeabilities in the TSw are determined by forcing the simulated and observed gas pressure signals at the upper and lower sensor locations in the TSw to have similar degrees of attenuation for Borehole SD-12. Borehole SD-12 was chosen because the distance between the two TSw sensors within this borehole is the largest among all the relevant boreholes. The degree of attenuation in the barometric signal through the TSw in SD-12, or the relative difference between the signals at the two sensor locations, is determined quantitatively by evaluating
 | [2] |
where N is the total number of calibration time points, P is the gas pressure, and subscripts u and b refer to the sensors in the upper and lower (bottom) portions of the TSw within Borehole SD-12. Obviously, if the gas signals from the two sensors are identical, F should be equal to zero. For the SD-12 gas-signal data, the F value is 2.01 x 10–3 kPa. In this study, fracture permeabilities to be determined should give F values similar to the value calculated from the data, such that the simulated and observed gas-pressure signals have similar degrees of attenuation.
Since the gas-pressure data from the TSw unit are relatively limited compared with the two upper units, TCw and PTn, the insignificant attenuation and time lag between the upper-most and lower-most sensors are used for calibration. Otherwise, fracture permeabilities for different model layers in the TSw unit could not be independently estimated reliably or uniquely. Note that the attenuation and time lag are determined by the overall hydraulic properties between the two sensors, rather than by properties in a single model layer (or subunit). Permeability values derived from small-scale air-k test data for each model layer are multiplied by a common factor, d. To calculate an F value for the d factor, we run forward simulations, generating gas pressures for Eq. [2]. The best estimate of the d factor is obtained by matching the calculated F value with that derived from the observations through changing fracture permeabilities for the TSw unit only. Figure 5 shows calculated F values as a function of factor d. The F value is quite sensitive to d.
View larger version (7K):
[in this window]
[in a new window]
|
Fig. 5. Estimated F ratio as a function of factor d. The F ratio is determined as the ratio of the calculated F value to the F value for the gas pressure signal.
|
|
Figure 6 shows a comparison between simulated and observed pneumatic pressures at Borehole SD-12, obtained through one-dimensional model inversion. Similarly excellent matches were also obtained for other boreholes. Fracture permeability values (as listed in the third column of Table 3) obtained from one-dimensional models are on average about one to two orders of magnitude higher than those inferred from single-hole air-injection tests performed in the TSW or upper geological units (e.g., Pan and Liu, 2004). The difference between permeability estimates from in situ pneumatic pressure data and the air-injection test results is caused mainly by scale effects. The test condition for the air-injection tests implies that permeability estimates are representative of values on the order of meters (corresponding to the injection interval). On the other hand, pneumatic signals observed in the depth of the Yucca Mountain unsaturated zone are likely representative of large-scale flow processes on the order of 10 to 100 m (corresponding to the vertical distance from the ground surface and sensor locations where pneumatic signals are observed). It has been noted that large-scale effective permeabilities are generally larger than smaller-scale ones (e.g., Neuman, 1994). An intuitive explanation for such scale-dependent behavior is that a large observation scale, in an average sense, corresponds to a larger opportunity to encounter more permeable zones or paths where observations are made, which considerably increases the value of the observed effective permeability.
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 6. Comparison between simulated (broken line) and observed pneumatic gas-pressure data at Borehole SD-12, obtained though one-dimensional model inversion (subunits of the Tpcp and Tpbt2 are located within the TCw and PTn, respectively, while Tptrn and Tptpll correspond to upper and lower portions of the TSw, respectively).
|
|
Three-Dimensional Model Calibration
The second step is to calibrate, by a trial-and-error approach, fracture permeability using three-dimensional pneumatic simulation results against measured subsurface pneumatic data. A forward, rather than inverse, modeling approach was chosen primarily because of the computational intensity required by automatic inversion of the large-scale three-dimensional gas flow model. To capture the details of periodic gas-pressure variations for peak and valley values, the maximum time step is set to be 13000 s. The results of these gas flow simulations are then compared with field-measured pneumatic data from several boreholes simultaneously to examine the results of the above one-dimensional models and to re-estimate fracture permeability in several TSw layers. This section focuses on the model calibration and analysis using these three-dimensional pneumatic simulation results.
Note that the above one-dimensional calibration was performed for fracture networks located in nonfault zones. A two-dimensional model is used for the fault fracture permeability calibration using iTOUGH2 (Ghezzehei and Liu, 2004). Because faults are relatively planar in geometry, flow in and around a fault zone can be approximately captured by the two-dimensional model that is aligned approximately parallel to the dip of the beds and parallel to the dip of the fault. The data from Borehole UZ-7a represent the most complete data set within a fault zone and were used for calibrating the fault properties. The methodology and the data set used here are identical to those of Ahlers et al. (1999). The estimated fault fracture permeability values are slightly higher than nonfault ones (Ghezzehei and Liu, 2004), and they are used here directly (without further adjustment) in the three-dimensional simulations.
The three-dimensional pneumatic model was calibrated against the field-measured pneumatic data from four boreholes, SD-7, SD-12, NRG-7a and UZ-7a, observed between December 1995 and June 1996 (Wu et al., 2003). The model calibration results indicated that modification of fractured rock properties, as estimated by one-dimensional inversion in the TSw layers, is necessary for matching field-observed gas pressures. In particular, it was found necessary to reduce the fracture permeability of the subunits within the TSw by a factor of 15, as well as for the PTn units by a factor from 1.8 to 21. The final calibration results are provided in Table 3 (Column 4). In the top unit of the TCw, however, no adjustments in fracture permeability from one-dimensional model inversions were made. This is because one-dimensional flow appears to provide a good approximation for gas flow through the top, shallow TCw unit. Figure 7 shows a comparison between the observed gas-pressure and simulation results, in which the curve labeled "noncalibrated" is plotted using the simulations with the one-dimensional model estimated fracture properties and "calibrated" using the three-dimensional model results, with TSw subunit fracture permeability reduced by a factor of 15. As shown in Fig. 7, three-dimensional model calibrated results significantly improved the model match of the observed gas-pressure data, whereas the simulations with noncalibrated or one-dimensional model fracture permeability overestimated gas pressure responses at the corresponding elevation.
View larger version (21K):
[in this window]
[in a new window]
|
Fig. 7. Comparison of simulated and observed gas pressures at Borehole SD-12 during a 60-d period, using simulation results with and without three-dimensional calibration.
|
|
The comparison in Fig. 7 indicates that the fracture permeability of three-dimensional gas flow through the UZ should be lower than the one-dimensional model estimates. As discussed above, the lower fracture permeability for the three-dimensional model is attributed to the one-dimensional model fracture permeability being estimated by considering one-dimensional vertical flow paths only. In a three-dimensional model, where there exist some high-permeability channels, such as faults, or zones with high fracture density, or varying thicknesses of different permeable layers, three-dimensional gas flow is able to find those high-permeability pathways with the least resistance. As a result, the fracture permeability of a three-dimensional model will in general be lower than that estimated from one-dimensional models.
Comparisons of the model simulation results with the field-measurement data for boreholes SD-7, SD-12, and UZ-7a are shown in Fig. 8, 9, and 10, respectively. The simulation results demonstrate a good match with measurement data in general. Except in the TSw unit of Borehole SD-7, the three-dimensional simulation predicts a slightly smaller or delayed amplitude signal than the observation data. Many comparisons between simulated gas pressures, with and without fracture-permeability modifications, and field measurements show that the calibrated three-dimensional model has improved consistently in matching observation data. Note that for Borehole SD-7, the calibrated fracture-permeability values of the TSw unit may be even lower for a better match. This difference between simulations and data in this borehole might be caused by how the nearby fault affects pneumatic signal propagation. In addition, slightly greater differences between simulated and observed gas pressures in the lower TSw unit may be caused by coarse spatial discreteization or the effect of heterogeneity encountered over larger travel distance or through longer flow paths with depth.
View larger version (37K):
[in this window]
[in a new window]
|
Fig. 8. Comparison of simulated and observed gas pressures at Borehole SD-7 over a 60-d period, with three-dimensional calibration.
|
|
View larger version (31K):
[in this window]
[in a new window]
|
Fig. 9. Comparison of simulated and observed gas pressures at Borehole SD-12 over a 60-d period, with three-dimensional calibration.
|
|
View larger version (34K):
[in this window]
[in a new window]
|
Fig. 10. Comparison of simulated and observed gas pressures at Borehole UZ-7a over a 60-d period, with three-dimensional calibration.
|
|
Note that when comparing simulated and observed gas pressures at different locations of the three boreholes (Fig. 8, 9, and 10), we find that simulated gas pressures and their patterns of variation are in general consistent with observed values. In particular, the simulations consistently reproduce increases and decreases of pressure resulting from changes in barometric pressure at the ground surface. Overall, a reduction by a factor of 15 for the TSw fracture permeability and reduction by a factor of 1.8 to PTn layers 21 provide a better fit to observed pneumatic data for all locations and time periods and provides further confidence in the model's capability to simulate three-dimensional gas flow behavior within the UZ.
|
CONCLUSIONS
|
|---|
We presented a three-dimensional modeling study of gas flow to estimate the large-scale fracture permeability of the Yucca Mountain unsaturated zone. Our modeling approach incorporates changes in subsurface pneumatic pressure (in response to barometric pressure changes on the land surface) into the UZ flow model developed for the site. Results from gas flow simulations are compared with the measured pneumatic data from underground boreholes for the purposes of estimating fracture properties. As a result of calibration, fracture permeabilities, initially estimated by small-scale air-injection testing and one-dimensional model inversion, are found to need adjustment. In particular, a reduction factor of 15 is needed for the TSw layer fracture permeability (relative to the one-dimensional inversion results) to obtain an overall good match between the three-dimensional model predictions and pneumatic data. The ability to match field pneumatic data observed from multiple sources, including pneumatic data for a long time period, indicates the reliability of the numerical model to describe air and water flow processes within the Yucca Mountain UZ system through better estimates of fracture flow properties.
The results of this study indicate that using field-measured pneumatic data in combination with numerical modeling analyses provides a practical and powerful technique for estimating flow properties of vadose zone formations. As shown in this work, periodic responses of subsurface gas-pressure signals to surface barometric-pressure changes, which are easy to measure, reveal invaluable information on gas mobility in unsaturated porous media. Many other data sources, such as moisture and heat flow data, prove to be insensitive to the fracture permeability of the unsaturated fractured-porous media under ambient, low-infiltration conditions at Yucca Mountain. In addition, this work demonstrates that multidimensional effects on model-estimated permeability are significant when determining fracture permeability in heterogeneous fractured media. These effects can be captured only by three-dimensional modeling analyses on relevant model scales.
This study presents one example of our current research efforts in characterizing flow and transport processes of the UZ system at Yucca Mountain. It is important to mention that certain limitations and uncertainties exist for the mountain-scale gas flow model. One significant limitation presented in this paper is that it does not incorporate small-scale heterogeneities within each stratigraphic units, other than layer-wise conceptualization. The model grid used is relatively coarse for investigating details of gas flow in different units. In addition, the effects of seasonal changes in surface temperatures on gas flow are not included in the analysis. Future modeling efforts may be needed if these issues are considered to be important. Nevertheless, note that pneumatic data are generally easy to measure and effective to use for estimating permeabilities within the vadose zone. Pneumatic data analysis using a numerical modeling approach will likely find wider application in future characterization of the vadose zone.
|
ACKNOWLEDGMENTS
|
|---|
The authors thank Stefan Finsterle and Dan Hawkes for their review of this paper, and thank the USGS scientists for providing pneumatic data and thank our colleagues, Lehua Pan and Diana Swantek, for their help in preparing the manuscript. We also thank the VZJ reviewers for their insightful, critical, and constructive reviews and suggestions for improving this manuscript. This work was supported by the Director, Office of Civilian Radioactive Waste Management, U.S. Department of Energy, through Memorandum Purchase Order EA9013MC5X between Bechtel SAIC Company, LLC, and the Ernest Orlando Lawrence Berkeley National Laboratory (Berkeley Lab). The support is provided to Berkeley Lab through the U.S. Department of Energy Contract no. DE-AC03-76SF00098.
|
REFERENCES
|
|---|
- Ahlers, C.F., S. Finsterle, and G.S. Bodvarsson. 1999. Characterization and prediction of subsurface pneumatic response at Yucca Mountain, Nevada. J. Contam. Hydrol. 38(1–3):47–68.[CrossRef][ISI]
- Bandurraga, T.M., and G.S. Bodvarsson. 1999. Calibrating hydrogeological properties for the three-dimensional site-scale unsaturated zone model of Yucca Mountain, Nevada. J. Contam. Hydrol. 38(1–3):25–46.[CrossRef][ISI]
- BSC. 2004a. Geologic framework model (GFM2000). MDL-NBS-GS-000002 REV 02. Bechtel SAIC Company, Las Vegas, NV.
- BSC. 2004b. Simulation of net infiltration for modern and potential future climate. MDL-NBS-HS-000023 REV 00. Bechtel SAIC Company, Las Vegas, NV.
- Buesch, D.C., R.W. Spengler, T.C. Moyer, and J.K. Geslin. 1995. Nomenclature and macroscopic identification of lithostratigraphic units of the Paintbrush group exposed at Yucca Mountain, Nevada. Rep. USGS OFR 94-469. USGS, Denver, CO.
- Finsterle, S. 1999. ITOUGH2 user's guide. LBNL-40040. Lawrence Berkeley Natl. Lab., Berkeley, CA.
- Ghezzehei, T.A., and H.H. Liu. 2004. Calibrated properties model. MDL-NBS-HS-000003 REV 02. Bechtel SAIC Company, Las Vegas, NV.
- Huang, K., Y.W. Tsang, and G.S. Bodvarsson. 1999. Simultaneous inversion of air-injection tests in fractured unsaturated tuff at Yucca Mountain. Water Resour. Res. 35:2375–2386.[CrossRef]
- Illman, W.A. 2005. Type curve analyses of pneumatic single-hole tests in unsaturated fractured tuff: Direct evidence for a porosity scale effect. Water Resour. Res. 41:W05018. doi:10.1029/2004WR003703.[CrossRef]
- Illman, W.A., and S.P. Neuman. 2000. Type-curve interpretation of multi-rate single-hole pneumatic injection tests in unsaturated fractured rock. Ground Water 38:899–911.[CrossRef][ISI]
- Illman, W.A., and S.P. Neuman. 2001. Type-curve interpretation of a cross-hole pneumatic test in unsaturated fractured tuff. Water Resour. Res. 37:583–604.[CrossRef]
- Illman, W.A., and S.P. Neuman. 2003. Steady-state analyses of cross-hole pneumatic injection tests in unsaturated fractured tuff. J. Hydrol. 281:36–54.[CrossRef]
- Illman, W.A., and D.M. Tartakovsky. 2005. Asymptotic analysis of three-dimensional pressure interference tests: A point source solution. Water Resour. Res. 41:W01002. doi:10.1029/2004WR003431.[CrossRef]
- LeCain, G.D. 1999. In-situ pneumatic testing of boreholes. In J.P. Rousseau et al. (ed.) Hydrogeology of the unsaturated zone, North Ramp Area of the Exploratory Studies Facility, Yucca Mountain, Nevada. USGS Water-Resources Invest. Rep. 98-4050. USGS, Denver, CO.
- Montazer, P., and W.E. Wilson. 1984. Conceptual hydrologic model of flow in the unsaturated zone, Yucca Mountain, Nevada. USGS Water-Resources Invest. Rep. 84-4345. USGS, Lakewood, CO.
- Neeper, D.A. 2002. Investigation of the vadose zone using barometric pressure cycles. J. Contam. Hydrol. 54:59–80.[CrossRef][ISI][Medline]
- Neuman, S.P. 1994. Generalized scaling of permeabilities: Validation and effect of support scale. Geophys. Res. Lett. 21:349–352.[CrossRef][ISI]
- Pan, L., J. Hinds, and C. Haukwa, Y.S. Wu, and G.S. Bodvarsson. 2000. WinGrider: An interactive grid generator for TOUGH2. Version 1.0. Users' manual. Earth Sciences Division, Lawrence Berkeley Natl. Lab., Berkeley, CA.
- Pan, L., and H.H. Liu. 2004. Analysis of hydrologic properties data. ANL-NBS-HS-000042 REV 00. Bechtel SAIC Company, Las Vegas, NV.
- Pruess, K., and T.N. Narasimhan. 1985. A practical method for modeling fluid and heat flow in fractured porous media. SPE J. 25:14–26.
- Pruess, K., C. Oldenburg, and G. Moridis. 1999. TOUGH2 user's guide. Version 2.0. LBNL-43134. Lawrence Berkeley Natl. Lab., Berkeley, CA.
- Rousseau, J.P., E.M. Kwicklis, and C. Gillies ed. 1999. Hydrogeology of the unsaturated zone, North Ramp area of the exploratory studies facility, Yucca Mountain, Nevada. USGS Water-Resources Invest. Rep. 98-4050. USGS, Denver, CO.
- Sass, J.H., A.H. Lachenbruch, W.W. Dudley, Jr., S.S. Priest, and R.J. Munroe. 1988. Temperature, thermal conductivity, and heat flow near Yucca Mountain, Nevada: Some tectonic and hydrologic implications. Open-File Rep. 87-649. USGS, Denver, CO.
- Scott, R.B., and J. Bonk. 1984. Preliminary geologic map of Yucca Mountain, Nye County, Nevada, with geologic sections. Rep. USGS OFR-84-494. USGS, Denver, CO.
- Unger, A., S. Finsterle, and G.S. Bodvarsson. 2004. Transport of radon gas into a tunnel at Yucca Mountain—Estimating large-scale fractured tuff hydraulic properties and implications for the operation of the ventilation system. J. Contam. Hydrol. 70:153–171.[CrossRef][ISI][Medline]
- Vesselinov, V.V., S.P. Neuman, and W.A. Illman. 2001a. Three-dimensional numerical inversion of pneumatic cross-hole tests in unsaturated fractured tuff: 1. Methodology and borehole effects. Water Resour. Res. 37:3001–3017.[CrossRef]
- Vesselinov, V.V., S.P. Neuman, and W.A. Illman. 2001b. Three-dimensional numerical inversion of pneumatic cross-hole tests in unsaturated fractured tuff: 2. Equivalent parameters, high-resolution stochastic imaging and scale effects. Water Resour. Res. 37:3019–3041.[CrossRef]
- Warren, J.E., and P.J. Root. 1963. The behavior of naturally fractured reservoirs. SPE J., Trans. AIME 228:245–255.
- Weeks, E.P. 1978. Field determination of vertical permeability to air in the unsaturated zone. USGS Prof. Paper 1051.
- Weeks, E.P. 1987. Effects of topography on gas flow in unsaturated fractured rock: Concepts and observations, flow and transport through unsaturated rock. p. 165–170. In D.D. Evens and T.J. Nicholson (ed.) Flow and transport in fractured rocks. 2nd ed. Geophysical Monogr. 42. AGU, Washington, DC.
- Wu, Y.S., C.F. Ahlers, P. Fraser, A. Simmons, and K. Pruess. 1996. Software qualification of selected TOUGH2 modules. Rep. LBL-39490; UC-800. Lawrence Berkeley Natl. Lab., Berkeley, CA.
- Wu, Y.S., C. Haukwa, and G.S. Bodvarsson. 1999. A site-scale model for fluid and heat flow in the unsaturated zone of Yucca Mountain, Nevada. J. Contam. Hydrol. 38(1–3):185–215.[CrossRef][ISI]
- Wu, Y.S., G. Lu, K. Zhang, G. Zhang, H.H. Liu, T. Xu, and E.L. Sonnenthal. 2003. UZ flow models and submodels. Res. Rep. (AMR), LBID-2480, MDL-NBS-HS-000006, REV01. Lawrence Berkeley Natl. Lab., BSC, CRWMS M&O, Las Vegas, NV.
- Wu, Y.S., G. Lu, K. Zhang, and G.S. Bodvarsson. 2004. A mountain-scale model for characterizing unsaturated flow and transport in fractured tuffs of Yucca Mountain. Available at www.vadosezonejournal.org. Vadose Zone J. 3:796–805.[Abstract/Free Full Text]
- Wu, Y.S., L. Pan, W. Zhang, and G.S. Bodvarsson. 2002. Characterization of flow and transport processes within the unsaturated zone of Yucca Mountain, Nevada, under current and future climates. J. Contam. Hydrol. 54(3–4):215–247.[CrossRef][ISI][Medline]
TOP
ABSTRACT
INTRODUCTION
