A new computational strategy for solving two‐phase flow in strongly heterogeneous poroelastic media of evolving scales

We develop a new computational methodology for solving two‐phase flow in highly heterogeneous porous media incorporating geomechanical coupling subject to uncertainty in the poromechanical parameters. Within the framework of a staggered‐in‐time coupling algorithm, the numerical method proposed herein relies on a Petrov–Galerkin postprocessing approach projected on the Raviart–Thomas space to compute the Darcy velocity of the mixture in conjunction with a locally conservative higher order finite volume discretization of the nonlinear transport equation for the saturation and an operator splitting procedure based on the difference in the time‐scales of transport and geomechanics to compute the effects of transient porosity upon saturation. Notable features of the numerical modeling proposed herein are the local conservation properties inherited by the discrete fluxes that are crucial to correctly capture the fingering patterns arising from the interaction between heterogeneity and nonlinear viscous coupling. Water flooding in a poroelastic formation subject to an overburden is simulated with the geology characterized by multiscale self‐similar permeability and Young modulus random fields with power‐law covariance structure. Statistical moments of the poromechanical unknowns are computed within the framework of a high‐resolution Monte Carlo method. Numerical results illustrate the necessity of adopting locally conservative schemes to obtain reliable predictions of secondary recovery and finger growth in strongly heterogeneous deformable reservoirs. Copyright © 2011 John Wiley & Sons, Ltd.     

Numerical modeling of the flow and transport of radionuclides in heterogeneous porous media

This paper is concerned with numerical methods for the modeling of flow and transport of contaminant in porous media. The numerical methods feature the mixed finite element method over triangles as a solver to the Darcy flow equation and a conservative finite volume scheme for the concentration equation. The convective term is approximated with a Godunov scheme over the dual finite volume mesh, whereas the diffusion–dispersion term is discretized by piecewise linear conforming triangular finite elements. It is shown that the scheme satisfies a discrete maximum principle. Numerical examples demonstrate the effectiveness of the methodology for a coupled system that includes an elliptic equation and a diffusion–convection–reaction equation arising when modeling flow and transport in heterogeneous porous media. The proposed scheme is robust, conservative, efficient, and stable, as confirmed by numerical simulations.      

A finite element formulation of gradient-based plasticity for porous media with C1 interpolation of internal variables

In this paper a new finite element formulation for numerical analysis of diffused and localized failure behavior of saturated and partially saturated gradient poroplastic materials is proposed. The new finite element includes interpolation functions of first order (C₁) for the internal variables field while classical C₀ interpolation functions for the kinematic fields and pore pressure. This finite element formulation is compatible with a thermodynamically consistent gradient poroplastic theory previously proposed by the authors. In this material theory the internal variables are the only ones of non-local character. To verify the numerical efficiency of the proposed finite element formulation, the non-local gradient poroplastic constitutive theory is combined with the modified Cam Clay model for partially saturated continua. Thereby, the volumetric strain of the solid skeleton and the plastic porosity are the internal variables of the constitutive theory. The numerical results in this paper demonstrate the capabilities of the proposed finite element formulation to capture diffuse and localized failure modes of boundary value problems of porous media, depending on the acting confining pressure and on the material saturation degree.     

Sequential implicit vertex approximate gradient discretization of incompressible two-phase Darcy flows with discontinuous capillary pressure

This work deals with sequential implicit schemes for incompressible and immiscible two-phase Darcy flows which are commonly used and well understood in the case of spatially homogeneous capillary pressure functions. To our knowledge, the stability of this type of splitting schemes solving sequentially a pressure equation followed by the saturation equation has not been investigated so far in the case of discontinuous capillary pressure curves at different rock type interfaces. It will be shown here to raise severe stability issues for which stabilization strategies are investigated in this work. To fix ideas, the spatial discretization is based on the Vertex Approximate Gradient (VAG) scheme accounting for unstructured polyhedral meshes combined with an Hybrid Upwinding (HU) of the transport term and an upwind positive approximation of the capillary and gravity fluxes. The sequential implicit schemes are built from the total velocity formulation of the two-phase flow model and only differ in the way the conservative VAG total velocity fluxes are approximated. The stability, accuracy and computational cost of the sequential implicit schemes studied in this work are tested on oil migration test cases in 1D, 2D and 3D basins with a large range of capillary pressure parameters for the drain and barrier rock types. It will be shown that usual splitting strategies fail to capture the right solutions for highly contrasted rock types and that it can be fixed by maintaining locally the pressure saturation coupling at different rock type interfaces in the definition of the conservative total velocity fluxes. The numerical investigation of the sequential schemes is also extended to the widely used finite volume Two-Point Flux Approximation spatial discretization.

     

