Stochastic simulations of steady state unsaturated flow in a three-layer, heterogeneous, dual continuum model of fractured rock

Journal of Hydrology 307 (2005) 17–37 www.elsevier.com/locate/jhydrol

Stochastic simulations of steady state unsaturated flow in a three-layer, heterogeneous, dual continuum model of fractured rock Walter A. Illmana,*,1, Debra L. Hughsonb,2 a

Department of Geoscience, The University of Iowa, 121 Trowbridge Hall, Iowa City, IA 52242-1379, USA b National Park Service, Barstow, CA 92311, USA Received 5 January 2004; revised 6 September 2004; accepted 13 September 2004

Abstract Unsaturated flow through fractured rocks is a concern in the siting and performance of waste disposal facilities such as the proposed geological repository at Yucca Mountain, Nevada. We simulated a small two-dimensional cross-section of welded volcanic tuff, representative of Yucca Mountain stratigraphy, using spatially correlated, randomly heterogeneous fracture permeability fields and homogeneous matrix permeability continua representing various degrees of welding. Ten realizations each of fracture permeability fields for four different variances (s2Z0.5, 1.0, 1.5, and 2.0) were generated by the direct Fourier transform method (Robin, M.J.L., Gutjahr, A.L., Sudicky, E.A., Wilson, J.L., 1993. Cross-correlated random field generation with the direct Fourier transform method, Water Resour. Res. 29(7)2385–2398) independently for the welded Tiva Canyon Tuff (TCw), the non-welded Paintbrush Tuff (PTn), and the welded Topopah Spring Tuff (TSw), the latter being the proposed primary repository horizon. Numerical simulations were run for steady state flow at three different uniform water flux boundary conditions. Boundary conditions along the sides were impermeable and the base was open to gas and liquid flow. Numerical simulations were performed using the dual-continuum, two-phase flow simulator METRA, which represents matrix and fractures as dual overlapping continua, where liquid flux between continua can be restricted by a uniform factor. Fracture– matrix interaction was modeled as being less restricted in the PTn as compared to the TCw/TSw. Heterogeneous fracture permeability fields generated strong preferential flow in the TCw/TSw fracture continuum and significant preferential flow in the uniformly permeable, PTn matrix continuum. Flow focusing led to a local increase in saturation, which in turn increased relative permeability to water along the preferential pathways, causing water to flow faster. The development of the preferential pathways reduced the wetted surface area for fracture–matrix interaction leaving a large volume of the rock isolated from the preferential flow pathways. Statistical analysis of water flux values in the three units showed that the magnitude of the ensemble variance, indicating preferential flow pathways, increased with both the variance of fracture permeability and the water flux boundary condition, reaching a plateau in the TCw/TSw fracture continuum units after flowing approximately 10 correlation scales of the vertical fracture permeability. Ensemble covariance of water flux normal to the layering revealed long range * Corresponding author. Fax: C1 319 335 1821. E-mail addresses: [email protected] (W.A. Illman), [email protected] (D.L. Hughson). 1 Secondary appointment at the Department of Civil and Environmental Engineering and also affiliated with IIHR-Hydroscience and Engineering, University of Iowa, Iowa City, IA, 52242. 2 Fax: C1 760 255 8809. 0022-1694/$ - see front matter q 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2004.09.015

18

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

correlations in all units, which is longer than the correlation length of the fracture permeability fields. This suggests that the water flux boundary condition influences preferential flow in unsaturated dual continua media. These results suggest that careful analysis of information on fracture permeability variability obtained from pneumatic and hydraulic tests is an important component in understanding deep percolation processes. q 2004 Elsevier B.V. All rights reserved. Keywords: Fractured materials; Infiltration; Numerical models; Preferential flow; Stochastic processes; Unsaturated zone

1. Introduction 1.1. Problem statement Predicting unsaturated flow under various geological conditions is a concern for many environmental and water resource matters. Unsaturated flow through fractured rocks is a concern in the siting and performance of the high-level waste disposal facility planned at Yucca Mountain, Nevada. It is also a concern at many low-level nuclear waste disposal sites, landfills, underground storage tanks, and mine tailings. Water flow in the vadose zone is also a concern to geoscientists studying speleothems for paleoclimate reconstruction. Therefore, understanding flow and transport in fractured vadose zones is crucial to management of these problems. Predicting unsaturated flow in fractured rocks, however, is complicated by the highly heterogeneous nature of fractured media, interactions between fractures and rock matrix, and the effect of localized infiltration. Wetting front instability and a poor understanding of constitutive relationships for unsaturated fractured rock are additional sources of uncertainty. There is at present no well-established field methodology for characterizing fluid flow and contaminant transport properties of unsaturated fractured rocks. Unsaturated flow and transport in fractured rocks have been under investigation at field sites around the world. Some examples of field studies relevant to unsaturated flow in fractured rocks include those in layered welded and nonwelded tuff (Bodvarsson et al., 1997; Wang et al., 1999; Bussod et al., 1999; Salve and Oldenburg, 2001; Trautz and Wang, 2002; Hu et al., 2001; Salve et al., 2002 ; Glass et al., 2002; Salve et al., 2003), basalt (Faybishenko et al., 2000), chalk (Nativ et al., 1995; Dahan et al., 1999), and partially welded tuff (Yao et al., 2002). Many new important insights have been gained through the conduct and observation of these

experiments. In particular, some of the experiments revealed the development of preferential flow pathways under unsaturated flow conditions in fractured media (Dahan et al., 1999; Faybishenko et al., 2000; Glass et al., 2002; Salve et al., 2002; Yao et al., 2002). Potential causes of preferential flow include: (1) the presence of fractures (and/or macropores, voids); (2) wetting front instability; (3) funnel flow caused by structural features (Pruess, 1999); and (4) film flow along fracture surfaces (Tokunaga and Wan, 1997; Tokunaga et al., 2000). Despite this recognition, conditions that lead to preferential flow are often neglected in predictive modeling studies because of the costs of extensive site characterization and prohibitive computational demands. Coarse numerical grids may result in volume-averaging of flow processes and consequent loss of small-scale preferential flow and transport behavior. That may ultimately affect the predictions of flow and transport. The US Department of Energy (DOE) is evaluating the suitability of Yucca Mountain, Nevada, approximately 135 km northwest of Las Vegas, Nevada, as a potential repository for high-level nuclear waste. The DOE’s safety case includes multiple barriers, with the waste being sealed in corrosion resistant canisters prior to emplacement in drifts. A significant component of the nuclear waste isolation relies on the integrity of the geologic barrier consisting of the thick, layered, unsaturated tuff at Yucca Mountain. The geologic formation consists of alternating strata of welded and non-welded ash flow and air fall tuffs with overlying alluvium. The hydrogeologic units have been categorized (Montazer and Wilson, 1984) based roughly on the degree of welding. The units from the soil surface include the Tiva Canyon welded tuff (TCw), the Paintbrush non-welded tuff (PTn), the Topopah Spring welded tuff (TSw), the Calico Hills non-welded tuff (CHn), and the Crater Flat undifferentiated (CFu) units. The potential repository is to be situated primarily within the 300 m thick TSw unit.

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

