A new lumped-parameter model for flow in unsaturated dual-porosity media

A new lumped-parameter approach to simulating unsaturated flow processes in dual-porosity media such as fractured rocks or aggregated soils is presented. Fluid flow between the fracture network and the matrix blocks is described by a non-linear equation that relates the imbibition rate of the local difference in liquid-phase pressure between the fractures and the matrix blocks. Unlike a Warren-Root-type equation, this equation is accurate in both the early and late time regimes. The fracture/matrix interflow equation has been incorporated into an existing unsaturated flow simulator, to serve as a source/sink term for fracture gridblocks. Flow processes are then simulated using only fracture gridblocks in the computational grid. This new lumped-parameter approach has been tested on two problems involving transient flow in fractured/porous media, and compared with simulations performed using explicit discretisation of the matrix blocks. The new procedure seems to accurately simulate flow processes in unsaturated fractured rocks, and typically requires an order of magnitude less computational time than do simulations using fully-discretised matrix blocks.     

The modified Reynolds equation for non-wetting fluid flow through a rough-walled rock fracture

The no-slip boundary condition has usually been assumed to hold for the Reynolds equations (local cubic law) for fluid flow through rough-walled fractures. However, its validity for non-wetting fluid flow, such as prevails in fractured oil reservoirs, has been questioned. A series of experiments with a rough-walled fracture with mean aperture of 760 μm finds a higher flow rate for non-wetting fluid than wetting fluid. A modified Reynolds equation with a slip boundary condition is derived for non-wetting fluid flow through rough-walled fractures. Comparison of the modified Reynolds equation predictions with experimental results confirms that slip was a plausible explanation for a higher flow rate. The amount by which the flow rate for non-wetting fluids exceeds that for wetting fluids is found to depend highly on, and increase with, the degree to which the flowing fluid was non-wetted to thin immobile films on the surfaces. Numerical studies using the modified Reynolds equation indicate that the flow rate of non-wetting fluid became higher than that of wetting fluid as the roughness of the fracture increases. As the aperture becomes smaller, the flow rate ratio of non-wetting fluid to wetting fluid becomes large, leading to the endpoint relative permeability for the non-wetting fluid to exceed 1. The experimental and numerical studies clearly show that as the aperture of the fracture became less than a few hundred microns, the modified Reynolds equation with slip boundary conditions provides a better model for flow of a non-wetting fluid through rough-walled fractures.     

Two-region non-Darcian flow toward a well in a confined aquifer

We have derived an analytical solution for two-region flow toward a well in a confined aquifer based on a linearization method. The two-region flow includes Izbash non-Darcian flow near the well and Darcian flow in the rest of the aquifer. The wellbore storage is also considered. The type curves in the non-Darcian and Darcian flow domains are obtained by a numerical Laplace inversion method incorporated in MATLAB programs. We have compared our results with the one-region Darcian flow model (Theis). Our solutions agree with those of Sen [Sen Z. Type curves for two-region well flow. J Hydr Eng 1988;114(12):1461–84] which were obtained using the Boltzmann transform at late times for fully turbulent flow, while some difference has been found at early and moderate times. We have defined a dimensionless non-Darcian hydraulic conductivity term which is shown to be a key parameter for analyzing the two-region flow. A smaller dimensionless non-Darcian hydraulic conductivity results in a larger drawdown in the non-Darcian flow region at late times. However, the dimensionless non-Darcian hydraulic conductivity does not affect the slope of the dimensionless drawdown versus the logarithmic dimensionless time in the non-Darcian flow region at late times. The dimensionless non-Darcian hydraulic conductivity does not affect the late time drawdown in the Darcian flow region.     

Solutions for Non‐Darcian Flow to an Extended Well in Fractured Rock