The use of laboratory experiments for the study of conservative solute transport in heterogeneous porous media

 Laboratory experiments on heterogeneous porous media (otherwise known as intermediate scale experiments, or ISEs) have been increasingly relied upon by hydrogeologists for the study of saturated and unsaturated groundwater systems. Among the many ongoing applications of ISEs is the study of fluid flow and the transport of conservative solutes in correlated permeability fields. Recent advances in ISE design have provided the capability of creating correlated permeability fields in the laboratory. This capability is important in the application of ISEs for the assessment of recent stochastic theories. In addition, pressure-transducer technology and visualization methods have provided the potential for ISEs to be used in characterizing the spatial distributions of both hydraulic head and local water velocity within correlated permeability fields. Finally, various methods are available for characterizing temporal variations in the spatial distribution (and, thereby, the spatial moments) of solute concentrations within ISEs. It is concluded, therefore, that recent developments in experimental techniques have provided an opportunity to use ISEs as important tools in the continuing study of fluid flow and the transport of conservative solutes in heterogeneous, saturated porous media. Received, December 1996 · Revised, July 1997 · Accepted, August 1997     

A fully implicit and consistent finite element framework for modeling reservoir compaction with large deformation and nonlinear flow model. Part I: theory and formulation

Reservoir depletion results in rock failure, wellbore instability, hydrocarbon production loss, oil sand production, and ground surface subsidence. Specifically, the compaction of carbonate reservoirs with soft rocks often induces large plastic deformation due to rock pore collapse. On the other hand, following the compaction of reservoirs and failure of rock formations, the porosity and permeability of formations will, in general, decrease. These bring a challenge for reservoir simulations because of high nonlinearity of coupled geomechanics and fluid flow fields. In this work, we present a fully implicit, fully coupled, and fully consistent finite element formulation for coupled geomechanics and fluid flow problems with finite deformation and nonlinear flow models. The Pelessone smooth cap plasticity model, an important material model to capture rock compaction behavior and a challenging material model for implicit numerical formulations, is incorporated in the proposed formulation. Furthermore, a stress-dependent permeability model is taken into account in the formulation. A co-rotational framework is adopted for finite deformation, and an implicit material integrator for cap plasticity models is consistently derived. Furthermore, the coupled field equations are consistently linearized including nonlinear flow models. The physical theories, nonlinear material and flow models, and numerical formulations are the focus of part I of this work. In part II, we verify the proposed numerical framework and demonstrate the performance of our numerical formulation using several numerical examples including a field reservoir with soft rocks undergoing serious compaction.     

Finite element methods for variable density flow and solute transport

Saltwater intrusion into coastal freshwater aquifers is an ongoing problem that will continue to impact coastal freshwater resources as coastal populations increase. To effectively model saltwater intrusion, the impacts of increased salt content on fluid density must be accounted for to properly model saltwater/freshwater transition zones and sharp interfaces. We present a model for variable density fluid flow and solute transport where a conforming finite element method discretization with a locally conservative velocity post-processing method is used for the flow model and the transport equation is discretized using a variational multiscale stabilized conforming finite element method. This formulation provides a consistent velocity and performs well even in advection-dominated problems that can occur in saltwater intrusion modeling. The physical model is presented as well as the formulation of the numerical model and solution methods. The model is tested against several 2-D and 3-D numerical and experimental benchmark problems, and the results are presented to verify the code.     

Finite Element Solvers for Biot’s Poroelasticity Equations in Porous Media

We study and compare five different combinations of finite element spaces for approximating the coupled flow and solid deformation system, so-called Biot’s equations. The permeability and porosity fields are heterogeneous and depend on solid displacement and fluid pressure. We provide detailed comparisons among the continuous Galerkin, discontinuous Galerkin, enriched Galerkin, and two types of mixed finite element methods. Several advantages and disadvantages for each of the above techniques are investigated by comparing local mass conservation properties, the accuracy of the flux approximation, number of degrees of freedom (DOF), and wall and CPU times. Three-field formulation methods with fluid velocity as an additional primary variable generally require a larger number of DOF, longer wall and CPU times, and a greater number of iterations in the linear solver in order to converge. The two-field formulation, a combination of continuous and enriched Galerkin function space, requires the fewest DOF among the methods that conserve local mass. Moreover, our results illustrate that three out of the five methods conserve local mass and produce similar flux approximations when conductivity alteration is included. These comparisons of the key performance indicators of different combinations of finite element methods can be utilized to choose the preferred method based on the required accuracy and the available computational resources.

     

