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.      

Analysis of subdiffusion in disordered and fractured media using a Grünwald-Letnikov fractional calculus model

The increasing applications of fractional calculus in simulating the anomalous transport behavior in disordered and fractured heterogeneous porous media has grown rapidly over the past decade. In the present study, a temporal fractional flux relationship is employed as a constitutive equation to relate the volumetric flow rate to the gradient of the pore pressure. The novelty of this paper entails interpreting the time fractional derivative operator in the flux relationship by the Grünwald-Letnikov (G-L) definition as opposed to the Caputo interpretation which has been widely considered. Subsequently, a numerical scheme based on the block-centered finite-difference discretization is formulated to handle the resulting non-linear fractional diffusion model. In addition, a linear stability analysis is successfully performed to establish the stability criterion of the developed numerical scheme. An expression for the modified incremental material balance index was derived to assess the effectiveness of the numerical discretization process. Finally, numerical experiments were performed to provide qualitative insights into the nature of pressure evolution in a hydrocarbon reservoir under the influence subdiffusion. In summary, the results establish that subdiffusion regime results in the development of higher pressure drop in the reservoir. This paper will provide a strong foundation for researchers interested in investigating anomalous diffusion phenomena in porous media.     

Micromechanical approach to unsaturated water flow in structured geomaterials by two-scale computations

This article presents a micromechanical approach to the problem of unsaturated water flow in heterogeneous porous media in transient conditions. The numerical formulation is based on the two-scale model obtained previously by periodic homogenization. It allows for a coupled solution of the non-linear flow equations at macroscopic and microscopic scales and takes into account the macroscopic anisotropy of the medium and the local non-equilibrium of the capillary pressure. The model was applied to simulate two-dimensional water infiltration at constant flux into an initially dry medium containing inclusions of square and rectangular shapes. The numerical results showed the influence of the inclusion–matrix conductivity ratio and the local geometry on the macroscopic behavior. The influence of the conductivity ratio manifested itself by the acceleration or retardation of the onset of the macroscopic water flux at the outlet, while the local geometry (anisotropy) significantly affected the macroscopic spatial distribution of the water flux. Such type of approach can be extended to simulate coupled phenomena (for example hydro-mechanical problems) with evolving local geometry.     

非饱和裂隙岩体渗流的统计模型研究方法

常规方法研究非饱和裂隙介质往往采用宏观连续体的概念,然而现场试验证明,非饱和裂隙岩体中的渗流具有相当的非均质性。天然裂隙岩体的非饱和渗流是发生在三维裂隙网络中的多场非等温流,要在不完整的信息基础上刻画和表现介质的不均匀性,通常借助随机模拟的手段。将三维裂隙系统近似为概化的二维非均匀多孔介质的二维平面裂隙。采用模拟退火法将取自各种来源的信息资料通过建立适当的目标函数汇集至模型之中,以此模拟裂隙的空间特征。     

饱和孔隙介质中非均匀震电平面波

采用描述地下流体饱和孔隙介质中弹性波与电磁场耦合现象的Pride方程组,利用Helmholtz分解,求解了非均匀震电耦合平面波解,推导了流体饱和孔隙介质中传播的非均匀震电平面波的能流表达式,讨论了非均匀震电平面波的能量特征.结果表明,非均匀震电平面波的平均能流在传播矢量和衰减矢量构成的平面内传播,伴随电磁场对纵波模式的能流没有贡献.当动电耦合系数为0时,震电平面波的能流表达式可退化为Biot弹性波的能流表达式。     

Analysis and comparison of mixed finite element and multipoint flux approximation methods for homogeneous media

We consider discretization on quadrilateral grids of an elliptic operator occurring, for example, in the pressure equation for porous-media flow. In a realistic setting – with non-orthogonal grid, and anisotropic, heterogeneous permeability – special discretization techniques are required. Mixed finite element (MFE) and multipoint flux approximation (MPFA) are two methods that can handle such situations. Previously, a framework for analytical comparison of MFE and MPFA in special cases has been suggested. A comparison of MFE and MPFA-O (one of two main variants of MPFA) for isotropic, homogeneous permeability on a uniformly distorted grid was also performed. In the current paper, we utilize the suggested framework in a slightly different manner to analyze and compare MFE, MPFA-O and MPFA-U (the second main variant of MPFA). We reconsider the case previously analyzed. We also consider the case of generally anisotropic, homogeneous permeability on an orthogonal grid.     

