# Vadose Zone Journal

- Soil Science Society of America

## Abstract

Simulations of flow and transport in variably saturated fractured rock generally assume equilibrium conditions between the fractures and the porous matrix, leading to predictions that are dominated by a diffusive process. Contrary to these predictions, an increasing body of evidence suggests that fracture-dominated flow, under nonequilibrium conditions between the fractures and porous matrix, occurs frequently in field and laboratory settings. Flow processes, such as fluid cascades and flow path switching, are often observed in laboratory experiments, but are generally not captured by diffusion-based conceptual and numerical models. Many of these processes are assumed to be averaged out at some representative elemental volume scale; however, anecdotal evidence from field experiments conducted at various scales of investigation suggest that this may not be the case. Comparison of experimental observations with numerical simulations illustrates at least two potential problems with standard equivalent continuum and discrete fracture conceptual models of unsaturated fractured and porous media flow. First, such models tend to overestimate the strength of interaction between the fracture and matrix domains. Second, model calibration may allow diffusion-based models to accurately reproduce experimental observations without providing a complete description of the physics governing the system. Failure to incorporate convective transport, reduced fracture–matrix interaction, and other sub-grid-scale processes in models of flow in fractured porous media may result in erroneous descriptions of system behavior.

- aniso-ECM, equivalent continuum model with anisotropic properties
- DFM, discrete fracture model
- ECM, equivalent continuum model
- INEEL, Idaho National Engineering and Environmental Laboratory
- iso-ECM, equivalent continuum model with isotropic properties

Flow and transport in unsaturated fractured rock has received considerable attention in the literature in the last decade, largely due to the need to better understand the long-term risk of radioactive and hazardous waste at western U.S. sites situated in areas with thick, fractured-rock vadose zones. Continuum-based conceptual models commonly used to describe flow and transport in these systems predict matrix-dominated flow, smooth wetting fronts, and idealized contaminant breakthrough curves (e.g., Wang and Narasimhan, 1985); however, system behavior observed in the field is often more complex (Fabryka-Martin et al., 1996; Flint et al., 2001; Finsterle et al., 2002; Yang et al., 1996; Faybishenko et al., 2001; Glass et al., 2002b). Field experiments conducted at the Idaho National Engineering and Environmental Laboratory (INEEL) have exhibited complex dynamical behaviors, such as flow path switching, sensitive dependence on initial conditions, failure to achieve steady state after long time periods (up to 20 d), and order of magnitude variations in infiltration rates under constant head boundary conditions. These behaviors are difficult to reproduce with continuum-based conceptual models and have been observed at a variety of scales, from about 1 m at the Hell's Half Acre site (Podgorney et al., 2000; Wood et al., 2000), through intermediate (≈10 m) at the Analog Site for Fractured Rock Characterization (Faybishenko et al., 1998, 2000), to a scale of approximately 100 m at the Large Scale Pumping and Infiltration Test site (Wood and Norrell, 1996). In analyzing these field experiments, Faybishenko et al. (2001) concluded that the fractured basalt vadose zone must be conceptualized at a hierarchy of scales, with no clear scaling relationships evident between different scales of observation.

Laboratory experiments conducted at the single fracture (e.g., Su et al., 1999; Glass and Nicholl, 1996; Nicholl et al., 1994) and fracture network scales (e.g., Wood et al., 2004; Glass et al., 2002a; LaViolette et al., 2003) have also shown behavior not well described with a continuum-based conceptual model. Fracture-dominated flow and dynamical behavior were predicted to exist on the basis of a consideration of capillary and gravitational forces interacting in a fracture network as early as 1995 (Glass et al., 1995). However, phenomena such as convergent flow with depth, flow path switching, and fluid cascades have only recently been clearly identified at the laboratory scale (Wood et al, 2004; LaViolette et al., 2003). The importance, or possibly even the existence, of these behaviors cannot be fully quantified at the field scale, although it may be inferred to be important from field tests at the INEEL and other sites. Incorporation of these behaviors into flow and transport models is a challenging task, but one that may be necessary to provide a complete representation of the behavior in natural fractured rock vadose zones (National Research Council, 2001).