Theoretical and numerical analyses of chemical‐dissolution front instability in fluid‐saturated porous rocks

The chemical‐dissolution front propagation problem exists ubiquitously in many scientific and engineering fields. To solve this problem, it is necessary to deal with a coupled system between porosity, pore‐fluid pressure and reactive chemical‐species transport in fluid‐saturated porous media. Because there was confusion between the average linear velocity and the Darcy velocity in the previous study, the governing equations and related solutions of the problem are re‐derived to correct this confusion in this paper. Owing to the morphological instability of a chemical‐dissolution front, a numerical procedure, which is a combination of the finite element and finite difference methods, is also proposed to solve this problem. In order to verify the proposed numerical procedure, a set of analytical solutions has been derived for a benchmark problem under a special condition where the ratio of the equilibrium concentration to the solid molar density of the concerned chemical species is very small. Not only can the derived analytical solutions be used to verify any numerical method before it is used to solve this kind of chemical‐dissolution front propagation problem but they can also be used to understand the fundamental mechanisms behind the morphological instability of a chemical‐dissolution front during its propagation within fluid‐saturated porous media. The related numerical examples have demonstrated the usefulness and applicability of the proposed numerical procedure for dealing with the chemical‐dissolution front instability problem within a fluid‐saturated porous medium. Copyright © 2007 John Wiley & Sons, Ltd.     

Conservative, shock-capturing transport methods with nonconservative velocity approximations

Conservative high-resolution, or shock-capturing, methods have become widely used for modeling transport equations described by conservation laws. In many geoscience applications, the transport equation is coupled to a conservation or continuity equation for a velocity field. Depending on how the velocity is approximated, the continuity equation may or may not be satisfied, either locally or globally. In this paper, we discuss the effect this has on a typical high resolution scheme, and propose a correction which accounts for the fact that the velocity may be nonconservative. We present several numerical examples and prove stability bounds and an a priori error estimate for the corrected method.     

A parallelizable method for two‐phase flows in naturally‐fractured reservoirs

A parallelizable, semi‐implicit numerical method is proposed for the study of naturally‐fractured reservoir systems. It has proved to be computationally efficient in producing accurate numerical solutions for the dual‐porosity model for immiscible, two‐phase flow in such reservoirs. The method combines hybridized mixed finite elements, a new version of the modified method of characteristics, a sophisticated operator‐splitting procedure for separating the pressure calculation in the fractures from that of the saturation, another operator splitting to handle the interaction of the matrix blocks and the fractures, and domain decomposition iterative procedures for both the pressure and the saturation. It permits moderately long time steps for the pressure and the saturation in the fractures and matrix blocks by using short, inexpensive microsteps to treat the transport portion of the saturation equation in the fractures. This paper is devoted to the formulation of the method and a discussion of numerical results for five‐spot and vertical cross‐section examples.     

Nonconforming Finite Volume Methods

We extend and generalize recently proposed finite volume methods using the framework of mixed finite element methods. Proposed discretizations are defined for tensor permeabilities and naturally produce a generalization of harmonic averaging, and are therefore well suited for heterogeneous and anisotropic media. They are locally mass conservative and work on extremely flexible distorted meshes. Flux variables can be excluded locally and the resulting discretizations for the pressures has the same stencils for Voronoi/PEBI grids as 2-point finite volume discretizations currently used in many simulators.     

Modelling Two-Phase Incompressible Flow in Porous Media Using Mixed Hybrid and Discontinuous Finite Elements

In this paper, we present a numerical model for simulating two-phase (oil–water and air–water) incompressible and immiscible flow in porous media. The mathematical model which is based on a fractional flow formulation is formed of two nonlinear partial differential equations: a mean pressure equation and a water saturation equation. These two equations can be solved in a sequential manner. Two numerical methods are used to discretize the equations of the two-phase flow model: mixed hybrid finite elements are used to treat the pressure equation, h-based Richards' equation and the diffusion term in the saturation equation, the advection term in the saturation equation is treated with the discontinuous finite elements. We propose a better way to calculate the nonlinear coefficients contained in our equations on each element of the discretized domain. In heterogeneous porous media, the saturation becomes discontinuous at the interface between two porous media. We show in this paper how to use the capillary pressure–saturation relationship in order to handle the saturation jump in the mixed hybrid finite element method. The two-phase flow simulator is verified against analytical solutions for some flow problems treated by other authors.     