We have investigated non‐Darcian flow to a vertical fracture represented as an extended well using a linearization procedure and a finite difference method in this study. Approximate analytical solutions have been obtained with and without the consideration of fracture storage based on the linearization procedure. A numerical solution for such a non‐Darcian flow case has also been obtained with a finite difference method. We have compared the numerical solution with the approximate analytical solutions obtained by the linearization method and the Boltzmann transform. The results indicate that the linearized solution agrees generally well with the numerical solution at late times, and underestimates the dimensionless drawdown at early times, no matter if the fracture storage is considered or not. When the fracture storage is excluded, the Boltzmann transform solution overestimates the dimensionless drawdown during the entire pumping period. The dimensionless drawdowns in the fracture with fracture storage for different values of dimensionless non‐Darcian hydraulic conductivity β approach the same asymptotic value at early times. A larger β value results in a smaller dimensionless drawdown in both the fracture and the aquifer when the fracture storage is included. The dimensionless drawdown is approximately proportional to the square root of the dimensionless time at late times.     

Weathering and erosion of fractured bedrock systems

We explore the contribution of fractures (joints) in controlling the rate of weathering advance for a low‐porosity rock by using methods of homogenization to create averaged weathering equations. The rate of advance of the weathering front can be expressed as the same rate observed in non‐fractured media (or in an individual block) divided by the volume fraction of non‐fractured blocks in the fractured parent material. In the model, the parent has fractures that are filled with a more porous material that contains only inert or completely weathered material. The low‐porosity rock weathers by reaction‐transport processes. As observed in field systems, the model shows that the weathering advance rate is greater for the fractured as compared to the analogous non‐fractured system because the volume fraction of blocks is < 1. The increase in advance rate is attributed both to the increase in weathered material that accompanies higher fracture density, and to the increase in exposure of surface of low‐porosity rock to reaction‐transport. For constant fracture aperture, the weathering advance rate increases when the fracture spacing decreases. Equations describing weathering advance rate are summarized in the ‘List of selected equations’. If erosion is imposed at a constant rate, the weathering systems with fracture‐bounded bedrock blocks attain a steady state. In the erosional transport‐limited regime, bedrock blocks no longer emerge at the air‐regolith boundary because they weather away. In the weathering‐limited (or kinetic) regime, blocks of various size become exhumed at the surface and the average size of these exposed blocks increases with the erosion rate. For convex hillslopes, the block size exposed at the surface increases downslope. This model can explain observations of exhumed rocks weathering in the Luquillo mountains of Puerto Rico. Published 2017. This article is a U.S. Government work and is in the public domain in the USA     

Steady state flow with hydraulic conductivity change around large diameter wells

Darcian flow law in aquifers assumes that the aquifer hydraulic conductivity is constant and the groundwater movement is due only to the piezometric level changes through hydraulic gradient. In practice, after the well development the aquifer just around the well has comparatively larger hydraulic conductivity and gradient. Patchy aquifer solutions in the literature consider sudden hydraulic conductivity changes with distance for the steady state flow. The change of transmissivity is demonstrated by the application of slope‐matching procedure to actual field data. It is the main purpose of this paper to derive simple analytical expressions for aquifer parameter evaluations with steadily decreasing hydraulic conductivity around the well. Spatial nonlinear hydraulic conductivity changes around a large‐diameter well within the depression cone of a confined aquifer are considered as exponentially decreasing functions of the radial distance. Copyright © 2010 John Wiley & Sons, Ltd.     

An equivalent medium model for wave simulation in fractured porous rocks

Seismic wave propagation in reservoir rocks is often strongly affected by fractures and micropores. Elastic properties of fractured reservoirs are studied using a fractured porous rock model, in which fractures are considered to be embedded in a homogeneous porous background. The paper presents an equivalent media model for fractured porous rocks. Fractures are described in a stress‐strain relationship in terms of fracture‐induced anisotropy. The equations of poroelasticity are used to describe the background porous matrix and the contents of the fractures are inserted into a matrix. Based on the fractured equivalent‐medium theory and Biot's equations of poroelasticity, two sets of porosity are considered in a constitutive equation. The porous matrix permeability and fracture permeability are analysed by using the continuum media seepage theory in equations of motion. We then design a fractured porous equivalent medium and derive the modified effective constants for low‐frequency elastic constants due to the presence of fractures. The expressions of elastic constants are concise and are directly related to the properties of the main porous matrix, the inserted fractures and the pore fluid. The phase velocity and attenuation of the fractured porous equivalent media are investigated based on this model. Numerical simulations are performed. We show that the fractures and pores strongly influence wave propagation, induce anisotropy and cause poroelastic behaviour in the wavefields. We observe that the presence of fractures gives rise to changes in phase velocity and attenuation, especially for the slow P‐wave in the direction parallel to the fracture plane.     

