Efficient simulation of geothermal processes in heterogeneous porous media based on the exponential Rosenbrock–Euler and Rosenbrock-type methods

Simulation of geothermal systems is challenging due to coupled physical processes in highly heterogeneous media. Combining the exponential Rosenbrock–Euler method and Rosenbrock-type methods with control-volume (two-point flux approximation) space discretizations leads to efficient numerical techniques for simulating geothermal systems. In terms of efficiency and accuracy, the exponential Rosenbrock–Euler time integrator has advantages over standard time-discretization schemes, which suffer from time-step restrictions or excessive numerical diffusion when advection processes are dominating. Based on linearization of the equation at each time step, we make use of matrix exponentials of the Jacobian from the spatial discretization, which provide the exact solution in time for the linearized equations. This is at the expense of computing the matrix exponentials of the stiff Jacobian matrix, together with propagating a linearized system. However, using a Krylov subspace or Léja points techniques make these computations efficient.The Rosenbrock-type methods use the appropriate rational functions of the Jacobian of the ODEs resulting from the spatial discretization. The parameters in these schemes are found in consistency with the required order of convergence in time. As a result, these schemes are A-stable and only a few linear systems are solved at each time step. The efficiency of the methods compared to standard time-discretization techniques are demonstrated in numerical examples.     

Grid convergence of variable-density flow simulations in discretely-fractured porous media

Grid convergence in space and time of variable-density flow in fractured-porous rock is systematically assessed. Convergence of the flow simulation is attained using both uniform and adaptive time-stepping. This contrasts to variable-density flow in unfractured porous media where grid convergence variable-density flow problems is almost never achieved. At high discretization levels, the number of fingers in fractured-porous rock is no longer influenced by spatial-temporal grid discretization, which is not the case in unfractured porous media. However, similar to unfractured porous media, the number of fingers in fractured-porous media varies at low discretization levels. Simulated convective pattern and penetration depth of the dense plume in fractured rock depend more on spatial discretization than on temporal discretization. The appropriate spatial-temporal grid is then used to examine some aspects of mixed convection in fractured-porous rock, characterized by the mixed convection number M. The critical mixed convection number M_c = 46 represents the transition between forced and free convection in fractured porous media, which is much higher than M_c = 1 in unfractured porous media. Thus, for mixed convective flow problems, the value of M_c is not a sufficient indicator to predict the convective mode (free convection-forced convection), and the presence of vertical fractures must be included in the prediction of convective flow modes.     

Application of the discontinuous spectral Galerkin method to groundwater flow

The discontinuous spectral Galerkin method uses a finite-element discretization of the groundwater flow domain with basis functions of arbitrary order in each element. The independent choice of the basis functions in each element permits discontinuities in transmissivity in the flow domain. This formulation is shown to be of high order accuracy and particularly suitable for accurately calculating the flow field in porous media. Simulations are presented in terms of streamlines in a bidimensional aquifer, and compared with the solution calculated with a standard finite-element method and a mixed finite-element method. Numerical simulations show that the discontinuous spectral Galerkin approximation is more efficient than the standard finite-element method (in computing fluxes and streamlines/pathlines) for a given accuracy, and it is more accurate on a given grid. On the other hand the mixed finite-element method ensures the continuity of the fluxes at the cell boundaries and it is particular efficient in representing complicated flow fields with few mesh points. Simulations show that the mixed finite-element method is superior to the discontinuous spectral Galerkin method producing accurate streamlines even if few computational nodes are used. The application of the discontinuous Galerkin method is thus of interest in groundwater problems only when high order and extremely accurate solutions are needed.     

Coupled groundwater flow and transport: 2. Thermohaline and 3D convection systems

This work continues the analysis of variable density flow in groundwater systems. It focuses on both thermohaline (double-diffusive) and three-dimensional (3D) buoyancy-driven convection processes. The finite-element method is utilized to tackle these complex non-linear problems in two and three dimensions. The preferred numerical approaches are discussed regarding appropriate basic formulations, balance-consistent discretization techniques for derivative quantitites, extension of the Boussinesq approximation, proper constraint conditions, time marching schemes, and computational strategies for solving large systems. Applications are presented for the thermohaline Elder and salt dome problem as well as for the 3D extension of the Elder problem with and without thermohaline effects and a 3D Bénard convection process. The simulations are performed by using the package FEFLOW. Conclusions are drawn with respect to numerical efforts and the appropriateness for practical needs.     