The welded units typically have lower matrix porosities and are highly fractured, while the bedded and non-welded tuffs tend to have higher matrix porosities and low fracture densities (Montazer and Wilson, 1984; Flint et al., 2001b). The distribution of permeability from the surface of Yucca Mountain to the proposed repository horizon is a major controlling factor in the spatial and temporal distribution of deep percolation flux at the repository horizon and drift seepage. There is significant uncertainty associated with the estimates of unsaturated flow rate and its spatial distribution in the mountain. 1.2. Conceptual and numerical models of unsaturated flow in fractured rocks at Yucca Mountain The conceptual model of unsaturated flow and transport through layered, fractured tuff developed over many years and worked on by a number of scientists at Yucca Mountain has evolved over the past two decades. The National Research Council defines a conceptual model as an evolving hypothesis identifying the important features, processes, and events controlling fluid flow and contaminant transport of consequence at a specific field site in the context of a recognized problem (National Research Council, 2001). Flint et al. (2001a) provide a comprehensive review of the evolution of conceptual models of groundwater flow in the unsaturated zone at Yucca Mountain. Most early models assumed negligible or very low net infiltration rates, flow predominantly in the matrix because of its higher capillary suction, extensive lateral diversion by the PTn, and virtually no fracture flow through the TSw, which results in extremely low vertical percolation rates and long travel times. These conceptual models have evolved to reflect the identification of bomb pulse 36Cl at the repository horizon, which indicates the existence of rapid transport through the unsaturated zone (Sweetkind et al., 1996; Fabryka-Martin et al., 1997; Wolfsberg et al., 2000). Elevated levels of the isotopes 36Cl and 99 Tc observed at the repository level in the Exploratory Studies Facility (ESF) (Fabryka-Martin et al., 1996, 1997) and tritium in the Calico Hills unit (Yang et al., 1996, 1998) indicate travel times to the repository horizon of 50 years or less. This conceptual model for rapid fracture flow is further supported by

19

secondary calcite deposits observed in fractures throughout the exploratory studies facility (ESF), indicating significant fracture flow (Paces et al., 1996). The current conceptual model of unsaturated flow at Yucca Mountain (Flint et al., 2001a) has relatively high infiltration rates (as much as 80 mm/yr at some locations), fracture dominated flow in the welded units (TCw and TSw), and vertical matrix-dominated flow in the PTn (Flint et al., 2003), all of which lead to much shorter travel times through the system. This most recent conceptual model also includes dampening of episodic infiltration pulses in the PTn so that steady-state conditions generally prevail in the lower units. However, despite decades of research at the site, there is still significant uncertainty in unsaturated flow processes at Yucca Mountain, although significant progress has been made in recent years. New insights on unsaturated flow processes have been gained from recent field experiments performed in individual stratigraphic units at Yucca Mountain (Wang and Bodvarsson, 2003). As there is presently no method to characterize deep percolation directly, mathematical models (primarily numerical methods) are used in conjunction with sensitivity studies of input parameters to analyze unsaturated flow and transport. Recent efforts to model unsaturated flow and transport at Yucca Mountain include those by Wu et al. (1999) (see also Bandarruga and Bodvarsson, 1999), Birkho¨lzer et al. (1999), Hinds et al. (1999), Finsterle (2000), McLaren et al. (2000), Wu et al. (2002), Doughty et al. (2002), Zhou et al. (2003), Tseng et al. (2003), Finsterle et al. (2003), Haukwa et al. (2003), Li and Tsang (2003), Liu et al. (2003) and Bodvarsson et al. (2003). In particular, Haukwa et al. (2003) investigated the effect of spatial variability in fracture permeability on liquid seepage and moisture distribution in the vicinity of a waste emplacement drift in the unsaturated zone of Yucca Mountain. They generated three realizations of 2D fracture permeability fields for the TSw. Li and Tsang (2003) calculated seepage rates for various drift-degradation scenarios and for different values of percolation flux in the subunits of TSw. They generated three separate 3D permeability fields based on statistics of fracture permeability obtained from Bechtel SAIC Company (2003). Finsterle et al. (2003) analyzed data from air

20

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

injection tests geostatistically and generated a single conditional permeability field that was used to analyze seepage into the drift through forward and inverse modeling. Bodvarsson et al. (2003) carried out 2 and 3D numerical modeling studies using a singlecontinuum approach to investigate the development of discrete-fracture flow paths and flow-focusing phenomena in five subunits of TSw. In this study, parameters were derived from Ahlers and Liu (2000) and in particular, the standard deviation of the logarithm of fracture permeability ranged between 0.55–0.66 (corresponding to variances of 0.30 and 0.44), while correlation lengths of 1 and 3 m were chosen for their analysis. Zhou et al. (2003) developed a 2D numerical model with grid blocks of 20!2 m to investigate the effects of multiscale heterogeneity of fracture permeability, matrix permeability and the a-value of the van Genuchten model for the matrix. Layer scale spatial variability of properties was obtained from inverse modeling of field data. 1.3. Study objectives The main objective of the study is to investigate how unsaturated flow takes place in deep, layered, fractured rocks under the assumption that the rock can be treated as a dual continuum system with a spatially variable fracture continuum. We constructed a stochastic model that incorporates spatial variability in fracture permeability within three separate layers as a random spatially correlated multivariate normal field and compared our results to the case in which the fracture permeability is treated as homogeneous for each layer. To study preferential flow in the units, we investigated the evolution of flux variances and covariances for a range of variances (s2Z0.5, 1.0, 1.5, and 2.0) in fracture permeability with depth for varying infiltration rates in the different units. Our motivation for constructing the numerical model was not to create a Yucca Mountain sitespecific model of unsaturated flow. Instead, we simulated unsaturated flow in a small section of the TCw, PTn, and TSw to study unsaturated flow across these units with randomly heterogeneous fracture permeability fields coupled to a uniform rock property field. It is well known that preferential flow can cause rapid fluid flow and solute transport but to our knowledge, no analytical or numerical study has

been undertaken to show the effects of nonstationary statistical variability in fracture permeability on unsaturated flow in thick, layered fractured rocks for two-phase, dual continuum systems.

2. Description of numerical model 2.1. Mass and energy transport (METRA) features Numerical simulations were conducted using a two-phase, nonisothermal, dual continuum simulator METRA (Lichtner et al., 2000). METRA is based on a fully implicit time step scheme, and a block-centered grid, employing an integral finite difference scheme suitable for an unstructured grid. It is analogous to the computer programs TOUGH2 (Pruess, 1991) and NUFT (Nitao, 1998). The primary field variables for the two-phase problem are the total gas pressure (pg), partial air pressure (pa), and gas saturation (sg). The code allows for the use of various expressions relating relative permeability to saturation to describe unsaturated flow properties, including the van Genuchten (1980) relationship used in this study. Spatial heterogeneity in many of the input parameters can be incorporated by METRA. Details of the program are described by Lichtner et al. (2000). 2.2. Setup of numerical experiments The numerical model considered flow of water and air in three geologic units, based on the dual continuum concept. Exchange of fluids between the two continua is governed by Darcy’s law and by the area of the matrix-fracture interface that is open to flow. Cross-sectional area for liquid, gas, and energy flux between the fracture and matrix continua in this model can be restricted independently by a spatially uniform factor between zero and one. The model includes the region above the PTn in the TCw to the TSw units. The thickness of PTn is representative of conditions at Yucca Mountain. The computational domain consisted of overlapping 99 by 99 elements of 1 m3 for the matrix and fracture continua. We decided to use cubic grid blocks for the numerical mesh of 1 m dimension, commensurate with the support volume of fractured rock properties estimated from single-hole pneumatic injection tests.

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

This amounts to a total of 19,602 elements and 48,609 connections. The top 36 m of the model domain consisted of a single layer of TCw, while the next 33 m comprised the PTn, and the last 30 m consisted of the underlying TSw. We modeled a region deep enough in the unsaturated zone such that evapotranspiration and near surface thermal processes could be neglected. Numerical simulations were conducted at three rates of water flux (qaZ12.5, 22.5, and 42.5 mm/yr) applied uniformly at the top boundary for all cases studied. These rates were consistent with the DOE infiltration model for Yucca Mountain (CRWMS M & O, 2000), which predicted mean annual infiltration ranging from 0 to 20 mm/yr for modern-day climate conditions and 1–50 mm/yr for projected future climate conditions. For all cases, 2.5 mm/yr of the total water flux was applied to the matrix continuum uniformly at the top boundary to avoid numerical difficulties. A mixed boundary condition was specified at the entire top boundary for both cases, in which the gas saturation, temperature and liquid flux were specified. Fluxes were set for both fractures and matrix at the top boundary. The left and right boundaries were specified as no flow boundaries, while at the bottom boundary, we specified constant field variables for all cases. Isothermal conditions were maintained by fixing the temperature at the top and bottom boundaries at 20 8C. 2.3. Model parameters The base case model consisted of uniform parameters for each layer taken from the unsaturated zone hydrology model developed for the Total System Performance Assessment-Viability Assessment (TSPA-VA) (CRWMS M & O, 1998). These model parameters are shown in Table 1. These are averaged values for the subunits of TCw, PTn, and TSw hydrogeologic units from the TSPA-VA and are comparable to recently updated values (Ahlers and Liu, 2000), obtained from the 1D calibration of site data assuming variable infiltration scenarios. Fracture permeability is treated as a distinct, stochastic continuum while the other rock properties are uniform for each layer. The matrix permeability is treated as uniform since spatially distributed data at

