Fracture-matrix Interaction Research Articles

Streamline simulation of counter-current imbibition in naturally fractured reservoirs

Describing fluid transport in naturally fractured reservoirs entails additional challenge because of the complicated physics arising from matrix–fracture interactions. In this paper, the streamline-based simulation is generalized to describe fluid transport in naturally fractured reservoirs through a dual-media approach. The fractures and matrix are treated as separate continua that are connected through a transfer function, as in conventional finite difference simulators for modeling fractured systems. The transfer functions that describe fluid exchange between the fracture and matrix system can be implemented easily within the framework of the current single-porosity streamline models. In particular, the streamline time of flight concept is utilized to develop a general dual-porosity, dual-permeability system of equations for water injection in naturally fractured reservoirs. The saturations equations are solved using an operator splitting approach that involves ‘convection’ along streamline followed ‘matrix–fracture’ exchange calculations on the grid. The proposed formulation reduces to the commonly used dual-porosity model when the flow in the matrix is considered negligible. For modeling matrix–fracture interactions, two different transfer functions are examined: an empirical transfer function (ETF) and a conventional transfer function (CTF). The ETF allows for analytical solution of the saturation equation for dual porosity systems and is used to validate numerical implementation. Results obtained using the CTF are compared with a commercial finite-difference simulator (ECLIPSE) for waterflooding in five-spot and nine-spot patterns. The streamline approach shows close agreement in terms of recovery histories and saturation profiles with a marked reduction in numerical dispersion and grid orientation effects. Finally, an examination of the scaling behavior of the computation time indicates that the streamline approach is likely to result in significant savings for large-scale field applications.

Unsaturated Flow through a Small Fracture-Matrix Network: Part 2. Uncertainty in Modeling Flow Processes

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

Unsaturated Flow through a Small Fracture–Matrix Network: Part 2. Uncertainty in Modeling Flow Processes

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

The fracture flow equation and its perturbation solution

This work derives the fracture flow equation from the two‐dimensional steady form of the Navier‐Stokes equation. Asymptotic solutions are obtained whereby the perturbation parameter is the ratio of the mean width over the length of the fracture segment. The perturbation expansion can handle arbitrary variation of the fracture walls as long as the dominant velocity is in the longitudinal direction. The effect of the matrix‐fracture interaction is also taken into account by allowing leakage through the fracture walls. The perturbation solution is used to obtain an estimate of the flow rate and the fracture transmissivity as well as the velocity and the pressure distribution in fractures of various geometries. The analysis covers eight different configurations of fracture geometry including linear and curvilinear variation as well as sinusoidal variation in the top and bottom walls with varying horizontal alignment and roughness wavelengths. The zero‐order solution yields the Reynolds lubrication approximation, and the higher‐order equations provide a correction term to the flow rate in terms of the roughness frequency and the Reynolds number. For sinusoidal and linear walls, the mathematical analysis shows that the zero‐order flow rate could be expressed in terms of the maximum to minimum width ratio. For equal widths at both ends of the fracture, the first‐order correction is zero. For sinusoidal fractures, the flow rate decreases with increasing Reynolds number and with increasing roughness amplitude and frequency. The effect of leakage is to create a nonuniform flow distribution in the fracture that deviates significantly from the flow rate estimate for impermeable walls. The derived flow expressions can provide a more reliable tool for flow and transport predictions in fractured domain.

A reaction-transport model for calcite precipitation and evaluation of infiltration fluxes in unsaturated fractured rock

The percolation flux in the unsaturated zone (UZ) is an important parameter addressed in site characterization and flow and transport modeling of the potential nuclear-waste repository at Yucca Mountain, NV, USA. The US Geological Survey (USGS) has documented hydrogenic calcite abundances in fractures and lithophysal cavities at Yucca Mountain to provide constraints on percolation fluxes in the UZ. The purpose of this study was to investigate the relationship between percolation flux and measured calcite abundances using reactive transport modeling. Our model considers the following essential factors affecting calcite precipitation: (1) infiltration, (2) the ambient geothermal gradient, (3) gaseous CO 2 diffusive transport and partitioning in liquid and gas phases, (4) fracture–matrix interaction for water flow and chemical constituents, and (5) water–rock interaction. Over a bounding range of 2–20 mm/year infiltration rate, the simulated calcite distributions capture the trend in calcite abundances measured in a deep borehole (WT-24) by the USGS. The calcite is found predominantly in fractures in the welded tuffs, which is also captured by the model simulations. Simulations showed that from about 2 to 6 mm/year, the amount of calcite precipitated in the welded Topopah Spring tuff is sensitive to the infiltration rate. This dependence decreases at higher infiltration rates owing to a modification of the geothermal gradient from the increased percolation flux. The model also confirms the conceptual model for higher percolation fluxes in the fractures compared to the matrix in the welded units, and the significant contribution of Ca from water–rock interaction. This study indicates that reactive transport modeling of calcite deposition can yield important constraints on the unsaturated zone infiltration-percolation flux and provide useful insight into processes such as fracture–matrix interaction as well as conditions and parameters controlling calcite deposition.