Probabilistic collocation and lagrangian sampling for advective tracer transport in randomly heterogeneous porous media

The Karhunen-Loeve (KL) decomposition and the polynomial chaos (PC) expansion are elegant and efficient tools for uncertainty propagation in porous media. Over recent years, KL/PC-based frameworks have successfully been applied in several contributions for the flow problem in the subsurface context. It was also shown, however, that the accurate solution of the transport problem with KL/PC techniques is more challenging. We propose a framework that utilizes KL/PC in combination with sparse Smolyak quadrature for the flow problem only. In a subsequent step, a Lagrangian sampling technique is used for transport. The flow field samples are calculated based on a PC expansion derived from the solutions at relatively few quadrature points. To increase the computational efficiency of the PC-based flow field sampling, a new reduction method is applied. For advection dominated transport scenarios, where a Lagrangian approach is applicable, the proposed PC/Monte Carlo method (PCMCM) is very efficient and avoids accuracy problems that arise when applying KL/PC techniques to both flow and transport. The applicability of PCMCM is demonstrated for transport simulations in multivariate Gaussian log-conductivity fields that are unconditional and conditional on conductivity measurements.     

Mixed finite element methods and higher order temporal approximations for variably saturated groundwater flow

Richards’ equation (RE) is commonly used to model flow in variably saturated porous media. However, its solution continues to be difficult for many conditions of practical interest. Among the various time discretizations applied to RE, the method of lines (MOL) has been used successfully to introduce robust, accurate, and efficient temporal approximations. At the same time, a mixed-hybrid finite element method combined with an adaptive, higher order time discretization has shown benefits over traditional, lower order temporal approximations for modeling single-phase groundwater flow in heterogeneous porous media. Here, we extend earlier work for single-phase flow and consider two mixed finite element methods that have been used previously to solve RE using lower order time discretizations with either fixed time steps or empirically based adaption. We formulate the two spatial discretizations within a MOL context for the pressure head form of RE as well as a fully mass-conservative version. We conduct several numerical experiments for both spatial discretizations with each formulation, and we compare the higher order, adaptive time discretization to a first-order approximation with formal error control and adaptive time step selection. Based on the numerical results, we evaluate the performance of the methods for robustness and efficiency.     

Nonlinear seismic response analysis of earth dams

The objective of this paper is to propose a general and efficient numerical procedure for analysing the dynamic response of geotechnical structures, which are considered as both nonlinear and two phase systems. In Section 2, the appropriate coupled dynamic field equations for the response of a two-phase soil system are briefly reviewed. The finite element spatial discretization of the field equations is described and time integration for the resulting nonlinear semi-discrete finite element equations is discussed. In Section 3, iterative techniques are examined for the solution of the global nonlinear system of finite element equations. A large amount of computational effort is expended in the iterative phase of the solution and so the iterative procedure used must be both reliable and efficient. The performance of three iterative procedure is examined: Newton Raphson, Modified Newton Raphson and Quasi-Newton methods, including BGFS and Broyden updates. Finally, in Section 4, the elasto-plastic earthquake response analysis of a two phase nonhomogeneous earth dam is presented. Extensive documentation exists¹ for the particular problem selected including recorded earthquake motions at the base and crest of the dam. The results of the numerical calculations are compared to the recorded response of the dam.     

A computational model for coupled multiphysics processes of CO2 sequestration in fractured porous media