The application of continuum-based numerical models to unsaturated flow in fractured rock systems relies on the assumption that computational grid blocks are sufficiently large to average out noncontinuum behavior (Mantoglou and Gelhar, 1987a, 1987b, 1987c). Numerical studies have shown that a properly calibrated continuum model can produce reasonable results even when the flow processes are discrete (Finsterle, 2000). When considering large-scale flow problems, however, it has been pointed out that volume averaging may be meaningless if the important components of flow occur as localized, preferential fracture flowpaths (Pruess, 1999). This was observed by Podgorney and Fairley (2003), who modeled flow in an experiment containing a single fracture embedded in a porous matrix. The investigators found that a continuum approach worked well at sufficiently low flow rates, but failed to capture the dynamics of the system at higher rates corresponding to fracture-dominated flow conditions. Kwicklis and Healy (1993) simulated a hypothetical fracture network embedded in a 5 by 5 m flow region, with the goal of examining the interaction between the fracture and matrix domains and to identify conditions where equilibrium between them exists. They were able to calculate equivalent permeability for the fracture network and illustrate how information from single fractures may be used to estimate flow properties at larger scales, but they expected flow behavior in a natural fracture network to be much more complicated. Liu et al. (2002) used a randomly generated fracture network containing thousands of fractures in a 10 by 10 m flow domain to examine the unsaturated flow of water around a drift. The authors stated that the current lack of understanding of processes occurring within a fracture network is largely due to the technical difficulties in making direct observations. While several examples of simulations for simple, single fracture experiments and idealized fracture networks can be found, very few examples can be found of models attempting to simulate data collected under controlled and highly instrumented experimental conditions. We suspect that discrepancies found in the literature between models and observations are the result of sub-grid-scale processes, which are not treated in the simulations. A number of authors (e.g., Liu et al., 1998) have suggested ad hoc corrections to numerical simulators as a means of overcoming such discrepancies; however, such approaches have not been tested against detailed experimental observations.

In an effort to further our understanding of how existing models can evolve to make better predictions of fluid flow in fracture networks, we present a comparison of flow in a well-characterized fracture network to behavior predicted using standard, continuum-based modeling approaches. Examination of the observed vs. predicted behaviors can provide insight into the physical processes captured by the continuum approach and identify noncontinuum processes that can exhibit control on the flow through a fractured rock network. In addition, coupling field observations from various scales of investigation with those from the laboratory experiments, in light of the discrepancies encountered in the numerical models, provides some insight into the dominant processes occurring within fracture networks across multiple scales of investigation.

We presented a set of physical experiments considering infiltration into a well-characterized, two-dimensional fracture–matrix system in a companion paper (Wood et al., 2004). In Part 2, we use that data set to demonstrate that employing a continuum-based approach and using standard simplifying assumptions can lead to erroneous model predictions at this scale.

## EXPERIMENTAL SET-UP

As presented in Part 1 of this paper (Wood et al., 2004), an analog fracture–matrix system was constructed in the laboratory using 12 limestone blocks stacked four blocks wide by three blocks high. Tap water equilibrated with the limestone material was supplied to the top of the fractures at a rate of 1 mL min^{−1} fracture^{−1}. Both quantitative and qualitative data were collected to compare with numerical simulations, including water inflow rate, water outflow rate, and breakthrough time of water at the bottom boundary. Time-lapse photography was used to document the evolution of the wetting front throughout the fracture network. Water exiting the system was conveyed to measurement bottles to provide spatial and temporal information on the outflow from the network.

Because the objective of the study was to test the ability of common numerical simulators to replicate the behavior of a well-characterized system, considerable effort was made to obtain accurate measurements of material properties wherever possible. Matrix properties of porosity and permeability were determined using standard laboratory characterization methods and testing in a UFA ultra-centrifuge (UFA Ventures, Inc., Richland, WA), with the saturation-dependent properties of the van Genuchten (1980) model determined by fitting the pressure–saturation data using the RETC computer code (van Genuchten et al., 1991). It is standard practice to use unsaturated properties based on drying curves, and this convention was followed by the present study to allow comparison with the results of other researchers. Repeated testing established the matrix properties with a high degree of precision, and demonstrated the limestone blocks to be surprisingly homogeneous in character.