On the physics of unstable infiltration, seepage, and gravity drainage in partially saturated tuffs

To improve understanding of the physics of dynamic instabilities in unsaturated flow processes within the Paintbrush nonwelded unit (PTn) and the middle nonlithophysal portion of the Topopah Spring welded tuff unit (TSw) of Yucca Mountain, we analyzed data from a series of infiltration tests carried out at two sites (Alcove 4 and Alcove 6) in the Exploratory Studies Facility (ESF), using analytical and empirical functions. The analysis of infiltration rates measured at both sites showed three temporal scales of infiltration rate: (1) a macro-scale trend of overall decreasing flow, (2) a meso-scale trend of fast and slow motion exhibiting three-stage variations of the flow rate (decreasing, increasing, and [again] decreasing flow rate, as observed in soils in the presence of entrapped air), and (3) micro-scale (high frequency) fluctuations. Infiltration tests in the nonwelded unit at Alcove 4 indicate that this unit may effectively dampen episodic fast infiltration events; however, well-known Kostyakov, Horton, and Philip equations do not satisfactorily describe the observed trends of the infiltration rate. Instead, a Weibull distribution model can most accurately describe experimentally determined time trends of the infiltration rate. Infiltration tests in highly permeable, fractured, welded tuff at Alcove 6 indicate that the infiltration rate exhibits pulsation, which may have been caused by multiple threshold effects and water–air redistribution between fractures and matrix. The empirical relationships between the extrinsic seepage from fractures, matrix imbibition, and gravity drainage versus the infiltration rate, as well as scaling and self-similarity for the leading edge of the water front are the hallmark of the nonlinear dynamic processes in water flow under episodic infiltration through fractured tuff. Based on the analysis of experimental data, we propose a conceptual model of a dynamic fracture flow and fracture–matrix interaction in fractured tuff, incorporating the time-dependent processes of water redistribution in the fracture–matrix system.

Evolution of the unsaturated zone testing at Yucca Mountain

The evaluation of the Yucca Mountain site has evolved from intensive surface-based investigations in the early 1980s to current focus on testing in underground drifts. Different periods of site characterization activities and prominent issues concerning the unsaturated zone (UZ) are summarized. Data collection activities have evolved from mapping of faults and fractures to estimation of percolation through tuff layers, and to quantification of seepage into drifts. Evaluation of discrete flow paths in drifts has led to fracture–matrix interaction and matrix diffusion tests over different scales. The effects of tuff interfaces and local faults are evaluated in fractured-welded and porous-nonwelded units. Mobilization of matrix water and redistribution of moisture are measured in thermal tests. Lessons learned from underground tests are used to focus on processes needed for additional quantification. Migration through the drift shadow zone and liquid flow through faults are two important issues that have evolved from current knowledge.

Reactive transport model for the ambient unsaturated hydrogeochemical system at Yucca Mountain, Nevada

To assist a technical review of a potential application for a geologic repository, a reactive transport model is presented for the ambient hydrogeochemical system at Yucca Mountain (YM). The model simulates two-phase, non-isothermal, advective and diffusive flow and transport through one-dimensional matrix and fracture continua (dual permeability) containing ten kinetically reactive hydrostratigraphic layers. We developed a thermodynamic interpretation of ambient groundwater chemistry at YM from limited analytical data and then combined it with insights from experimental and site-specific data to craft a conceptual model describing thermodynamic and kinetic relationships in the mineral/glass–water–gas system. The model was calibrated by adjusting uncertain thermodynamic and kinetic parameters to reflect observed trends at YM. The following calibration criteria were used to ensure that model predictions are consistent with hydrochemical and petrologic data from YM: (1) calculated multicomponent matrix pore water compositions do not vary significantly from reinterpreted analytical data; (2) simulated variations in matrix silica concentrations with depth are bounded by the observed analytical range; and (3) feldspars and glass dissolve while clays and calcite precipitate. The largest changes in calculated species concentrations occur in response to glass dissolution. Simulated matrix groundwater compositions at the depth of the potential repository are largely inherited from percolating waters rather than controlled by in situ chemical reactions. Predictions about groundwater flow pathways in the ambient YM system are sensitive to assumptions regarding percolation flux and fracture–matrix interactions. Confidence in complex reactive transport models of the thermally perturbed YM repository system depends on successful representation of ambient system conditions.