Comparison of Numerical Formulations for Two-phase Flow in Porous Media

Numerical approximation based on different forms of the governing partial differential equation can lead to significantly different results for two-phase flow in porous media. Selecting the proper primary variables is a critical step in efficiently modeling the highly nonlinear problem of multiphase subsurface flow. A comparison of various forms of numerical approximations for two-phase flow equations is performed in this work. Three forms of equations including the pressure-based, mixed pressure–saturation and modified pressure–saturation are examined. Each of these three highly nonlinear formulations is approximated using finite difference method and is linearized using both Picard and Newton–Raphson linearization approaches. Model simulations for several test cases demonstrate that pressure based form provides better results compared to the pressure–saturation approach in terms of CPU_time and the number of iterations. The modification of pressure–saturation approach improves accuracy of the results. Also it is shown that the Newton–Raphson linearization approach performed better in comparison to the Picard iteration linearization approach with the exception for in the pressure–saturation form.     

Solving the Buckley–Leverett equation with gravity in a heterogeneous porous medium

Immiscible two‐phase flow in porous media can be described by the fractional flow model. If capillary forces are neglected, then the saturation equation is a non‐linear hyperbolic conservation law, known as the Buckley–Leverett equation. This equation can be numerically solved by the method of Godunov, in which the saturation is computed from the solution of Riemann problems at cell interfaces. At a discontinuity of permeability this solution has to be constructed from two flux functions. In order to determine a unique solution an entropy inequality is needed. In this article an entropy inequality is derived from a regularisation procedure, where the physical capillary pressure term is added to the Buckley‐Leverett equation. This entropy inequality determines unique solutions of Riemann problems for all initial conditions. It leads to a simple recipe for the computation of interface fluxes for the method of Godunov. This revised version was published online in July 2006 with corrections to the Cover Date.     

A comparison of multiscale methods for elliptic problems in porous media flow

We review and perform comparison studies for three recent multiscale methods for solving elliptic problems in porous media flow; the multiscale mixed finite-element method, the numerical subgrid upscaling method, and the multiscale finite-volume method. These methods are based on a hierarchical strategy, where the global flow equations are solved on a coarsened mesh only. However, for each method, the discrete formulation of the partial differential equations on the coarse mesh is designed in a particular fashion to account for the impact of heterogeneous subgrid structures of the porous medium. The three multiscale methods produce solutions that are mass conservative on the underlying fine mesh. The methods may therefore be viewed as efficient, approximate fine-scale solvers, i.e., as an inexpensive alternative to solving the elliptic problem on the fine mesh. In addition, the methods may be utilized as an alternative to upscaling, as they generate mass-conservative solutions on the coarse mesh. We therefore choose to also compare the multiscale methods with a state-of-the-art upscaling method – the adaptive local–global upscaling method, which may be viewed as a multiscale method when coupled with a mass-conservative downscaling procedure. We investigate the properties of all four methods through a series of numerical experiments designed to reveal differences with regard to accuracy and robustness. The numerical experiments reveal particular problems with some of the methods, and these will be discussed in detail along with possible solutions. Next, we comment on implementational aspects and perform a simple analysis and comparison of the computational costs associated with each of the methods. Finally, we apply the three multiscale methods to a dynamic two-phase flow case and demonstrate that high efficiency and accurate results can be obtained when the subgrid computations are made part of a preprocessing step and not updated, or updated infrequently, throughout the simulation. The research is funded by the Research Council of Norway under grant nos. 152732 and 158908.     

The representer method for two-phase flow in porous media

The representer method is applied to a one-dimensional two-phase flow model in porous media; capillary pressure and gravity are neglected. The Euler–Lagrange equations must be linearized, and one such linearization is presented here. The representer method is applied to the linear system iteratively until convergence, though a rigorous proof of convergence is out of reach. The linearization chosen is easy to calculate but does not converge for certain weights; however, a simple damping restores convergence at the cost of extra iterations. Numerical experiments are performed that illustrate the method, and quick comparison to the ensemble Kalman smoother is made. This research was supported by NSF grant EIA-0121523.     

A space-time discontinuous Galerkin method applied to single-phase flow in porous media

