Conceptual Models of Flow and Transport in the Fractured Vadose Zone
A significant number of subsurface environmental problems involve fluid flow and solute transport in the fractured vadose zone. In this report, the vadose zone refers to that part of the subsurface from land surface to the lowest seasonal water table elevation. The vadose zone may be composed of consolidated rock and/or unconsolidated granular material, including soils. Contamination of the vadose zone can result from many sources, including chemical spills, leaky underground storage tanks, leachate from waste disposal sites and mine tailings, and application of agricultural chemicals. Another major environmental concern is the potential for long-term migration of radionuclides from low- and high-level nuclear waste disposal facilities. The presence of fractures and other channel-like openings in the vadose zone poses a particularly significant problem, because such features are potential avenues for rapid transport of chemicals from contamination sources to the water table.
A key component in assessing contamination hazards and designing remedial actions is the development of flow and transport models to approximately represent real systems. Present-day models often use computer programs to simulate flow and transport processes, and such models are applied throughout the scientific, engineering, and regulatory arena. When carefully developed and supported by field data, models can be effective tools for understanding complex phenomena and for making informed predictions for a variety of assumed future scenarios. However, model results are always subject to some degree of uncertainty due to limitations in field data and incomplete knowledge of natural processes.
Thus, when models form the basis for decision-making, uncertainty will be an inescapable component of environmental management and regulation. A key consideration in the modeling process is whether or not the model has undergone sufficient development and testing to address the problem being analyzed in a sufficiently meaningful manner.
The underpinning of any vadose zone fluid transport model is the conceptualization of (1) the relevant processes, (2) the structure of the subsurface, and (3) the potential events or scenarios that impact the behavior of the modeled system. These conceptualizations together form a “conceptual model,” and it is such conceptual models that are the focus of this Panel's study. The evolution of a conceptual model, from the initial formulation through the subsequent revisions and refinements, is the crux of the model development process. Conceptual models of fully saturated flow in granular aquifers have progressed to a relatively mature state due to a long history of model development and application. By contrast, conceptual models of partially saturated flow and transport in fractured vadose zone environments are poorly developed and largely untested. In cases where multiple competing conceptual models could lead to drastically different conclusions, strategies for evaluating these models must be based on sound technical criteria. The need to develop such strategies and criteria was a key reason for the appointment of this Panel.
The purpose of this study is to investigate the processes through which conceptual models of flow and transport in the fractured vadose zone are developed, tested, refined, and reviewed. The Panel convened a two-day workshop in March 1999, during which a large group of specialists (see Appendix A) from the hydrogeologic, geochemical, soil science, and related fields discussed the current state of knowledge, lessons learned from field investigations, and needs for future research. A series of invited presentations provided a basis for much of the discussion at the workshop. Individually authored papers based on these presentations are presented as Chapter 2, Chapter 3, Chapter 4, Chapter 5, Chapter 6, Chapter 7, Chapter 8, Chapter 9, Chapter 10 and Chapter 11.
The Panel was charged with preparing a consensus report on the development and testing of conceptual models for fluid flow and transport in the fractured vadose zone. The Panel's conclusions and recommendations were based in large part on the workshop presentations and discussions. The report is intended to:
provide information on contemporary philosophies, approaches, and techniques for conceptual model building;
provide guidance to regulatory agencies on the review process for conceptual models developed for site licensing;
bring together knowledge and experiences from related disciplines so that technical communities can benefit from advances in related fields; and
identify future research needed to further the technical basis for developing and evaluating conceptual models of flow and transport in the fractured vadose zone.
The Panel devoted a major portion of its study to an analysis of fluid flow models in the fractured vadose zone because (1) understanding flow is a prerequisite for understanding transport, (2) questions regarding the nature of fast pathways are a major concern, and (3) recent studies, especially investigations at the potential site for a nuclear waste repository at Yucca Mountain, have accumulated a significant body of knowledge applicable to fractured vadose zone models. In its consideration of transport in the vadose zone, the Panel focused on the application of environmental tracers (such as tritium and chlorine-36) because they provide integrated responses that are difficult to determine by point measurements of fluid potential or moisture content. The Panel briefly reviewed approaches for modeling transport of conservative solutes, but the scope of the study did not include reactive solutes or water-rock interactions.
The Panel report is composed of three main sections. First, we discuss general considerations applicable to the development and testing of conceptual models. Second, we summarize the current state of knowledge of flow and transport processes in the fractured vadose zone. Third, we present our conclusions and recommendations. Appended to the Panel report are the invited papers based on presentations at the workshop. These papers served as the starting point for workshop discussions, and form much of the background material used in preparation of this report. Although the papers have undergone peer review, their inclusion in this report does not constitute any specific endorsement of their contents, either by the Panel or by the National Research Council.
DEVELOPMENT AND TESTING OF CONCEPTUAL MODELS
Definition of Conceptual Model
Models representing natural systems are often viewed as composed of two components, a conceptual model and a mathematical model. In general terms, a conceptual model is qualitative and expressed by ideas, words, and figures. A mathematical model is quantitative and expressed as mathematical equations. The two are closely related. In essence, the mathematical model results from translating the conceptual model into a well-posed mathematical problem that can be solved.
Various definitions of conceptual models can be found in the scientific and technical literature. These definitions are generally consistent in their fundamental meaning, and differ mainly in scope, detail, and context. The statement of the conceptual model often reflects the key questions and unknowns to be investigated. For example, Anderson and Woessner (1992) give the following definition for the purpose of modeling groundwater flow in aquifers:
A conceptual model is a pictorial representation of the groundwater flow system, frequently in the form of a block diagram or a cross section.
This relatively brief statement is suitable for modeling groundwater flow because the process is well understood, and the principal concern is the delineation of subsurface structure, such as stratigraphy, and the spatial distribution of hydraulic properties, such as hydraulic conductivity. By contrast, Tsang (1991) gives the following definition in the context of modeling a broad range of coupled physical-chemical processes:
A site-specific conceptual model consists of three main components: structure, processes, and boundary and “initial” conditions. “Structure” refers to the geometric structure of the system, such as stratigraphy, faults, heterogeneity, fracture density and lengths, and other geometric and geologic characteristics. “Processes” are physical and chemical phenomena such as buoyancy flow, colloidal transport, matrix diffusion, and dissolution and precipitation. “Boundary conditions” are constant or time-dependent conditions imposed on the boundaries of the model domain. “Initial conditions” are the physical and chemical conditions over the model domain at a particular instant of time. This is usually taken at the initial instant of time, though in general it can be any specified point in time.
The more detailed statement reflects the complexity of the modeled processes, some of which may be poorly understood.
In the present study, we define a conceptual model as follows: A conceptual model is 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. This definition stresses several ideas. A conceptual model is a hypothesis because it must be tested for internal consistency and for its ability to represent the real system in a meaningful way. The hypothesis evolves (is revised and refined) during testing and as new information is gathered. Although a conceptual model is by necessity a simplification of the real system, the degree of simplification (or conversely, the amount of complexity retained in the model) should be commensurate with the problem being addressed. In order to present a clear and easily understandable definition, we have used the term ‘events' rather than the mathematically explicit terms ‘initial conditions' and ‘boundary conditions' that have been used in earlier definitions. The context in which the model is developed constrains the range of applicability of the model.
The Modeling Process
We refer to the modeling process as an iterative sequence of actions that includes (1) identifying a site-specific problem; (2) conceptualizing important features, processes, and events; (3) implementing a quantitative description; (4) collecting and assimilating field data that are used to calibrate the model and evaluate its predictive capabilities; and (5) developing predictions that are used to resolve the identified problem. This process is illustrated by the flow chart in Figure 1-1.
The modeling process usually begins with some initial perception of a field problem, expressed as questions to be answered or decisions to be made (as indicated in the upper left of Figure 1-1). Available site-specific data, related experience, and generic scientific knowledge are then combined to identify the factors (features, processes, and events) that are important for the identified problem. This typically involves preliminary calculations using generic parameters to understand the relative importance of various aspects of the problem. The result of this assessment is the initial conceptual model.
The formulation of this initial conceptual model is arguably the most important step in the modeling process. The conceptual model is the foundation of the mathematical model, and strongly influences the type of computer code to be used and the design and priority of site characterization activities. Philips (this report, Chapter 9) describes a case study where processes that were neglected in the initial conceptual model were later found to be of fundamental importance, leading to a substantial revision of the conceptual model. Such revisions, possibly with additional iterations, may be expected in complex environments where often poorly understood physical, chemical, and biological processes interact. Therefore, the conceptual modeling process should consider a broad range of reasonable alternative hypotheses. It is also important to employ a variety of different types of data. For example, an initial conceptual model that is based on regional water budgets, general hydraulic properties, and sparse environmental tracer data is more likely to include relevant and significant processes than if only general hydraulic properties were considered. The conceptual model is the foundation upon which all aspects of the modeling process are constructed, and this should be clearly understood by all members of a project team, from data gatherers to modelers to managers. In any situation where a particular field problem requires a linkage between site characterization and performance assessment, the development of a conceptual model and the continued testing and assessment of the model should be an explicitly required component of efforts to address the problem.
The next step in the modeling process is the development of a quantitative description of the conceptual model in terms of mathematical equations to be solved—the mathematical model. Except for simple problems, mathematical equations are usually solved by use of a computer code, also called a numerical model. A newly developed computer code must be tested to assure that it correctly solves the equations of the mathematical model. This “verification” process typically involves comparisons of the numerical results with known solutions for related simple configurations. Because the verification process of a complex code usually involves piecemeal testing of the individual transport processes represented by the model, without simultaneous testing of all or most features of a code, it is rarely possible to be assured that a code is fully verified. In addition, even a well-established computer code should be used carefully by knowledgeable analysts to avoid potential problems such as numerical instability.
A mathematical model always includes a number of parameters, such as the distribution of hydraulic properties within the modeled region. At the early stages of the modeling process, the values of some of these parameters are unknown. The process of estimating these parameters is known as model calibration or parameter estimation. The calibration process will typically include sensitivity analyses, designed to elucidate how changes in parameters can influence the simulation results. Calibration procedures based on statistical/stochastic formulations can also provide quantitative measures of uncertainty in the parameter estimates. Such information is extremely valuable for guiding the collection of additional field data to refine the model. As shown in Figure 1-1, the calibration step can involve significant feedback to the conceptual model.
Model calibration is often achieved by adjusting model parameters, either manually or by automated methods, until the model simulated results agree with field data to an acceptable level. As the model is a simplification of the real system, a perfect fit between the simulated results and field data is unlikely to be achieved. A serious lack of fit indicates that the conceptual model should be reexamined. On the other hand, a good fit does not necessarily prove that the conceptual model is adequate to address the issues in question. Another complicating factor may be that a measured parameter may not be the same as a modeled parameter, e.g., it is difficult to determine what the measured moisture tension at a porous tip actually represents for a situation where the moisture originates as film-flow in a fracture. Additional model testing, which creates another feedback loop to the conceptual model element in Figure 1-1, is essential to gain confidence that the model provides meaningful results that can be used for decision-making and problem resolution. Model testing is discussed in further detail in the next section.
It is important to recognize that model predictions require assumptions about future events or scenarios. The importance of this “event” component of the conceptual model is made clear by postaudit studies of groundwater models (e.g., Anderson and Woessner, 1992, p. 288-293). These studies suggest that, when model predictions were compared to actual outcomes, a major reason for inaccurate prediction was that the assumed future scenarios (e.g., future stresses on the system) did not occur. This finding points to the need for developing model predictions for many different possible future scenarios.
Another important aspect of model prediction is uncertainty, which arises from many factors, such as simplification of the real system, limited field data, measurement errors, and multiple conceptualizations or interpretations that cannot be resolved by existing data. Even if the modeling process has undergone several cycles of revision, testing, and gathering of additional field data, prediction uncertainty is reduced but not entirely eliminated. Decisions based on model predictions must take this into account by approaches such as analysis of risks, worst-case scenario, or incorporation of safety margins. If the conceptual model
uncertainty is too great, the project goals may not be achievable within a reasonable time or at acceptable cost. Alternatively, it may be necessary to redefine the original problem to be resolved.
The formulation of a conceptual model is inherently subjective in that it relies on limited available site-specific observations and data, as well as experience and insights developed through work on similar sites and/or related problems. A conceptualization is susceptible to biases arising from the disciplinary background and experience of the analyst, and/or by differing perceptions of the problem as influenced by external social and political forces. Although model calibration represents one level of testing by requiring the model to reproduce one set of field data, the model may be applied to field conditions that are significantly different than the conditions under which calibration data were collected. For these reasons, it is important to test the predictive capabilities of a model.
Due to the large variability in the objectives of modeling projects, it is not possible to provide a prescriptive step-by-step procedure for model testing. The amount of effort devoted to model testing will strongly depend on the question to be resolved, the scope of the investigation, available resources, and the consequence of an inadequate or inappropriate model. In this section, we discuss some general issues that should be considered in conceptual model testing. We assume at the outset that the computer code used for simulation is already verified, in other words, that the program logic and numerical algorithm are correctly implemented, and the results are free of computational errors.
A traditional procedure of model testing is to use the calibrated model to simulate a set of field data that was not used during the calibration process. For example, if historical data exist for the site, a portion of the data may be used for calibration and the remaining data used for testing. Alternatively, the model is used to predict the result of a field experiment, and the experiment is then carried out as a check against model predictions. The test is successful if the model simulation agrees with the test data, within reasonable limits, without the need to further adjust the model parameters. A model that passes one or a series of such tests would have demonstrated a certain level of predictive capability. If the model fails the test, then the test data are used to further calibrate the model and additional field data are needed for a new round of testing.
Although the above test procedure is conceptually appearing, it may not always be feasible. In most cases, field data are limited and all the data must be used for calibration, leaving none for testing. Therefore, a broader view of model testing is often necessary. One approach is to evaluate how well the conceptual model represents the real system in terms of features, processes, and events. This
involves both examining the model assumptions and evaluating alternative hypotheses. The underlying rationale is that a model that has undergone such evaluations can be used with an increased degree of confidence.
As an example of evaluating a “feature” component of the conceptual model, consider the question of whether a geologic stratum that is represented in the model as a continuous layer underlying the entire model region may in fact be discontinuous and absent at certain locations. The first step in evaluating this alternative hypothesis is to consider its consequence. This can be done by modifying the model according to the alternative hypothesis (discontinuous layer), recalibrating the model, and developing new predictions. If the modified model with changed parameters cannot match the data, this is an indication that the alternative hypothesis is inconsistent with available field data. If recalibration is successful but the new predictions are similar to the old predictions, this is an indication that the alternative hypothesis is of little consequence (within the context of the question to be resolved). In both cases, the alternative hypothesis may not warrant additional consideration. However, if the recalibrated model leads to new predictions that are significantly different from the old predictions, then further investigations are needed. The continuity of the layer in question may be evaluated by geophysical methods to image the subsurface, by test borings at strategic locations, or by using an understanding of the depositional environment to infer the likelihood that the layer may be continuous or discontinuous.
Whether or not a model can be “validated” is a topic that has seen substantial debate in recent years. Certain authors (e.g., Konikow and Bredehoeft, 1992; Oreskes et al., 1994) have presented the opinion that models cannot be validated, because the truth of scientific theory can never be proven. Describing a model as “validated” implies, especially to a nontechnical audience, a level of correctness and certainty that is unattainable. Other authors (e.g., see discussion by Jarvis and Larsson, this report, Chapter 6)argue that a “validated model” can be taken to mean that a model is acceptable for its intended use. Currently, there is a lack of consensus on the definition of “validation.”
From an operational point of view, model testing and evaluation can be viewed as activities designed to establish the credibility of a model. In this regard, peer review is an important part of the modeling process. Review by a group of objective, independent, and respected experts can utilize knowledge and opinions beyond those of the study team. Maintaining a free flow of information (e.g., field data, model results) can also add a significant measure of credibility, because the model is open to scrutiny by concerned parties (e.g., regulators, license applicants, government agencies, public citizen groups, and the general scientific community). In the final analysis, whether or not a model is acceptable to concerned parties depends on whether or not these parties have confidence that the model predictions provide meaningful input for decisionmaking and problem resolution.
FLOW AND TRANSPORT IN THE FRACTURED VADOSE ZONE
Nature of Problem
A major issue in the investigation of the vadose zone is estimating fluid flux and solute travel time from land surface to the water table. Field observations suggest that infiltrating water and solute may not necessarily advance downward as a uniform infiltration front, as would be simulated by a model of infiltration into a homogeneous medium. Instead, fluid and solute may travel at widely different fluxes and velocities through different parts of the vadose zone. Field evidence for uneven water movement includes the detection of environmental tracers indicating that relatively young water has moved to significant depths in arid regions, where infiltration rates are expected to be small. For example, at the Exploratory Studies Facility, Yucca Mountain, Nevada, water with a bomb-pulse chlorine-36 signature (that is, water less than 50 years old) was found several hundred meters below the land surface (Fabryka-Martin et al., 1998; see also discussion by Phillips, this report, Chapter 9), whereas a simple model of infiltration through a homogeneous vadose zone would predict travel times of hundreds to thousands of years to reach such a depth.
In this study, we focus on five major issues that cause difficulties for estimating fluid flux and travel times through the fractured vadose zone. First, flow in unsaturated fractures may occur as either capillary flow or film flow. Recent studies have also reported intermittent flow behaviors that are not considered by classical theory. The field conditions under which these flow mechanisms occur (either simultaneously or one predominating over the other) are poorly understood. Second, preferential flow can occur in the vadose zone as a result of heterogeneities and/or flow instability. The factors controlling preferential flow are difficult to characterize and quantify. Third, as water moves through variably saturated fractures, a portion of the flow is imbibed into the rock matrix. Although the degree of imbibition exerts a strong influence on fluid movement, the nature of fracture-matrix interaction is not well known. Fourth, solute transport in fractured rocks can exhibit complex behaviors that are difficult to interpret. Development of conceptual models may be more difficult for transport than for flow. Fifth, in a number of important cases, the interpretation of environmental tracers has lead to conclusions that seemingly contradict the initial conclusions based on classical hydrodynamic analysis. Resolution of this apparent conflict is a necessary requirement before robust conceptual models can be developed. In the remainder of this section, we present a summary of the current understanding of these five issues.
Capillary Flow, Film Flow, and Intermittent Behaviors
The basic notion of capillarity, as developed in classical theory of unsaturated flow, is that fluid is held under tension in pore space within an unsaturated
porous medium. In drier conditions (lower saturation, larger magnitude of negative pressure head), fluid occupies the smaller interstices between solid grains. With increasingly wet conditions (higher saturation, negative pressure head approaching zero), the fluid occupies increasingly larger pores. By drawing an analogy between an unsaturated porous medium and an unsaturated fracture, an understanding of fluid within interstitial space can be applied to fluid within a fracture aperture. A basic assumption is that the fluid is in contact with both sides of the fracture wall. This analogy suggests that in drier conditions, the fluid is held within small apertures, and that with increasingly wetter conditions, fluid will occupy increasingly larger apertures. The term “capillary flow” refers to flow, either in interstitial pore space or within fractures, under these conditions.
Figure 1-2 illustrates the current conceptualization of capillary flow in a fracture subject to increasing levels of saturation. Figure 1-2a illustrates a fracture that is essentially dry, and there is no fluid flow in the fracture plane. In Figure 1-2b, the shaded areas illustrate islands of water around the contact regions. At this saturation, there is still no flow in the fracture plane because the wetted areas are disconnected. In Figure 1-2c, saturation has increased to a point where the wetted areas are connected to form contiguous pathways within the fracture, thus allowing fluid flow along the fracture plane. The rate of fluid flow is controlled by the local fracture aperture of the wetted areas. In Figure 1-2d, the fracture is close to complete saturation (except for a trapped air bubble), and the flow rate along the fracture plane is higher than in Figure 1-2c.
In the case of an air-filled fracture in the vadose zone (Figure 1-3), “film flow” occurs as a thin film of fluid flowing down the fracture wall. Unlike capillary flow, in which the fluid contacts both fracture walls, the film of fluid contacts only one fracture wall, with an air phase between itself and the opposing fracture wall. The film of fluid flows under the force of gravity and is not affected by capillary forces within the fracture. Consequently, film flow can occur in small-aperture as well as large-aperture fractures, flow rate is not directly controlled by the width of the fracture aperture, and the onset of flow does not require a contiguous fluid-filled pathway within the fracture. Dragila and Wheatcraft (this report, Chapter 7) suggest that film flow is likely to occur in fractures larger than approximately 1 mm, because less energy is required to transport fluid that contacts one fracture wall compared with the transport of fluid that contacts both fracture walls.
Film flow can exhibit behaviors that are not expected in capillary flow. In laboratory experiments on a sample block of Bishop Tuff, Tokunaga and Wan (1997) demonstrated that the average fluid film velocity could be about 1,000 times faster compared to the average pore water velocity if the sample is saturated and subjected to a unit hydraulic gradient. Analysis by Dragila and Wheatcraft (this report, Chapter 7) suggests that the free-surface film flow may behave in a chaotic manner and may develop solitary waves, that is, water traveling in “lumps” over a thin film substrate. These solitary waves can carry mass (as
opposed to just energy), and travel at a speed about twice the average film speed. Solitary waves may also lead to sporadic flow behavior. If a solitary wave grows in amplitude or enters a region of smaller fracture aperture, the film may contact the opposite fracture wall. This produces a situation in which capillary force becomes important, and results in a temporary change in flow rate.
In both capillary flow and film flow situations, a portion of the infiltrating fluid can travel at velocities significantly higher than the remainder of the fluid. Thus, estimation of solute travel time based on the average advance of an infiltration front may not accurately represent the fastest solute travel. In capillary flow, heterogeneity and flow instability can cause large variations in velocity, resulting in “preferential flow” (discussed at greater length below). The combination of film flow in fractures and capillary flow in the matrix may also increase transport rates. Dragila and Wheatcraft (this report, Chapter 7) hypothesize a scenario in which fractures that are not connected with each other can nonetheless enhance infiltration rates. During an infiltration event, a fracture in the vadose zone not exposed to the surface may generate seepage by creating a capillary barrier against
the unsaturated porous matrix. If the fracture is inclined towards the vertical, the seepage may generate a free-surface film that will flow at very rapid speeds to the lower terminus of the fracture, where the fluid is again absorbed into the matrix. By this mechanism, the fracture becomes a short circuit to the porous matrix flow. A series of fortuitously placed fractures could substantially increase the transport rate.
At the field scale, the approach for modeling film flow (in a fracture network) is unclear. Current studies of film flow have been conducted primarily at the laboratory scale (in a single fracture). In Chapter 7, Dragila and Wheatcraft model film flow by a set of ordinary differential equations that is a simplified form of the momentum and conservation equations of fluid mechanics. By contrast, Tokunaga and Wan (1997) characterize film flow by measuring how “surface transmissivity” and average film thickness vary with (negative) pressure head. The resultant curves appear analogous to the hydraulic conductivity and retention curves used in traditional models of capillary flow (discussed in next section). Nonetheless, it remains an open issue whether a film flow model can be combined with the capillary flow model at the field scale, or whether the two types of flow require a fundamentally different approach.
There is also a diversity of opinions concerning whether film flow plays a significant role in infiltration. Because fluid films can travel at high velocities, Tokunaga and Wan (1997) state that “film flow can be an important mechanism contributing to fast flow in unsaturated fractures and macropores.” However, based on order-of-magnitude calculations for conditions applicable to Yucca Mountain, Pruess (1999) suggests that “rates of film flow will be small.” Resolution of this and associated issues relating to film flow will have to await future research.
Recent laboratory experiments using a transparent epoxy replica of a rock fracture (Su et al., 1999) have demonstrated intermittent flow behaviors that are not considered by classical theory. Observations of water invasion into an initially dry fracture replica indicated that flow paths in the fracture consisted of broad, water-filled regions, known as “capillary islands,” connected by thin threads of water, known as “rivulets.” Even though inflow to the fracture replica was kept constant, intermittent flow persisted in the form of cyclic snapping and re-forming of rivulets.
Figure 1-4 illustrates the intermittent flow cycle observed by Su et al. (1999). The cycle begins just after the rivulet has snapped, and water accumulates in a small capillary island just upstream of the snap point (Figure 1-4b). The capillary island grows in size (Figure 1-4c) and then migrates down the fracture replica, leaving behind a new rivulet (Figure 1-4d). Finally, the capillary island drains towards the bottom outlet (Figure 1-4e), and the rivulet eventually snaps, starting a new cycle.
Su et al. (1999) hypothesized that the type of intermittent flow observed in their experiments evolves from liquid flowing through a sequence of small to
large to small apertures. In Chapter 8, Doe analyzes the physics of capillary islands, based on the work by Furmidge (1962) on sliding of liquid drops on solid surfaces. These studies have enhanced the conceptual understanding of intermittent flow in unsaturated fractures. How to incorporate these small-scale mechanisms into field-scale models remains a difficult challenge.
While the simplest conceptualization of infiltration in the vadose zone is that of a uniform, downward-advancing wetting front in a homogeneous medium, field observations indicate that infiltration can be non-uniform (e.g., Kung, 1990a, see also discussion by Hendrickx and Flury, this report, Chapter 5). The term “preferential flow” is often used to describe non-uniform (or uneven) water flow characterized by widely different local-scale velocities. The consequence of preferential flow is that a portion of the infiltrating water can be concentrated to move along certain pathways at rates that are significantly faster than the rest of the infiltrating water. With the occurrence of preferential flow, contaminants carried by the infiltrating water can reach a given depth in less time than predicted by calculations assuming a uniform wetting front.
As pointed out by Hendrickx and Flury (this report, Chapter 5), the term “preferential flow,” in itself, does not distinguish among the causes of the non-uniform flow. Nonetheless, preferential flow is generally recognized to arise from three factors (Figure 1-5): (a) the presence of macropores (including fractures), (b) the development of flow instability, and (c) “funneling” of flow due to the presence of sloping layers that redirect the downward water movement. More sophisticated geological characterization of any particular site, with a well developed understanding of inhomogeneity, will contribute to the development of models that appropriately describe preferential flow.
Macropores refer to void space whose characteristic dimensions are significantly larger than the characteristic pore size in the rest of the medium. Soil macropores include decayed root channels, earthworm burrows, gopher holes, and drying cracks in fine-textured (clay) soils. In structured soils, where the primary grains are clustered into aggregates, the interaggregate pore space is another macropore example. Rock fractures may be also considered as a kind of macropore. Under certain conditions, the macropores allow water to move rapidly through the subsurface while bypassing the smaller pore spaces. The term “macropore flow” is often used to describe this bypassing process.
Although the physics of macropore flow in near-surface soils and in deep subsurface fractured rock are similar, macropore properties may be dynamically altered to different extents and in different ways. In near-surface soils, processes such as shrink-swell, freeze-thaw, biological activity (leading to earthworm holes and root channels), and physical manipulation (e.g., plowing of an agricultural field) can dynamically alter the preferential flow pathways. Because these processes
and their dynamics tend to decrease with depth, they are less important in the deeper subsurface, where fractured rocks predominate. Nonetheless, counter to common assumptions, rock fractures can also be dynamically altered. Berkowitz et al. (this report, Chapter 4) indicate that in some situations, unsaturated flow pathways within fractured rock can be altered by chemical dissolution and precipitation.
Preferential flow may also occur in the absence of macropores, in the form of unstable flow in which an initially subhorizontal, downward-moving wetting front breaks into “fingers.” This type of flow pattern is sometimes referred to as “fingering.” Wetting front instability can occur in a layered soil profile, as the front moves from a fine-textured layer into coarse-textured layer (Hill and Parlange, 1972). Instability and the resultant fingering can also occur in seemingly homogeneous soils, due to water repellency and/or air entrapment, and in fracture planes, as demonstrated by the laboratory experiments of Nicholl et al. (1994).
Another cause of preferential flow is lateral redirection of water through an interbedded soil profile that consists of sloping layers. Where a sloping layer acts as a barrier to flow, the downward infiltrating water will be redirected laterally above that layer. As illustrated in Figure 1.5c, the redirected flow is focused into a narrow column that can move downward at a faster rate than a broad wetting front in a homogeneous medium. The term “funnel flow” is often used for this process (Kung, 1990b). The same mechanism can occur in a subvertical fracture containing a subhorizontal obstacle, as illustrated by the computer simulation of Pruess (1999). In this simulation (Figure 1-6), the wetting front travels down a uniform fracture to a depth of 100 m in 456 days. However, when the fracture contains a subhorizontal obstacle, the redirected flow is funneled into a narrow column that reaches the 100-m depth in significantly less time (203-113 days).
A variety of approaches have been developed for modeling preferential flow. For field application, it is difficult to distinguish among the various mechanisms that cause preferential flow, and these mechanisms often act together to reinforce each other. Many models of preferential flow employ a combination of “fast flow” and “slow flow” components. Figure 1-7 (modified after Altman et al., 1996) summarizes a series of increasingly complex models that have been used to simulate preferential flow. This figure was originally prepared in the context of flow through fractured rocks, but may be viewed in the broader context of preferential flow (that is, fractures may be replaced by macropores, and matrix may be replaced by micropores).
Figure 1-7a illustrates the traditional porous medium model of capillary flow in the vadose zone. The governing equation is Richards' equation, which is derived from a combination of Darcy's law (for unsaturated flow) and the principle of mass conservation. The arrow in the schematic representation indicates fluid flow through the porous medium. The medium is characterized by two functional relationships: (1) the hydraulic conductivity curve, which gives the relation between (unsaturated) hydraulic conductivity and pressure head, and (2) the retention curve, which gives the relation between water content and pressure head. On the right side of Figure 1-7, the hydraulic conductivity curve is schematically drawn as a single curve. However, the relation is known to be hysteretic, which means that the curve during wetting is different than the curve during drying.
(unsaturated) hydraulic conductivity and pressure head, and (2) the retention curve, which gives the relation between water content and pressure head. On the right side of Figure 1-7, the hydraulic conductivity curve is schematically drawn as a single curve. However, the relation is known to be hysteretic, which means that the curve during wetting is different than the curve during drying.
Figure 1-7b illustrates a simple extension of the traditional porous medium model to account for fracture flow. This is achieved by using a composite hydraulic conductivity function that represents both matrix and fractures. The arrow in the schematic representation indicates fluid flow through the composite (matrix and fracture) system. The underlying assumption is that locally, the pressure head in the fracture is equal to the pressure head in the matrix. At low saturation (greater magnitude of negative pressure head), flow is assumed to occur only in the matrix, and so the conductivity function represents only the matrix. At high saturation (pressure head approaching zero), flow is assumed to occur in both the
matrix and the fractures, and the conductivity function represents both components. Near full saturation, the fracture conductivity may be significantly higher than the matrix conductivity. Thus, the composite hydraulic conductivity curve may have a double hump appearance. By using such an unconventional hydraulic conductivity curve, a traditional model of unsaturated flow can be adapted for preferential flow. However, this approach may not always be applicable. Mohanty et al. (1997) successfully used this approach to simulate preferential flow in a flood-irrigated agricultural field. However, Flint et al. (this report, Chapter 2) noted that the assumptions behind this approach might not hold for arid environments such as Yucca Mountain.
In the dual-porosity model (Figure 1-7c), the fracture and the matrix are separately represented by two interacting continua. For saturated conditions, this approach is widely used in modeling flow in fractured porous media. The assumption is that flow takes place only through the fracture network, as indicated by the flow arrow in the fracture continuum. Fluid exchange may occur between fracture and matrix, as indicated by the two smaller arrows between the fracture and matrix continua. However, there is no flow through the matrix blocks. For unsaturated conditions, however, this approach is seldom used because it is generally too restrictive. One exception is modeling of solute transport under steady-state flow where solute exchange may occur between the flowing water (in the fractures) and the immobile water (in the matrix) by diffusion.
The dual-permeability model (Figure 1-7d) is an extension of the dual-porosity model to allow for flow through the matrix blocks, as indicated by the additional flow arrow in the matrix continuum. This type of model is widely used for simulating preferential flow. The model involves two water retention functions, one for the fracture network and one for the matrix, and two hydraulic conductivity functions: Kf(hf) for the fracture network, and Km(hm) for the matrix, where hf and hm are the pressure heads in the fracture network and matrix, respectively. The flow between the fracture and matrix, denoted by Γw, is often expressed as:
Γw = αw (hf − hm), (1.1)
where αw is a transfer coefficient (Gerke and van Genuchten, 1993b). The value of αw exerts a strong control on flow through the fracture-matrix system, and is discussed in greater detail in the next section on fracture-matrix interaction.
In a discrete fracture model (Figure 1-7e and Figure 1-7f), fractures are explicitly represented in the model. This approach was first developed for saturated domains, where the matrix can be assumed to be impermeable (Figure 1-7e). To simulate saturated flow, discrete fracture models require data on the geometry and transmissivity of individual fractures. Such data are almost always far from complete for a given field site. In this regard, discrete fracture models share the data burden of any model that attempts to capture the detailed heterogeneity of the flow system. Field applications of discrete fracture models typically employ
(1) major, flow-controlling features determined by field characterization, (2) stochastically generated fracture networks based on statistics of fracture geometry and transmissivity, and (3) geologic understanding of the fracture origin and growth process. This topic was discussed in detail by the Committee on Fracture Characterization and Fluid Flow (NRC, 1996, Chapter 6).
When applied to the fractured vadose zone, discrete fracture models require the ability to simulate flow in the matrix as well as in the fractures (Figure 1-7f). Because fractures and matrix occur together in the same domain, the flow field can be highly complex, as indicated by the numerous arrows in the schematic representation. Fractures may be represented by two-dimensional elements embedded within a three-dimensional mesh representing the matrix. In addition to geometry data, each fracture in the model must be assigned a hydraulic conductivity curve and a retention curve. Because the data requirements are extremely demanding, unsaturated discrete fracture models have primarily been applied to simple ideal systems for conceptual understanding (e.g., Wang and Narasimhan, 1985), to laboratory experiments (e.g., Kwicklis et al., 1998), or to hypothetical field settings with orthogonal fracture networks (e.g., Therrien and Sudicky, 1996).
Discrete fracture models are also valuable for understanding the effects of fracture geometry. As discussed by Doe (this report, Chapter 8), these geometric effects include (1) diversion, (2) focusing, and (3) discontinuity in the fracture network. Diversion occurs when fractures have a strongly preferred orientation, thus imparting anisotropy to the rock. Fluid flow deviates from the direction of hydraulic-head gradient towards the preferred orientation. Focusing refers to the localization of seepage due to heterogeneities within the flow system. According to Pruess (1999), focusing can be induced by fracture intersections and terminations, as well as by asperity contacts in individual fractures. Discontinuity in the fracture network may result in the absence of a through-going preferential pathway. Thus, fluid may flow rapidly in fractures, but will slowly bleed into the matrix at fracture terminations that do not connect to other fractures.
The models illustrated in Figure 1-7 all assume applicability of Richards' equation, and hence of Darcy's law. This assumption may not be strictly correct for the fracture system. Single or dual continuum models (Figure 1-7b, Figure 1-7c, and Figure 1-7d) assume that the fracture network can be replaced by a porous medium. Analysis by Kwicklis and Healy (1993) suggests that a continuum representation of a fracture network may be more suitable for some ranges in pressure head than for other ranges. In discrete fracture models, an individual fracture is typically treated essentially as a thin layer of porous medium. However, preferential flow within a fracture (e.g., Figure 1-6) can focus transport to a narrow region. Unless each fracture is represented by a fine mesh in the model, such mechanisms will not be explicitly simulated. Furthermore, at the small-scale level, intermittent flow behaviors observed in recent visualization experiments (e.g., Su et al., 1999) are not characterized by Darcy's law.
For practical applications, however, the real issue may not necessarily be the validity of Darcy's law as such, but whether Darcy's law, even if formally invalid, can still provide an adequate and useful description of the preferential flow process. Alternative descriptions of the flow regime in fractures include Manning's equation for turbulent overland flow (Chen and Wagenet, 1992), kinematic wave theory (e.g., Germann and Beven, 1985), and simple gravity-flow models (e.g., Jarvis et al., 1991, see also Jarvis and Larsson, this report, Chapter 6). These alternatives, however, also involve simplifying approximations. It is unclear whether or not these alternatives would provide substantive enhancements for practical applications.
The extent to which fracture-matrix interaction controls unsaturated flow depends on a number of factors. Upon the initiation of fracture flow (e.g., after a precipitation event), water in the fracture may be imbibed into the rock matrix (Figure 1-8). Factors that may create “strong” fracture-matrix interaction include:
high matrix suction, large contact area between water and fracture wall, and absence of fracture coatings that impede matrix imbibition. Under these conditions, water in the fracture will be quickly absorbed into the matrix, and fracture flow cannot be sustained. By contrast, factors that create “weak” fracture-matrix interaction include: low matrix suction, small contact area between water and fracture wall, and presence of fracture coatings. Under these conditions, water in the fracture is not readily absorbed into the matrix, and fracture flow may occur to significant depths.
Although fracture-matrix interaction can, in principle, be analyzed by a detailed simulation of water flow in a fracture and propagation of the imbibition front in the matrix, current models of flow in the fractured vadose zone generally treat fracture-matrix interaction in a simplified manner. Whether or not the simplifications are justified generally depends on site-specific conditions. For example, the equivalent continuum model using a composite hydraulic conductivity function (Figure 1-7b) assumes that locally, the pressure head in the fracture is equal to the pressure head in the matrix. Under this simplification, fracture flow will not occur until the matrix is fully saturated. Although this simplification may be appropriate for certain field conditions (e.g., Mohanty et al., 1997), there are other situations where it cannot be applied (e.g., it may not be suitable for drier environments; see Flint et al., this report, Chapter 2).
In dual porosity and dual permeability models (Figure 1-7c and Figure 1-7d), fracture-matrix interaction is represented by equation (1.1), where the transfer parameter αw describes the extent of fracture-matrix flow. Many investigators, starting with Warren and Root (1963), have studied the nature of αw. Typically, αw depends on the matrix hydraulic conductivity and average fracture spacing (or equivalently, the average size of the matrix blocks). It is commonly assumed that fracture-matrix flow can occur over the entire fracture surface. However, recent findings from experiments at Yucca Mountain suggest that the above approach for calculating αw significantly overestimates the degree of fracture-matrix interaction (see discussion by Flint et al., this report, Chapter 2). A relatively large value for αw implies little to no fracture flow. By contrast, evidence indicating significant fracture flow at Yucca Mountain includes: detection of an environmental tracer (bomb-pulse chlorine-36) at depth, an estimated infiltration rate that is greater than the saturated matrix hydraulic conductivity of welded tuff units, and geochemical nonequilibrium between perched water and pore water in the rock matrix. These findings suggest that fracture-matrix interaction is relatively weak at Yucca Mountain.
To avoid overestimating the degree of fracture-matrix interaction at Yucca Mountain, Ho (1997) applied a reduction factor to decrease the fracture-matrix interface area. This was based on an assumption that if preferential flow occurs in the fracture plane, then only a fraction of the fracture plane is wetted, and the connection area between fracture and matrix is reduced. In addition, Liu et al. (1998) proposed the “active fracture model,” which hypothesized that only a
portion of the fractures in a connected network contribute to water flow, while other fractures are simply bypassed. This would further reduce the connection area between fractures and matrix. These factors led Ho (1997) to reduce the fracture-matrix interaction by four orders of magnitude in order to better model infiltration at Yucca Mountain and simulate the fast transport of bomb-pulse nuclides.
Fracture-matrix interaction can also be reduced by the presence of fracture surface coatings that impede matrix imbibition (Thoma et al., 1992). To account for fracture coatings, Gerke and van Genuchten (1993b) characterized the fracture-matrix interface by a hydraulic conductivity function, where is defined as some type of average (e.g., arithmetic, geometric) of pressure heads in the fracture and in the matrix. This approach can also be applied to structured soils. Soil aggregate can have a higher local bulk density (and hence lower conductivity) near its surface than in the aggregate center, due to deposition of organic matter, fine-texture mineral particles, or various oxides and hydroxides on the aggregate exteriors or macropore walls. Analysis by Gerke and van Genuchten (1993a) suggests that when Ka is roughly equal to Km (the conductivity of the soil matrix), pressure head in the fractures should be nearly in equilibrium with pressure head in the matrix (except for very large blocks), and this will not result in significant fracture flow. To allow fracture flow over a significant depth, Ka must be smaller than Km by several orders of magnitude, or the matrix block sizes must be relatively large. These findings emphasize the importance of fracture coatings for controlling fracture flow.
For field-scale models, the parameters describing fracture-matrix interactions will probably have to be determined by model calibration or inversion (see Bodvarsson et al., this report, Chapter 11), as there are no methods to determine these parameters by direct measurements. Although calibration may provide a good match between simulated results and existing field data, use of the model for different stresses or climate conditions may require extrapolation to conditions beyond those used in the calibration. For example, as flux increases, the fracture-matrix interaction factor may change in two ways: by an increase in the contact area between the fracture and the matrix, or by an increase in the number of fractures actively carrying water. How to accommodate such changes remains a challenging topic for future research on fracture-matrix interaction.
Modeling solute transport in the fractured vadose zone adds several layers of complexity compared with modeling fluid flow. The following discussion is limited to modeling of nonreactive solute of sufficiently low concentration that density effects can be neglected. In the traditional approach to modeling solute transport, the fluid flow problem is first solved to determine the distribution of fluid velocity. The transport of solute by this velocity field is known as advection,
which is one of two transport mechanisms recognized in classical transport theory. The second transport mechanism is hydrodynamic dispersion, which describes the spreading and mixing of solute due to (1) microscopic and local-scale variations in velocity that are not explicitly described by the above velocity field, and (2) molecular diffusion. In the classical approach, hydrodynamic dispersion is represented as a Fickian (diffusion-like) process. The resulting equation governing solute transport is known as the advection-dispersion equation.
For modeling solute transport in saturated fractured rocks, a simple extension to the classical model is applied to incorporate the effects of solute diffusion between mobile and immobile fluids. In a number of tracer experiments in saturated fractured rocks, the breakthrough curve (plot of concentration versus time) exhibits a long tail, or skewness towards later times (e.g., Novakowski et al., 1995). This feature is commonly interpreted to result from the solute moving through the fracture network while diffusing into the immobile fluid within the matrix. The conceptual flow model is that of the dual-porosity model illustrated by Figure 1-7c. The long tail in the breakthrough curve is explained by the later diffusion of solute from the matrix back into the fractures. This “matrix diffusion” can be considered as one example of solute exchange between mobile and immobile fluids. In addition to the matrix, regions of immobile fluid may exist within an individual fracture (if flow through the fracture is channeled) and in “dead-end” fractures of a network. Solute exchange between mobile and immobile fluid is commonly represented by a diffusion equation (e.g., Tang et al., 1981). Models of this type typically assume steady-state flow.
For the fractured vadose zone, solute transport models typically assume that fluid flow occurs in the matrix blocks and in the fractures. The conceptual flow model is that of the dual permeability model illustrated by Figure 1-7d. Gerke and van Genuchten (1993a) use two advection-dispersion equations, coupled by a solute transfer term, to describe solute transport in the fracture and in the matrix. Solute transfer between the two domains can occur both by fluid flow (advection) and by diffusion. The latter is often approximated as a first-order process; that is, the diffusive flux is proportional to the difference between solute concentration in the fracture and in the matrix. For certain field settings, the advective mass transfer can act in a direction opposite to the diffusive mass transfer. For example, Gerke and van Genuchten (1993a) simulated the infiltration of solute-free water into a fractured vadose zone that contained initial solute in the matrix blocks. As the solute-free water infiltrated down the fractures, a portion was imbibed into the matrix by capillary forces. At the same time, the solute in the matrix tended to diffuse back into the fractures because of the concentration difference between the two domains. The result was a highly complex response in the fracture-matrix system.
Further refinements of solute transport modeling have led to the development of multiregion models (e.g., Gwo et al., 1995; Hutson and Wagenet, 1995). In these models, the fractured rock or macroporous soil is represented by more
than two domains (or regions). These domains may be delineated in a somewhat arbitrary manner, or may be associated with geological features such as different types of fractures and matrix blocks. For example, Gwo et al. (1995) used a three-region model, composed of macropores, mesopores, and micropores, to simulate flow and transport through a laboratory column of soil from a forested watershed. In Chapter 3, Jardine et al. conceptualize the fractured weathered shale at the Oak Ridge National Laboratory to be composed of primary fractures, secondary fractures, and soil matrix. They also describe various experiments to examine the interaction and mass transfer between the different domains. In general, multiregion models are more flexible than dual- or single-domain models. However, this flexibility comes at the cost of additional model parameters. Calibration of multiregion models requires a significant amount of field data.
Solute transport can also be simulated in a discrete fracture model as shown in Figure 1-7f. Substantial insight can be gained by using such a model for examining the complex interactions between solutes in fractures and in the matrix. For example, Therrien and Sudicky (1996) simulated three-dimensional transport from a hypothetical waste facility in a fractured, low-permeability stratum overlying an aquifer. Their results show that where vertical fractures are sparse, the contaminant plume can become discontinuous, making it difficult to interpret the data from a field-monitoring program. Another value of discrete fracture models is for assessing the adequacy of continuum models such as dual-permeability or multidomain models. For a realistic field application, however, a discrete-fracture model may be problematic due to the large data requirement.
As noted by Phillips (this report, Chapter 9), environmental tracers have not been widely applied in hydrologic investigations of fractured rocks, even though they offer a powerful and direct method for assessing solute transport effects. Often, the approach to subsurface investigation is first to characterize fluid flow processes, and then to evaluate the additional processes that affect solute transport. While this approach can lead to a detailed understanding of flow processes, key solute transport mechanisms can be overlooked. In some cases, interpretation of environmental tracers has lead to drastic revisions of a conceptual model based initially solely on hydrodynamic analysis. For example, Phillips (this report, Chapter 9) points out that measurements of tritium in the fractured chalk of southern England (see Foster, 1975, and Foster and Smith-Carrington, 1980) revealed that the process of matrix diffusion was active, and this in turn explained an apparent contradiction between the rapid hydraulic response of the aquifer to recharge versus the highly attenuated movement of dissolved pollutants.
Environmental tracers consist of solutes that have a widespread or global occurrence from both natural and anthropogenic sources in precipitation, and have been entering the subsurface over large spatial and typically long temporal
scales (tens to millions of years). The large space and time scales associated with these tracers result in a signal that is naturally integrated. When samples for environmental tracer analyses are collected from a point in the subsurface, the results represent an integration of upstream transport mechanisms. This is in contrast to most hydraulic measurements where point sampling yields a point measurement.
While numerous environmental tracers have been used in subsurface hydrology, the existence of both air and water phases in the vadose zone is likely to complicate the use of some tracers, especially those that utilize dissolved gases. Environmental tracers that are likely to be useful for directly investigating fluid flow in fractured vadose zones include tritium, stable isotopes of oxygen and hydrogen in the water molecule, halides, and chlorine-36; see Phillips (this report, Chapter 9) for a brief discussion of these tracers. Other tracers that may be useful for indirectly examining fluid flow include carbon-14, uranium-series nuclides, strontium isotopes, boron-11, silicon-32, and iodine-129; see Cook and Herczeg (2000) for a discussion of these tracers. In addition, dissolved gas tracers such as helium-3, chlorofluorocarbons (CFCs, synthetic gases that were released into the atmosphere beginning in the early 1940s), and krypton-85 can provide a measure of groundwater age in the saturated zone. These age estimates may also be useful for indirectly characterizing fractured vadose zones.
Characterizing the fractured vadose zone can be accomplished by measuring the distribution of environmental tracers in the subsurface. In the absence of preferential flow, or in stratigraphic intervals of fractured rocks where matrix flow dominates, the vertical profile of tritium in the vadose zone can be used to estimate infiltration rates by identifying the depths of tritium peaks corresponding to atmospheric testing of nuclear weapons. However, if the majority of recharge at a site occurs via fracture flow, a detailed investigation of matrix processes in the vadose zone may be misleading. Furthermore, if the infiltration is localized by a few high porosity features or faults or other heterogeneities, an investigation at a random point in the map may miss the dominant vertical flow paths. In these cases, the vadose zone can be evaluated by sampling at the water table, that is, where the infiltrating fluid and solute enter into the underlying saturated zone. While the details of flow and transport processes will tend to be less well resolved by this approach, the integration (averaging) of processes may be enhanced.
The following three examples illustrate the importance of sampling from both the vadose and the saturated zones. In the first example, Wood and Sanford (1995) found significant differences in the chloride concentrations in pore waters in the unsaturated zone, compared with mean chloride concentrations in underlying groundwater in Texas. These differences were interpreted to mean that the majority of recharge occurred along preferential (macropore) pathways, which were not detected by sampling in the vadose zone. In the second example, Solomon et al. (1995) obtained groundwater ages as a function of depth below the
water table, and then calculated vertical fluid velocities and recharge rates through the vadose zone on Cape Cod, Massachusetts. Groundwater ages were zero at the water table because the dating method (tritium/helium-3) utilized a dissolved gas, which exchanged into pore air in the vadose zone. Solomon et al. (1995) also delineated the depth of the tritium bomb peak that occurred in 1963. The difference between the total travel time of water as delineated by the tritium bomb peak (30 years), and the tritium/helium-3 age of the water since it was at the water table (16 years), provided a measure of the travel time in the vadose zone (14 years). In the third example, Busenberg et al. (1993) found that concentrations of CFCs dissolved in groundwater near the water table were much larger than expected based on measurements of the CFC content of pore air immediately above the water table. They concluded that preferential flow was not equilibrating with the vadose zone atmosphere and was responsible for significant recharge to the Snake River Plain aquifer.
A mass balance of environmental tracers within the vadose zone can be useful for evaluating regional-scale solute and fluid fluxes, without necessarily understanding the details of fracture flow, fracture-matrix interaction, etc. (e.g., Cook et al., 1994). This approach assumes that a mass balance can be adequately formulated using subsurface measurements (i.e., that spatial variations in the tracer distribution can be adequately sampled), and the tracer input is sufficiently well known. It can be argued that over long time scales, diffusion processes tend to reduce small-scale spatial variation in tracer concentrations, making a mass balance feasible. However, this approach has not been extensively tested in fractured vadose zones.
Environmental tracers may also be useful for evaluating and testing specific processes that are included in a conceptual model. For example, Desaulniers et al. (1981) evaluated the influence of fractures on solute transport in a clay-rich till in Southern Ontario. They found that profiles of oxygen isotopes and chloride were best explained by molecular diffusion, with minimal advective transport in fractures. Because of their usefulness, environmental tracers should be considered a primary method for investigating fractured vadose zones and for formulating and testing conceptual models, and they should be included in the field investigation strategy from the very beginning of site characterization.
CONCLUSIONS AND RECOMMENDATIONS
This report discusses the process through which conceptual models of flow and transport in the fractured vadose zone are developed, tested, refined, and reviewed. A conceptual model is defined 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. The conclusions presented below are grouped according to the two major topics addressed in this report: (1) general considerations during
the development and testing of conceptual models, and (2) flow and transport in the fractured vadose zone. These conclusions are followed by the Panel's suggestions for research activities that will contribute to the conceptual modeling process.
Development of the conceptual model is the most important part of the modeling process. The conceptual model is the foundation of the quantitative, mathematical representation of the field site (i.e., the mathematical model), which in turn is the basis for the computer code used for simulation. Given a sufficiently robust conceptual model, different mathematical formulations will probably produce similar results. By contrast, an inappropriate conceptual model can easily lead to predictions that are orders of magnitude in error.
The context in which a conceptual model is developed constrains the range of its applicability. A conceptual model is by necessity a simplification of the real system, but the degree of simplification must be commensurate with the problem being addressed. Thus, a conceptual model developed for addressing one type of problem may not be adequate for another type of problem. For example, a conceptual model developed for estimating recharge flux may not be adequate for estimating contaminant travel time from land surface to the water table.
It is important to recognize that model predictions require assumptions about future events or scenarios, and are subject to uncertainty. Quantitative assessment of prediction uncertainty should be an essential part of model prediction. A suite of predictions for a range of different assumptions and future scenarios is more useful than a single prediction. Uncertainty in model predictions can be partially quantified by sensitivity analysis (how uncertainties in estimated model parameters affect model predictions), by using a statistical-or stochastic-based calibration procedure, or by formulating the mathematical model in a probabilistic framework. However, it is difficult to quantitatively assess the possibility that the conceptual model might not adequately represent the major features and processes in the real system.
Testing and refinement of the conceptual model are critical parts of the modeling process. The initial conceptual model is developed based upon limited field data and is susceptible to biases reflected by the disciplinary background and experience of the analyst. Therefore, site investigation should not be designed solely to support the initial conceptual model. Reasonable alternative conceptualizations and hypotheses should be developed and evaluated. In some cases, the early part of a study might involve multiple conceptual models until alternatives are eliminated by field results.
Although model calibration does provide a certain level of model testing, a good fit to the calibration data does not necessarily prove that the
model is adequate to address the issues in question. The significance of a good fit to calibration data generally depends on the nature of the data. A model that matches different types of calibration data (e.g., heads and fluxes) collected under different field conditions (e.g., at different water contents) is likely to be more robust than a model that matches a limited range of calibration data. However, if the model cannot be calibrated to match the calibration data, this is an indication that the conceptualization should be re-examined.
Checking model simulation results against field data (that were not used for calibration) is one, but not the exclusive, approach to model testing. In some cases, all field data are needed for calibration, and none are left for further testing. However, this does not mean that the model cannot be used for prediction. A broader view of model testing is to develop greater confidence that the model provides a good representation of the real system. This can be achieved by strengthening the justifications for model assumptions, and by evaluating alternative hypotheses.
From an operational perspective, the goal of model testing is to establish the credibility of the model. A credible model is essential if it is to gain acceptance by parties involved in decision-making or problem resolution. In addition to testing and evaluation, the credibility of a model can be enhanced by peer review undertaken by an independent panel of experts, and by maintaining an open flow of information so that the model is available for scrutiny by concerned parties.
There exists a body of field evidence indicating that infiltration through fractured rocks and structured soils does not always occur as a wetting front advancing at a uniform rate. Large variations in fluid velocity (i.e., preferential flow) may be caused by (a) the presence of macropores and fractures, (b) flow instability, or (c) funneling effects. Thus, model simulation based upon a uniform wetting front advancing down a homogeneous medium may provide erroneous estimates of flux and travel times through the vadose zone. Sophisticated characterization of geological inhomogeneity within the vadose zone increases the likelihood that non-uniform flow can be appropriately modeled.
The current state of knowledge is not adequate to determine which processes are likely to control unsaturated flow and transport at a given field site. Laboratory and theoretical analyses demonstrate that film flow in fractures can transport fluid and solute at rates substantially higher than transport by capillary flow. However, at the field scale, the significance of film flow and the modeling approach are topics of controversy. The field environments in which film flow plays a significant role during infiltration are poorly understood. In addition, it is unclear whether film flow can be incorporated into traditional
models of capillary flow by defining effective curves for hydraulic conductivity and retention, or whether it requires a fundamentally different set of governing equations to describe the dynamics of a water film.
Although not identical, structured soils and fractured rocks exhibit many similarities in flow and transport processes. Macropores and aggregates in structured soils are respectively analogous to fractures and matrix blocks in rock. However, soil studies are typically conducted in the shallow subsurface, where macropores may be dynamically altered to a greater degree than rock fractures at greater depths. Nonetheless, knowledge gained from study of one medium may be useful for the other. Communication between workers in the soil science field and in the fractured rock field will be of benefit to both groups.
Models of varying complexity have been developed for preferential flow, but their adequacy for field-scale application requires further testing. The approach in many current models is to avoid explicitly simulating the mechanisms that cause preferential flow. Instead, the model is implemented to simulate fast and slow flow, by use of a composite hydraulic conductivity curve or by dual-permeability domains. Such approaches have been successfully applied to laboratory and small-scale experiments. However, further testing is needed to examine whether these models are adequate for field-scale application over a broad range of field conditions. This issue is of particular concern in the fractured vadose zone because of the inherently nonlinear nature of processes involved. As flow conditions change, different flow and transport mechanisms, not represented in the model, may become important, leading to large errors in predictions.
The interaction between fracture and matrix exerts a strong control on fluid and solute movement. However, the strength of this interaction in the field is not well known. The simplified representation of this interaction in current models also requires further evaluation. Factors controlling fracture-matrix interaction include the density of water-transmitting fractures, the amount of wetted area on the fracture surface, the hydraulic conductivity of the matrix, and hydraulic conductivity at the fracture-matrix interface. Current models lump these factors into a transfer coefficient that is determined by model calibration rather than by direct measurements. Whether or not this approach can adequately simulate flow under a range of field conditions requires further evaluation.
Solute transport in the fractured vadose zone can exhibit complex behavior due to the large variations in fluid velocity, and the interplay of advective and diffusive transport between fractures and matrix. Better understanding of such systems is required in order to effectively analyze complex responses. Solute transport models are more complex than flow models, and can involve multiple regions to represent the diversity of macropore and micropore sizes. To apply these models, greater guidance is needed on how to delineate different pore regions, and how to determine the parameters that characterize solute exchange between pore regions.
Environmental tracers should be included in field investigation strategies from the very beginning of a site characterization program. In a number of studies, geochemistry and environmental tracer data have led to substantial revisions of the conceptual models initially developed based upon hydrodynamic analysis. These experiences emphasize the need for better integration of geochemistry and environmental tracers early in the model development process.
Flow and transport in the fractured vadose zone have been and will continue to be an active area of research in both the soil science and subsurface hydrology disciplines. The research recommended in this report is not meant to be inclusive. Instead, the list below reflects topics that deserve greater attention so that conceptual models of flow and transport can be improved, and to address the issues identified in this report.
Fundamental research to understand flow and transport processes in unsaturated fractures should continue. Better understanding of fundamental processes will improve the model representation of the real system. Particular emphasis should be placed on understanding mechanisms that cause non-uniform (preferential) flow, film flow, and intermittent behavior.
Research is needed to understand the spatial variability in vadose zone properties, and to develop upscaling methods. Spatial variability is a key cause of model uncertainty, because the subsurface cannot be exhaustively sampled. Furthermore, vadose zone properties are typically determined by small-scale laboratory measurements. To use these small-scale measurements, upscaling methods are needed to derive field-scale flow and transport properties needed in models. Such upscaling methods should be based on a thorough understanding of small-scale processes, together with an understanding of how these interact and contribute to large-scale phenomena.
There is a need for comprehensive field experiments in several fractured vadose zone geologic environments. These experiments should be designed to understand the controlling processes (capillary flow, film flow, and intermittent behavior) for a broad range of field conditions, to evaluate methods of parameter upscaling, and to test alternative conceptual models.
Current models should be evaluated for their adequacy for simulating flow and transport in the presence of fingering, flow instability, and funneling. One approach is to construct a model with detailed representation of small-scale heterogeneities based on high-resolution field or synthetic data, so that the processes causing preferential flow are explicitly simulated. The model results can then be compared to results from simpler models that do not explicitly simulate preferential flow. Of particular importance is the evaluation of transfer coefficients to represent fluid and solute exchange between fracture and matrix.
There is a need to develop quantitative assessment of prediction uncertainty for models of flow and transport in the fractured vadose zone. Meaningful quantification of uncertainty should be considered an integral part of any modeling endeavor, as it establishes confidence bands on predictions given the current state of knowledge about the system. If prediction uncertainties are realistically quantified, postaudit studies can be carried out in a systematic hypothesis-testing framework, which can provide a great deal of insight about the predictive capabilities of the model.
Research should be undertaken to develop improved techniques for geochemical sampling from the fractured vadose zone. Current sampling technology is very limited, especially for sampling at depth. Destructive core sampling followed by fluid extraction is a viable approach for sampling matrix water in the unsaturated zone. Sampling fluids directly from fractures remains problematic. Improved sampling techniques will facilitate the use of environmental tracers and geochemical data for conceptual model building.
Altman, S. J., B. W. Arnold, R. W. Barnard, G. E. Barr, C. K. Ho, S. A. McKenna, and R. R. Eaton, 1996. Flow Calculations for Yucca Mountain Groundwater Travel Time (GWTT-95) . Report SAND96-0819. Albuquerque, N. Mex.: Sandia National Laboratories. 170 pp.
Anderson, M. P., and W. W. Woessner, 1992. Applied Groundwater Modeling. New York: Academic Press. 381 pp.
Busenberg, E., E. P. Weeks, L. N. Plummer, and R. C. Bartholemay, 1993. Age Dating Ground Water by Use of Chlorofluorocarbons (CCl3F and CCl2F2), and Distribution of Chlorofluorocarbons in the Unsaturated Zone, Snake River Plain Aquifer, Idaho National Engineering Laboratory, Idaho. U.S. Geological Survey Water-Resources Investigations Report 93-4054 . Reston, Va.: U.S. Geological Survey. 47 pp.
Chen, C., and R. J. Wagenet, 1992. Simulation of water and chemicals in macropore soils. I. Representation of the equivalent macropore influence and its effect on soil-water flow. Journal of Hydrology 130: 105-126.
Cook, P., and A. L. Herczeg, 2000. Environmental Tracers in Subsurface Hydrology. Boston: Kluwer Academic Publishers. 529 pp.
Cook, P. G., Jolly, I. D., Leaney, F. W., Walker, G. R., Allan, G. L., Fifield, L. K., and Allison, G. B., 1994. Unsaturated zone tritium and chlorine-36 profiles from southern Australia: their use as tracers of soil water movement. Water Resources Research 30(6): 1709-1719.
Desaulniers, D.E., J.A. Cherry, and P. Fritz, 1981. Origin, age and movement of pore water in argillaceous quaternary deposits at four sites in southwestern Ontario. Journal of Hydrology 50: 231-257.
Fabryka-Martin, J. T., A. V. Wolfsberg, S. S. Levy, J. L. Roach, S. T. Winters, L. E. Wolfsberg, D. Elmore, and P. Sharma, 1998. Distribution of Fast Hydrologic Paths in the Unsaturated Zone at Yucca Mountain. Paper presented at 8th Annual International High-Level Radioactive Waste Management Conference, Las Vegas, Nev. American Nuclear Society, La Grange Park, Ill. p. 264-268.
Foster, S. S. D., 1975. The chalk groundwater tritium anomaly, a possible explanation. Journal of Hydrology 25: 159-165.
Foster, S. S. D., and A. Smith-Carrington, 1980. The interpretation of tritium in the chalk unsaturated zone. Journal of Hydrology 46: 343-364.
Furmidge, C. G. L, 1962. Studies at phase interfaces, 1. The sliding of liquid drops on solid surfaces and a theory for spray retention. Journal of Colloidal Science 17: 309-324
Gerke, H. H., and M. Th. van Genuchten, 1993a. A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. Water Resources Research 29(2): 305-319.
Gerke, H. H., and M. Th. van Genuchten, 1993b. Evaluation of a first-order water transfer term for variably saturated dual-porosity flow models. Water Resources Research 29(4): 1225-1238.
Germann, P., and K. J. Beven, 1985. Kinematic wave approximation to infiltration into soils with sorbing macropores. Water Resources Research 21(7): 990-996.
Glass, R. J., and J-Y. Parlange, 1989. Wetting front instability. 1. Theoretical discussion and dimensional analysis. Water Resources Research 25(6): 1187-1194.
Gwo, J. P., P. M. Jardine, G. V. Wilson, and G. T. Yeh, 1995. A multiple-pore-region concept to modeling mass transfer in subsurface media. Journal of Hydrology 164: 217-237.
Hill, D. E., and J-Y. Parlange, 1972. Wetting front instability in layered soils. Proceedings Soil Science Society of America 36: 697-702.
Ho, C. K., 1997. Models of Fracture-Matrix Interactions During Multiphase Heat and Mass Flow in Unsaturated Fractured Porous Media. Paper presented at 6th Symposium on Multiphase Transport in Porous Media, ASME International Mechanical Engineering Congress and Exposition , American Society of Mechanical Engineering, Dallas, Texas.
Hutson, J. L., and R. J. Wagenet, 1995. A multiregion model describing water flow and solute transport in heterogeneous soils. Soil Society of America Journal 59: 743-751.
Jarvis, N. J., P-E. Jansson, P. E. Dik, and I. Messing, 1991. Modelling water and solute transport in macroporous soils. I. Model description and sensitivity analysis. Journal of Soil Science 42: 59-70.
Konikow, L. F., and J. D. Bredehoeft, 1992. Ground-water models cannot be validated. Advances in Water Resources 15(1): 75-83.
Kung, K-J. S., 1990a. Preferential flow in a sandy vadose zone. 1. Field observations. Geoderma 46: 51-58.
Kung, K-J. S., 1990b. Preferential flow in a sandy vadose zone. 2. Mechanism and implications. Geoderma 46: 59-71.
Kwicklis, E. M., and R. W. Healy, 1993. Numerical investigation of steady liquid water flow in a variably saturated fracture network. Water Resources Research 29(12): 4091-4102.
Kwicklis, E. M., F. Thamir, R. W. Healy, and D. Hampson, 1998. Numerical Simulation of Air- and Water-Flow Experiments in a Block of Variably Saturated, Fractured Tuff from Yucca Mountain, Nevada . U.S. Geological Survey Water-Resources Investigations Report 97-4274 . Denver, Colo.: U.S. Geological Survey. 64 pp.
Liu, H. H., C. Doughty, and G. S. Bodvarsson, 1998. An active fracture model for unsaturated flow and transport in fractured rocks. Water Resources Research 34(10): 2633-2646.
Mohanty, B. P., R. S. Bowman, J. M. H. Hendrickx, and M. Th. van Genuchten, 1997. New piecewise-continuous hydraulic functions for modeling preferential flow in an intermittentflood-irrigated field. Water Resources Research 33(9): 2049-2063.
National Research Council (NRC), 1996. Rock Fractures and Fluid Flow: Contemporary Understanding and Applications . Washington, D.C.: National Academy Press. 551 pp.
Nicholl, M. J., R. J. Glass, and S. W. Wheatcraft, 1994. Gravity-driven infiltration instability in initially dry nonhorizontal fractures. Water Resources Research 30(9): 2533-2546.
Novakowski, K. S., P. A. Lapcevic, J. Voralek, and G. Bickerton, 1995. Preliminary interpretation of tracer experiments conducted in a discrete rock fracture under conditions of natural flow. Geophysical Research Letters 22(11): 1417-1420.
Oreskes, N., K. Shrader-Frechette, and K. Belitz, 1994. Verification, validation, and confirmation of numerical models in the earth sciences. Science 263: 641-646.
Peters, R. R., and E. A. Klavetter, 1988. A continuum model for water movement in an unsaturated fractured rock mass. Water Resources Research 24(3): 416-430.
Pruess, K., 1999. A mechanistic model for water seepage through thick unsaturated zones in fractured rocks of low matrix permeability. Water Resources Research 34(4): 1039-1051.
Solomon, D. K., R. J. Poreda, P. G. Cook, and A. Hunt, 1995. Site characterization using 3He/3He ground water ages, Cape Cod, Mass. Ground Water 33(6): 988-996.
Su, G. W., J. T. Geller, K. Pruess, and F. Wen, 1999. Experimental studies of water seepage and intermittent flow in unsaturated, rough-walled fractures. Water Resources Research 35(4): 1019-1037.
Tang, D. H., E. O. Frind, and E. A. Sudicky, 1981. Contaminant transport in fractured porous media: Analytical solution for a single fracture. Water Resources Research 17(3): 555-564.
Therrien, R., and E. A. Sudicky, 1996. Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media. Journal of Contaminant Hydrology 23: 1-44.
Thoma, S. G., D. P. Gallegos, and D. M. Smith, 1992. Impact of fracture coating on fracture/matrix flow interactions in unsaturated porous media. Water Resources Research 28(5): 1357-1367.
Tokunaga, T. K., and J. Wan, 1997. Water film flow along fracture surfaces of porous rock. Water Resources Research 33(6): 1287-1295.
Tsang, C. F., 1991. The modeling process and model validation. Ground Water 29(6): 825-831.
Wang, J. S. Y., and T. N. Narasimhan, 1985. Hydrologic mechanisms governing fluid flow in a partially saturated, fractured, porous medium. Water Resources Research 21(12): 1861-1874.
Warren, J. E., and P. J. Root, 1963. The behavior of naturally fractured reservoirs. Society of Petroleum Engineers Journal 3(5): 245-255.
Wood, W. W., and W. E. Sanford, 1995. Chemical and isotopic methods for quantifying groundwater recharge in a regional, semiarid environment. Ground Water 33(3): 458-468.