Finite element modelling of transport of organic chemicals through soils

Aquifer contamination by organic chemicals in subsurface flow through soils due to leaking underground storage tanks filled with organic fluids is an important groundwater pollution problem. The problem involves transport of a chemical pollutant through soils via flow of three immiscible fluid phases: namely air, water and an organic fluid. In this paper, assuming the air phase is under constant atmospheric pressure, the flow field is described by two coupled equations for the water and the organic fluid flow taking interphase mass transfer into account. The transport equations for the contaminant in all the three phases are derived and assuming partition equilibrium coefficients, a single convective – dispersive mass transport equation is obtained. A finite element formulation corresponding to the coupled differential equations governing flow and mass transport in the three fluid phase porous medium system with constant air phase pressure is presented. Relevant constitutive relationships for fluid conductivities and saturations as function of fluid pressures lead to non-linear material coefficients in the formulation. A general time-integration scheme and iteration by a modified Picard method to handle the non-linear properties are used to solve the resulting finite element equations. Laboratory tests were conducted on a soil column initially saturated with water and displaced by p-cymene (a benzene-derivative hydrocarbon) under constant pressure. The same experimental procedure is simulated by the finite element programme to observe the numerical model behaviour and compare the results with those obtained in the tests. The numerical data agreed well with the observed outflow data, and thus validating the formulation. A hypothetical field case involving leakage of organic fluid in a buried underground storage tank and the subsequent transport of an organic compound (benzene) is analysed and the nature of the plume spread is discussed.     

Relationships among some locally conservative discretization methods which handle discontinuous coefficients

This paper presents the relationships between some numerical methods suitable for a heterogeneous elliptic equation with application to reservoir simulation. The methods discussed are the classical mixed finite element method (MFEM), the control-volume mixed finite element method (CVMFEM), the support operators method (SOM), the enhanced cell-centered finite difference method (ECCFDM), and the multi-point flux-approximation (MPFA) control-volume method. These methods are all locally mass conservative, and handle general irregular grids with anisotropic and heterogeneous discontinuous permeability. In addition to this, the methods have a common weak continuity in the pressure across the edges, which in some cases corresponds to Lagrange multipliers. It seems that this last property is an essential common quality for these methods. T.F. Russell: Partially supported by the National Science Foundation Grant Nos. DMS-0084438 and DMS-0222300.     

Analytical solutions for advective–dispersive solute transport in double‐layered finite porous media

Analytical solutions for advection and dispersion of a conservative solute in a one‐dimensional double‐layered finite porous media are presented. The solutions are applicable to five scenarios that have various combinations of fixed concentration, fixed flux and zero concentration gradient conditions at the inlet and outlet boundaries that provide a wide number of options. Arbitrary initial solute concentration distributions throughout the media can be considered via explicit formulations or numerical integration. The analytical solutions presented have been verified against numerical solutions from a finite‐element‐based approach and an existing closed‐form solution for double‐layered media with an excellent correlation being found in both cases. A practical application pertaining to advective transport induced by consolidation of underlying sediment layers on contaminant movement within a capped contaminated sediment system is presented. Comparison of the calculated concentrations and fluxes with alternative approaches clearly illustrates the need to consider advection processes. Consideration of the different features of contaminant transport due to varying pore‐water velocity fields in primary consolidation and secondary consolidation stages is achieved via the use of non‐uniform initial concentration distributions within the proposed analytical solutions. Copyright © 2010 John Wiley & Sons, Ltd.     

Using global node-based velocity in random walk particle tracking in variably saturated porous media: application to contaminant leaching from road constructions