In this paper, a computational model for the simulation of coupled hydromechanical and electrokinetic flow in fractured porous media is introduced. Particular emphasis is placed on modeling CO₂ flow in a deformed, fractured geological formation and the associated electrokinetic flow. The governing field equations are derived based on the averaging theory and the double porosity model. They are solved numerically with a mixed discretization scheme, formulated on the basis of the standard Galerkin finite element method, the extended finite element method, the level-set method and the Petrov–Galerkin method. The standard Galerkin method is utilized to discretize the equilibrium and the diffusive dominant field equations, and the extended finite element method, together with the level-set method and the Petrov–Galerkin method, are utilized to discretize the advective dominant field equations. The level-set method is employed to trace the CO₂ plume front, and the extended finite element method is employed to model the high gradient in the saturation field front. The proposed mixed discretization scheme leads to a convergent system, giving a stable and effectively mesh-independent model. The accuracy and computational efficiency of the proposed model is evaluated by verification and numerical examples. Effects of the fracture spacing on the CO₂ flow and the streaming potential are discussed.     

Petrophysical characterization of porous media starting from micro-tomographic images

A quasi-static scheme based on pore space spatial statistics is presented to simulate pore-scale two-phase capillary-dominant displacement processes. The algorithm is coupled with computational fluid dynamics in order to evaluate saturation functions. Wettability heterogeneity in partial and fractional/mixed-wet media is implemented using a contact angle map. The simulation process is pixel-wised and performed directly on binary images. Bypassing and snap-off are tackled as non-wetting phase trapping mechanisms. Post-processing results include residual saturations, effective permeability and capillary pressure curves for drainage and imbibition scenarios. The primary advantages of the proposed workflow are eliminating pore space skeletisation/ discretization, superior time efficiency and minimal numerical drawbacks when compared to other direct or network-based simulation techniques.     

A parallel laboratory for simulation and visualization of seismic wavefields

The use of parallel computers makes simulation of elastic waves feasible throughout large structures by means of recent advances in domain decomposition methods. We introduce a competitive parallel algorithm for the propagation of elastic waves in complex heterogeneous media using finite-element discretization. This parallel method, called the multiblock method , performs more efficiently than classical domain decomposition techniques based on substructuration, such as the Schur complement technique. It reduces considerably the amount of communication amongst processors because the interface problem between subdomains is solved by taking advantage of Huygens' principle for wave propagation. We provide some numerical examples and detailed studies on the efficiency and performance of the algorithm, proving that it is competitive and less costly, from the computational viewpoint, than algorithms based on the Schur technique.     

Application of implicit sub-time stepping to simulate flow and transport in fractured porous media

In general, the accuracy of numerical simulations is determined by spatial and temporal discretization levels. In fractured porous media, the time step size is a key factor in controlling the solution accuracy for a given spatial discretization. If the time step size is restricted by the relatively rapid responses in the fracture domain to maintain an acceptable level of accuracy in the entire simulation domain, the matrix tends to be temporally over-discretized. Implicit sub-time stepping applies smaller sub-time steps only to the sub-domain where the accuracy requirements are less tolerant and is most suitable for problems where the response is high in only a small portion of the domain, such as within and near the fractures in fractured porous media. It is demonstrated with illustrative examples that implicit sub-time stepping can significantly improve the simulation efficiency with minimal loss in accuracy when simulating flow and transport in fractured porous media. The methodology is successfully applied to density-dependent flow and transport simulations in a Canadian Shield environment, where the flow and transport is dominated by discrete, highly conductive fracture zones.     

华北北部三次地震前的应变释放加速特征

强地震发生前,区域地震活动往往出现应变释放加速(AMR)现象。本文以华北地区的1976年唐山7.8级地震、1998年张北6.2级地震、2006年文安5.1级地震作为目标地震,分别分析三次主震前4~8年期间、40~120 km半径范围内的AMR现象及其特征。结果显示,唐山7.8级地震与张北6.2级地震之前均不同程度出现AMR异常,但文安5.1级地震之前并未出现。另外,本研究发现80 km半径可能是研究区寻找AMR异常有效、合理的空间扫描尺度。基于这三个震例研究的认识,本文进一步针对晋冀蒙交界地区和首都圈地区进行分析计算,结果识别出首都圈地区至目前存在AMR异常地区,这可作为判定华北北部未来强震危险性的参考。     

Appropriate vertical discretization of Richards' equation for two‐dimensional watershed‐scale modelling