Fracture aperture between the blocks was measured at intervals of approximately 2.5 cm in both the vertical and horizontal fractures. The measured apertures were variable with values ranging from 3.0 × 10^{−5} to 6.3 × 10^{−4} m for vertical fractures and 2.5 × 10^{−5} to 3.2 × 10^{−3} m for horizontal fractures. Figure 1 shows a histogram of the fracture apertures and summary statistics for Tests 1 through 4. All of the aperture distributions were positively skewed, with the smallest apertures dominating the aperture distributions. Fracture properties were calculated using the mean value of the vertical and horizontal apertures.

## NUMERICAL MODELING

The numerical simulations of flow used three conceptual models: (i) an equivalent continuum model (ECM) with isotropic properties (iso-ECM), (ii) an ECM with anisotropic properties (aniso-ECM), and (iii) a discrete fracture model (DFM) with isotropic matrix properties and anisotropic fracture properties. Hydraulic property sets for the three models were estimated for the experimental apparatus using commonly accepted techniques. Matrix properties determined from laboratory testing were used without modification; however, fracture properties are notoriously difficult to characterize through laboratory testing, and were therefore assumed equivalent to parallel plates, with apertures equal to the mean of the measured apertures. Permeability was calculated using the well-known solution to the Navier–Stokes equations (cubic law solution), with the fractures represented numerically in the DFM simulations as porous material with a porosity of nearly unity.

The van Genuchten (1980) model was used to describe the unsaturated capillary pressure, and is given by1where *m*, *n*, and *P*_{o} are fitting parameters. *S** is defined by2and *m* is given as3In the above, *S*_{w} is the water phase saturation, *S*_{r} is the residual (irreducible) water phase saturation, and *S*_{s} is the water phase saturation at zero capillary pressure. *P*_{o} is defined as4where the α parameter, defined as one over the pressure differential at which air will displace water (the air-entry pressure), was estimated for the fractures from capillary theory as the interfacial tension divided by one-half the mean fracture aperture (assuming a perfectly water-wet matrix). No accepted method exists for making a priori estimations of the van Genuchten *m* parameter for a fracture; however, Reitsma and Kueper (1994) reported an experimentally determined value of approximately of *m* = 0.8 for a natural limestone fracture. Therefore we chose an initial value of 0.8 for the simulations, and sensitivity studies were conducted for a range of values between 0.2 and 0.9.

Properties estimated as described above were used directly in the DFM simulations. Parameter adjustments for the DFM simulations were largely confined to fracture properties, since the level of uncertainty for matrix properties was considered small in comparison. For the ECM simulations, single-continuum properties were calculated based on arithmetic and harmonic averaging of fracture and matrix domain properties. The small volume occupied by fractures resulted in a single continuum property set dominated by matrix permeability in the iso-ECM; for the aniso-ECM, vertical and horizontal permeabilities differed by a factor of slightly more than 3. Because matrix properties dominated hydraulic properties regardless of the averaging scheme used, matrix values of van Genuchten *m* and α were used directly for both single continuum models. Selected values of parameters used in the simulations are shown in Table 1.

The laboratory experiments were simulated with the Richards equation (i.e., multiphase diffusive flow, where the gas phase is considered “passive”) using the numerical code TOUGH2 (Pruess et al., 1999a). The simulations neglected the effects of hysteresis, which is a common practice in vadose zone modeling.

## RESULTS AND DISCUSSION

Most standard hydrogeology textbooks (e.g., Bouwer, 1978; Freeze and Cherry, 1979) discuss the use of harmonic and arithmetic averaging for estimating effective hydraulic conductivity in layered heterogeneous systems. Strictly speaking, effective properties derived using these methods are valid for steady-state, one-dimensional flow either parallel or perpendicular to layering. For two-dimensional steady flow in a heterogeneous medium with a lognormal conductivity distribution, the geometric mean of local conductivities provides a good estimate of effective conductivity. More generally, it has been shown that the effective conductivity of a heterogeneous medium is bounded by the arithmetic and harmonic means of local conductivities, regardless of dimensionality or the statistical distribution of conductivity for the case of uniform flow. Unfortunately, no method currently exists for calculating effective conductivities in a nonuniform flow field (Matheron, 1967; de Marsily, 1986). Despite these well-known limitations on methods for estimating effective hydraulic properties, lack of viable alternatives often dictates their use by modelers working in fractured, or otherwise heterogeneous, porous media systems.