21

Table 1 Dual continuum parameters for model simulations Symbol

TCw

PTn

TSw

km (m2) fm am (PaK1) mm Srm kf(m2) ff af (PaK1) mf Srf Xfm (m2)

9.68!10K18 0.10 9.83!10K7 0.33 0.23 1.48!10K12 1.85!10K4 1.93!10K4 0.49 0.01 4.90!10K4

1.29!10K14 0.38 6.60!10K5 0.33 0.12 1.26!10K13 6.69!10K5 1.49!10K3 0.45 0.01 3.10!10K1

8.32!10K17 0.09 1.65!10K5 0.26 0.08 3.24!10K12 1.09!10K4 6.66!10K5 0.49 0.01 5.76!10K5

km, matrix permeability; fm, matrix porosity; am, van Genuchten aparameter for matrix; mm, exponent m in the van Genuchten equation for matrix; Srm, residual saturation for matrix; kf,fracture permeability; ff, fracture porosity; af, van Genuchten a-parameter for fracture; mf, exponent m in the van Genuchten equation for fracture; Srf, residual saturation for fracture; Xfm, fracture–matrix connection area.

Yucca Mountain on matrix permeabilites are not available at the same scale as the fracture permeability estimates derived from single-hole pneumatic injection tests. In all three units, matrix permeability is significantly lower than fracture permeability, while the matrix porosity is significantly larger than fracture porosity. 2.4. Treatment of fracture permeability as a random field The treatment of fracture permeability as random and other input parameters as uniform is based on the geostatistical analysis of interstitial, pneumatic, hydraulic, and thermal property data of unsaturated fractured tuffs at the Yucca Mountain analogue Apache Leap Research Site (ALRS) in central Arizona, USA (Rasmussen et al., 1990). Analyzed data include 3-m scale estimates of permeability, laboratory (core-scale) estimates of permeability, fracture density, effective porosity, alpha-value of the van Genuchten function, Klinkenberg coefficient (Klinkenberg, 1941), water content, thermal conductivity, and dry rock specific heat obtained from Rasmussen et al. (1990). This and other studies (Illman et al., 1998; Chen et al., 2000) at the site have shown that fracture permeability varies the most spatially, while the other variables

22

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

are considerably less variable. We relied on data collected at the ALRS as similar geostatistical analysis of rock parameters collected at 1–3 m scales have not yet been published. Air injection tests conducted at the ALRS (Guzman et al., 1996; Illman et al., 1998; Illman and Neuman, 2001, 2003) show that those tests primarily measure the properties of the fracture system in partially welded tuffs. Among the findings and conclusions at the site (Chen et al., 2000), the most relevant to Yucca Mountain are (a) the pneumatic pressure behavior of fractured tuff at the site is amenable to analysis by treating the rock as a continuum on scales ranging from meters to tens of meters, and (b) this continuum is representative primarily of interconnected fractures. Both the single- and cross-hole test results have proven to be virtually free of skin effect (Illman and Neuman 2000, 2001), implying that they represent rock conditions unperturbed by the presence of boreholes. We expect similar conditions at Yucca Mountain. 2.5. Estimates of fracture permeability variability and spatial correlation At Yucca Mountain, fracture permeabilities obtained from single-hole pneumatic injection tests (LeCain, 1997, 1998; Wang et al., 1998) have been found to vary over four orders of magnitude in the Exploratory Studies Facility (ESF). These data, all collected in the subunits of TSw, indicate that heterogeneity of fracture permeability can be quite large within a stratum. More recently, Finsterle et al. (2003) published results from a geostatistical analysis of fracture permeabilities obtained from single-hole pneumatic injection tests conducted in the ESF and found that the mean base-10 logarithm of permeability (log10 k) equaled K12.13 with a standard deviation of 0.80 from 78 data points in Niche 3107. Corresponding values for Niche 4788 were log10 kZK11.79 with a standard deviation of 0.84 based on 63 values. Computed experimental variograms revealed a weak spatial correlation with correlation lengths ranging between 0.6 and 1.3 m. According to Finsterle et al. (2003), permeability obtained from single-hole injection tests located outside of the footprint of the niche were removed from the data set. They reported that the removed data set belonged

to a separate population from an area of relatively undisturbed rock of lower permeability. The variance and the correlation scale of the log10 k for the fracture continuum for our random field simulations were chosen based on the values reported in Finsterle et al. (2003). We assumed an exponential correlation structure of fracture permeability for convenience. We also obtained estimates of variances from geostatistical analysis of single- (Guzman et al., 1996) and cross-hole (Illman et al., 1998; Illman and Neuman, 2001, 2003) pneumatic injection test data by Chen et al., (2000) and Vesselinov et al. (2001) of ALRS for the generation of spatially correlated, random fracture permeability fields. Estimates of variance obtained from the geostatistical analysis of single-hole estimates of fracture permeability at the ALRS (s2Z0.87 from single-hole data and s2Z0.45 from kriged field) are comparable to estimates from Yucca Mountain (s2Z0.64–0.71). Information about the spatial variability of flow parameters such as fracture permeability is commonly obtained through geostatistical inference from small scale measurements. Estimates of fracture permeability variance from the numerical inversion of cross-hole test data are also useful because the estimates are based on larger scale tests sampling a larger volume of the rock that takes into account the connectivity of fractures better than single-hole estimates and are thus more representative of the numerical simulations that we consider here. Our higher estimates of variance (s2Z1.5 and 2.0) are based on the geostatistical inversion of cross-hole pneumatic injection tests by Vesselinov et al. (2001). 2.6. Generation of fracture permeability fields Spatially correlated random fields were generated using a direct Fourier transform method (Robin et al., 1993). Anisotropy in correlation lengths of fracture permeability were not considered since geostatistical analyses have not revealed such directional effects (Illman et al., 1998; Chen et al., 2000) at the ALRS. Random fields were generated for layered, heterogeneous media using the mean of the log10 of fracture permeability provided in Table 1 and four different variances for the simulations (0.5, 1.0, 1.5, and 2.0) to simulate varying degrees of heterogeneity in fracture permeability. The correlation length of the fracture

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

23

Fig. 1. Random fields of base 10 logarithm of fracture permeability (log10 k).

permeability field was set at 2 m. Ten realizations of random fields were generated for each value of variance bringing the total number of random fields to 40. Random fields were generated independently with different statistics for each layer and combined for the model runs. According to Journel and Huijbregts (1978), a two step approach for random field simulation is not only more consistent with the underlying physical phenomenon, but it also avoids violating the stationarity and statistical homogeneity requirements underlying most continuous random function models. Therefore, we focus on the reproduction of major heterogeneities, in our case

considering the stratification of tuff and differences in the hydrological properties of individual units. Fig. 1a–d show one realization of log10 fracture permeability, that was varied from s2Z0.5–2.0. No attempt was made to explicitly model discrete features including fractures, joints, and faults.

3. Results The steady state distribution of saturation in the fracture and matrix continua, respectively, with qaZ42.5 mm/yr for s2Z1.0 is shown in Fig. 2a and b. Plotted alongside are contour plots of the flux

24

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

Fig. 2. Steady state distribution of: (a) fracture saturation; (b) matrix saturation; (c) fracture water flux; and (d) matrix water flux with homogeneous fracture permeability for each layer. The applied flux at the top boundary is qaZ42.5 mm/yr.