Comment on “Determination of particle transfer in random walk particle methods for fractured porous media” by H. H. Liu et al.

Comment on “Determination of particle transfer in random walk particle methods for fractured porous media” by H. H. Liu et al.

Experimental and analytical study of multidimensional imbibition in fractured porous media

Capillary imbibition is an important mechanism during water injection and aquifer influx in fractured porous media. Better understanding of matrix–fracture interaction and imbibition in general is needed to model effectively these processes. Using an X-ray computerized tomography (CT) scanner and a novel CT-compatible core holder, a series of experiments to study air and oil expulsion from rock samples by capillary imbibition of water in a three-dimensional (3-D) geometry were conducted. The air–water system was useful because a relatively large number of experiments could be conducted to delineate physical processes. Different injection rates and fracture apertures were utilized. Two different fracture flow regimes were identified. The “filling fracture” regime shows a plane source that grows in length due to relatively slow water flow through fractures. In the second regime, the “instantly filled fracture” regime, the time to fill the fracture is much less than the imbibition time and the imbibition performance scales as the square root of time. In the former regime, the mass of water imbibed scales linearly with time. A new analytical model is proposed for filling fractures incorporating implicit matrix/fracture coupling. Good agreement is found between experiments and calculation. This analytic coupling was obtained by solving the saturation diffusion equation with appropriate initial and boundary conditions. The solution provides the location of the wetting phase front in the fracture and the saturation distribution in the matrix. The solution is analogous to that obtained by Marx and Langenheim [Trans. AIME 216 (1959) 312] for the areal extent of an equivalent heated zone in thermal recovery methods. Analogous terms among flow and heat transfer in porous media were found and are also presented.

Migration of a water pulse through fractured porous media

Contaminant transport from waste-disposal sites is strongly affected by the presence of fractures and the degree of fracture–matrix interaction. Characterization of potential contaminant plumes at such sites is difficult, both experimentally and numerically. Simulations of water flow through fractured rock were performed to examine the penetration depth of a large pulse of water entering such a system. Construction water traced with lithium bromide was released during the excavation of a tunnel at Yucca Mountain, Nevada, which is located in an unsaturated fractured tuff formation. Modeling of construction-water migration is qualitatively compared with bromide-to-chloride ratio (Br/Cl) data for pore-water salts extracted from drillcores. The influences of local heterogeneities in the fracture network and variations in hydrogeologic parameters were examined by sensitivity analyses and Monte Carlo simulations. The simulation results are qualitatively consistent with the observed Br/Cl signals, although these data may only indicate a minimum penetration depth, and water may have migrated farther through the fracture network.

Scaling capillary imbibition during static thermal and dynamic fracture flow conditions

To incorporate the effect of temperature, previously reported scaling formulations were modified considering the wettability, change in oil properties, and thermal expansion. The results of static imbibition experiments conducted using different types of heavy-oil samples at different temperatures up to 90 °C were used to verify the formulation. The modified scaling formulations responded well if the thermal expansion effect is included, except in high-temperature crude oil imbibition cases. The results obtained through water injection experiments in artificially fractured core samples with matrix permeabilities ranging from 300 to 0.075 mD were used to verify the numerical model of core scale displacements. Then the convergence constant in Aranofsky's abstract recovery–time relationship was correlated to an imbibition group consisting of the injection rate and matrix properties. Three different formulations were presented and tested. Good agreement with the experimental observations could be achieved when the wettability and permeability of the matrix and the flow rate in fracture is included in the formulations. When the matrix size is included in the regression analysis, not all correlation types yielded a good agreement with the experimental data. This study will provide an insight into the scaling of the laboratory scale experiments to the reservoir conditions and constitute a base to define matrix fracture interaction in numerical models dealing with fractured reservoir modeling.

Development of an Integrated Reservoir Model for a Naturally Fractured Volcanic Reservoir in China