Precise and efficient numerical simulation of transport processes in subsurface systems is a prerequisite for many site investigation or remediation studies. Random walk particle tracking (RWPT) methods have been introduced in the past to overcome numerical difficulties when simulating propagation processes in porous media such as advection-dominated mass transport. Crucial for the precision of RWPT methods is the accuracy of the numerically calculated ground water velocity field. In this paper, a global node-based method for velocity calculation is used, which was originally proposed by Yeh (Water Resour Res 7:1216–1225, 1981). This method is improved in three ways: (1) extension to unstructured grids, (2) significant enhancement of computational efficiency, and (3) extension to saturated (groundwater) as well as unsaturated systems (soil water). The novel RWPT method is tested with numerical benchmark examples from the literature and used in two field scale applications of contaminant transport in saturated and unsaturated ground water. To evaluate advective transport of the model, the accuracy of the velocity field is demonstrated by comparing several published results of particle pathlines or streamlines. Given the chosen test problem, the global node-based velocity estimation is found to be as accurate as the CK method (Cordes and Kinzelbach in Water Resour Res 28(11):2903–2911, 1992) but less accurate than the mixed or mixed-hybrid finite element methods for flow in highly heterogeneous media. To evaluate advective–diffusive transport, a transport problem studied by Hassan and Mohamed (J Hydrol 275(3–4):242–260, 2003) is investigated here and evaluated using different numbers of particles. The results indicate that the number of particles required for the given problem is decreased using the proposed method by about two orders of magnitude without losing accuracy of the concentration contours as compared to the published numbers.     

Comparing global and local implementations of nonlinear complementary problems for the modeling of multi-component two-phase flow with phase change phenomena

Compositional multiphase flow is considered to be one of the fundamental physical processes in the field of water resources research. The strong nonlinearity and discontinuity emerging from phase transition phenomena pose a serious challenge for numerical modeling. Recently, Lauser et al. (Adv Water Resour 34(8):957–966, 2011) have proposed a numerical scheme, namely the nonlinear complementary problem (NCP), to handle this strong nonlinearity. In this work, the NCP is implemented at both local and global levels of a finite element algorithm. In the former case, the NCP is integrated into the local thermodynamic equilibrium calculation, while in the latter one, it is formulated as one of the governing equations. The two different formulations have been investigated through three well-established benchmarks and analyzed for their efficiency and robustness. It is found that both globally and locally implemented NCP formulations are numerically more efficient and robust in comparison with traditional primary variable switching approach. In homogeneous media, the globally implemented NCP formulation leads to an approximately 20% faster simulation compared to the local NCP. This is because a nested Newton iteration for the local phase state identification can be avoided, and thus, the overall computational resources are saved accordingly. However, for problems involving strongly heterogeneous media, the locally integrated NCP formulation suppresses numerical oscillations and delivers more accurate and robust results, especially at the phase boundary.     

Numerical evaluation of multicomponent cation exchange reactive transport in physically and geochemically heterogeneous porous media

Experimental evidence and stochastic studies strongly show that the transport of reactive solutes in porous media is significantly influenced by heterogeneities in hydraulic conductivity, porosity, and sorption parameters. In this paper, we present Monte Carlo numerical simulations of multicomponent reactive transport involving competitive cation exchange reactions in a two-dimensional vertical physically and geochemically heterogeneous medium. Log hydraulic conductivity, log K, and log cation exchange capacity (log CEC) are assumed to be random Gaussian functions with spherical semivariograms. Random realizations of log K and log CEC are used as input data for the numerical simulation of multicomponent reactive transport with CORE^2D, a general purpose reactive transport code. Longitudinal features of the fronts of reactive and conservative species are computed from the temporal and spatial moments of depth-averaged concentrations. Monte Carlo simulations show that: (1) the displacement of reactive fronts increases with increasing variance of log K, while it decreases with the variance of log CEC; (2) second-order spatial moments increase with increasing variances of log K and log CEC; (3) uncertainties in the mean arrival time are largest (smallest) for negatively (positively) correlated log K and Log CEC; (4) cations undergoing competitive cation exchange exhibit different apparent velocities and retardation factors due to both physical and geochemical heterogeneities; and (5) the correlation between log K and log CEC affects significantly apparent cation retardation factors in heterogeneous aquifers.     

A semi-implicit method for incompressible three-phase flow in porous media

In this paper, we present a semi-implicit method for the incompressible three-phase flow equations in two dimensions. In particular, a high-order discontinuous Galerkin spatial discretization is coupled with a backward Euler discretization in time. We consider a pressure-saturation formulation, decouple the pressure and saturation equations, and solve them sequentially while still keeping each equation implicit in its respective unknown. We present several numerical examples on both homogeneous and heterogeneous media, with varying permeability and porosity. Our results demonstrate the robustness of the scheme. In particular, no slope limiters are required and a relatively large time step may be taken.