magnitudes in the fracture (Fig. 2c) and matrix (Fig. 2d) continua. Streamtraces were generated using the software TECPLOT (Amtec Engineering, 2001) and are included to indicate flow paths in the fracture and matrix continua. A predictor-corrector integration scheme was used to trace particles through a steady-state vector field generated by METRA. Saturation and water flux contours are plotted in a similar manner (Figs. 4 and 5) for the other cases. Results show that water saturation is more-or-less uniform in the TCw/TSw fracture and matrix continua (Fig. 2a and b) except at unit contacts and at the bottom boundary, where the effects of the boundary

condition can be seen. Saturation is, in general, low in the fracture continuum in all units but is also low in the PTn matrix continuum. Saturation decreases in the PTn fracture continuum as the water moves progressively into the PTn matrix due to capillary effects. There is a thin zone of saturation buildup at the TCw/ PTn contact in the fracture continuum, due to the formation of a permeability barrier. There is also a saturation buildup at the PTn/TSw contact in the fracture continuum, due to a formation of a capillary barrier. Fig. 2c and d show that water fluxes are highest in the TCw and TSw fracture continua. Water enters

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

the PTn fracture continuum upon reaching the TCw/ PTn contact and converts to matrix flow reflecting the strong influence of the matrix on flow within the PTn unit. Streamtraces (Fig. 2c and d) indicate gravity dominated flow with a straight, downward pathway. Water fluxes are highest within the PTn matrix despite the lower saturations (Fig. 2d). As expected, capillary pressure is also uniform within the TCw and TSw fracture and matrix continua (Figs. 3a and b) but higher in TSw beneath the PTn. Capillary pressure is lower in the upper portion of the PTn fracture continuum and in the TSw matrix continuum.

25

Numerical simulations were also conducted with random fracture permeability fields by applying water uniformly at the top boundary at a rate of 12.5, 22.5, and 42.5 mm/yr to allow direct comparison with results from the base case model runs. The steady state distribution of capillary pressure in the fracture continuum (Fig. 3c) for a single realization of fracture permeability corresponding to s2 Z1.0 with qaZ42.5 mm/yr reveals a large variation in capillary pressure. In particular, the lower values of capillary pressure (corresponding to higher saturations) in the PTn fracture continuum extends downward further at xZ4, 16, 46, and 84 m in comparison to the base case.

Fig. 3. Steady state distribution of: (a) capillary pressure in the fracture continuum for the base case; (b) capillary pressure in the matrix continuum for the base case; (c) capillary pressure in the fracture continuum with log10 k of s2Z1.0 for realization nine; and (d) capillary pressure in the matrix continuum with log10 k of s2Z1.0 for realization 9. The applied flux at the top boundary is qaZ42.5 mm/yr.

26

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

Fig. 4. Steady state distribution of: (a) fracture saturation; (b) matrix saturation; (c) fracture water flux; and (d) matrix water flux with log10 k of s2Z1.0 for realization 9. The applied flux at the top boundary is qaZ42.5 mm/yr.

Similar behavior is seen in all other realizations not shown here. Fig. 4 shows the steady state distribution of fracture saturation (Fig. 4a), matrix saturation (Fig. 4b), fracture water flux (Fig. 4c), and matrix water flux (Fig. 4d) for a single realization of fracture permeability corresponding to s 2Z1.0 with qaZ42.5 mm/yr. Examination of Fig. 4a shows that the distribution of saturation in the fracture continuum is highly variable in the nonwelded and welded units. The saturation is also highest locally, where permeability in the fracture continuum is low and at

the TCw/PTn boundary in the fracture continuum, where the contrast in permeability causes a permeability barrier. As water moves progressively from the fracture to the matrix continuum in the PTn unit, the variability in saturation decreases with depth in the PTn fracture continuum. The saturation distribution in the PTn fracture continuum is markedly smoother, due to the diffusive nature of the matrix continuum. The highest saturations in the PTn fracture continuum correspond to low capillary pressures (Fig. 3c) and are due to the high input fluxes at the bottom of the TCw unit in the fracture

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

continuum. In contrast, higher values of saturations in PTn matrix continuum are found, where the saturations in the PTn fracture continuum are lower. Fig. 4c and d show the development of preferential flow paths (defined here as a flowpath that carries more water than the applied rate) in all three units of the fracture continuum and in the PTn matrix continuum. Narrow preferential flow paths develop within the fracture continuum almost immediately beneath the upper boundary despite the uniform application of water and without explicitly building in high permeability pathways or discrete features that represent fractures or fracture zones. The bulk of water flows through these narrow pathways, isolating a large volume of the fracture and matrix continua from flow. Preferential flow paths also develop in the PTn matrix continuum despite the uniform matrix permeability. The flow paths are broader in comparison to flow paths in the TCw and TSw fracture continua reflecting the strong capillary effects in the PTn matrix. Results from the cases in which the applied water flux was set lower (12.5 and 22.5 mm/yr) also show the development of preferential flow paths. The random field simulations revealed that varying the water flux rate does not affect the frequency of flow channels under steady state conditions. The degree of heterogeneity (i.e. the magnitude of fracture permeability variance) has a large effect on the development of preferential flow paths, the magnitude and direction of water flow and its corresponding distribution. Lowering the degree of heterogeneity to s2Z0.5 causes a more uniform distribution of saturation and capillary pressure in the fracture continuum and wider regions of higher saturations in the PTn matrix continuum. Preferential flow paths develop in both the fracture and matrix continua, with lowering of flow focusing. In contrast, as the degree of heterogeneity is increased to s2Z1.5 and s2Z2.0 (Fig. 5a–d), the degree of flow focusing increases further and causes the convergence of flow paths, leaving a large volume of the rock isolated from preferential flow pathways (Fig. 5c) particularly in the fracture continua of all three units. In fact, some of the flow paths terminate because water moves completely from the fracture to matrix continuum, as fewer number of flow paths carry a larger volume of water.

27

The steady state distribution of water flux for a single realization (realization nine) in the fracture continuum of the TCw unit at zZ20.5 m from the top of the boundary is shown in Fig. 6, where ‘MF’ stands for medium (22.5 mm/yr) flux. Here, ‘v05’ stands for s2Z0.5, ‘v1’ for s2Z1.0, ‘v15’ for s2Z1.5, and ‘v2’ for s2Z2.0. ‘Homog’ stands for the base case with homogeneous fracture permeability for each layer. The large variability in water flux is evident from this plot despite the uniform application of water at the top boundary. We see from the plots that in all cases water flow rates in preferential flow pathways can be locally very high (approaching ten times the applied flux) in the fracture continuum for the s2Z2.0 case. For example, water flux in the preferential flowpath 12 m from the left boundary is 208 mm/yr which is slightly higher than ten times the applied flux. For this case, water flux in the other preferential flow paths range between 21.3 and 208 mm/yr. However, for the regions outside of the preferential flow paths, water fluxes are comparable for varying degrees in heterogeneity (s2Z0.5–2.0) and in general lower than the water fluxes of the base case. We also found that water fluxes in the preferential flow paths are correspondingly lower for the lower variance cases. Fig. 7 shows the steady state distribution of water flux for a single realization (realization nine) in the fracture continuum of the TSw unit at zZ80.5 m from the top of the boundary. We note the differences in the magnitude of water flux and distribution of preferential flow paths in the TSw compared to the TCw fracture continua. Plots of water flux made for other realizations not shown here exhibit similar behavior. This suggests that one may not be able to extrapolate percolation flux estimates obtained in the shallow subsurface to depths reaching the repository level in the TSw unit. This is especially true when stratification of rocks with varying hydrological properties can cause the redistribution of water.

4. Statistical analysis of results To analyze the results from the numerical simulations from all realizations (total of 120 model runs), we computed the ensemble mean, variance, and covariance of water flux along each row of the model grid. Fig. 8 shows the mean water flux in the fracture

28

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

Fig. 5. Steady state distribution of: (a) fracture saturation; (b) matrix saturation; (c) fracture water flux; and (d) matrix water flux with log10 k of s2Z2.0 for realization 9. The applied flux at the top boundary is qaZ42.5 mm/yr.