A number of watershed‐scale hydrological models include Richards' equation (RE) solutions, but the literature is sparse on information as to the appropriate application of RE at the watershed scale. In most published applications of RE in distributed watershed‐scale hydrological modelling, coarse vertical resolutions are used to decrease the computational burden. Compared to point‐ or field‐scale studies, application at the watershed scale is complicated by diverse runoff production mechanisms, groundwater effects on runoff production, runon phenomena and heterogeneous watershed characteristics. An essential element of the numerical solution of RE is that the solution converges as the spatial resolution increases. Spatial convergence studies can be used to identify the proper resolution that accurately describes the solution with maximum computational efficiency, when using physically realistic parameter values. In this study, spatial convergence studies are conducted using the two‐dimensional, distributed‐parameter, gridded surface subsurface hydrological analysis (GSSHA) model, which solves RE to simulate vadose zone fluxes. Tests to determine if the required discretization is strongly a function of dominant runoff production mechanism are conducted using data from two very different watersheds, the Hortonian Goodwin Creek Experimental Watershed and the non‐Hortonian Muddy Brook watershed. Total infiltration, stream flow and evapotranspiration for the entire simulation period are used to compute comparison statistics. The influences of upper and lower boundary conditions on the solution accuracy are also explored. Results indicate that to simulate hydrological fluxes accurately at both watersheds small vertical cell sizes, of the order of 1 cm, are required near the soil surface, but not throughout the soil column. The appropriate choice of approximations for calculating the near soil‐surface unsaturated hydraulic conductivity can yield modest increases in the required cell size. Results for both watersheds are quite similar, even though the soils and runoff production mechanisms differ greatly between the two catchments. Copyright © 2003 John Wiley & Sons, Ltd.     

Adaptive finite element simulation of Stokes flow in porous media

The Stokes problem describes flow of an incompressible constant-viscosity fluid when the Reynolds number is small so that inertial and transient-time effects are negligible. The numerical solution of the Stokes problem requires special care, since classical finite element discretization schemes, such as piecewise linear interpolation for both the velocity and the pressure, fail to perform. Even when an appropriate scheme is adopted, the grid must be selected so that the error is as small as possible. Much of the challenge in solving Stokes problems is how to account for complex geometry and to capture important features such as flow separation. This paper applies adaptive mesh techniques, using a posteriori error estimates, in the finite element solution of the Stokes equations that model flow at pore scales. Different selected numerical test cases associated with various porous geometrics are presented and discussed to demonstrate the accuracy and efficiency of our methodology.     

On the predictivity of pore-scale simulations: Estimating uncertainties with multilevel Monte Carlo

A fast method with tunable accuracy is proposed to estimate errors and uncertainties in pore-scale and Digital Rock Physics (DRP) problems. The overall predictivity of these studies can be, in fact, hindered by many factors including sample heterogeneity, computational and imaging limitations, model inadequacy and not perfectly known physical parameters. The typical objective of pore-scale studies is the estimation of macroscopic effective parameters such as permeability, effective diffusivity and hydrodynamic dispersion. However, these are often non-deterministic quantities (i.e., results obtained for specific pore-scale sample and setup are not totally reproducible by another “equivalent” sample and setup). The stochastic nature can arise due to the multi-scale heterogeneity, the computational and experimental limitations in considering large samples, and the complexity of the physical models. These approximations, in fact, introduce an error that, being dependent on a large number of complex factors, can be modeled as random. We propose a general simulation tool, based on multilevel Monte Carlo, that can reduce drastically the computational cost needed for computing accurate statistics of effective parameters and other quantities of interest, under any of these random errors. This is, to our knowledge, the first attempt to include Uncertainty Quantification (UQ) in pore-scale physics and simulation. The method can also provide estimates of the discretization error and it is tested on three-dimensional transport problems in heterogeneous materials, where the sampling procedure is done by generation algorithms able to reproduce realistic consolidated and unconsolidated random sphere and ellipsoid packings and arrangements. A totally automatic workflow is developed in an open-source code [1], that include rigid body physics and random packing algorithms, unstructured mesh discretization, finite volume solvers, extrapolation and post-processing techniques. The proposed method can be efficiently used in many porous media applications for problems such as stochastic homogenization/upscaling, propagation of uncertainty from microscopic fluid and rock properties to macro-scale parameters, robust estimation of Representative Elementary Volume size for arbitrary physics.     