Summary This paper presents the development of an integrated, multidiscipline reservoir model for dynamic flow simulation and performance prediction of a geologically complex, naturally fractured volcanic reservoir in the Shang 741 Block of the Shengli field in China. A static geological model integrates lithological information, petrophysics, fracture analysis, and stochastic fracture network modeling with Formation MicroImage (FMI) log data and advanced 3D seismic interpretations. Effective fracture permeability, fracture-matrix interaction, reservoir compartmentalization, and flow transmissibility of conductive faults are obtained by matching various dynamic data. As a result of synergy and multiple iterations among various disciplines, a history-matched dynamic reservoirsimulation model capable of future performance prediction for optimum asset management is constructed.

Tritium and 36Cl as constraints on fast fracture flow and percolation flux in the unsaturated zone at Yucca Mountain

An analysis of tritium and 36Cl data collected at Yucca Mountain, Nevada suggests that fracture flow may occur at high velocities through the thick unsaturated zone. The mechanisms and extent of this “fast flow” in fractures at Yucca Mountain are investigated with data analysis, mixing models and several one-dimensional modeling scenarios. The model results and data analysis provide evidence substantiating the weeps model [Gauthier, J.H., Wilson, M.L., Lauffer, F.C., 1992. Proceedings of the Third Annual International High-level Radioactive Waste Management Conference, vol. 1, Las Vegas, NV. American Nuclear Society, La Grange Park, IL, pp. 891–989] and suggest that fast flow in fractures with minimal fracture–matrix interaction may comprise a substantial proportion of the total infiltration through Yucca Mountain. Mixing calculations suggest that bomb-pulse tritium measurements, in general, represent the tail end of travel times for thermonuclear-test-era (bomb-pulse) infiltration. The data analysis shows that bomb-pulse tritium and 36Cl measurements are correlated with discrete features such as horizontal fractures and areas where lateral flow may occur. The results presented here imply that fast flow in fractures may be ubiquitous at Yucca Mountain, occurring when transient infiltration (storms) generates flow in the connected fracture network.

Field tracer-transport tests in unsaturated fractured tuff

This paper presents the results of a field investigation in the unsaturated, fractured welded tuff within the Exploratory Studies Facility (ESF) at Yucca Mountain, NV. This investigation included a series of tests during which tracer-laced water was released into a high-permeability zone within a horizontal injection borehole. The tracer concentration was monitored in the seepage collected in an excavated slot about 1.6 m below the borehole. Results showed significant variability in the hydrologic response of fractures and the matrix. Analyses of the breakthrough curves suggest that flow and transport pathways are dynamic, rather than fixed, and related to liquid-release rates. Under high release rates, fractures acted as the predominant flow pathways, with limited fracture–matrix interaction. Under low release rates, fracture flow was comparatively less dominant, with a noticeable contribution from matrix flow. Observations of tracer concentrations rebounding in seepage water, following an interruption of flow, provided evidence of mass exchange between the fast-flowing fractures and slow- or non-flowing regions. The tests also showed the applicability of fluorinated benzoate tracers in situations where multiple tracers of similar physical properties are warranted.

Modeling multiphase non-isothermal fluid flow and reactive geochemical transport in variably saturated fractured rocks: 2. Applications to supergene copper enrichment and hydrothermal flows

Reactive fluid flow and geochemical transport in unsaturated fractured rocks has been of increasing interest to investigators in the areas of geo- and environmental- sciences. To test geochemical hypotheses based on petrologic observation and to predict geochemical reactions that occur through a complex dynamic interplay of physical and chemical processes, we use the methods presented in a companion paper (part 1, this issue p. 16–33) to investigate two problems: (1) supergene copper enrichment in unsaturated-saturated media and (2) predicted effects of thermohydrology on geochemistry during the Drift Scale Heater Test at the Yucca Mountain potential nuclear waste repository, Nevada. Through these two examples we address the importance of the following issues on geochemical processes: (1) participation of gas phase in transport and reaction, (2) interactions between fractures and rock matrix for water and chemical constituents, (3) heat effects on fluid flow and reaction properties and processes. In the supergene enrichment system, oxygen gas diffusion from the land surface through fractured rock promotes the alteration of the primary sulfide minerals and the subsequent deposition of secondary minerals. Modeling of the large-scale heater test shows effects of fracture-matrix interaction, heat-driven vaporizing fluid flow, and CO₂ degassing on mineral alteration patterns. The two examples also serve as a demonstration of our methods for reactive transport in variably saturated fractured rocks.

Numerical simulation of non-isothermal multiphase tracer transport in heterogeneous fractured porous media