continuum for all simulations (realizations 1 through 10 for all fluxes and variances) along the z-axis where ‘LF’ and ‘HF’ stand for low (12.5 mm/yr) and high (42.5 mm/yr) fluxes, respectively. It is evident from Fig. 8 that the mean water flux in the TCw fracture continuum stays constant, but decreases exponentially as it enters the PTn fracture continuum and approaches a much smaller value compared to the applied flux within 15 m of the TCw/ PTn contact. The increase in the mean water flux is abrupt at the PTn/TSw contact and increases gradually to a stable value in the fracture continuum for all

cases. A smaller percentage of the water flux is in the TSw fracture continuum compared to the TCw fracture continuum. The variance of water flux along each row of the model grid through the fracture and matrix continua represents the degree in preferential flow. Low variance in water flux implies its uniformity and increasing variance suggests convergence of flow causing larger variability in fluxes forming preferential flow paths. Examination of Fig. 9 shows that the computed variance is low near the top boundary indicating the influence of the uniform boundary

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

29

Fig. 6. Water flux in the fracture continuum of TCw at zZ20.5 m from top boundary. Here, ‘v05’ stands for s2Z0.5, ‘v1’ for s2Z1.0, ‘v15’ for s2Z1.5, and ‘v2’ for s2Z2.0, where ‘MF’ stands for medium (22.5 mm/yr) flux. ‘Homog’ stands for the base case with homogeneous fracture permeability for each layer.

Fig. 7. Water flux in the fracture continuum of TSw at zZ80.5 m from top boundary.

30

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

Fig. 8. Mean water flux in the fracture continuum for all simulations (realizations 1-10) of along the z-axis, where ‘LF’ and ‘HF’ stand for low (12.5 mm/yr) and high (42.5 mm/yr) fluxes, respectively.

Fig. 9. Variance of water flux in the fracture continuum for all simulations (realizations 1–10) along the z-axis.

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

31

Fig. 10. Covariance computed normal to each layer of TCw, PTn, and TSw units for log10 k with s2Z1.0 for realizations 1 through 10. The applied flux at the top boundary is qaZ42.5 mm/yr. Symbols are computed values of covariances while the solid and dashed curves are fit using nonlinear regression.

condition. As water flows downward, the variance increases with both the variance of fracture permeability and the infiltration boundary condition, reaching a plateau in each of the TCw and TSw fracture continua after flowing about 10 correlation lengths (20 m) in the fracture permeability. The variance decreases in the PTn fracture continuum exponentially as the water imbibes into the matrix from the fracture. Therefore, flow becomes uniform in the fracture continuum. As the water re-emerges into the TSw, the variance again increases signifying the redevelopment of preferential flow paths in the TSw fracture continuum. We also computed the ensemble covariances of vertical water flux in the same direction as the mean and variance, but calculated separately for each layer. Fig. 10 shows results from such a computation for 10 realizations with a log10 k variance of 1.0 and an applied flux rate of 42.5 mm/yr. Symbols represent computed covariances. The solid and dashed curves are model fits of the computed covariance values using the nonlinear regression routine PEST (Doherty et al.,

1994) in which an exponential covariance model was fit. It reveals that the covariances are highest at the TCw/PTn boundary and then decay rapidly with depth. Some oscillations are visible for covariances in the TCw and TSw while the covariance of the PTn decays monotonically due to the strong capillary forces smoothing the water flux variability. Fig. 11 shows covariances in the TSw with varying degrees of heterogeneity for the applied flux rate of 42.5 mm/yr. Examination of the two cases (Figs. 10 and 11) reveals that a range in correlation may be reached for all cases approximately at 12.5 m. Similar behavior is seen in the low flux case (12.5 mm/yr). Table 2 summarizes the covariance of vertical water flux at zZ0 m (R(0)), correlation length (l) estimated using PEST, their 95% confidence intervals, and depth estimates at which the covariance reaches zero calculated based on the fitted model. We note that the confidence limits provide only an indication of parameter uncertainty. They rely on a linearity assumption which may not extend as far in parameter space as the confidence limits themselves

32

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

Fig. 11. Covariance computed normal to TSw unit for log10 k with s2Z0.5, 1.0, 1.5, and 2.0 for realizations 1 through 10. The applied flux at the top boundary is qaZ42.5 mm/yr. Symbols are computed values of covariances while the solid and dashed curves are fit using nonlinear regression.

(Doherty et al., 1994). The correlation length of vertical flux is approximately twice that of the correlation length of the logarithm of fracture permeability (2 m) in the welded units. In all cases, the depth at which the covariance reaches zero varies from 18.3 to 33.0 m, is highest for the high flux case, and increases slightly with the variance.

the site. However, it should be possible to qualitatively compare our numerical results against a number of field experiments in order to assess the validity of our model results. 5.1. Visual field evidence of preferential flow in Topopah Spring welded tuff Glass et al. (2002) conducted an infiltration experiment in nonlithophysal Topopah Spring welded

5. Comparison of simulation results with field observations The numerical experiments that we conducted were not designed to reproduce any experiments conducted at Yucca Mountain but were intended to illustrate the general behavior of unsaturated flow in heterogeneous fractured rock. However, there are a number of field experiments conducted within individual layers at Yucca Mountain (Wang and Bodvarsson, 2003), some of which can be compared with these results. Many of the field experiments were conducted at 1–10 m scales in drifts and niches to investigate specific issues such as seepage into drifts and unsaturated flow in the PTn and so were conducted at different locations throughout

Table 2 Covariances at zZ0 m (R(0)), correlation lengths (l) estimating using PEST, their 95% confidence intervals, and the calculated depth at which the correlation reaches zero based on the fitted model Case

R(0)

l

z (m) when R(z)Z0

TCw-10-HF PTn-10-HF TSw-05-HF TSw-10-HF TSw-15-HF TSw-20-HF TSw-05-LF TSw-10-LF TSw-15-LF TSw-20-LF

450.8G55.6 1210.1G38.0 497.7G31.0 963.6G58.6 1386.1G88.8 1766.9G121.2 37.0G1.9 73.8G4.0 108.3G6.5 139.8G9.3

4.0G0.8 2.4G0.1 4.2G0.4 4.1G0.4 4.0G0.4 4.0G0.4 4.1G0.3 4.0G0.4 3.9G0.4 3.9G0.4

36.7 23.8 38.9 40.5 41.2 41.7 27.3 29.3 30.3 31.1

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

tuff at Fran Ridge (several kilometers away from Yucca Mountain) to visualize the subsurface movement of dyed water. Their principal objective was to map the fracture network and liquid structure beneath an infiltration event and to capture the gross dynamics of phase structure development within the network, during both infiltration and drainage (Glass et al., 2002). The experiment lasted over a 46 min period in which the infiltration and drainage events were followed by electrical resistance tomography (ERT) along 2D transects beneath the infiltration surface. The ERT data showed rapid infiltration of water and its drainage in a fraction of the rock monitored indicating that preferential flow was taking place. The rock was then excavated to a depth of w5 m and mapped in detail looking for traces of the dye. This mapping revealed strong connectivity in both the vertical and horizontal planes. Some of the notable features reported by Glass et al. (2002) include (1) fragmented flow along a single fracture; (2) concentration of flow at fracture intersections; and (3) bypassing of fractures, despite obvious physical connection to the flow field. In particular, they found pervasive staining near the infiltration surface to more fragmented and intricate structure at depth bypassing a large fraction of the fracture network despite its connectivity. These field observations agree with our numerical results qualitatively. In particular, the calculated water flux variance (Fig. 9) revealed strong channelization of water flux in the TCw and TSw fracture continua. The stabilization of flow field indicated by a constant water flux variance in the numerical experiments was not observed in this field experiment. The observed penetration of dye into the matrix was found to be negligible across the entire excavated region confirming that flow was constrained to the fracture network. This observation supports our model results in the TCw and TSw fracture continua. In addition, the observed dye tracer structure shows that our model generally represents the behavior of preferential flow. We acknowledge that our model obviously cannot capture the complex dye structure seen, as many mechanisms including gravitational instability and hysteresis in the saturation-capillary pressure curve were excluded from the numerical model. Despite the obvious limitations of our numerical model, we do see a resemblance of our model results to the visual

33