A space-time discontinuous Galerkin finite element method is proposed and applied to a convection-dominant single-phase flow problem in porous media. The numerical scheme is based on a coupled space-time finite element discretization allowing for discontinuous approximations in space and in time. The continuities on the element interfaces are weakly enforced by the flux treatments, so that no extra penalty factor has to be determined. The resulting space-time formulation possesses the advantage of capturing the steep concentration front with sharp gradients efficiently. The stability and reliability of the proposed approach is demonstrated by numerical experiments. The author is grateful to the DFG (German Science Foundation—Deutsche Forschungsgemeinschaft) for the financial support under the grant number Di 430/4-2.     

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.     

Numerical simulation of NAPL flow in the subsurface

A three-dimensional, three-phase numerical model is presented for simulating the movement of immiscible fluids, including nonaqueous-phase liquids (NAPLs), through porous media. The model is designed to simulate soil flume experiments and for practical application to a wide variety of contamination scenarios involving light or dense NAPLs in heterogeneous subsurface systems. The model is derived for the three-phase flow of water, NAPL, and air in porous media. The basic governing equations are based upon the mass conservation of the constitutents within the phases. The descretization chosen to transform the governing equations into the approximating equations, although logically regular, is very general. The approximating equations are a set of simultaneous coupled nonlinear equations which are solved by the Newton-Raphson method. The linear system solutions needed for the Newton-Raphson method are obtained using a matrix of preconditioner/accelerator iterative methods. Because of the special way the governing equations are implemented, the model is capable of simulating many of the phenomena considered necessary for the sucessful simulation of field problems including entry pressure phenomena, entrapment, and preferential flow paths. The model is verified by comparing it with several exact analytic test solutions and three soil flume experiments involving the introduction and movement of light nonaqueous-phase liquid (LNAPL) or dense nonaqueous-phase liquid (DNAPL) in heterogeneous sand containing a watertable. This revised version was published online in August 2006 with corrections to the Cover Date.     

A parallel block preconditioner for large-scale poroelasticity with highly heterogeneous material parameters

Large-scale simulations of coupled flow in deformable porous media require iterative methods for solving the systems of linear algebraic equations. Construction of efficient iterative methods is particularly challenging in problems with large jumps in material properties, which is often the case in realistic geological applications, such as basin evolution at regional scales. The success of iterative methods for such problems depends strongly on finding effective preconditioners with good parallel scaling properties, which is the topic of the present paper. We present a parallel preconditioner for Biot’s equations of coupled elasticity and fluid flow in porous media. The preconditioner is based on an approximation of the exact inverse of the two-by-two block system arising from a finite element discretisation. The approximation relies on a highly scalable approximation of the global Schur complement of the coefficient matrix, combined with generally available state-of-the-art multilevel preconditioners for the individual blocks. This preconditioner is shown to be robust on problems with highly heterogeneous material parameters. We investigate the weak and strong parallel scaling of this preconditioner on up to 512 processors and demonstrate its ability on a realistic basin-scale problem in poroelasticity with over eight million tetrahedral elements.     

Exponential time integrators for stochastic partial differential equations in 3D reservoir simulation

The transport of chemically reactive solutes (e.g. surfactants, CO₂ or dissolved minerals) is of fundamental importance to a wide range of applications in oil and gas reservoirs such as enhanced oil recovery and mineral scale formation. In this work, we investigate exponential time integrators, in conjunction with an upwind weighted finite volume discretisation in space, for the efficient and accurate simulation of advection–dispersion processes including non-linear chemical reactions in highly heterogeneous 3D oil reservoirs. We model sub-grid fluctuations in transport velocities and uncertainty in the reaction term by writing the advection–dispersion–reaction equation as a stochastic partial differential equation with multiplicative noise. The exponential integrators are based on the variation of constants solution and solve the linear system exactly. While this is at the expense of computing the exponential of the stiff matrix representing the finite volume discretisation, the use of real Léja point or the Krylov subspace technique to approximate the exponential makes these methods competitive compared to standard finite difference-based time integrators. For the deterministic system, we investigate two exponential time integrators, the second-order accurate exponential Euler midpoint (EEM) scheme and exponential time differencing of order one (ETD1). All our numerical examples demonstrate that our methods can compete in terms of efficiency and accuracy compared with standard first-order semi-implicit time integrators when solving (stochastic) partial differential equations that model mixing and chemical reactions in 3D heterogeneous porous media. Our results suggest that exponential time integrators such as the ETD1 and EEM schemes could be applied to typical 3D reservoir models comprising tens to hundreds of thousands unknowns.     

