Heterogeneous Permeability Fields Research Articles

Effect of permeability variations on solute transport in highly heterogeneous porous media

The effect of aquifer heterogeneity on flow and solute transport in two-dimensional isotropic porous media was analyzed using the Monte Carlo method. The two-dimensional logarithmic permeability (ln K) was assumed to be a non-stationary random field with its increments being a truncated fractional Lévy motion (fLm). The permeability fields were generated using the modified successive random additions (SRA) algorithm code SRA3DC [ http://www.iamg.org/CGEditor/index.htm]. The velocity and concentration fields were computed respectively for two-dimensional flow and transport with a pulse input using the finite difference codes of MODFLOW 2000 and MT3DMS. Two fLm control parameters, namely the width parameter ( C) and the Lévy index ( α), were varied systematically to examine their effect on the resulting permeability, flow velocity and concentration fields. We also computed the first- and second-spatial moments, the dilution index, as well as the breakthrough curves at different control planes with the corresponding concentration fields. In addition, the derived breakthrough curves were fitted using the continuous time random walk (CTRW) and the traditional advection-dispersion equation (ADE). Results indicated that larger C and smaller α both led to more heterogeneous permeability and velocity fields. The Lévy-stable distribution of increments in ln K resulted in a Lévy-stable distribution of increments in logarithm of the velocity (ln v). Both larger C and smaller α created sharper leading edges and wider tailing edges of solute plumes. Furthermore, a relatively larger amount of solute still remained in the domain after a relatively longer time transport for smaller α values. The dilution indices were smaller than unity and increased as C increased and α decreased. The solute plume and its second-spatial moments increased as C increased and α decreased, while the first-spatial moments of the solute plume were independent of C and α values. The longitudinal macrodispersivity was scale-dependent and increased as a power law function of time. Increasing C and decreasing α both resulted in an increase in longitudinal macrodispersivity. The transport in such highly heterogeneous media was slightly non-Gaussian with its derived breakthrough curves being slightly better fitted by the CTRW than the ADE, especially in the early arrivals and late-time tails.

Non-Fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media

We investigate the temporal scaling properties of mixing in heterogeneous permeability fields with variances ranging from very small ( σ ln K 2 = 0.01 ) to very large ( σ ln K 2 = 9 ) . We quantify mixing by the scalar dissipation rate, which we estimate over a large range of temporal scales. For an initial pulse line injection, we find that moderate and strong heterogeneity induce anomalous temporal scaling of the scalar dissipation rate, which we call non-Fickian mixing. This effect is particularly relevant for upscaling reactive transport as it implies a non-Fickian scaling of reactive transport. Although spreading and mixing are intimately coupled, we find that their scaling properties are not directly related in general. In the non-Fickian mixing regime, the temporal scaling of the scalar dissipation rate depends on the complex spatial distribution of the concentration field that generates transverse mixing. For times larger than the characteristic diffusion time associated with one permeability field correlation length, the heterogeneity of concentration in the plume is attenuated and progressively erased by diffusion. Thus, at large times, the temporal scaling of mixing and spreading can be related through a simple analytical expression.

Experimental study of the relation between the permeability of kaolinite and its deformation at micro and macro scale

We present an experimental study of relations between hydraulic properties of clayey rocks and their deformation at bulk and micro scale. The experimentation is based on triaxial deformation of kaolinitic core samples, on permeability measurements, and on several methods of analysis of the microstructure. To produce the samples, remoulded kaolinite material is consolidated by axial loading up to 20 MPa in an oedometric cell. Then, in a triaxial cell, an undrained stress path is applied in order to induce shearing damage in the sample. The permeability of the sample is measured by fluid flow across the sample. Before and after these hydro-mechanical tests, the evolution and the homogeneity of the pore space and moreover the texture of the material are analysed by X-ray tomography, environmental scanning electron microscopy and thermic neutron diffraction. Contrary to expectations, the permeability at the sample scale does not vary significantly, although quite large deformation is reached. The final deformation is heterogeneous, with well-developed shear bands. It appears that the development of compacting and dilating zones generates a heterogeneous permeability field, which results in a negligible variation of the bulk permeability at the sample scale.

Stochastic simulation of uranium migration at the Hanford 300 Area

This work focuses on the quantification of groundwater flow and subsequent U(VI) transport uncertainty due to heterogeneity in the sediment permeability at the Hanford 300 Area. U(VI) migration at the site is simulated with multiple realizations of stochastically-generated high resolution permeability fields and comparisons are made of cumulative water and U(VI) flux to the Columbia River. The massively parallel reactive flow and transport code PFLOTRAN is employed utilizing 40,960 processor cores on DOE's petascale Jaguar supercomputer to simultaneously execute 10 transient, variably-saturated groundwater flow and U(VI) transport simulations within 3D heterogeneous permeability fields using the code's multi-realization simulation capability. Simulation results demonstrate that the cumulative U(VI) flux to the Columbia River is less responsive to fine scale heterogeneity in permeability and more sensitive to the distribution of permeability within the river hyporheic zone and mean permeability of larger-scale geologic structures at the site.

Simulation of the subsurface mobility of carbon nanoparticles at the field scale

The production and use of nanomaterials will inevitably lead to their disposal in the natural environment. To assess the risk that these materials pose to human and ecosystem health an understanding of their mobility and ultimate fate is essential. To date, however, relatively little research has been conducted on the fate of nanoparticles in subsurface systems. In this study the subsurface mobility of two carbon nanoparticles: nano-fullerenes (nC 60) and multi-walled carbon nanotubes (MWCNTs) is assessed. A two-dimensional finite element model was used to simulate the movement of these nanoparticles under a range of hydrologic and geological conditions, including a heterogeneous permeability field. The numerical model is based on colloid filtration theory (CFT) with a maximum retention capacity term. For the conditions evaluated the carbon nanotubes are much more mobile than nC 60 due to the smaller collector efficiency associated with carbon nanotubes. However, the mobility of nC 60 increased significantly when a maximum retention capacity term was included in the model. Model results also demonstrate that, for the systems examined, nanoparticles were predicted to be less mobile in heterogeneous systems compared to the homogeneous systems with the same average hydraulic properties.

Transport behavior of coupled continuous-time random walks

The origin of anomalous or non-Fickian transport in disordered media is the broad spectrum of transition rates intrinsic to these systems. A system that contains within it heterogeneities over multiple length scales is geological formations. The continuous time random walk (CTRW) framework, which has been demonstrated to be an effective means to model non-Fickian transport features in these systems and to have predictive capacities, has at its core this full spectrum represented as a joint probability density psi(s,t) of random space time displacements (s,t) . Transport in a random fracture network (RFN) has been calculated with a coupled psi(s,t) and has subsequently been shown to be approximated well by a decoupled form psi(s,t)=F(s)psi(t) . The latter form has been used extensively to model non-Fickian transport in conjunction with a velocity distribution Phi(xi),xi identical with 1v, where v is the velocity magnitude. The power-law behavior of psi(t) proportional to (-1-beta), which determines non-Fickian transport, derives from the large xi dependence of Phi(xi) . In this study we use numerical CTRW simulations to explore the expanded transport phenomena derived from a coupled psi(s,t) . Specifically, we introduce the features of a power-law dependence in the s distribution with different Phi(xi) distributions (including a constant v) coupled by t=s(xi) . Unlike Lévy flights in this coupled scenario the spatial moments of the plumes are well defined. The shapes of the plumes depend on the entire Phi(xi) distribution, i.e., both small and large xi dependence; there is a competition between long displacements (which depend on the small xi dependence) and large time events (which depend on a power law for large xi). These features give rise to an enhanced range of transport behavior with a broader scope of applications, e.g., to correlated migrations in a RFN and in heterogeneous permeability fields. The approximation to the decoupled case is investigated as a function of the nature of the s distribution.

Operator-Based Multiscale Method for Compressible Flow

Summary Multiscale methods have been developed for accurate and efficient numerical solution of flow problems in large-scale heterogeneous reservoirs. A scalable and extendible Operator-Based Multiscale Method (OBMM) is described here. OBMM is cast as a general algebraic framework. It is natural and convenient to incorporate more physics in OBMM for multiscale computation. In OBMM, two operators are constructed: prolongation and restriction. The prolongation operator is constructed by assembling the multiscale basis functions. The specific form of the restriction operator depends on the coarse-scale discretization formulation (e.g., finite-volume or finite-element). The coarse-scale pressure equation is obtained algebraically by applying the prolongation and restriction operators to the fine-scale flow equations. Solving the coarse-scale equation results in a high-quality coarse-scale pressure. The fine-scale pressure can be reconstructed by applying the prolongation operator to the coarse-scale pressure. A conservative fine-scale velocity field is then reconstructed to solve the transport (saturation) equation. We describe the OBMM approach for multiscale modeling of compressible multiphase flow. We show that extension from incompressible to compressible flows is straightforward. No special treatment for compressibility is required. The efficiency of multi-scale formulations over standard fine-scale methods is retained by OBMM. The accuracy of OBMM is demonstrated using several numerical examples including a challenging depletion problem in a strongly heterogeneous permeability field (SPE 10).

A multilevel multiscale mimetic (M 3) method for two-phase flows in porous media

We describe a multilevel multiscale mimetic (M 3) method for solving two-phase flow (water and oil) in a heterogeneous reservoir. The governing equations are the elliptic equation for the reservoir pressure and the hyperbolic equation for the water saturation. On each time step, we first solve the pressure equation and then use the computed flux in an explicit upwind finite volume method to update the saturation. To reduce the computational cost, the pressure equation is solved on a much coarser grid than the saturation equation. The coarse-grid pressure discretization captures the influence of multiple scales via the subgrid modeling technique for single-phase flow recently proposed in [Yu. A. Kuznetsov. Mixed finite element method for diffusion equations on polygonal meshes with mixed cells. J. Numer. Math., 14 (4) (2006) 305–315; V. Gvozdev. discretization of the diffusion and Maxwell equations on polyhedral meshes. Technical Report Ph.D. Thesis, University of Houston, 2007; Yu. Kuznetsov. Mixed finite element methods on polyhedral meshes for diffusion equations, in: Computational Modeling with PDEs in Science and Engineering, Springer-Verlag, Berlin, in press]. We extend significantly the applicability of this technique by developing a new robust and efficient method for estimating the flux coarsening parameters. Specifically, with this advance the M 3 method can handle full permeability tensors and general coarsening strategies, which may generate polygonal meshes on the coarse grid. These problem dependent coarsening parameters also play a critical role in the interpolation of the flux, and hence, in the advection of saturation for two-phase flow. Numerical experiments for two-phase flow in highly heterogeneous permeability fields, including layer 68 of the SPE Tenth Comparative Solution Project, demonstrate that the M 3 method retains good accuracy for high coarsening factors in both directions, up to 64 for the considered models. Moreover, we demonstrate that with a simple and efficient temporal updating strategy for the coarsening parameters, we achieve accuracy comparable to the fine-scale solution, but at a fraction of the cost.

Phase separation and convection in heterogeneous porous media: Implications for seafloor hydrothermal systems

We examine experimentally the effect of phase separation of a binary fluid on convection and heat transfer in porous domains with heterogeneous permeability fields. In mid‐ocean ridge hydrothermal systems, a dense brine phase and a coexisting vapor phase form during the phase separation of seawater at supercritical conditions; this process strongly influences the chemistry of vents and the convective transfer of heat. Furthermore, convection and phase separation are likely to occur in oceanic crust that is highly heterogeneous in nature. Because of high temperatures and pressures required for phase separation in saline systems, a binary fluid [H2O‐2‐butoxyethanol], with a lower consolute point of 48.5°C, was used as a proxy for seawater. A saturated pseudo‐two‐dimensional porous domain containing high permeability [k] zones embedded in a lower permeability matrix was heated from below and cooled from above. As with homogeneous systems, chemical stratification develops in heterogeneous domains, with dense fluid accumulating in low‐velocity regions at the bottom of the porous matrix. Furthermore, in both types of porous media, the efficiency of thermal transport is reduced relative to single‐phase systems. Heterogeneities in the permeability field, however, can act to amplify these effects; in systems with vertical high‐k zones, the thickness of the dense‐phase layer increases in the low‐k regions, and thermal transport is suppressed even more. Also, the experiments show that regions of both high and low permeabilities can act as a store for brines and facilitate the formation of much thicker brine layers than are otherwise predicted by models based on homogeneous oceanic crust.

Characterization of NAPL Source Zone Architecture and Dissolution Kinetics in Heterogeneous Porous Media Using Magnetic Resonance Imaging

A direct visualization method using magnetic resonance imaging (MRI) was developed to characterize sand grain size distribution, nonaqueous phase liquid (NAPL) source zone architecture, and aqueous flowpaths in a three-dimensional (3-D) flowcell (26.5 cm x 10.5 cm x 10.5 cm) packed with a heterogeneous distribution of five different sand fractions. All images were acquired at a resolution of 0.1875 cm x 0.1875 cm x 0.225 cm. A 1H image of pore water resolved the heterogeneous permeability field; grain size differences as small as 0.1 mm could be distinguished. A time series of 1H images of water doped with the paramagnetic tracer MnCl2 were acquired and used to obtain voxel-scale breakthrough curves. Water preferentially flowed through coarse sands before NAPL release. After NAPL release, the flow bypassed NAPLzones, and bypassing was more evident for high NAPL saturation zones. A time series of 19F images of NAPL were acquired and used to determine voxel-scale NAPL saturation (Sn) during dissolution. Results show that 93% of NAPL mass was in the coarsest sand, most NAPL was trapped as pools and not as residual ganglia, NAPL saturation increased with depth, and the NAPL dissolution front moved vertically from the top to the bottom of the flowcell during the first 170 pore volumes of waterflushed. NAPL component effluent concentrations initially increased due to the development of flow in zones with decreasing NAPL saturation. Flowpath images suggest that this occurs as NAPL transitions from pools (Sn > 0.15) to residual ganglia. The results highlight the importance of flow bypassing and provide the opportunity to develop more accurate NAPL dissolution models.

Characterization of the velocity field organization in heterogeneous media by conditional correlation

The purpose of the present work is to quantify the correlation structure of simulated velocity fields in heterogeneous permeability fields and to discuss how to represent it in upscaled transport models. We investigate the velocity field correlation structure for multinormal log permeability fields. The simulated velocity distributions are analyzed in a Lagrangian framework, i.e., along the particles' paths. To quantify the different spatial organization of low‐ and high‐velocity zones, we condition the estimated velocity correlation length and time on the initial particle velocity. The velocity correlation length is found to increase with the initial particle velocity, following a power law. Such an effect is likely due to the channeling of high‐velocity zones, which implies that particles keep memory of their initial velocity over longer distances for high initial velocities than for low initial velocities. Two distinct regimes are identified for the velocity correlation time. For low initial particle velocity the correlation time is controlled by the large time needed to escape from the low‐velocity zones. For high initial particle velocity it is controlled by the large time needed for particles to sample the whole velocity field, in particular low‐velocity zones. One of the consequences of these results is that for such velocity fields the nonlinear dependence of both the correlation length and time on the particle initial velocity restricts the use of spatial or temporal Markovian assumptions for modeling velocity transitions in effective transport models.

NUMERICAL SIMULATION OF THE IMPACT OF WATER-AIR FRONTS ON RADIONUCLIDE PLUMES IN HETEROGENEOUS MEDIA

The goal of this paper is to investigate the interaction of water-air fronts with radionuclide plumes in unsaturated heterogeneous porous media. This problem is modeled by a system of equations that describes both water-air flow and radionuclide transport. The water-air flow problem is solved numerically by a mixed finite element combined with a non-oscillatory central difference scheme. For the radionuclide transport equation we use the Modified Method of Characteristics (MMOC).We present results of numerical simulations for heterogeneous permeability fields taking into account sorption effects.

Unstructured grid optimization for improved monotonicity of discrete solutions of elliptic equations with highly anisotropic coefficients

Multipoint flux approximation (MPFA) techniques are commonly applied for discretizing the porous media flow equations within the context of finite volume numerical procedures. Although these methods can be applied to heterogeneous, anisotropic systems on generally unstructured grids, the inverse of the resulting linear operator can suffer from a loss of monotonicity at high permeability anisotropy ratios, resulting in spurious oscillations of the pressure solution. The purpose of this paper is to develop a method for optimizing unstructured grids in two and three dimensions such that the monotonicity behavior of the MPFA technique is significantly improved. The method employs anisotropic triangulation and can be readily combined with permeability upscaling procedures. Results are presented for a variety of examples and the technique is shown to perform well on problems involving complex grid point distributions and heterogeneous permeability fields with strong anisotropy ratios (of O(100)). Oscillation free pressure solutions are achieved in all cases considered, even for examples in which the original grid shows large oscillations.

Homogenization and effective parameters for the Henry problem in heterogeneous formations

General upscaling of density‐dependent flow is investigated for the classical two‐dimensional Henry problem of saltwater intrusion in coastal aquifers. Combined theoretical and numerical results are here presented. Effective flow and transport parameters for saline intrusions could be derived for statistically isotropic and anisotropic heterogeneous permeability fields by use of homogenization theory, applying also to the preasymptotic regime. Our numerical results indicate that heterogeneities in permeability affect foremost the transient evolution of saltwater intrusion, whereas the steady state saltwater distribution is less sensitive to spatially varying permeabilities and longitudinal dispersion. For the isotropic media the effective permeability is found to correspond to the geometric mean, as under conditions without fluid density contrast. The extension of results to anisotropic media requires nontrivial corrections for the effective permeability, which include the variance of log(k) and the directional correlation lengths. The appropriate dispersion coefficients for the problem correspond to the local dispersion coefficients, rather than macrodispersion coefficients. These results are discussed in light of the specific flow configuration posed in the Henry problem.

Stochastic simulations of steady state unsaturated flow in a three-layer, heterogeneous, dual continuum model of fractured rock

Unsaturated flow through fractured rocks is a concern in the siting and performance of waste disposal facilities such as the proposed geological repository at Yucca Mountain, Nevada. We simulated a small two-dimensional cross-section of welded volcanic tuff, representative of Yucca Mountain stratigraphy, using spatially correlated, randomly heterogeneous fracture permeability fields and homogeneous matrix permeability continua representing various degrees of welding. Ten realizations each of fracture permeability fields for four different variances ( σ 2=0.5, 1.0, 1.5, and 2.0) were generated by the direct Fourier transform method (Robin, M.J.L., Gutjahr, A.L., Sudicky, E.A., Wilson, J.L., 1993. Cross-correlated random field generation with the direct Fourier transform method, Water Resour. Res. 29(7)2385–2398) independently for the welded Tiva Canyon Tuff (TCw), the non-welded Paintbrush Tuff (PTn), and the welded Topopah Spring Tuff (TSw), the latter being the proposed primary repository horizon. Numerical simulations were run for steady state flow at three different uniform water flux boundary conditions. Boundary conditions along the sides were impermeable and the base was open to gas and liquid flow. Numerical simulations were performed using the dual-continuum, two-phase flow simulator METRA, which represents matrix and fractures as dual overlapping continua, where liquid flux between continua can be restricted by a uniform factor. Fracture–matrix interaction was modeled as being less restricted in the PTn as compared to the TCw/TSw. Heterogeneous fracture permeability fields generated strong preferential flow in the TCw/TSw fracture continuum and significant preferential flow in the uniformly permeable, PTn matrix continuum. Flow focusing led to a local increase in saturation, which in turn increased relative permeability to water along the preferential pathways, causing water to flow faster. The development of the preferential pathways reduced the wetted surface area for fracture–matrix interaction leaving a large volume of the rock isolated from the preferential flow pathways. Statistical analysis of water flux values in the three units showed that the magnitude of the ensemble variance, indicating preferential flow pathways, increased with both the variance of fracture permeability and the water flux boundary condition, reaching a plateau in the TCw/TSw fracture continuum units after flowing approximately 10 correlation scales of the vertical fracture permeability. Ensemble covariance of water flux normal to the layering revealed long range correlations in all units, which is longer than the correlation length of the fracture permeability fields. This suggests that the water flux boundary condition influences preferential flow in unsaturated dual continua media. These results suggest that careful analysis of information on fracture permeability variability obtained from pneumatic and hydraulic tests is an important component in understanding deep percolation processes.

Combining finite element and finite volume methods for efficient multiphase flow simulations in highly heterogeneous and structurally complex geologic media

AbstractThe permeability of the Earth's crust commonly varies over many orders of magnitude. Flow velocity can range over several orders of magnitude in structures of interest that vary in scale from centimeters to kilometers. To accurately and efficiently model multiphase flow in geologic media, we introduce a fully conservative node‐centered finite volume method coupled with a Galerkin finite element method on an unstructured triangular grid with a complementary finite volume subgrid. The effectiveness of this approach is demonstrated by comparison with traditional solution methods and by multiphase flow simulations for heterogeneous permeability fields including complex geometries that produce transport parameters and lengths scales varying over four orders of magnitude.

Streamline‐based dual‐porosity simulation of reactive transport and flow in fractured reservoirs

We present a streamline‐based formulation to model two processes: transport of a tracer undergoing rate‐limited sorption and two‐phase (water/oil or air/water) transport in fractured systems, using a dual‐porosity approach. We show that these two processes can be simulated using mathematically equivalent formulations. In both cases the system conceptually has two components: a flowing fraction connected to stagnant regions with transfer between the two domains. Streamlines capture movement through the flowing fraction. Fluid transfer between flowing and stagnant regions enters as a source/sink term in the one‐dimensional transport equations along a streamline. To model flow and transport in fractured systems, we develop a new formulation for the transfer function that matches experimental imbibition data. Then we illustrate the streamline approach with synthetic reservoir problems. We use a finely gridded (over one million grid blocks) three‐dimensional domain with a highly heterogeneous permeability field to study both fracture flow and tracer transport. We find breakthrough curves that are consistent with anomalous transport described by an exponent that characterizes the longtime tail of the transit time distribution. For fracture flow we demonstrate that the speed of fluid advance in the fractures is controlled by the imbibition rate. The run times for the simulations scale approximately linearly with system size, making the method appropriate for the simulation of large numerical models.

The averaging of gravity currents in porous media

We explore the problem of a moving free surface in a water-saturated porous medium that has either a homogeneous or a periodically heterogeneous permeability field. We identify scaling relations and derive similarity solutions for the homogeneous, constant coefficient case in both a Cartesian and an axisymmetric, radial coordinate system. We utilize these similarity scalings to identify half-height slumping time scales as a rough guide for field groundwater cleanup strategies involving injected brines. We derive averaged solutions using homogenization for a vertically periodic, a horizontally periodic, and a two-dimensional periodic case—the solution of which requires solving a cell problem. Using effective coefficients, we connect the first two of these homogenized solutions to the similarity scaling solution derived for the homogeneous case. By simplifying to a thin limit, retaining variations of the porous media in the horizontal direction, we derive a homogenization solution in agreement with the general horizontally layered solution and an expression for the leading-order correction. Finally, we implement two numerical solution approaches and show that self-similar scaling and agreement with leading-order averaging emerge in finite time, and demonstrate the accuracy and convergence rate of the leading order correction for both the interior and the boundary of the domain.

A Generalized Convection-Diffusion Model for Subgrid Transport in Porous Media

Models for transport in heterogeneous subsurface formations usually require some type of treatment for subgrid effects. In this work, a generalized nonlinear convection-diffusion model for subgrid transport in two-dimensional systems is developed and applied. The model, although somewhat heuristic, is motivated by previous findings within both the stochastic and the deterministic frameworks. The numerical calculation of the diffusive and convective subgrid terms is described. The model is applied to several example cases involving heterogeneous permeability fields and different global boundary conditions. Both linear and nonlinear fine scale flux functions are considered. Coarse scale results for oil cut (fraction of oil in the produced fluid) and the global saturation fieldgenerated using the new subgrid model are shown to be in consistently better agreement with reference fine scale solutions than are coarse scale results using standard subgrid treatments. Extensions of the method required to treat more realistic subsurface systems are discussed.

Determination of anisotropy of spatial correlation structure in a three‐dimensional permeability field accompanied by shallow faults

Structural planes, such as fractures, bedding planes, schistosity, and so on, induce heterogeneity of permeability field in geologic bodies. One of the ways to demonstrate the heterogeneity is to geostatistically infer the spatial correlation structure. This paper presents a newly developed geostatistical technique for estimating the three‐dimensional heterogeneous permeability field. In the presented method, the plane‐directional anisotropy of spatial correlation structure of the three‐dimensional permeability field is determined by estimating semivariograms in multiplane directions, using limited in situ measurement data of basement rock permeability. A case study of a dam site shows that the spatial correlation of permeability in the direction parallel to the fault planes is greater than that perpendicular to it. Direction of the largest anisotropy of the spatial correlation of the permeability field is concordant with the strike of the major faults.