images provided in Nicholl and Glass (2002) and Glass et al. (2002). 5.2. Preferential flow in the Topopah Springs unit (TSw) Evidence for dominance in fracture flow at Yucca Mountain has largely been indirect field observations such as, poor attenuation of barometric pressure response in the tuff (Rousseau et al., 1997), young ages of perched waters (Yang et al., 1996, 1998), water associated mineral deposits on fracture surfaces (Paces et al., 1996) and bomb pulse 36Cl found in the ESF (Fabryka-Martin et al., 1996, 1997). A recent study by Salve et al. (2002) showed strong direct field evidence for fracture flow in the TSw. They presented results of their field experiments conducted at the scale of 1–3 m in order to investigate fracture flow and fracture/matrix interaction in the TSw. They conducted multiple releases of tracer-laden water in low and high permeability zones, where the permeabilities were determined through air injection tests at the 0.3 m scale. Saturation and water potential measurements were made in three boreholes drilled beneath the injection borehole. Slots were excavated beneath the injection and monitoring boreholes to quantify the volumes and rates of water flux. Results from these borehole infiltration experiments under constant head conditions showed intake rates varying from 0.5 to 119 ml/min indicating its large variability from one injection interval to the next. In addition, localized arrival of water at the slots suggests flow of water in a few well connected preferential pathways. Such variation is also evident in our model when we examine the variability in water fluxes on Figs. 4–7. 5.3. Unsaturated flow processes in the Paintbrush Nonwelded Tuff (PTn) At Yucca Mountain, the presence of the PTn in the unsaturated zone has led to the conceptualization of deep percolation becoming attenuated and laterally diverted once it reaches the PTn. Recent field experiments in Alcove 4 of ESF (Salve et al., 2003) provided important insights into resolving this issue. Salve et al. (2003) conducted a series of borehole infiltration experiments under constant head

34

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

conditions into the matrix and along minor subvertical normal faults. Saturation and water potentials were monitored in adjacent boreholes, while seepage was monitored in an excavation beneath the injection and monitoring boreholes. The field tests showed that a small portion of PTn matrix can transmit water quickly through discrete flow paths while the adjoining bulk matrix has much lower permeability. It was shown that a relatively dry PTn matrix can dampen infiltration events, however, a fault can convey a pulse of water over large distances if the matrix is wet (Salve et al., 2003). Comparable behavior is seen in our numerical simulations in which a preferential flow path develops in the PTn matrix continuum. In addition, all of our simulations for different variances and applied flux showed a strong dampening effect in PTn despite the presence of a spatially variable, yet correlated random fracture permeability field.

6. Discussion and concluding remarks Random field simulations using heterogeneous fracture permeability revealed the development of preferential pathways and focusing of flow, which can have a significant effect on unsaturated flow and transport processes in layered, fractured rocks. The development of preferential flow paths was evident in the TCw and TSw fracture continua and in the PTn matrix continuum. Simulations revealed the development of preferential pathways increasing water saturation locally. This local increase in saturation caused an increase in relative permeability for water along the pathways and reduced the wetted surface area for fracture–matrix interaction. Focusing of water flux increased this interaction locally, but much of the rock volume remained isolated from preferential flow paths. Thus, many fractures and large areas of the matrix far from the flowing fractures remained inactive despite the uniform application of water at the top boundary. Such reduction in fracture– matrix interaction reduced the ability of the matrix to attenuate water flow in unsaturated fractured rock model. Field studies of flow in unsaturated fractured chalk (Dahan et al., 1999) and in unsaturated fractured tuff (Glass et al., 2002; Salve et al., 2002) showed a similar behavior in which a small volume of rock allowed percolation of a large volume of water.

The results from our simulations suggest that the detection of preferential flow pathways will be difficult using site characterization techniques that have a large sampling volume as water flow takes place through narrow channels. A recent infiltration experiment at the ALRS (Yao et al., 2002) has shown that the neutron probe was found to be ineffective in detecting changes in water content, while tensiometers were able to detect the arrival of the wetting front in a number of monitoring stations. Electrical resistance tomography, however, as shown by Glass et al. (2002) during an infiltration experiment can detect anomalies in resistivity corresponding to the rapid migration of infiltrating water. Simulation results presented here showed an increase in the degree of focusing and decrease in the tortuosity of the flow paths as a function of variance in fracture heterogeneity. Random field simulations with larger variance values (s2Z1.5 and 2.0) showed the convergence of flow paths further reducing the fracture–matrix interaction. This is because the convergence of flow paths causes higher levels of saturation in a preferential flow channel, which reduces capillary pressure. Reduction of capillary pressure causes gravitational effects on unsaturated flow to become more important than capillary effects. Other factors such as anisotropy in fracture permeability correlation lengths, may cause convergence of flow paths and preferential pathways to be less tortuous, but were not investigated because this anisotropy was not observed in geostatistical analysis of air injection test data at the ALRS (Illman et al., 1998; Chen et al., 2000). Liu et al. (1998) developed an active fracture model in an attempt to represent the effects of preferential (rivulet or fingering type) flow in fractures. They modified the van Genuchten (1980) capillary pressure-saturation and relative permeability relationships in manner that allows only a portion of the fracture to be active. In this model, only a portion of connected fractures are active in conducting water. Inactive fractures are considered to be part of the matrix continuum when calculating liquid water flow between fractures and matrix. Parameters for the active-fracture model (particularly their gamma parameter) are inferred through inverse modeling using the site scale unsaturated zone model and fieldscale data collected at Yucca Mountain. Our random

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

field simulations show that a significant portion of the water flux can bypass a large volume of rock similar to the results obtained from the model of Liu et al. (1998), without the inclusion of parameters that represent active fractures. Calculation of water flux covariances revealed long range correlations of water flux in the vertical direction despite the short correlation length of fracture permeability. The long range correlation of fluxes implies that a uniform water flux boundary condition may not be applicable at the top of each stratum and this assumption may underestimate the distribution and quantity of water flux in deep unsaturated zones consisting of layered, fractured geologic media. The present model does not incorporate many complexities including geological features such as sloping beds, large faults, and offsets in strata. We have also neglected spatial variability in other parameters including matrix permeability and have not considered gravity instability or hysteresis in moisture-capillary pressure relationship in this numerical model. Consideration of those factors and mechanisms should be attempted in the future for an improved understanding of unsaturated flow processes in deep, stratified fractured media. Despite those simplifications, our results show that the inclusion of heterogeneity in modeling unsaturated flow processes is important and use of effective homogeneous parameter may lead to erroneous conclusions regarding the distribution of unsaturated flow through thick, fractured vadose zones. Conceptual and numerical models of unsaturated flow should consider all available site data accounting for spatial variability. In particular, our findings stress the importance of determining geostatistical parameters of fracture permeabilities through the spatial analysis of data obtained from methods such as single- and cross-hole pneumatic injection tests.

Acknowledgements The first author was supported in part by the 2002 and 2003 Old Gold Fellowships from the University of Iowa, as well as by funding from the National Science Foundation (NSF) and US Department of Defense’s Strategic Environmental Research and Development Program (SERDP). We thank Quanlin

35

Zhou and an anonymous reviewer for their constructive reviews of this manuscript.