Even in ubiquitously fractured rock, matrix volume far outweighs that occupied by fractures; this is reflected in the effective property estimates for the ECM (Table 1), which are largely or wholly representative of matrix properties. Although the ECM simulation results plotted in Fig. 2b do not demonstrate any of the key behaviors of the laboratory experiments (i.e., flow path focusing or wetting front in the fractures significantly leading the wetting front in the matrix), the model predictions conform to expectations for a system in which matrix and fractures are in capillary equilibrium. Given conditions of fracture–matrix capillary equilibrium, the relatively low capillary suction of the fractures causes matrix properties to dominate system response, as predicted by the conceptual model of Wang and Narasimhan (1985). Fracture–matrix capillary equilibrium is a reasonable assumption for low, constant or slowly changing fluxes. Conversely, it has been shown that the assumptions underlying the ECM representation of fractured porous media are violated by fluxes significantly in excess of a “specific flux” (m^{3} s m^{−1} of fracture length) that is a function of matrix porosity, saturation, diffusivity, and the time scale of response of the matrix (Nitao and Buscheck, 1991; Nitao et al., 1993). The calculated specific critical flux for the experimental domain is on the order of 3 mL min^{−1} fracture^{−1}, a value that is slightly larger than that used in the laboratory experiment, yet flow in the fractures can still be expected to dominate the system (Podgorney and Fairley, 2003).

Improved experiment–simulation agreement might be expected for the DFM simulations, because the DFM does not assume capillary equilibrium between the fracture and matrix domains. Although the conceptual model underlying the DFM provides a better description of the experimental system, a number of difficulties complicate its application in a practical setting. The use of the parallel plate solution of the Navier–Stokes equations to approximate fracture permeability is common practice among modelers, in spite of the fact that this model has been shown to be generally inapplicable to natural fractures (Zimmerman and Bodvarsson, 1997). Although the more general Reynolds Lubrication equation may be theoretically valid for some natural fractures, it is more cumbersome than the parallel plate model, and is therefore rarely used. Additionally, the narrow limits on hydraulic gradient make its application to natural systems questionable (Zimmerman and Bodvarsson, 1997).

Although estimation of fracture permeability is a familiar problem in the development of discrete fracture models, sensitivity studies in the DFM simulations showed that the amount of interaction between fracture and matrix domains was by far the most important factor for obtaining agreement between simulation results and experimental observations. It is clear that the uncalibrated DFM simulations significantly overestimated the strength of coupling between the two domains. This resulted in poor agreement between observed and predicted moisture distributions, wetting front behavior (Fig. 2c), breakthrough time, and cumulative flux at the lower boundary of the model, and the predicted values differed only slightly from ECM predictions (Fig. 3) .

The poor representation of the experimental system offered by the DFM is not surprising for a number of reasons. Numerical models of flow in fractured rock appear to overestimate fracture–matrix interaction, generating unrealistic representations of diffuse wetting fronts, matrix-dominated flow regimes, and unreliable travel time predictions. Liu et al. (1998) suggested reducing the numerical contact area between the fracture and matrix domains by a factor proportional to fracture saturation raised to a power. Observations of channelized flow in laboratory experiments support the hypothesis that the fracture–matrix interaction area is a relatively small fraction of the geometric area under unsaturated conditions (Su et al., 1999). Unfortunately, the Liu et al. (1998) model requires a parameter that is system specific, and can only be obtained by means of inverse modeling.