Nonlinear Flow in Karst Formations

The variation of effective hydraulic conductivity as a function of specific discharge in several 0.2-m and 0.3-m cubes of Key Largo Limestone was investigated. The experimental results closely match the Forchheimer equation. Defining the pore-size length scale in terms of Forchheimer parameters, it is demonstrated that significant deviations from Darcian flow will occur when the Reynolds number exceeds 0.11. A particular threshold model previously proposed for use in karstic formations does not show strong agreement with the data near the onset of nonlinear flow.     

Static and dynamic conceptual model of a complexly fractured crystalline rock aquifer

A conceptual model of anisotropic and dynamic permeability is developed from hydrogeologic and hydromechanical characterization of a foliated, complexly fractured, crystalline rock aquifer at Gates Pond, Berlin, Massachusetts. Methods of investigation include aquifer‐pumping tests, long‐term hydrologic monitoring, fracture characterization, downhole heat‐pulse flow meter measurements, in situ extensometer testing, and earth tide analysis. A static conceptual model is developed from observations of depth‐dependent and anisotropic permeability that effectively compartmentalizes the aquifer as a function of foliation intensity. Superimposed on the static model is dynamic permeability as a function of hydraulic head in which transient bulk aquifer transmissivity is proportional to changes in hydraulic head due to hydromechanical coupling. The dynamic permeability concept is built on observations that fracture aperture changes as a function of hydraulic head, as measured during in situ extensometer testing of individual fractures, and observed changes in bulk aquifer transmissivity as determined from earth tides during seasonal changes in hydraulic head, with higher transmissivity during periods of high hydraulic head, and lower transmissivity during periods of relatively lower hydraulic head. A final conceptual model is presented that captures both the static and dynamic properties of the aquifer. The workflow presented here demonstrates development of a conceptual framework for building numerical models of complexly fractured, foliated, crystalline rock aquifers that includes both a static model to describe the spatial distribution of permeability as a function of fracture type and foliation intensity and a dynamic model that describes how hydromechanical coupling impacts permeability magnitude as a function of hydraulic head fluctuation. This model captures important geologic controls on permeability magnitude, anisotropy, and transience and therefor offers potentially more reliable history matching and forecasts of different water management strategies, such as resource evaluation, well placement, permeability prediction, and evaluating remediation strategies.     

Experimental study of the effect of roughness and Reynolds number on fluid flow in rough‐walled single fractures: a check of local cubic law

Local cubic law (LCL) is one of the most commonly applied physical laws for flow in single fractures (SF) and fractured media. The foundation of LCL is Darcian flow. This experimental study examines if LCL is valid for flow in a single rough fracture and how the fracture roughness and Reynolds number (Re) affect flow. Similar to the Moody diagram for flow in pipes, a diagram for flow in a single rough fracture has been generated to relate the friction coefficient with Re and the roughness. Under the experimental condition of this study, flow appears to be substantially different from Darcian flow. The flow law of q ∝ eⁿJ^m appears to be valid for describing the flow scheme where q, e, and J are the unit width flux, the average aperture, and the hydraulic gradient. The value of the power index m is found to be around 0·83 ～ 0·98, less than what has been used in Darcian flow (m = 1). The power index n is around 11·2 and 13·0, much greater than the n value used in the LCL (n = 3), and it increases with the average velocity. The Moody type of diagram shows that the friction factor for flow in SFs is influenced by Re and the roughness. It decreases with Re when Re is small, and becomes less sensitive to Re when Re is large enough. It also increases with the roughness. Copyright © 2010 John Wiley & Sons, Ltd.     

An efficient numerical model for multicomponent compressible flow in fractured porous media