References Ahlers, C.F., Liu, H.H., 2000. Calibrated properties model, MDLNBS-HS-000003, Lawrence Berkeley National Laboratory, Berkeley, CA. Amtec Engineering, 2001. Reference Manual for TECPLOT 9.0. Bandurruga, T.M., Bodvarsson, G.S., 1999. Calibrating hydrogeologic parameters for the 3-D site-scale unsaturated zone model of Yucca Mountain, Nevada. J. Cont. Hydrol. 38 (1–3), 25–46. Birkho¨lzer, J.T., Li, G., Tsang, C.F., Tsang, Y.W., 1999. Modeling studies and analysis of seepage into drifts at Yucca Mountain. J. Cont. Hydrol. 38 (1–3), 349–384. Bodvarsson, G.S., Bandurraga, T.M., Wu, Y.S., 1997. The SiteScale Unsaturated Zone Model of Yucca Mountain, Nevada, for the Viability Assessment. Yucca Mountain Project Milestone Report SP24UFM4, Lawrence Berkeley National Laboratory Report LBNL-40376, Berkeley, CA. Bodvarsson, G.S., Wu, Y-S., Zhang, K., 2003. Development of discrete flow paths in unsaturated fractures at Yucca Mountain. J. Cont. Hydrol. 62-63, 23–42. BSC (Bechtel SAIC), Seepage Calibration Model and Seepage Testing Data, MDL-NBS-HS-000004 REV 02. Las Vegas, Nevada: Bechtel SAIC Company, ACC: DOC.20030408.0004, 2003. Bussod, G.T.Y., Turin, H.J., Robinson, B.A., 1999. Unsaturated zone tracer testing at Busted Butte, Nevada: phase 1 results and model predictions. Geol. Soc. Am. Abstr. Programs 31 (7), A-491. Chen, G., Illman, W.A., Thompson, D.L., Vesselinov, V.V., Neuman, S.P., 2000. Geostatistical, type curve and inverse analyses of pneumatic injection tests in unsaturated fractured tuffs at the Apache Leap Research Site near Superior, Arizona,. Dynamics of Fluids in Fractured Rocks In: Faybishenko, B. et al. (Ed.), Geophysical Monograph, 122. AGU, Washington, DC, pp. 73–98. Civilian Radioactive Waste Management System Management and Operating Contractor, 1998. Chapter 2 Total System Performance Assessment-Viability Assessment (TSPA-VA) Analyses Technical Basis Document, Unsaturated Zone Hydrology Model, B00000000-01717-4301-00002 REV 00, Las Vegas, NV, Civilian Radioactive Waste Management System Management and Operating Contractor. Civilian Radioactive Waste Management System Management and Operating Contractor, 2000. Simulation of Net Infiltration for Modern and Potential Future Climates, ANL-NBS-GS-000032, Revision 00, Las Vegas, NV. Dahan, O., Nativ, R., Adar, E., Berkowitz, B., Ronen, Z., 1999. Field observation of flow in a fracture intersecting unsaturated chalk. Water Resour. Res. 35 (11), 3315–3326. Doherty, J., Brebber, L., Whyte, P., 1994. PEST: Model Independent Parameter Estimation. Watermark Computing, Brisbane, Australia.

36

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37

Doughty, C., Salve, R., Wang, J.S.Y., 2002. Liquid-release tests in unsaturated fractured welded tuffs: II. Numerical modeling. J. Hydrol. 256, 80–105. Fabryka-Martin, J.T., Turin, H.J., Wolfsberg, A.V., Brenner, D., Dixon, P.R., Musgrave, J.A., 1996. Summary Report of Chlorine36 Studies, LA-CST-TIP-96-003, YMP Milestone Report 3782M, Los Alamos National Laboratory, Los Alamos, NM. Fabryka-Martin, J.T., Flint, A.L., Sweetkind, D.S., Wolfsberg, A.V., Levy, S.S., Roemer, G.J.C., Roach, J.L., Wolfsberg, L.E., Duff, M.C., 1997. Evaluation of flow and transport models of Yucca Mountain, based on chlorine-36 and chloride studies for FY97, Los Alamos National Laboratory, Yucca Mountain Project Milestone Report SP2224M3. Faybishenko, B., Doughty, C., Steiger, M., Long, J.C.S., Wood, T.R., Jacobsen, J.S., Lore, J., Zawislanski, P.T., 2000. Conceptual model of the geometry and physics of water flow in a fractured basalt vadose zone. Water Resour. Res. 36 (12), 3499–3520. Finsterle, S., 2000. Using the continuum approach to model unsaturated flow in fractured rock. Water Resour. Res. 36 (8), 2055–2066. Finsterle, S., Ahlers, C.F., Trautz, R.C., Cook, P.J., 2003. Inverse and predictive modeling of seepage into underground openings. J. Cont. Hydrol 62-63, 89–109. Flint, A.L., Flint, L.E., Bodvarsson, G.S., Kwicklis, E.M., FabrykaMartin, J.M., 2001a. Evolution of the conceptual model of unsaturated zone hydrology at Yucca Mountain, Nevada. J. Hydrol. 247 (1–2), 1–30. Flint, A.L., Flint, L.E., Kwicklis, E.M., Bodvarsson, G.S., FabrykaMartin, J.M., 2001b. Hydrology of Yucca Mountain, Nevada. Rev. Geophys. 39 (4), 447–470. Flint, L.E., Flint, A.L., Selker, J.S., 2003. Influence of transitional volcanic strata on lateral diversion at Yucca Mountain, Nevada. Water Resour. Res. 39 (4), 1084 (doi:10.1029/2002WR001503). Glass, R.J., Nicholl, M.J., Ramirez, A.L., Daily, W.D., 2002. Liquid phase structure within an unsaturated fracture network beneath a surface infiltration event: field experiment. Water Resour. Res. 38 (10), 1199. Guzman, A.G., Geddis, A.M., Henrich, M.J., Lohrstorfer, C.F., Neuman, S.P., 1996. Summary of air permeability data from single-hole injection tests in unsaturated fractured tuffs at the apache leap research site: results of steady-state test interpretation, NUREG/CR-6360. US Nucl. Regul. Comm., Washington, DC 1996;. Haukwa, C.B., Tsang, Y.W., Wu, Y.-S., Bodvarsson, G.S., 2003. Effect of heterogeneity in fracture permeability on the potential for liquid seepage into a heated emplacement drift of the potential repository. J. Cont. Hydrol. 62–63, 509–527. Hinds, J.J., Ge, S., Fridrich, C.J., 1999. Numerical modeling of perched water under Yucca Mountain, Nevada. Ground water 37 (4), 498–504. Hu, Q., Salve, R., Stringfellow, W.T., Wang, J.S.Y., 2001. Field tracer-transport tests in unsaturated fractured tuff. J. Cont. Hydrol. 51, 1–12. Illman, W.A., Neuman, S.P., 2000. Type-curve interpretation of multi-rate single-hole pneumatic injection tests in unsaturated fractured rock. Ground water 38 (6), 899–911.

Illman, W.A., Neuman, S.P., 2001. Type-curve interpretation of a cross-hole pneumatic test in unsaturated fractured tuff. Water Resour. Res. 37 (3), 583–604. Illman, W.A., Neuman, S.P., 2003. Steady-state analyses of crosshole pneumatic injection tests in unsaturated fractured tuff. J. Hydrol. 281, 36–54. Illman, W.A., Thompson, D.L., Vesselinov, V.V., Neuman, S.P., 1998. Single-Hole and Cross-Hole Pneumatic Tests in Unsaturated Fractured Tuffs at the Apache Leap Research Site: Phenomenology, Spatial Variability, Connectivity and Scale, Rep. NUREG/CR-5559, prepared for US Nuclear Regulatory Commission, Washington, DC, September. Journel, A.G., Huijbregts, Ch.J., 1978. Mining Geostatistics. Academic Press, London. 600 pp. Klinkenberg, L.J., 1941. The permeability of porous media to liquids and gases. Am. Pet. Inst., Drill. Prod. Practice 1941;, 200–213. LeCain, G.D., 1997. Air-injection testing in vertical boreholes in welded and nonwelded tuff Yucca Mountain, Nevada, US Geological Survey Water-Resources Investigations Report 96-4262. LeCain, G.D., 1998. Results from Air-Injection and Tracer Testing in the Upper Tiva Canyon, Bow Ridge Fault, and Upper Paintbrush Contact Alcoves of the Exploratory Studies Facility, August 1994 through July 1996, Yucca Mountain, Nevada, Water-Resources Investigations Report 98-4058. Li, G., Tsang, C.-F., 2003. Seepage into drifts with mechanical degradation. J. Cont. Hydrol. 62–63, 157–172. Lichtner, P.C., Seth, M.S., Painter, S., 2000. MULTIFLO user’s manual, MULTIFLO version 1.2, two-phase nonisothermal coupled thermal-hydrologic-chemical flow simulator, center for nuclear waste regulatory analyses, San Antonio, Texas, Rev. 2 change 0. Liu, H.-H., Doughty, C., Bodvarsson, G.S., 1998. An active fracture model for unsaturated flow and transport in fractured rocks. Water Resour. Res. 34 (10), 2633–2646. Liu, H.-H., Haukwa, C.B., Ahlers, C.F., Bodvarsson, G.S., Flint, A.L., Guertal, W.B., 2003. Modeling flow and transport in unsaturated rock: an evaluation of the continuum approach. J. Cont. Hydrol. 62-63, 173–188. McLaren, R.G., Forsyth, P.A., Sudicky, E.A., VanderKwaak, J.E., Schwartz, F.W., Kessler, J.H., 2000. Flow and transport in fractured tuff at Yucca Mountain: numerical experiments on fast preferential flow mechanisms. J. Cont. Hydrol. 43 (3–4), 211–238. Montazer, P., Wilson, W.E., 1984. Conceptual hydrologic model of flow in the unsaturated zone, Yucca Mountain, Nevada. Water Resources Investigations Report 84-4355, US Geological Survey, Denver, CO. National Research Council, 2001. Conceptual Models of Flow and Transport in the Fractured Vadose Zone. National Academy Press, Washington DC.. 374 pp. Nativ, R., Adar, E., Dahan, O., Geyh, M., 1995. Water recharge and solute transport through the vadose zone of fractured chalk under desert conditions. Water Resour. Res. 31 (2), 253–262. Nicholl, M.J., Glass, R.J., 2002. Field investigation of flow processes associated with infiltration into an initially dry