In addition to the interaction area, the ratio of air-entry pressures for the fracture and matrix domains is an important factor in determining fracture–matrix interaction in DFM simulations. This is illustrated in Fig. 2c, which shows the results of simulations in which matrix capillary suction was severely limited. These simulations showed improved agreement with experimental observations, but at the expense of a physically unrealistic parameter set for the porous matrix. The difficulty of accurately representing capillary pressure in a fractured medium is increased by uncertainty regarding appropriate models relating saturation and capillary pressure in natural fractures. A number of well-known and generally accepted models relating moisture content with capillary pressure and relative permeability are available for soils and other porous media. It is unclear, however, to what extent these relationships are applicable to flow in fractures. In our study, this uncertainty is compounded because the saturation–capillary pressure relationship is characterized using an “effective” fracture air-entry pressure that overly simplifies a complex situation and contributes to poor model performance.

As mentioned above, no accepted method exists for making a priori estimations of the van Genuchten (1980) *m* parameter, although reports in the literature for a natural limestone fracture with a hydraulic aperture ranging from 5.9 × 10^{−5} to 1.3 × 10^{−4} m give a value of approximately of 0.8 (Reitsma and Kueper, 1994). An initial value of 0.8 was therefore chosen for the simulations, and sensitivity studies were conducted for a range of values between 0.2 and 0.9. As shown on Fig. 4 , as the *m* value is decreased from 0.8 to 0.2, the time required to reach a pseudo-steady-state condition in the simulation is increased. The lower *m* values tended to match the early time better while higher values did a better job overall matching the experimental cumulative flux out of the system. A “calibrated” DFM model (see Fig. 2 and 3) was obtained by limiting the maximum capillary pressure in the matrix to approximately 1.1 × 10^{3} Pa and using a fracture *m* value of 0.8.

Many of the problems discussed in the previous paragraphs are assumed to be minimized or overcome at some representative elemental volume scale; however, comparison between observations and simulations indicate there are advective processes at work in the experiment that fall outside the range of behaviors reproducible with a diffusive model. Similarly, observations from large-scale field experiments, such as intermittent detection of water and contaminants at vadose zone monitoring wells during ponded head infiltration tests (Wood and Norrell, 1996), may be attributable to the types of flow behaviors observed in the laboratory experiments.

It is important to point out that the discrepancies between the numerical modeling results and the experimental data in no way reflect poorly on the code chosen. Any diffusion-based numerical simulation would have produced similar results and produced the same discrepancies. There are numerous codes available for modeling unsaturated systems (e.g., FEHM, Zyvoloski et al., 1995; NUFT, Nitao, 1996; TETRAD, Vinsome and Shook, 1993). We chose the TOUGH2 code (Pruess et al., 1999a) because of personal preference and our familiarity with the code.

## SUMMARY AND CONCLUSIONS

We applied common continuum-based numerical modeling approaches to experimental results of flow through a surrogate fracture network within a porous matrix. Having a well-characterized experimental dataset allowed for the formulation of the numerical simulations based on measured property values for the porous matrix and calculated values for the fractures. Three basic conceptual models were formulated and applied to the experimental data: isotropic ECM, anisotropic ECM, and DFM. The initial, or uncalibrated, results from all of the numerical models compared poorly against the experimental data, with the numerical results showing considerably more fracture matrix interaction and significantly later arrival of the wetting front at the bottom boundary of the experiment. Calibration of the DFM model to match the experimental data resulted in an unrealistic property set for the matrix domain.

The discrepancies between the observed and modeled behavior may be attributed to the fact that some of the processes observed in the experiments are not captured by the continuum-based approach. Dynamical behaviors such as flow path switching and fluid cascades are advective processes and are important mechanisms in the experimental dataset. The question that remains is whether these advective processes exist and are important at the field scale. Anecdotal evidence from field experiments at various scales and locations suggests that this may indeed be the case. Incorporation of these processes into conceptual and numerical models may lead to better predictions of water flow and contaminant transport through unsaturated fractured vadose zones.

## Acknowledgments

The authors would like to thank M.J. Nicholl, R.J. Glass, G.M. Shook, R.J. Lenhard, and two anonymous reviewers for their many helpful comments and suggestions. This work was supported by the Environmental Systems Research Program under contract number DE-AC07-99ID13727 from the Office of Environmental Management, Department of Energy to the Idaho National Engineering and Environmental Laboratory.

- Received May 10, 2003.