An efficient and accurate numerical model for multicomponent compressible single-phase flow in fractured media is presented. The discrete-fracture approach is used to model the fractures where the fracture entities are described explicitly in the computational domain. We use the concept of cross flow equilibrium in the fractures. This will allow large matrix elements in the neighborhood of the fractures and considerable speed up of the algorithm. We use an implicit finite volume (FV) scheme to solve the species mass balance equation in the fractures. This step avoids the use of Courant–Freidricks–Levy (CFL) condition and contributes to significant speed up of the code. The hybrid mixed finite element method (MFE) is used to solve for the velocity in both the matrix and the fractures coupled with the discontinuous Galerkin (DG) method to solve the species transport equations in the matrix. Four numerical examples are presented to demonstrate the robustness and efficiency of the proposed model. We show that the combination of the fracture cross-flow equilibrium and the implicit composition calculation in the fractures increase the computational speed 20–130 times in 2D. In 3D, one may expect even a higher computational efficiency.     

Near‐Well Nonlinear Flow Identified by Various Displacement Well Response Testing

The impact of nonlinear flow phenomena on well response tests is still not completely understood today. With the present paper, a set of 10 well response tests is investigated. The tests were conducted in a fractured Devonian limestone formation close to the western national border of Germany. The test design incorporates a packer as well as different solid cylinders to initiate a series of slug‐injection and slug‐withdrawal tests by various initial displacements. Nonlinear response characteristics were observed in the course of the tests, which cannot be explained by tubing‐controlled flow inside the cased well. The analysis shows that the nonlinear response characteristics are likely to be either formation controlled according to non‐Darcian flow developing in a high‐conductivity fracture compartment of the tested limestone formation or a consequence of a severe well inefficiency caused by some sort of screen clogging. This inference is obtained from analyzing the data by a nonlinear well response test model, which differentiates between wellbore internal hydraulic head losses and a generalized rate‐dependent skin effect accounting for nonlinear flow processes in the vicinity of the well. The potential of identifying near‐well nonlinear flow by various displacement well response testing may indicate this methodology to be a valuable complement to modern high‐resolution borehole imaging techniques used when characterizing fractured reservoirs and the tightness of fractured reservoir cap rocks.     

A practical model for fluid flow in discrete-fracture porous media by using the numerical manifold method

In this study, a numerical manifold method (NMM) model is developed to analyze flow in porous media with discrete fractures in a non-conforming mesh. This new model is based on a two-cover-mesh system with a uniform triangular mathematical mesh and boundary/fracture-divided physical covers, where local independent cover functions are defined. The overlapping parts of the physical covers are elements where the global approximation is defined by the weighted average of the physical cover functions. The mesh is generated by a tree-cutting algorithm. A new model that does not introduce additional degrees of freedom (DOF) for fractures was developed for fluid flow in fractures. The fracture surfaces that belong to different physical covers are used to represent fracture flow in the direction of the fractures. In the direction normal to the fractures, the fracture surfaces are regarded as Dirichlet boundaries to exchange fluxes with the rock matrix. Furthermore, fractures that intersect with Dirichlet or Neumann boundaries are considered. Simulation examples are designed to verify the efficiency of the tree-cutting algorithm, the calculation's independency from the mesh orientation, and accuracy when modeling porous media that contain fractures with multiple intersections and different orientations. The simulation results show good agreement with available analytical solutions. Finally, the model is applied to cases that involve nine intersecting fractures and a complex network of 100 fractures, both of which achieve reasonable results. The new model is very practical for modeling flow in fractured porous media, even for a geometrically complex fracture network with large hydraulic conductivity contrasts between fractures and the matrix.     

Non-Fickian mass transport in fractured porous media

The paper provides an introduction to fundamental concepts of mathematical modeling of mass transport in fractured porous heterogeneous rocks. Keeping aside many important factors that can affect mass transport in subsurface, our main concern is the multi-scale character of the rock formation, which is constituted by porous domains dissected by the network of fractures. Taking into account the well-documented fact that porous rocks can be considered as a fractal medium and assuming that sizes of pores vary significantly (i.e. have different characteristic scales), the fractional-order differential equations that model the anomalous diffusive mass transport in such type of domains are derived and justified analytically. Analytical solutions of some particular problems of anomalous diffusion in the fractal media of various geometries are obtained. Extending this approach to more complex situation when diffusion is accompanied by advection, solute transport in a fractured porous medium is modeled by the advection-dispersion equation with fractional time derivative. In the case of confined fractured porous aquifer, accounting for anomalous non-Fickian diffusion in the surrounding rock mass, the adopted approach leads to introduction of an additional fractional time derivative in the equation for solute transport. The closed-form solutions for concentrations in the aquifer and surrounding rocks are obtained for the arbitrary time-dependent source of contamination located in the inlet of the aquifer. Based on these solutions, different regimes of contamination of the aquifers with different physical properties can be readily modeled and analyzed.     

