4
Conceptual and Numerical Model Development
Site models are necessary to build a conceptual picture of a site’s geology and hydrology, to interpret field experiments, to estimate where and how fast contaminants may migrate, and to aid design of remediation processes and forecast performance. This chapter outlines an approach for fractured aquifer model development and emphasizes a conceptualization of the fractured rock environment that captures important features and processes of the geologic, hydrologic, and geomechanical environments and can be used to inform engineering and resource management decisions. This report will refer to such conceptual models as hydrostructural models. A site conceptual model developed early in site management and engineering can guide characterization, design, and choice of detailed numerical modeling approach appropriate for a given site, project goals, and available resources. An appropriate range of conceptual models underlies more detailed modeling efforts. Relatively simple calculation methods are available to estimate the bounds of fluid and contaminant behavior and uncertainties for even the most complex fractured rock site. This type of calculation is hereafter referred to as a scoping calculation.
Hydrostructural conceptual model development has been used by some in the European high-level radioactive waste community (e.g., Blomqvist et al., 1998) to characterize more significant fractures as deterministically as feasible. A motivation to use hydrostructural models for fractured rock is the general observation that a few discrete fractures tend to dominate flow and transport (see Box 3.1). The art of hydrostructural modeling includes distinguishing between these significant fractures and “background” fractures, which can be characterized adequately using, for example, stochastic approaches. This construct for conceptual model development is emphasized in this report because it encourages a more interdisciplinary approach to decision making related to engineering, management, or remediation of the fractured rock environment. The term emphasizes that a new way of thinking about engineering in the fractured rock environment is necessary. There remains, however, the significant challenge of how hydrostructural models and site characterization methods can contribute to the deterministic identification and interconnection of key permeable fractures. Lacking their identification, practice depends on stochastic models, which are often less satisfactory.
Conceptual models may be as simple as a pictorial representation of features, processes, and events, or can be detailed to the point at which parameters for those features, processes, and events are assigned theoretical values that can be used and tested in numerical models. This chapter describes how to conceptualize and quantify parameters of the fractured rock site conceptual model and the utility of doing so. The discussion concludes with an outline of methods used to construct and analyze numerical models to ensure there are quantitative connections between observations and numerical models. Simple approaches and quantifying both conceptual and parameter uncertainty are emphasized. This chapter focuses on approaches to numerical modeling rather than specific programming packages. Several commercial and non-commercial numerical models are available. Chapter 5 describes methods to collect data to support detailed numerical models.
DEFINING AND DEVELOPING HYDROSTRUCTURAL MODELS
Hydrostructural models for fractured rock, like numerical models discussed later, are approximations of reality used to describe the geologic and hydraulic characteristics of a site. Hydrostructural models are site specific, and usually include some combination of discrete features (i.e., fractures, faults, karst, dikes, brecciated formations, and permeable bed boundaries) and rock matrix blocks (defined by the rock material between discrete features). Without an appropriate underlying hydrostructural model as a conceptual framework, analysis or numerical modeling will not yield a thorough understanding of discrete fracture geometry and its effect on flow, transport, and storage for even the smallest project.
The term “upscaling” is used to describe the process of going from a detailed, high-resolution model to a coarse model, for example, through averaging (e.g., Rubin and Gómez-Hernández, 1990; Durlofsky, 1992; Ouenes and Hartley, 2000). Upscaling is often necessary to allow numerical computation at larger scales with fewer computational resources (e.g., Durlofsky, 2002). It is required, generally, to translate from a hydrostructural to a numerical model. Hydrostructural models, therefore, can be simple, but need to represent the full range of geologic features, characteristics, and geometries that are important hydraulically to the processes of interest, so that those processes, features, and behaviors are not lost during upscaling.
Not all hydrostructural model elements need to be quantified, and models should be informed but not limited by existing data. Sparse, biased, or censored data may result in an erroneous conceptual model. A vertical borehole, for example, may not capture hydraulically important vertical fractures known to exist at a site. It is therefore better to include those vertical fractures in the hydrostructural model, even if their parameters are unknown, rather than omitting them because site-specific data are not (and may never be) available. Hydrostructural models are not constrained by the same assumptions as specific numerical models (e.g., gridding and grid resolution, equivalent continuum concepts, or simplified porosities).
Hydrostructural Model Templates
It has been concluded, based on limited empirical data, that analysts often have difficulty selecting an appropriate conceptual model and that initial models chosen are often found to be incorrect (Bredehoeft, 2005). Hydrostructural conceptualization begins with the simplest visualization of site geology and grows in complexity as information is gathered until transport and remediation processes can be incorporated into the model. A dolomite site may be conceptualized, for example, by first assuming a set of parallel bedding planes and two sets of fractures perpendicular to bedding. Detail is added to the model to reflect solution features likely found in dolomite. More details are added as fracture and solution feature properties are characterized. The hydrostructural model should eventually describe the roles in and relative importance of fractures and rock matrix in flow (permeability; see Box 4.1), storage (porosity, storativity), and connectivity (discrete pathways). Table 4.1 provides generalized hydrostructural model examples for some geologic settings that may be good starting points for conceptualization. The maturing hydrostructural model can inform characterization strategies (see Chapter 5), data from which can validate and inform the hydrostructural model (see Chapter 7 for information about the use of the observational method to update models).
Representing a fractured rock as a system of discrete fractures and rock matrix is critical to conceptualizing contaminant transport and remediation. The conceptualized division is not strictly
consistent with physical reality, however, because fractures can occur at scales from microns to kilometers. The scale of the contaminant plume being studied, for example, may determine which fractures and fracture characteristics (e.g., porosities and fracture permeabilities) are relevant to transport or storage, and therefore whether they should be conceptualized as part of the effective rock matrix or fractures. The full range of storage and flow porosities needs to be considered and assigned in light of the information needed for forecast of contaminant transport and design of a remediation process.
TABLE 4.1 Sample Hydrostructural Model Templates for Generic Geologic Settings^{a}
Geologic Setting | Starting Hydrostructural Model Template |
---|---|
Fractured crystalline rock (e.g., granite, diorite) | Very low (< 1%) rock matrix porosity, with fractures at a continuous range of scales from mm to 10s of km providing. Fractures occur in one or more sets, with both lateral and vertical spatial variation. Rock mass porosity and permeability based on smaller fractures and rock matrix. Essentially no flow and very little storage in matrix. (Photo of Sentinel Ridge. Courtesy of L. Lau.) |
Fractured carbonate rocks (e.g., limestone, dolomite, oil shales) | Mechanical unit bedding with fractures perpendicular to bedding. Because bedding can be distorted by folding or faulting, fractures will be oriented in a local coordinate system consistent with bedding. Fractures will be bound by mechanical units. There are frequently solution enhanced features (karst), including thin brecciated beds. Superposed on this may be fracture sets related to faulting (in a global coordinate system) and folding (in a bed-local coordinate system). Rock matrix porosity and permeability can be significant. Photo of part of the Supai Group in the Grand Canyon. Courtesy of the U.S. Geological Survey.) |
Argillaceous rocks (e.g., shales) | Similar to fractured carbonates, with both bed-local coordinate system fractures, and superposed tectonic and structural fractures. Rock matrix permeability is generally low enough to inhibit flow such that smaller fractures may form a major portion of the rock mass permeability and porosity. Photo of the Geneseo/Burket gas shale at Taughannock Falls State Park, Trumansburg, New York. Natural hydraulic factures emerge from the top of the black gas shale and penetrate upward into a gray shale. Courtesy of T. Engelder, The Pennsylvania State University.) |
Extrusive Igneous rocks (e.g., tuff, basalt) | Unique in their complexity. Fracture patterns can be systematic (particularly in basalt) but can exist in layering of varying porosities made complex by cooling processes, multiple eruption and faulting events, igneous intrusions, brittleness, and interbedding of paleo-alluvial and ash-flow materials. Photo of pillow basalts at Mary’s Peak, Oregon. Courtesy of R. Keller.) |
^{a} This table does not represent all geologic environments and should not be considered a comprehensive tool.
Microstructural Models of Fractures
Fractures are often not defined by a single pair of parallel flat surfaces as they are often modelled. Because diffusion and surface sorption rates can be controlled by the ratio of flowing volumes to fracture surface area, diffusion and sorption rates are a function of the geometry and topology of surfaces within a fracture. Individual fractures may therefore need to be characterized at smaller scales—even down to the micro level: a microstructural model. A microstructural model allows conceptualization for multiple fracture surfaces, infillings, coatings, and altered rock zones that make up a fracture. It is these characteristics that determine the geochemistry of retention and the relationships between mechanical, filtering, transport, and hydraulic apertures. See Box 4.2 for an example of a microstructural model.
Hydrostructural Conceptualization of Faults
Details about faults that control fluid flow and contaminant fate and transport can also be incorporated into hydrostructural models. Only a limited number of fractures, some of which may be faults, need to be considered explicitly in hydrostructural modeling because some fractures do not conduct fluid or accommodate storage or transport. For the purpose of conceptualization, fractures can be grouped as
- Low transmissivity fractures that add to pore volume but do not conduct flow. Fluid behavior in these fractures is similar to fluid behavior in the surrounding rock matrix, and the fractures can therefore be conceptualized as part of the effective rock mass. Transmissivity considered “low” depends on the application—for radioactive waste, this value might be 1 × 10^{–8} m^{2}/s, whereas it might be 1 × 10^{–7} m^{2}/s for oil.
- Higher transmissivity fractures that provide significant flow and reactive surface area.
- Significant conductive and flow barrier features defined explicitly to describe flow pathways and rock matrix storage.
The fault itself can be modeled as a fault core and a damage zone. Each may be conductive or not, and the damage zone can be small or large compared to the fault core. Elements incorporated into the hydrostructural model include
- Fault core gouge—ground material found between fault walls formed as result of mechanical erosion resulting from fault displacement. Fault core gauge can be several meters wide and create barriers to flow across a fault.
- Fault core breccia and other high-permeability and high-porosity infill materials that can provide major storage for contaminants or key flow pathways in the fault core.
- Hangingwall and footwall disturbed zones. These zones can have significant permeability and porosity parallel to faults, even if the fault itself is a flow barrier (low permeability) perpendicular to the fault due to core gouge.
Figure 4.1 illustrates the range of fault-zone architectures and associated permeable structures possible in a hydrostructural model. Each corner of the figure represents an end member along a continuum of possible scenarios. When the damage zone is small, and the fault core is also a small fraction of the local volume, the fault itself can be a conductive feature, with fluids traveling along the fault. If the damage zone is a region that contains many fractures, and the fault core is relatively small, then both the fault core and damage zone contribute to flow in a distributed conduit. If the fault core is large and acts as a flow barrier, and the damage zone is also large, then the conduit created by the fault can allow enhanced flow along the fault but block flow across it. The combination of a small damage zone with a fault core barrier can create a barrier to flow transport across the fault and not enhance flow along the fault.
These end member examples represent very different settings for flow and exchange of contaminants with nearby zones. Failure to account for these elements in the hydrostructural model can limit the usefulness of the models for design and forecast of contaminant transport and mitigation.
Rock Matrix Conceptualization
There are at least two timescales for remediation: one associated with groundwater residence times and another associated with the diffusion of contaminant in the immobile porosity of the rock. This is true for single-porosity porous (non-fractured) and multiple-porosity fractured rock. Blocks of rock matrix may be of different porosity than the fractures surrounding them, and the shape and size of the matrix blocks affect diffusive, advective, geochemical, and other types of exchanges between fracture and matrix. The oil industry conceptualizes interactions between rock matrix and fractures in terms of three alternative rock matrix block aspect ratios: slab, matchstick, and sugar cube (or “sphere”; see Figure 4.2). Each of these model blocks exhibits distinct contaminant transfer behavior between the advective (fracture) porosity and the immobile and storage (rock matrix) porosity, affecting both hydraulic pressures and geochemical equilibrium. The timescales of transfer depend on the characteristic size defined by the minimum rock block dimension (Kazemi and Gilman, 1993; Carrera et al., 1998). In remediation, contaminants in matchstick-shaped blocks are accessed quickly, while those in a matrix with sugar cube dimensions are accessed more slowly. The shape of the tail of the contaminant plumes is therefore influenced by matrix block shape, and block size also influences contaminant storage.
Matrix blocks are commonly modeled as sugar cubes, but better models result if variability is accounted for. The distribution of rock block shapes in the hydrostructural model can be calculated from estimates of the storage and flow porosities and the estimated spatial distribution of contaminants.
QUANTIFYING THE HYDROSTRUCTURAL MODEL
Site characterization data will inform but not completely constrain the hydrostructural model. The next stage of model development is to add parameters that describe each of the fracture and matrix elements.
Data Integration
Generic parameterizations applied early in hydrologic modeling avoid common bias in the models. Values derived from field tests are not always representative across a site. Every field test or measurement has some spatial scale of investigation (see Chapter 5), and results of those tests should be applied to other length scales only with caution. For example, a single hydraulic test might yield a specific local transmissivity, but this value may or may not be representative of transmissivity on larger scales or in adjacent areas. Such measurements can do little to constrain the range of values that may be present at the site scale. Initial parameter value estimates can be based on conditions observed at similar sites.
Not all elements of a site hydrostructural model can be parameterized deterministically. Successful parameterization depends on the ability to integrate multiple kinds of information and lines of evidence from multiple sources (see Chapter 5 for discussion of data integration). Better predictions and practices are possible if both parametric uncertainties and uncertainties related to the appropriateness of the hydrostructural conceptualization are considered. Table 4.2 lists the types of data that could be integrated to determine various hydrostructural elements associated with fractures and rock matrix. Hydrostructural elements related to boundary conditions do not appear on the table because boundary conditions do not conveniently fit into a table of this type. The next sections describe more fully the parameterization of fractures, rock matrix, and boundary conditions.
TABLE 4.2 Data Integration for Parameterization of Hydrostructural Models
Method of Parameter Determination | Hydrostructural Element Related to Fractures | Hydrostructural Element Related to Rock Matrix |
---|---|---|
Geologic (including core, outcrops, and lab testing) | Spatial pattern sets and orientation distributions Size, shape, roughness, and infilling Microstructural model Aperture |
Shape Permeability and porosity Geochemical properties |
Geophysical (including fracture image, wireline, and surface seismic methods) | Spatial pattern Transmissivity distributions Orientation distributions Anisotropy and heterogeneity Major structures (i.e., large faults, fold crests and hinge zones, dikes, large karst features) |
Elastic properties |
Hydrodynamic (including well- and multi-well-based tests, long-term pressure and temperature monitoring, and hydrophysical flow logging) | Transmissivity, storativity, aperture measures, porosity measures Functional relationships among hydrodynamic properties Functional relationships between hydrodynamic property and geometric properties |
Hydraulic conductivity Porosities measures |
Geomechanical | Strength Shear and normal stiffness In situ stress |
Strength Deformability In situ stress |
Thermal | Thermal deformability Functional relationships between hydrodynamic properties and temperature Heat capacity |
Thermal deformability Functional relationships between hydrodynamic property and temperature Heat capacity |
Geochemical | Porosity measures Diffusion measures Sorption/retention measures |
Porosity measures Diffusion measures Sorption/retention measures Fluid/rock interactions: dissolution/precipitation reactions |
Fracture Parameterization
Fracture parameterization within hydrostructural models is most useful for many applications if based on site microstructural models, spatial patterns, hydrodynamic properties, and geochemical and biologic properties (e.g., retention, natural attenuation, and those related to mitigation). Table 4.3 provides a list of various fracture properties and their measures. Following are examples of how various fracture elements could be parameterized.
TABLE 4.3 Fracture Measures
SOURCE: Courtesy of W.S. Dershowitz.
Hydrodynamic Properties
Hydrodynamic properties that dominate contaminant migration rates need to be parameterized. These include transmissivity, aperture measures, storativity, and hydraulic diffusivity (transmissivity/storativity). Transmissivity, storativity, and hydraulic diffusivity should be quantified and applied considering the correlations between them and the importance of those parameters for understanding flow and transport behavior. Unless there is evidence of little variability in fracture geometry and in hydrodynamic properties, assume there will be spatial and statistical variability of these within and among fractures.
Fracture Intensity
As discussed in previous chapters, fracture intensity is the number of fractures in a unit volume. The number of fractures indicates matrix diffusion rates associated with surface fractures. Table 4.3 lists measures of fracture intensity.
Fracture Roughness and Infilling
Fracture roughness and infilling are often neglected fracture characteristics. Under any given set of boundary conditions, contaminant travel time, retention, and reactive surface areas depend on whether flow is effectively one dimensional (e.g., a pipe), two dimensional (e.g., a plate), three dimensional (e.g., a volume), or some intermediate (fractional) dimension. The nature of flow and transport dimension can be influenced strongly by roughness and mineralization, particularly at
fracture intersections, and can also be a function of both fluid pressure and in situ stress magnitude and orientation.
Relative Permeability Curves
Relative permeability curves (also referred to as characteristic curves) describe permeability to individual fluids as a function of the phase saturation. They are used to parameterize multiphase flow and transport behavior. Distinct and appropriate relative permeability curves are dependent on multiple geologic and geomechanical properties (e.g., fracture roughness, infilling, and persistence). Multiphase relationships need to be described frequently through correlation and probabilistic relationships with other fracture characteristics. Caution is necessary, however, when using relative permeability relationships because the use of such pressure-saturation relationships in fractures poses challenges. These relationships were developed for continuum representation of porous media and are not necessarily directly relevant to discrete features such as fractures.
Rock Matrix Parameterization
Rock matrix parameterization is dependent on determining which fractures should be considered part of the effective rock mass and the matrix block geometries for fracture-rock matrix interaction calculations. Properties such as permeability (and its function of phase saturation), total porosity, and pore size distributions for intact rock can be characterized generally through laboratory and other testing methods (see Chapter 5). The model needs to define matrix porosity and pore size distributions to evaluate mobile and immobile contaminant transport, and mitigation and remediation activity designs. As discussed in Chapter 3, pore size is important in filtering and biological processes. Anisotropy of rock matrix parameters should always be assumed and at least approximately quantified based on available geologic and geophysical data.
Boundary Condition Parameterization
Hydrostructural model boundary conditions describe hydrogeologic features such as water levels, faults, end of aquifers, and recharge or discharge zones that affect local hydraulic conditions. Boundary conditions describe the characteristics assigned to the boundaries or edges of a numerical model that represent regional hydraulic behaviors, as well as the parameters assigned at the beginning of a numerical model test to represent background conditions. Flow, pressure, and concentration boundary conditions are important parts of a hydrostructural model—particularly as they help define dominant discrete features. Standard techniques for defining boundary conditions in hydrogeology (Bear, 1972) also are applicable to fractured rock. Boundary conditions, however, can be discontinuous. For example, a fracture connected to a fault may provide a constant head boundary condition, while an adjacent fracture connected to a pumping well may provide a flow rate boundary condition. The details of the hydrostructual model and its connectivity thus create a spatially varying hydraulic distance to defined boundary conditions. Sensitivity studies can be used to quantify the effect of alternative boundary condition definitions. In many cases, larger discrete features may themselves be the boundary conditions—for example, a sealing fault may be a no-flow boundary condition.
Boundary condition development is complex, and there is much literature on the topic (e.g., Anderson and Woessner, 1992; Barnett et al., 2012). Initial and boundary conditions—including
some combination of pressure, flow, temperature, geochemical, and biological boundary states—affect model predictions and design simulations strongly. Consequently, the level of effort for establishing appropriate boundary conditions should be comparable to that used for assessing material properties.
Boundary conditions to be considered include
- large, conductive faults that act like constant pressure boundaries;
- flow barrier faults that act like no-flow boundary conditions;
- recharge and discharge boundaries identified from geochemical signatures (e.g., highest total dissolved solids at discharge);
- topographic surface water divides (used in continuum models, but are not frequently groundwater divides in fractured and compartmentalized rock systems); and
- very low fracture intensity zones that may be modeled as no-fracture-flow boundaries.
SCOPING CALCULATIONS TO ASSESS AND REFINE MODELS AND UNCERTAINTIES
Analyses of contaminant mitigation and remediation engineering solutions need not be complicated—even at complex fractured rock sites. In some cases, a simpler analysis is more helpful than the results of a more complex regional groundwater model. It is worthwhile and cost effective to use simple scoping calculations to determine the level of analysis necessary for a given engineering problem. At sites where there is evidence that a single discrete feature dominates transport, an approximate analysis to identify and quantify the pathway flow and storage properties can forecast transport accurately. At other sites, the most important process may be the exchange of contaminant between mobile fracture porosities and immobile matrix porosities. In those cases, an analytical solution of mass exchange might provide an excellent indication of the timescales of site remediation. Many important questions related to fractured rock characterization and monitoring can be expressed in terms of the probability that a given array of boreholes will intersect features of a given geometry. Equations of stochastic geometry and geometric probability (e.g., Chiu et al., 2013) can answer such questions without groundwater flow or transport modeling.
A good analysis approach is to determine first the questions to answer, then to develop efficient and effective ways to answer them using the site hydrostructural model and knowledge about available simple and complex analysis tools. Simplified analyses can determine the most sensitive assumptions and parameters in the system so that the most resources can be devoted to quantify them. Other simple analyses can determine which porosities, geometric issues, and processes to evaluate. These include geometric issues such as probabilities of intersection and the probability of connections between two points in space (i.e., between plume and discharge locations). Simplified analyses can also indicate which hydrostructural and microstructural model features may be upscaled and the consequences of upscaling. Simple analyses can be used to determine the range of results expected from more complex and expensive modeling, and whether the resulting reduced uncertainties justify the additional modeling costs.
ANALYSIS TOOLS TO INFORM MODELING
Numerous commercial and public analysis tools are available to facilitate fractured rock modeling, but it is beyond the scope of this report to assess all of these tools. A list of publicly available
analytical tools has been produced by the Integrated Groundwater Modeling Center at the Colorado School of Mines.^{1} The choice of tool warrants careful consideration by the practitioner, and it behooves the engineer to understand the modeling approaches incorporated into software design. Geometric calculation tools should be used to answer geometry-related questions (e.g., determining matrix block diffusion distances, or the probability that a discrete fracture pathway exists between a contaminant source and compliance boundary). Solute transport calculation tools should include the ability to analyze matrix diffusion processes, surface sorption, and rock matrix sorption without representing fractured rock as a single porosity system. Well-test interpretation tools should determine the flow dimension effects of fracture network connectivity.
Correlations and Parameterizations
The assumption of independence between fracture geometric properties (e.g., fracture size, orientation, and depth) and hydrodynamic properties (e.g., transmissivity, storativity, aperture) can have significant impact on forcasts made with such models. Independence should only be assumed where there is clear evidence that geometry and hydrodynamic properties are independent. More research is needed to better define the functional relationships between fracture geometric and hydrodyanmic properties. In the absense of site-specific evidence to the contrary, it is better to assume a correlation than to assume there is no functional relationship. Klimczak et al. (2010) and Hjerne et al. (2010) provide example correlations between fracture size and flow and transport properties.
It is also important to emphasize that fracture intensity should not be expected to show clear correlation to effective rock mass permeability (see, for example, research at the Apache Leap Research Site in Arizona; Chen et al., 2000). This is due to the wide variability in fracture transmissivity (generally many orders of magnitude); a few high-conductivity fractures contribute more to rock mass hydraulic response than a much larger intensity of less transmissive fractures. Hypothesis testing approaches to determine whether such functional relationships can be expressed in terms of in situ stress and geomechanical properties can be useful. Available correlations allow estimation of fracture transport parameters (e.g., transport aperture, transmissivity, and storativity) from geologic and geometric measurements (see Figure 4.3).
Pathway Geometry and Topology
Many fractured rock contaminant hydrogeology issues can be thought of as topological (connectivity) rather than as hydrodynamic problems. As such, they can be addressed through analysis of geometry. Geometry and topology can answer questions such as the probability that a given monitoring well configuration will detect contaminants, about whether an injection/withdrawal pair will have sufficient reactive surface area to address a certain volume of contaminated rock, or whether contaminant flowing through a fracture network in an aquitard layer will transport contaminant to an aquifer.
___________________
^{1} See http://igwmc.mines.edu/software/category_list.html#STA (accessed December 11, 2013).
Scoping calculations based on graph theory analysis of the geometry of conductive features and flow barriers (e.g., Bodin et al., 2007) can address problems of topology and the probability of connections among features (i.e., Sahimi, 1994). If there is concern about the probability of rapid breakthrough at a well 100 meters from a contaminated location through a single high-transmissivity fracture, then the fracture size (radius) distribution can be used to calculate the probability of breakthrough for a single fracture of that size or greater. Such approaches allow transport prediction and remediation and monitoring design to be carried out with simpler, one- or two-dimensional methods.
Immobile Zone Interactions
Mass transfer rates and storage volumes for mobile and immobile porosities can be estimated with one-dimensional approximate calculation approaches. The approaches can therefore be used to better estimate where contaminants are likely stored, to estimate the timescales of various remediation processes, and to assess whether more complex and resource-intensive modeling and analysis is warranted (e.g., single porosity/single permeability, single permeability/dual porosity, dual porosity/dual permeability, or discrete fracture network [DFN] approaches).
Sensitivity to Coupled Mechanisms
Coupled processes can have a significant effect on flow and transport. This has been a particularly important issue in the petroleum industry, where fracture closure in response to decreased reservoir pressure with production has been found to reduce flow rates and oil production (e.g., Hansford and Fisher, 2009). Baseline scoping calculations need to account for such coupling—at least in a generic way—where significant changes over time in temperatures, pressures, or chemistry occur. Failure to include these correlations in analyses can degrade the ability to forcast flow and transport response over time, increasing both error and uncertainty.
However, coupled process calculations can complicate analyses, and they require more time or effort available for a particular application. In such cases, as a minimum, the effect of ignoring these couplings needs to be estimated. Better analysis will result when simplified scoping calculations are used to estimate the sensitivity of contaminant transport, retention, and remediation in fractured rock to coupled thermal, hydrodynamic, mechanical, biological, and geochemical processes. Estimates of the change in storativity and transmissivity as a function of stress can be used to explore fracture stiffness-aperture-transmissivity relationships. Equations such as those provided by Jiang et al. (2009) can be used. Revised storativity and stiffness values can then be applied in flow and contaminant transport calculations to estimate responses of flow and transport to changes in effective stress magnitude and orientation.
Type-Curve Derivative Approaches
The theoretical response of models to constant-pressure and other tests can be calculated, and the expression of the results can be used for practical applications. For example, calculation results can be generalized on dimensionless graphs (e.g., dimensionless pressure versus dimensionless time) for type curve analysis (e.g., Gringarten, 1987; Barker, 1988) and used to assess fracture network connectivity and dual-porosity behaviors in a simple, practical manner. Derivative approaches using type curves to analyze dynamic flow and transport data (e.g., Black et al., 1987; Doe, 1991; Acuna et al., 1995; Illman and Neuman, 2000, 2001) are powerful and underused tools. They can be used to examine fracture network connectivity and boundary condition issues that, in turn, control the flow system, interaction between mobile and immobile porosities, and boundary conditions such as flow barrier faults and connected high-porosity aquifers.
When dynamic flow or transport information is available, type-curve derivative analysis can provide significant insights. Derivative analysis and multi-rate type curve analysis both can be used to better understand flow geometry and to help constrain fracture hydrodynamic properties (transmissivity and storativity), boundary conditions, and dual-porosity or dual-permeability behavior. This approach can be used to determine when contaminant transport can be modeled as a single porosity, pipe (linear) flow process, a dual porosity (fracture and matrix) system, or a radial (confined aquifer; Doe, 1991) system. The relative importance of geometry (e.g., discrete flow channels) and matrix interactions can be assessed by reviewing hydrodynamic data using flow dimension and derivative approaches. Better forecast of contaminant transport, storage, monitoring, and remediation will be possible with firm understanding of dimension of flow and the effect of flow boundaries.
Simplified Transport Approaches
There are several multiple-immobile-zone solute transport approaches that provide quick estimates of solute transport, even in complex multiple-porosity systems and can therefore be used in sensitivity studies to examine assumptions, to prioritize data collection, and to scope methods for mitigation. These include Laplace Transform Galerkin methods (Sudicky and McLaren, 1992), Lagrangian methods (Cvetkovic and Dagan, 1994), and continuous time random walk analysis approaches (Bijeljic et al., 2011).
TYPES OF NUMERICAL MODELS
Discrete fracture pathways and mass transfer between mobile (generally fractures) and immobile (generally matrix) porosities can be modeled with most hydrogeologic modeling codes. The level of accuracy and effort required for modeling will depend on how well the selected numerical approach can be parameterized and discretized, given financial, technical, and schedule constraints (Selroos et al., 2002). In general, the more explicitly the fractures and rock matrix blocks are represented, the easier it will be to relate the hydrogeologic model to the underlying geology, and less dramatic simplification will be required. Modeling approaches that explicitly represent fractured rock systems include the DFN approach (Long, 1983; Robinson, 1984; Dershowitz, 1985), the hybrid equivalent porous media (EPM)/DFN approach (Neuman, 1987; Bordas, 2005; Dershowitz, 2006; Illman, 2014), and channel network approaches (Watanabe et al., 1997; Bruines, 2003). Reviews of these models and concepts can be found in the literature and include Evans et al. (1987), Haneberg et al. (1999), Faybinshenko et al. (2000), Berkowitz (2002), Selroos et al. (2002), Dershowitz et al. (2004), Ijiri et al. (2009), and Illman (2014). Neuman (2005) provides a broad overview of trends, prospects, and challenges when quantifying flow and transport in the fractured rock environment. Table 4.4 lists numerical modeling approaches for fractured rock. DFN models are available through commercial, government, and academic sources, but are not used as commonly as continuum models even though they offer advantages for modeling fractured rock system.
There is considerable room to improve the ability of even state-of-the-art numerical models to explicitly describe flow and transport in fractured low permeability systems. The Used Fuel Disposition Campaign Disposal Research and Development Roadmap (Used Fuel Disposition Campaign, 2012) highlighted several ways to characterize and model features in a natural system to determine suitability for waste disposal. That document called for improved numerical modeling tools that represent fractures as discrete features, such as DFN and hybrid DFN/EPM multiple-porosity methods. Single-porosity EPM-type models are only applicable for fractured rock systems when the consequences of severe simplification of the system have been addressed. Loosely coupled approaches to examine geochemical, thermal, geomechanical, and biological interactions can be combined with explicit DFN approaches for scoping calculations before resorting to EPM approaches.
Other effective tools are widely available but not widely used, including model-independent parameter estimation and uncertainty analysis (e.g., Sandia National Laboratories, 2010; Austria et al., 2015) and calibration-constrained Monte Carlo analysis of highly parameterized models using subspace techniques (e.g., Tonkin and Doherty, 2009). Given the number of choices, it is critical that the modeler is familiar enough with the range of available approaches to ensure that simplifications made, given project resource constraints and available data, are those that provide the most insight, the lowest uncertainty bounds, and the best possible forecasts, and therefore the most defensible decisions.
Tools for model optimization, parameter estimation, sensitivity analysis, and uncertainty quantification frameworks can be used to improve the predictive power and quantification of uncertainty from all hydrogeologic models. Some of these have been adapted specifically for fractured rock (e.g., Vesselinov et al., 2001; Reed and Hetland, 2002; Finsterle, 2004; Doherty et al., 2010; Vesselinov and Harp, 2012).
TABLE 4.4 Approaches to Numerical Modeling
Numerical Model Approach | Model Description | Sample References |
---|---|---|
Discrete Fracture Network (DFN) | Fractures are explicitly represented, usually as meshes of 2D elements (polygons) in 3D space or intersecting pipe elements in 2D space to approximate the discrete pathways and heterogeneous connectivity. Heterogeneity on fracture planes can be represented by assigning varying material properties to elements. The rock matrix between fractures may be represented with an EPM approach, using analytical solutions or 1D approximations. Efficiency of this approach arises from the ability to model a 3D volume with 2D elements. This is balanced by computational demands where large numbers of fractures are included and by the complex geometry of fracture intersections. | Long et al., 1982 Dershowitz et al., 1998, 2004 Jackson et al., 2000 Hartley et al., 2006 |
Equivalent Porous Media (EPM) | Fractured rock represented as a mesh containing cells (elements) with flow and transport properties defined by a continuity flow equation (Bear, 1972). Each element in the mesh connects to its neighbors so flow is continuous. Transport pathways are determined by flow field, not discrete features. The effect of fractures on flow and transport may be represented using anisotropic elements and by varying element properties (e.g., using stochastic or conditioned approaches). When gridded finely enough, variation of material properties can be similar to that of a DFN model and can provide the heterogeneous connectivity seen in fractured rock. EPM model geometry is the same as that of elements used (i.e., 2D EPM models use 2D elements; 3D EPM models use 3D elements). Multiple techniques can be used (e.g., “non-neighbor connections”) to approximate connectivity. EPM approaches may discretize both fractures and rock matrix into a single mesh, or may implement separate interacting meshes for fracture and matrix permeabilities (e.g., dual continuum). | Bear, 1972 |
Channel Network (CN) | A DFN approach in which 1D elements are used to implement the fracture geometry. The rock matrix may be implemented as a CN model or using an EPM approach. This can improve computational efficiency, but can cause problems in converting fracture networks to pipe networks. The approach is attractive where fracture intersections or solution-enhanced permeability are a significant factor in defining flow and transport pathways. CN approaches defined along transport pathways can be computationally efficient for modeling reactive transport (sorption, reactions) and mass transport between rock matrix and channels via idealized channel/matrix interaction models. | Moreno and Neretnieks, 1993 Eker and Akin, 2006 |
Stochastic Continuum (SC) | Assigns permeability and porosity to the elements of a DFN, EPM, or CN model using, for example, statistical, geostatistical, or geological approaches. Most DFN models and many EPM models incorporate at least some element of SC to account for uncertainty and variability. SC approaches are useful particularly where a single realization of the random fields will not represent a specific site but a well-defined suite of realizations can display a range of behavior representative expected at a specific site. | Neuman, 1987 Blessent et al., 2009 |
Hybrid Discrete Fracture Network/Equivalent Porous Media (DFN/EPM | The hybrid approach generally breaks the model volume into multiple domains in which the DFN approach is used in some volumes and the EPM approach is used in others. These nested models require a numerical strategy to communicate flow and transport between DFN and EPM domains. | Dershowitz, 2006 Bonneau et al., 2013 |
Numerical Model Approach | Model Description | Sample References |
---|---|---|
Dual Continuum/Discrete Fracture Network/Equivalent Porous Media (DC/DFN/EPM) | This approach is similar to the hybrid DFN/EPM, but a dual continuum approach is used in the continuum portion of the model. | Pruess, 1985 Tatomir et al., 2011 |
UPSCALING AND MODEL SIMPLIFICATION
The upscaling process simplifies model properties for application in larger-scaled models. Even single-porosity EPM approaches are implicitly upscaled from the actual geology. The properties of grains and pores are simplified to parameters such as “hydraulic conductivity.” In fractured rock models, fracture plane properties need to be simplified. There are limitations to upscaling in any application that may affect the result of analyses.
By upscaling some properties, numerical models can focus on other hydraulically important hydrostructural properties and features. When using an EPM, consideration must be given to how the selected parameters relate to the geometric and hydrodynamic properties of the underlying fractures and rock matrix, and how simplification to a single-porosity EPM model might avoid misleading results. At a minimum, the sensitivity of model results could be evaluated within the range of model parameters and boundary conditions consistent with the hydrostructural model.
Powerful, easy-to-use, and simple approaches are available to upscale fracture hydrostructural models and DFN to heterogeneous and anisotropic EPM models. These make it possible to respect the spatial pattern, heterogeneity, connectivity, and fracture/rock matrix interaction of a realistic hydrostructural conceptual model in a more conventional EPM model. Examples include the Oda fabric tensor (Will et al., 2005; Wang et al., 2008; Matthäi and Nick, 2009) and numerical permeameters (Lu and Kwicklis, 2012). These methods allow large-scale, standardized EPM approaches to preserve some of the discrete pathway and mobile-immobile zone mass transfer behavior.
The upscaling process for fractured rock requires definition of a “cut-off” scale between fractures that are considered explicitly and those that are included within the effective rock matrix porosity and permeability. Careful consideration is necessary of the cutoff scales between fractures upscaled to rock mass porosity and anisotropic permeability, fractures represented as a stochastic population, and fractures represented deterministically.
Hydrostructural models for many fractured rock sites include heterogeneously connected fractures and interacting blocks of rock matrix with high storage porosity. The important rock connectivity (discrete pathway) and rock matrix storage features need to be preserved during upscaling. If this is not possible, then supplemental analyses need to account for discrete pathway and dual porosity effects. Upscaling approaches are frequently better at retaining spatial heterogeneity and anisotropic permeability behaviors than at preserving underlying discrete pathways and connectivities. Mobile-immobile zone interactions, discrete flow pathways, and flow barrier features need to be preserved.
The shift of focus from properties of specific fractures to an equivalent continuum hydraulic conductivity, storativity, and porosity is inevitably a process of upscaling. Upscaled continuum approaches generally assume that flow streamlines correspond to transport pathways. This assumption can result in severe miscalculation of groundwater travel time, effective plume dispersion, and long-tail retention. In many cases, discrete pathways need to be considered separately.
This can be achieved in some of those cases by supplementing upscaled models with simplified scoping calculations for discrete pathways.
Another potentially inaccurate but common assumption is that hydraulic tests in fractured rock can directly define hydraulic properties such as hydraulic conductivity and porosity that in turn can be used as grid elements properties in a continuum model. Hydraulic tests apply hydraulic loads to specific discrete fractures and fracture network that happen to intersect the well. Better results are more likely if fracture and rock matrix hydraulic properties are derived from the well test. Properties can then be extrapolated to the population of fractures not intersected by the well, and then used to upscale from the fracture network to continuum grid hydraulic properties.
The process of translating a hydrostructural model to a numerical model often requires gridding. Grid assumptions and the details of upscaling affect the results of upscaled numerical modeling. Hydraulic and transport tests on fractures and matrix blocks can provide information that can be applied with reasonable confidence to other fractures and matrix blocks of the same scale. However, if the numerical model grid resolution is too large to resolve model elements at the scale of the hydraulic tests, then sort of averaging of properties is required to represent those elements numerically. Significant error in flow and transport forecasting can result if limitations imposed by model gridding and upscaling are ignored. Simulation of processes influenced by smaller-scale features (e.g., capillary driven flow, filtering, coupled processes, and reactive surface areas) can also be represented inaccurately or lost. Sensitivity analyses of alternate upscaling and gridding methods can quantify those errors and limitations, allowing decisions about the reasonableness of errors introduced by specific methods. If sensitivity studies are not feasible, as a minimum, then model results need to be used with an explicit caveat that the model error and uncertainty have not been quantified.
NUMERICAL MODEL ANALYSIS
Using the simplest model that matches data and observations for fractured rock systems can lead to incorrect results. Measurements often are not available in the quantities, distributions, or scales needed, and data may not reflect important processes that occur. Formal numerical model analysis applied to fractured rock systems quantitatively relates the hydrostructural conceptualization, the equations used to represent the conceptualization, the algorithms used to solve the equations in a numerical model, the model input parameters, and model outputs. It is important not only that numerical models adequately solve the equations that represent the hydrostructural conceptualization, but also that tools that quantify model uncertainties are applied so that decisions informed by numerical models can be informed by quantified uncertainties. Even when extensive data, observations, and experimental results are available, modeling approaches need to balance the value of more detailed numerical models (e.g., of multiple coupled processes) against the uncertainty reduction.
Numerical Testing of Alternate Conceptual Models
Mathematical methods are often designed for parameters with ranges (minimum, maximum) or distributions (normal, log-normal). When change in model output with respect to changes in input are continuous (i.e., derivative of outputs to inputs are continuous), mathematical tools that take advantage of that continuity can be applied. However, in some situations, alternative conceptual models need to be used that result in discrete, discontinuous changes in model behavior. For
example, there may be insufficient data to determine whether the fracture intensity is above or below the percolation threshold. If the percolation threshold is not reached, then the system does not have long-range connectivity of intersecting fractures and flow does not occur. Above the percolation threshold, there is long-range connectivity, and a pressure gradient will drive flow through permeable fractures. The system essentially goes through a phase change with distinctly different responses (no flow versus flow). In another example, subsurface stratigraphy may not be constrained well enough to know whether a fault offsets a highly conductive unit enough to cut flow through that unit.
Scenarios such as these can be conceptualized, and their consistency with data and model results can be tested (e.g., via pump tests with observation wells on both sides of the fault). In the best cases, alternate conceptualizations can be eliminated based on comparison of the data that is sensitive to the alternative numerical model output. In any case, however, even when a conceptualization is rejected, it is important to document and retain analyses that lead to acceptance or rejection of alternate conceptual models to inform future efforts. It is also important to communicate with sponsors, clients, regulators, and stakeholders how to visualize and represent such scenarios. More information on testing of conceptual models can be found in the literature (e.g., Neuman, 2003; Neuman and Wieranga, 2003).
Model Calibration
A productive approach to develop reasonable model inputs is to identify plausible input ranges based on hydrostructural model conceptualizations, run forward models within an automated parameter estimation framework, identify model output that is most reasonably and quantitatively consistent with observations, and identify the input parameter combinations that yield the most reasonable outputs. However, uninformed use of this approach can lead to incorrect interpretations. Many combinations of input parameters may yield output that is a good fit to observed conditions. A good fit, therefore, does not necessarily ensure the model correctly represents the flow and transport conditions. Data used typically for model calibration include heads, flow rates, and geochemical measurements under natural and pumped conditions. Examples of hydrologic model calibrations include study of potential sites for high-level nuclear waste repositories such as Yucca Mountain (Zyvoloski, 2003) and Äspö (Mazurek et al., 2003).
Sensitivity Analysis
In cases where numerical model output changes continuously with variations in model input values, quantitative measures of model sensitivity can be calculated by what is, essentially, computing the partial derivatives of model output with respect to input. This provides a quantitative measure of sensitivity of output to input. However, complex numerical models with large numbers of degrees of freedom often make it impractical to carry out an analysis of each of the systematic changes made to model input. Doing so would involve a full forward run of the numerical model for each input parameter change—an impractical brute-force approach that is inappropriate for all but the simplest of models. Instead, the sensitivity can be quantified by using sampling methods (e.g., Monte Carlo, screening methods) to statistically estimate the response of the multi-dimensional space of model behavior. This approach samples changes in inputs from a known or defined distribution and characterizes output variability in terms of observed output distributions. Again,
in models with many input parameters, sampling the high-dimensional space can be computationally expensive, but methods and tools to do this efficiently are available (e.g., Yeh, 1986).
As described earlier, quantifying model uncertainty can have considerable value. Analysis of initial scoping calculations often reveals ways to simplify numerical models when, for example, some numerical model inputs have little effect on output. Quantitatively analyzing the sensitivity—or in this case, the insensitivity—of a model can be powerful because it can reveal what additional data may or may not constrain a model parameter and decrease model uncertainty. This type of insight is more valuable than that resulting from a single numerical model calculation.
Quantitative sensitivity analysis can help to quantify correlated parameters and inform on ways to simplify models, as well as the decisions related to the types of observations and data needed to constrain the model. Sensitivity analysis is affected by issues of non-uniqueness. For example, many alternative combinations of fracture geometry and hydrodynamic properties may explain equally well the measured responses. Additionally, parameters studied in sensitivity analysis may not be those directly affecting observations. For example, many hydrodynamic responses are a function of hydraulic diffusivity (the ratio of transmissivity to storativity), such that sensitivity analysis of transmissivity alone is a poorly posed problem.
Uncertainty Quantification
Model uncertainty can arise from (a) simplifications necessary to implement a hydrostructural model in a numerical model, (b) limitations of the understanding or implementation of physical-chemical-biological processes in the model, (c) errors in the numerical model implementation, and (d) limitations in the match between measurements and model results. Quantifying these errors and uncertainties ensures appropriate application of numerical models and the conclusions drawn from them.
Just as examining the sensitivity of a model to input parameters can guide subsequent model development, quantifying the uncertainties that arise during the modeling process can guide interpretation of model output. Numerical models and the equations they solve are inevitably imperfect representations of the real world that may be biased or inadequate and result in different categories of uncertainty. Quantitative measures need to be developed to determine the likelihood of remediation success (e.g., that remediation strategies reduce concentrations to below maximum contaminant level; that the mass fraction of a contaminant is accessible given a specific remediation method; that the amount of time for contaminant to reach a compliance boundary is reasonable; and that uncertainty can be reduced through collection of additional, potentially expensive, data).
Structural uncertainties can be large and are often the most difficult to quantify and estimate when modeling fractured rock. Parameter (numerical model input) uncertainty is always present and often can be quantified. Because model inputs are usually based on limited observations, they must be interpolated and extrapolated. Formal (statistical data analysis) or heuristic (expert opinion) methods can be used to place bounds and characterize the possible parameter distribution (normal, log normal). This can be done through statistical analysis of the data if enough data are available. When data are limited, expert opinion with sound, documented explanations for the choices can be applied.
Uncertainty from all sources propagates forward through a numerical model in forward uncertainty quantification. Uncertainties in model outputs can be stated in terms of uncertainty in model inputs. In all but the simplest numerical models, the expense (e.g., time, computational intensity)
to evaluate a single forward model makes it prohibitive to sample the full high-dimensional parameter space. Sampling methods (e.g., Monte Carlo, Markov chain Monte Carlo, adaptive sampling) need to be used judiciously to reduce the total number of forward solutions of the numerical model constructed to represent the physical system.
ANALYSIS AND RESOURCES
Numerical models can improve decision processes by forecasting system performance. The processes of data collection, hydrostructural conceptualization, numerical model setup and execution, and model analysis is not a rigid and prescriptive one-pass process. A successful analysis process requires continuous reassessment of data, model assumptions, numerical model validity, and analysis of model output relative to data and observations (see Chapter 7). An application of some of the most advanced methods and practices to siting an underground nuclear waste repository in fractured rock is at the Forsmark site in Sweden. Decades of developing methodologies and algorithms to implement numerical models are described by Hartley and Joyce (2013).
Numerical models and associated software implementations will continue to grow in response to need, but it is often a slow process to bring the latest technologies into mainstream practice. Time and resource availability are often inconsistent with the complex and resource intensive workflows required to implement the latest technologies, and the interdisciplinary nature of the analysis of flow and transport in fractured rock will always make that analysis challenging. Many analysis tools are not widely available or are difficult to use and, as such, may not be accessible to everyone. The use of oversimplified or inappropriate numerical models, however, can lead to inadequate or incorrect results. Hydrogeologists that receive proper training that allows them to develop appropriate conceptual models and then apply inexpensive initial scoping calculations to bound system behavior will be able to make more informed decisions regarding data needs and more complex modeling.