On the aquifer's integrated balance equations

A general methodology is presented for describing transport phenomena in porous media at a macroscopic level. Then, these macroscopic balance equations are integrated (or averaged) along the vertical for confined, leaky and phreatic aquifers.The results are employed to derive (averaged) aquifer equations for the flow of water and of a solute (hydrodynamic dispersion). It is shown that in all cases, the resulting equation is identical to that derived on the basis of an assumption of horizontal flow (the Dupuit assumption).Macrodispersion, occurring at the aquifer level, is discussed and appropriate coefficients are proposed.     

A mixed-dimensional finite volume method for two-phase flow in fractured porous media

We present a vertex-centered finite volume method for the fully coupled, fully implicit discretization of two-phase flow in fractured porous media. Fractures are discretely modeled as lower dimensional elements. The method works on unstructured, locally refined grids and on parallel computers with distributed memory. An implicit time discretization is employed and the nonlinear systems of equations are solved with a parallel Newton-multigrid method. Results from two-dimensional and three-dimensional simulations are presented.     

The development of stochastic space transformation and diagrammatic perturbation techniques in subsurface hydrology

This paper develops concepts and methods to study stochastic hydrologic models. Problems regarding the application of the existing stochastic approaches in the study of groundwater flow are acknowledged, and an attempt is made to develop efficient means for their solution. These problems include: the spatial multi-dimensionality of the differential equation models governing transport-type phenomena; physically unrealistic assumptions and approximations and the inadequacy of the ordinary perturbation techniques. Multi-dimensionality creates serious mathematical and technical difficulties in the stochastic analysis of groundwater flow, due to the need for large mesh sizes and the poorly conditioned matrices arising from numerical approximations. An alternative to the purely computational approach is to simplify the complex partial differential equations analytically. This can be achieved efficiently by means of a space transformation approach, which transforms the original multi-dimensional problem to a much simpler unidimensional space. The space transformation method is applied to stochastic partial differential equations whose coefficients are random functions of space and/or time. Such equations constitute an integral part of groundwater flow and solute transport. Ordinary perturbation methods for studying stochastic flow equations are in many cases physically inadequate and may lead to questionable approximations of the actual flow. To address these problems, a perturbation analysis based on Feynman-diagram expansions is proposed in this paper. This approach incorporates important information on spatial variability and fulfills essential physical requirements, both important advantages over ordinary hydrologic perturbation techniques. Moreover, the diagram-expansion approach reduces the original stochastic flow problem to a closed set of equations for the mean and the covariance function.     

The development of stochastic space transformation and diagrammatic perturbation techniques in subsurface hydrology

This paper develops concepts and methods to study stochastic hydrologic models. Problems regarding the application of the existing stochastic approaches in the study of groundwater flow are acknowledged, and an attempt is made to develop efficient means for their solution. These problems include: the spatial multi-dimensionality of the differential equation models governing transport-type phenomena; physically unrealistic assumptions and approximations and the inadequacy of the ordinary perturbation techniques. Multi-dimensionality creates serious mathematical and technical difficulties in the stochastic analysis of groundwater flow, due to the need for large mesh sizes and the poorly conditioned matrices arising from numerical approximations. An alternative to the purely computational approach is to simplify the complex partial differential equations analytically. This can be achieved efficiently by means of a space transformation approach, which transforms the original multi-dimensional problem to a much simpler unidimensional space. The space transformation method is applied to stochastic partial differential equations whose coefficients are random functions of space and/or time. Such equations constitute an integral part of groundwater flow and solute transport. Ordinary perturbation methods for studying stochastic flow equations are in many cases physically inadequate and may lead to questionable approximations of the actual flow. To address these problems, a perturbation analysis based on Feynman-diagram expansions is proposed in this paper. This approach incorporates important information on spatial variability and fulfills essential physical requirements, both important advantages over ordinary hydrologic perturbation techniques. Moreover, the diagram-expansion approach reduces the original stochastic flow problem to a closed set of equations for the mean and the covariance function.