We have developed a new numerical method for modeling tracer and radionuclide transport through heterogeneous fractured rocks in a non-isothermal multiphase system. The model formulation incorporates a full hydrodynamic dispersion tensor, based on three-dimensional velocity fields with a regular or irregular grid in a heterogeneous geological medium. Two different weighting schemes are proposed for spatial interpolation of three-dimensional velocity fields and concentration gradients to evaluate the mass flux by dispersion and diffusion of a tracer or radionuclide. The fracture–matrix interactions are handled using a dual-continua approach, such as the double- or multiple-porosity, or dual-permeability. This new method has been implemented into a multidimensional numerical code to simulate processes of tracer or radionuclide transport in non-isothermal, three-dimensional, multiphase, porous/fractured subsurface systems. We have verified this new transport-modeling approach by comparing its results to the analytical solution of a parallel-fracture transport problem. In addition, we use published laboratory results and the particle-tracking scheme to further validate the model. Finally, we present two examples of field applications to demonstrate the use of the proposed methodology for modeling tracer and radionuclide transport in unsaturated fractured rocks.

Calibrating hydrogeologic parameters for the 3-D site-scale unsaturated zone model of Yucca Mountain, Nevada

An important issue in the evaluation of the unsaturated zone (UZ) at Yucca Mountain, Nevada is the calibration of the parameters used in the 3-D site-scale numerical flow model. The hydrogeologic parameters are calibrated using an inversion code (ITOUGH2) to fit measured core sample and in situ data according to statistical criteria. The available data include saturations, water potentials, pneumatic pressure records, and pore water age estimates based on geochemical analyses. Five 1-D columns are extracted from the 3-D site-scale model and used in an inversion that matches data from the five boreholes simultaneously. This results in layer-averaged fracture and matrix permeabilities and van Genuchten fitting parameters. The columns are then refined and individual inversions are performed with each 1-D model to develop a spatially varying parameter set. A number of sensitivity studies are performed to evaluate the effects of modeling assumptions on the calibrated parameters. Several fracture-matrix interaction formulations are also investigated during the study. The results show that the calibrated parameter set concluded to be most consistent with the conceptual model provides a good fit to the ambient data. The study also shows that the choice of the fracture-matrix interaction formulation has a significant effect on the calibrated parameter set.

Density-dependent solute transport in discretely-fractured geologic media: is prediction possible?

The development of a dense solute plume in a fractured geologic medium can be highly irregular due to both the complexity of the fracture network as well as the presence of convection cells that may arise as a result of the density contrast between the invading solute and the ambient groundwater. A two-dimensional numerical model has been developed here to investigate density-dependent groundwater flow and solute transport in geologic materials that contain discrete fractures in order to examine some of the complex forms into which plumes can evolve, particularly with regard to fracture–matrix interactions. Results from simulations which involve parallel vertical fractures show that the evolution of the solute plume is affected by the development of convection cells in the porous matrix blocks between the vertical fractures. In a geologic medium containing a network of regularly spaced horizontal and vertical fractures, complex migration pathways can develop that are unexpected even though the geometry and interconnectivity of the fractures are known a priori. Downward solute migration can occur in some vertical fractures, while upward migration of less dense fluid can occur in others with transient circulation patterns developing in the intervening porous matrix. Because of the inherent uncertainty associated with fracture delineation, and because of the irregular nature of unstable dense plumes, deterministic prediction of dense-plume migration pathways and travel times in fractured geologic media will be subject to considerable uncertainty.

Permeability porosity relationships from numerical simulations of fluid flow

Numerical calculations of permeability are obtained from a lattice Boltzmann simulation of flow in simplified 2D porous media over a range of solid fractions. The evolution equation governing the flow behaviour incorporates the effect of porous medium geometry through the definition of solid density ns, a real number, at each node of the simulation grid. The results obtained for homogeneous media are compared to commonly used theoretical and empirical relationships relating rock properties to permeability. Behaviour consistent with a Kozeny‐Carman type relationship between porosity ϕ and permeability k is obtained for low to intermediate solid fractions. At high solid fractions the rapid decrease in k is consistent with a percolation process giving a power‐law relationship for ϕ and k. Both the critical porosity and power‐law exponent are in agreement with quoted values for the lattice geometry used. A comparison of the results for homogeneous media with k values, obtained by embedding a spanning planar fracture into the matrix, illustrates the importance of matrix‐fracture flow interactions. The results for this case are consistent with experimental observations and illustrate the difficulties involved in using simplified assumptions to predict permeability from porosity in fractured porous rock.