A geomorphology‐based integrated stream–aquifer interaction model for semi‐gauged catchments

Many of the existing stream–aquifer interaction models available in the literature are very complex with limited applicability in semi‐gauged and ungauged catchments. In this study, to estimate the influent and effluent subsurface water fluxes under limited geo‐hydrometeorological data availability conditions, a simple stream–aquifer interaction model, namely, the variable parameter McCarthy–Muskingum (VPMM) hillslope‐storage Boussinesq (hsB) model, has been developed. This novel model couples the VPMM streamflow transport with the hsB groundwater flow transport modules in online mode. In this integrated model, the surface water–groundwater flux exchange process is modelled by the Darcian approach with the variable hydraulic heads between the river stage and groundwater table accounting for the rainfall forcing. Considering the exchange fluxes in the hyporheic zone and lateral overland flow contribution, this approach is field tested in a typical 48‐km stretch of the Brahmani River in eastern India to simulate the streamflow and its depth with the minimum Nash–Sutcliffe efficiency of 94% and 88%; the maximum root mean square error of 134 m³/s and 0.35 m; and the minimum index of agreement of 98% and 97%, respectively. This modelling approach could be very well utilized in data‐scarce world‐river basins to estimate the stream–aquifer exchange flux due to rainfall forcings.     

A general analytical solution for flow to a single horizontal well by Fourier and Laplace transforms

The objective of this paper is to present an analytical solution for describing the head distribution in an unconfined aquifer with a single pumping horizontal well parallel to a fully penetrating stream. The Laplace-domain solution is developed by applying Fourier sine, Fourier and Laplace transforms to the governing equation as well as the associated initial and boundary conditions. The time-domain solution is obtained after taking the inverse Laplace transform along with the Bromwich integral method and inverse Fourier and Fourier sine transforms. The upper boundary condition of the aquifer is represented by the free surface equation in which the second-order slope terms are neglected. Based on the solution and Darcy’s law, the equation representing the stream depletion rate is then derived. The solution can simulate head distributions in an aquifer infinitely extending in horizontal direction if the well is located far away from the stream. In addition, the solution can also simulate head distributions in confined aquifers if specific yield is set zero. It is shown that the solution can be applied practically to evaluate flow to a horizontal well.     

Matrix–fracture transfer shape factor for modeling flow of a compressible fluid in dual-porosity media

The matrix–fracture transfer shape factor is one of the important parameters in the modeling of fluid flow in fractured porous media using a dual-porosity concept. Warren and Root [36] introduced the dual-porosity concept and suggested a relation for the shape factor. There is no general relationship for determining the shape factor for a single-phase flow of slightly compressible fluids. Therefore, different studies reported different values for this parameter, as an input into the flow models. Several investigations have been reported on the shape factor for slightly compressible fluids. However, the case of compressible fluids has not been investigated in the past. The focus of this study is, therefore, to find the shape factor for the single-phase flow of compressible fluids (gases) in fractured porous media. In this study, a model for the determination of the shape factor for compressible fluids is presented; and, the solution of nonlinear gas diffusivity equation is used to derive the shape factor. The integral method and the method of moments are used to solve the nonlinear governing equation by considering the pressure dependency of the viscosity and isothermal compressibility of the fluid. The approximate semi-analytical model for the shape factor presented in this study is verified using single-porosity, fine-grid, numerical simulations. The dependency of the shape factor on the gas specific gravity, pressure and temperature are also investigated. The theoretical analysis presented improves our understanding of fluid flow in fractured porous media. In addition, the developed matrix–fracture transfer shape factor can be used as an input for modeling flow of compressible fluids in dual-porosity systems, such as naturally fractured gas reservoirs, coalbed methane reservoirs and fractured tight gas reservoirs.     

Unified fractional differential approach for transient interporosity flow in naturally fractured media