非均质介质中重非水相污染物运移受泄漏速率影响数值分析

重非水相污染物(DNAPL)在地下介质中运移和分布受多种因素控制,包括DNAPL本身的物理化学性质,土的性质,泄漏条件等等。由于介质的非均质性,使得多相流运移行为更为复杂。基于地下水随机理论构建渗透率随机场,采用蒙特卡罗方法探讨泄漏速率对非均质饱和介质中DNAPL运移的影响。数值结果表明,在泄漏总量一定的情况下,泄漏速率越低,介质非均质性对DNAPL运移的影响程度越高。反之,DNAPL的渗漏速率越高,小尺度地层的非均质性影响越低。由于DNAPL运移过程中在垂直方向受重力的影响,污染羽在空间上的质心位置(一阶矩)以及展布范围(二阶矩)在垂直方向上的变异程度要高于水平方向。     

A discretization and multigrid solver for a Darcy–Stokes system of three dimensional vuggy porous media

We develop a finite element discretization and multigrid solver for a Darcy–Stokes system of three-dimensional vuggy porous media, i.e., porous media with cavities. The finite element method uses low-order mixed finite elements in the Darcy and Stokes domains and special transition elements near the Darcy–Stokes interface to allow for tangential discontinuities implied by the Beavers–Joseph boundary condition. We design a multigrid method to solve the resulting saddle point linear system. The intertwining of the Darcy and Stokes subdomains makes the resulting matrix highly ill-conditioned. The velocity field is very irregular, and its discontinuous tangential component at the Darcy–Stokes interface makes it difficult to define intergrid transfer operators. Our definition is based on mass conservation and the analysis of the orders of magnitude of the solution. The coarser grid equations are defined using the Galerkin method. A new smoother of Uzawa type is developed based on taking an optimal step in a good search direction. Our algorithm has a measured convergence factor independent of the size of the system, at least when there are no disconnected vugs. We study the macroscopic effective permeability of a vuggy medium, showing that the influence of vug orientation; shape; and, most importantly, interconnectivity determine the macroscopic flow properties of the medium. This work was supported by the U.S. National Science Foundation under grants DMS-0074310 and DMS-0417431.     

Multidimensional upstream weighting for multiphase transport in porous media

Truly multidimensional methods for hyperbolic equations use flow-based information to determine the computational stencil, as opposed to applying one-dimensional methods dimension by dimension. By doing this, the numerical errors are less correlated with the underlying computational grid. This can be important for reducing bias in flow problems that are inherently unstable at simulation scale, such as in certain porous media problems. In this work, a monotone, multi-D framework for multiphase flow and transport in porous media is developed. A local coupling of the fluxes is introduced through the use of interaction regions, resulting in a compact stencil. A relaxed volume formulation of the coupled hyperbolic–elliptic system is used that allows for nonzero residuals in the pressure equation to be handled robustly. This formulation ensures nonnegative masses and saturations (volume fractions) that sum to one (Acs et al., SPE J 25(4):543–553, 1985). Though the focus of the paper is on immiscible flow, an extension of the methods to a class of more general scalar hyperbolic equations is also presented. Several test problems demonstrate that the truly multi-D schemes reduce biasing due to the computational grid.     

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     

Capillary hysteresis and gravity segregation in two phase flow through porous media

We study the gravity driven flow of two fluid phases in a one dimensional homogeneous porous column when history dependence of the pressure difference between the phases (capillary pressure) is taken into account. In the hyperbolic limit, solutions of such systems satisfy the Buckley-Leverett equation with a non-monotone flux function. However, solutions for the hysteretic case do not converge to the classical solutions in the hyperbolic limit in a wide range of situations. In particular, with Riemann data as initial condition, stationary shocks become possible in addition to classical components such as shocks, rarefaction waves and constant states. We derive an admissibility criterion for the stationary shocks and outline all admissible shocks. Depending on the capillary pressure functions, flux function and the Riemann data, two cases are identified a priori for which the solution consists of a stationary shock. In the first case, the shock remains at the point where the initial condition is discontinuous. In the second case, the solution is frozen in time in at least one semi-infinite half. The predictions are verified using numerical results.