W.A. Illman, D.L. Hughson / Journal of Hydrology 307 (2005) 17–37 fracture network at Fran Ridge, Yucca Mountain, Nevada: a photo essay and data summary, SAND2002-1369, 75 pp., Sandia Natl. Lab., Albuquerque, N.M. Nitao, J.J., 1998. User’s manual for the USNT Module of the NUFT code, Version 2.0 (NP-Phase, NC-Component, Thermal) UCRL-MA-1306533, 82 p. Paces, J.B., Newmark, L.A., Marshall, B.D., Whelan, J.F., Peterman, Z.E., 1996. Ages and origins of subsurface secondary minerals in the Exploratory Studies Facility (ESF), US Geological Survey Milestone Report 3GQH450M. Pruess, K., 1991. TOUGH2-A general purpose numerical simulator for multiphase fluid and heat flow, lawrence berkeley laboratory report LBL-29400 UC-251. Pruess, K., 1999. A mechanistic model for water seepage through thick unsaturated zones in fractured rocks of low matrix permeability. Water Resour. Res. 35 (4), 1039–1051. Rasmussen, T.C., Evans, D.D., Sheets, P.J., Blanford, J.H., 1990. Unsaturated Fractured Rock Characterization Methods and Data Sets at the Apache Leap Tuff Site, NUREG/CR-5596, US Nucl. Regul. Comm., Washington DC. Robin, M.J.L., Gutjahr, A.L., Sudicky, E.A., Wilson, J.L., 1993. Cross-correlated random field generation with the direct Fourier transform method. Water Resour. Res. 29 (7), 2385–2398. Rousseau, J.P., Loskot, C.L., Thamir, F., Lu, N., 1997. Results of borehole monitoring in the unsaturated zone within the main drift area of the exploratory studies facility, Yucca Mountain, Nevada, US Geological Survey, Administrative Report, Level 3 Milestone, SPH22M3. Salve, R., Oldenburg, C.M., 2001. Water flow within a fault in altered nonwelded tuff. Water Resour. Res. 37 (12), 3043–3056. Salve, R., Wang, J.S.Y., Doughty, C., 2002. Liquid-release tests in unsaturated fractured welded tuffs: I. Field investigations. J. Hydrol. 256, 60–79. Salve, R., Oldenburg, C.M., Wang, J.S.Y., 2003. Fault–matrix interactions in nonwelded tuff of the Paintbrush Group at Yucca Mountain. J. Cont. Hydrol. 62–63, 269–286. Sweetkind, D.S., Verbeek, E.R., Geslin, J.K., Moyer, T.C., 1996. Fracture character of the Paintbrush Tuff nonwelded hydrologic unit, Yucca Mountain, Nevada. USGS Administrative Report. US Geological Survey, Denver, CO. Tokunaga, T.K., Wan, J., 1997. Water film flow along fracture surfaces of porous rock. Water Resour. Res. 33 (6), 1287–1295. Tokunaga, T.K., Wan, J., Sutton, S.R., 2000. Transient film flow on rough fracture surfaces. Water Resour. Res. 36 (7), 1737–1746. Trautz, R.C., Wang, J.S.Y., 2002. Seepage into an underground opening constructed in unsaturated fractured rock under evaporative conditions. Water Resour. Res. 38 (10), 1188 (doi:10.1029/2001WR000690). Tseng, P.H., Soll, W.E., Gable, C.W., Turin, H.J., Bussod, G.Y., 2003. Modeling flow and transport at the Busted Butte field test site. J. Cont. Hydrol. 62-63, 303–318.

37

van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898. Vesselinov, V.V., Neuman, S.P., Illman, W.A., 2001. Threedimensional numerical inversion of pneumatic cross-hole tests in unsaturated fractured tuff: 2. Equivalent parameters, highresolution stochastic imaging and scale effects. Water Resour. Res. 37 (12), 3019–3042. Wang, J.S.Y., Bodvarsson, G.S., 2003. Evolution of the unsaturated zone testing at Yucca Mountain. J. Cont. Hydrol. 62-63, 337–360. Wang, J.S.Y., Cook, P.J., Trautz, R.C., Salve, R., James, A.L., Finsterle, S., Tokunaga, T.K., Solbau, R., Clyde, J., Flint, A.L., Flint, L.E., 1998. Field testing and observation of flow paths in niches: phase 1 status report of the drift seepage test and niche moisture study, level 4 milestone SPC314M4 for WBS 1.2.3.3.1.2.4, Yucca Mountain Site Characterization Project. Wang, J.S.Y., Trautz, R.C., Cook, P.J., Finsterle, S., James, A.L., Birkholzer, J., 1999. Field tests and model analyses of seepage into drift. J. Cont. Hydrol. 38 (1–3), 323–347. Wolfsberg, A.V., Campbell, K.S., Fabryka-Martin, J.T., 2000. Use of chlorine-36 and chlorine data to evaluate fracture flow and transport models at Yucca Mountain, In: Faybishenko, B. et al. (Ed.), Dynamics of Fluids in Fractured Rocks Geophysical Monograph, 122. AGU, Washington, DC, pp. 349–362. Wu, Y.S., Haukwa, C., Bodvarsson, G.S., 1999. A site-scale model for fluid and heat flow in the unsaturated zone of Yucca Mountain, Nevada. J. Cont. Hydrol. 38 (1-3), 185–215. Wu, Y.S., Zhang, W., Pan, L., Hinds, J., Bodvarsson, G.S., 2002. Modeling capillary barriers in unsaturated fractured rock. Water Resour. Res. 38 (11), 1253 (doi: 10.1029/2001WR000852). Yang, I.C., Rattray, G.W., Yu, P., 1996. Interpretations of chemical and isotopic data from boreholes in the unsaturated zone at Yucca Mountain, Nevada, US Geological Survey WaterResources Investigation Report 96-4058. US Geological Survey, Denver, CO. Yang, I.C., Yu, P., Rattray, G.W., Thorstenson, DC, 1998. Hydrogeochemical investigations and geochemical modeling in characterizing the unsaturated zone at Yucca Mountain, Nevada, US Geological Survey Water-Resources Investigation Report 98-4132. Yao, T., Mai, C.J., Wierenga, P.J., Neuman, S.P., Graham, A.R., 2002. Results of Ponding Infiltration and Tracer Experiments at the Apache Leap Research Site, Arizona, unpublished final study report prepared for US Nuclear Regulatory Commission, Washington, DC. Zhou, Q., Liu, H.H., Bodvarsson, G.S., Oldenburg, C.M., 2003. Flow and transport in unsaturated fractured rock: effects of multiscale heterogeneity of hydrogeologic properties. J. Cont. Hydrol. 60, 1–30.