A unified approach to modeling flows of slightly compressible fluids through naturally fractured media is presented. The unified fractional differential model is derived by combining the flow at micro scale for matrix blocks and macro scale for fractures, using the transient interporosity flow behavior at the interface between matrix blocks and fractures. The derived model is able to unify existing transient interporosity flow models formulated for different shapes of matrix blocks in any medium dimensions. The model is formulated in the form of a fractional order partial differential equation that involves Caputo derivative of order 1/2 with respect to time. Explicit solutions for the unified model are derived for different axisymmetrical spatial domains using Hankel or Hankel–Weber finite or infinite transforms. Comparisons between the predictions of the unified model and those obtained from existing transient interporosity flow models for matrix blocks in the form of slabs, spheres and cylinders are presented. It is shown that the unified fractional derivative model leads to solutions that are very close to those of transient interporosity flow models for fracture-dominant and transitional fracture-to-matrix dominant flow regimes. An analysis of the results of the unified model reveals that the pressure varies linearly with the logarithm of time for different flow regimes, with half slope for the transitional fracture-to-matrix dominant flow regime vs. the fracture and matrix dominant flow regimes. In addition, a new re-scaling that involves the characteristic length in the form of matrix block volume to surface area ratio is derived for the transient interporosity flow models for matrix blocks of different shapes. It is shown that the re-scaled transient interporosity flow models are governed by two dimensionless parameters Θ and Λ compared to only one dimensionless parameter Θ for the unified model. It is shown that the solutions of the transient interporosity flow models for different shapes of matrix blocks are almost identical for the re-scaled variables. Furthermore, the driving parameters for solution behavior are identified based on asymptotic approximations for different flow regimes. It is found that the matrix diffusion and the matrix area-to-volume ratio affect the solution behavior only for the transitional fracture-to-matrix dominant flow regime, that the capacitance ratio affects the solution behavior only for transitional and matrix dominant flow regimes and that the fracture diffusion is involved in all three flow regimes. Similar identification of the driving parameters is also presented in the re-scaled case.     

Understanding the Geometry of Connected Fracture Flow with Multiperiod Oscillatory Hydraulic Tests

An understanding of the spatial and hydraulic properties of fast preferential flow pathways in the subsurface is necessary in applications ranging from contaminant fate and transport modeling to design of energy extraction systems. One method for the characterization of fracture properties over interwellbore scales is Multiperiod Oscillatory Hydraulic (MOH) testing, in which the aquifer response to oscillatory pressure stimulations is observed. MOH tests were conducted on isolated intervals of wells in siliciclastic and carbonate aquifers in southern Wisconsin. The goal was to characterize the spatial properties of discrete fractures over interwellbore scales. MOH tests were conducted on two discrete fractured intervals intersecting two boreholes at one field site, and a nest of three piezometers at another field site. Fracture diffusivity estimates were obtained using analytical solutions that relate diffusivity to observed phase lag and amplitude decay. In addition, MOH tests were used to investigate the spatial extent of flow using different conceptual models of fracture geometry. Results indicated that fracture geometry at both field sites can be approximated by permeable two‐dimensional fracture planes, oriented near‐horizontally at one site, and near‐vertically at the other. The technique used on MOH field data to characterize fracture geometry shows promise in revealing fracture network characteristics important to groundwater flow and transport.     

Temporal Changes in VOC Discharge to Surface Water from a Fractured Rock Aquifer During Well Installation and Operation, Greenville, South Carolina

Analysis of the vapor in passive vapor samplers retrieved from a streambed in fractured rock terrain implied that volatile organic carbon (VOC) discharge from ground water to surface water substantially increased following installation of a contaminant recovery well using air rotary drilling. The air rotary technique forced air into the aquifer near the stream. The injection produced an upward hydraulic gradient that appears to have transported water and contaminants from deeper parts of the aquifer through fractures into shallow parts of the aquifer. Once in the shallow flow regime, the contamination was transported to the stream, where it discharged during the next several weeks following well installation. After the recovery well was activated and began continuously pumping contaminated ground water to a treatment facility, the VOC concentrations in the stream bottom passive vapor samplers decreased to below detectable concentrations, suggesting that the withdrawal had captured the contaminated ground water that previously had discharged to the stream.