Accurate local upscaling with variable compact multipoint transmissibility calculations

We propose a new single-phase local upscaling method that uses spatially varying multipoint transmissibility calculations. The method is demonstrated on two-dimensional Cartesian and adaptive Cartesian grids. For each cell face in the coarse upscaled grid, we create a local fine grid region surrounding the face on which we solve two generic local flow problems. The multipoint stencils used to calculate the fluxes across coarse grid cell faces involve the six neighboring pressure values. They are required to honor the two generic flow problems. The remaining degrees of freedom are used to maximize compactness and to ensure that the flux approximation is as close as possible to being two-point. The resulting multipoint flux approximations are spatially varying (a subset of the six neighbors is adaptively chosen) and reduce to two-point expressions in cases without full-tensor anisotropy. Numerical tests show that the method significantly improves upscaling accuracy as compared to commonly used local methods and also compares favorably with a local–global upscaling method.     

A local discontinuous Galerkin scheme for Darcy flow with internal jumps

We present a new version of the local discontinuous Galerkin method which is capable of dealing with jump conditions along a submanifold Γ_LG (i.e., Henry’s Law) in instationary Darcy flow. Our analysis accounts for a spatially and temporally varying, non-linear permeability tensor in all estimates which is also allowed to have a jump at Γ_LG and gives a convergence order result for the primary and the flux unknowns. In addition to this, different approximation spaces for the primary and the flux unknowns are investigated. The results imply that the most efficient choice is to choose the degree of the approximation space for the flux unknowns one less than that of the primary unknown. The only stabilization in the proposed scheme is represented by a penalty term in the primary unknown.     

Short Note: Effective Hydraulic Conductivity and Transmissivity for Heterogeneous Aquifers

Regional scale models of groundwater flow and transport often employ domain discretizations with grid blocks larger than typical scales of field data. For heterogeneous formations, this difference in scales is often handled by using effective (upscaled) parameters. We investigate the problem of upscaling hydraulic conductivity and transmissivity from a small scale of measurement to a larger scale of grid blocks. Transmissivity statistics is expressed in terms of statistics of hydraulic conductivity, and expressions for the effective (upscaled) hydraulic conductivity K _eff and transmissivity T _eff for steady state flow in confined heterogeneous aquifers are derived by means of stochastic averaging and perturbation analysis. These expressions reveal that the commonly used relation T _eff = BK _eff, where B is the confined aquifer thickness, is not generally valid.     

A monotone nonlinear finite volume method for diffusion equations and multiphase flows

We present a new nonlinear monotone finite volume method for diffusion equation and its application to two-phase flow model. We consider full anisotropic discontinuous diffusion or permeability tensors on conformal polyhedral meshes. The approximation of the diffusive flux uses the nonlinear two-point stencil which provides the conventional seven-point stencil for the discrete diffusion operator on cubic meshes. We show that the quality of the discrete flux in a reservoir simulator has great effect on the front behavior and the water breakthrough time. We compare two two-point flux approximations (TPFA), the proposed nonlinear TPFA and the conventional linear TPFA, and multipoint flux approximation (MPFA). The new nonlinear scheme has a number of important advantages over the traditional linear discretizations. Compared to the linear TPFA, the nonlinear TPFA demonstrates low sensitivity to grid distortions and provides appropriate approximation in case of full anisotropic permeability tensor. For nonorthogonal grids or full anisotropic permeability tensors, the conventional linear TPFA provides no approximation, while the nonlinear flux is still first-order accurate. The computational work for the new method is higher than the one for the conventional TPFA, yet it is rather competitive. Compared to MPFA, the new scheme provides sparser algebraic systems and thus is less computational expensive. Moreover, it is monotone which means that the discrete solution preserves the nonnegativity of the differential solution.     

Nonlinear two-point flux approximation for modeling full-tensor effects in subsurface flow simulations

Subsurface flow models can exhibit strong full-tensor anisotropy due to either permeability or grid nonorthogonality effects. Upscaling procedures, for example, generate full-tensor effects on the coarse scale even for cases in which the underlying fine-scale permeability is isotropic. A multipoint flux approximation (MPFA) is often needed to accurately simulate flow for such systems. In this paper, we present and apply a different approach, nonlinear two-point flux approximation (NTPFA), for modeling systems with full-tensor effects. In NTPFA, transmissibility (which provides interblock connections) is determined from reference global flux and pressure fields for a specific flow problem. These fields can be generated using either fully resolved or approximate global simulations. The use of fully resolved simulations leads to an NTPFA method that corresponds to global upscaling procedures, while the use of approximate simulations gives a method corresponding to recently developed local–global techniques. For both approaches, NTPFA algorithms applicable to both single-scale full-tensor permeability systems and two-scale systems are described. A unified framework is introduced, which enables single-scale and two-scale problems to be viewed in a consistent manner. Extensive numerical results demonstrate that the global and local–global NTPFA techniques provide accurate flow predictions over wide parameter ranges for both single-scale and two-scale systems, though the global procedure is more accurate overall. The applicability of NTPFA to the simulation of two-phase flow in upscaled models is also demonstrated.     

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.     

Interior boundary-aligned unstructured grid generation and cell-centered versus vertex-centered CVD-MPFA performance

Grid generation for reservoir simulation must honor classical key constraints and be boundary aligned such that control-volume boundaries are aligned with geological features such as layers, shale barriers, fractures, faults, pinch-outs, and multilateral wells. An unstructured grid generation procedure is proposed that automates control-volume and/or control point boundary alignment and yields a PEBI-mesh both with respect to primal and dual (essentially PEBI) cells. In order to honor geological features in the primal configuration, we introduce the idea of protection circles, and to generate a dual-cell feature based grid, we construct halos around key geological features. The grids generated are employed to study comparative performance of cell-centred versus cell-vertex control-volume distributed multi-point flux approximation (CVD-MPFA) finite-volume formulations using equivalent degrees of freedom. The formulation of CVD-MPFA schemes in cell-centred and cell-vertex modes is analogous and requires switching control volume from primal to dual or vice versa together with appropriate data structures and boundary conditions. The relative benefits of both types of approximation, i.e., cell-centred versus vertex-centred, are made clear in terms of flow resolution and degrees of freedom required.     

Efficient generation of conditional simulations by chebyshev matrix polynomial approximations to the symmetric square root of the covariance matrix

Consider the problem of generating a realization y₁ of a Gaussian random field on a dense grid of points ₁ conditioned on field observations y₂ collected on a sparse grid of points ₂. An approach to this is to generate first an unconditional realization y over the grid =_{1 2}, and then to produce y₁ by conditioning y on the data y₂. As standard methods for generating y, such as the turning bands, spectral or Cholesky approaches can have various limitations, it has been proposed by M. W. Davis to generate realizations from a matrix polynomial approximations to the square root of the covariance matrix. In this paper we describe how to generate a direct approximation to the conditional realization y₁, on ₁ using a variant of Davis' approach based on approximation by Chebyshev polynomials. The resulting algorithm is simple to implement, numerically stable, and bounds on the approximation error are readily available. Furthermore we show that the conditional realization y₁ can be generated directly with a lower order polynomial than the unconditional realization y, and that further reductions can be achieved by exploiting a nugget effect if one is present. A pseudocode version of the algorithm is provided that can be implemented using the fast Fourier transform if the field is stationary and the grid ₁ is rectangular. Finally, numerical illustrations are given of the algorithm's performance in generating various 2-D realizations of conditional processes on large sampling grids.     

Numerical modeling of two-phase binary fluid mixing using mixed finite elements

Diffusion coefficients of dense gases in liquids can be measured by considering two-phase binary nonequilibrium fluid mixing in a closed cell with a fixed volume. This process is based on convection and diffusion in each phase. Numerical simulation of the mixing often requires accurate algorithms. In this paper, we design two efficient numerical methods for simulating the mixing of two-phase binary fluids in one-dimensional, highly permeable media. Mathematical model for isothermal compositional two-phase flow in porous media is established based on Darcy’s law, material balance, local thermodynamic equilibrium for the phases, and diffusion across the phases. The time-lag and operator-splitting techniques are used to decompose each convection–diffusion equation into two steps: diffusion step and convection step. The Mixed finite element (MFE) method is used for diffusion equation because it can achieve a high-order and stable approximation of both the scalar variable and the diffusive fluxes across grid–cell interfaces. We employ the characteristic finite element method with moving mesh to track the liquid–gas interface. Based on the above schemes, we propose two methods: single-domain and two-domain methods. The main difference between two methods is that the two-domain method utilizes the assumption of sharp interface between two fluid phases, while the single-domain method allows fractional saturation level. Two-domain method treats the gas domain and the liquid domain separately. Because liquid–gas interface moves with time, the two-domain method needs work with a moving mesh. On the other hand, the single-domain method allows the use of a fixed mesh. We derive the formulas to compute the diffusive flux for MFE in both methods. The single-domain method is extended to multiple dimensions. Numerical results indicate that both methods can accurately describe the evolution of the pressure and liquid level.     

Modified Markov Chain Monte Carlo Method for Dynamic Data Integration Using Streamline Approach

In this paper, the Markov Chain Monte Carlo (MCMC) approach is used for sampling of the permeability field conditioned on the dynamic data. The novelty of the approach consists of using an approximation of the dynamic data based on streamline computations. The simulations using the streamline approach allows us to obtain analytical approximations in the small neighborhood of the previously computed dynamic data. Using this approximation, we employ a two-stage MCMC approach. In the first stage, the approximation of the dynamic data is used to modify the instrumental proposal distribution. The obtained chain correctly samples from the posterior distribution; the modified Markov chain converges to a steady state corresponding to the posterior distribution. Moreover, this approximation increases the acceptance rate, and reduces the computational time required for MCMC sampling. Numerical results are presented.     

Natural analogues: a potential approach for developing reliable monitoring methods to understand subsurface CO2 migration processes

One possible way of mitigating carbon dioxide (CO₂) emissions from fossil fuel combustion is using carbon dioxide capture and storage (CCS) technology. However, public perception concerning CO₂ storage in the geosphere is generally negative, being particularly motivated by perceived leakage risks. Therefore, a main issue when attempting to gain public acceptance is ensuring provision of appropriate monitoring practices, aimed at providing health, safety and environmental risk assessment, so that potential risks from CO₂ storage are minimized. Naturally occurring CO₂ deposits provide unique natural analogues for evaluating and validating methods used for the detection and monitoring of CO₂ spreading and degassing into the atmosphere. Geological and hydrological structures of the Cheb Basin (NW Bohemia, Czech Republic) represent such a natural analogue for investigating CO₂ leakage and offer a perfect location at which to verify monitoring tools used for direct investigation of processes along preferential migration paths. This shallow basin dating from the Tertiary age is characterized by up to 300?m thick Neogene sediment deposits and several tectonically active faults. The objectives of this paper are to introduce the CO₂ analogues concept to present the Eger Rift as a suitable location for a natural CO₂ analogue site and to demonstrate to what extent such an analogue site should be used (with a case study). The case study presents the results obtained from a joint application of geoelectrical measurements in combination with soil CO₂ concentration and flux determination methods, for the detection and characterization of natural CO₂ releases at gas seeps (as part of a hierarchic monitoring concept). To highlight discharge-controlling structural near surface features was the initial motivation for the application of geoelectrical measurements. Soil-gas concentration and flux measurement techniques are relatively simple to employ and are valuable methods that can be used to monitor seeping CO₂ along preferential pathways. Joint interpretation of both approaches yields a first insight into fluid paths and reveals that the thickness and permeability of site-specific near surface sedimentary deposits have a great influence upon the spatial distribution of the CO₂ degassing pattern at surface level.     

Nested gridding and streamline-based simulation for fast reservoir performance prediction

Detailed reservoir models routinely contain 10⁶–10⁸ grid blocks. These models often cannot be used directly in a reservoir simulation because of the time and memory required for solving the pressure grid on the fine grid. We propose a nested gridding technique that efficiently obtains an approximate solution for the pressure field. The domain is divided into a series of coarse blocks, each containing several fine cells. Effective mobilities are computed for each coarse grid block and the pressure is then found on the coarse scale. The pressure field within each coarse block is computed using flux boundary conditions obtained from the coarse pressure solution. Streamline-based simulation is used to move saturations forward in time. We test the method for a series of example waterflood problems and demonstrate that the method can give accurate estimates of oil production for large 3D models significantly faster than direct simulation using streamlines on the fine grid, making the method overall approximately up to 1,000 times faster than direct conventional simulation.     

A method for automatic history matching of a compositional reservoir simulator with multipoint flux approximation

A method for history matching of an in-house petroleum reservoir compositional simulator with multipoint flux approximation is presented. This method is used for the estimation of unknown reservoir parameters, such as permeability and porosity, based on production data and inverted seismic data. The limited-memory Broyden–Fletcher–Goldfarb–Shanno method is employed for minimization of the objective function, which represents the difference between simulated and observed data. In this work, we present the key features of the algorithm for calculations of the gradients of the objective function based on adjoint variables. The test example shows that the method is applicable to cases with anisotropic permeability fields, multipoint flux approximation, and arbitrary fluid compositions.     

Fractal properties of multi-order folding as a tool for exploration of low-grade banded iron ores in the Krivoy Rog basin (Ukraine)

The low-grade Palaeoproterozoic stratabound banded iron ores of the Krivoy Rog basin (Ukraine) underwent strong tectonometamorphic deformation into superimposed folds of several orders, with amplitudes from centimetres to hundreds of metres. The across-strike sections of bed surfaces defining the low-grade ore bodies resemble self-similar fractal curves; hence, a fractal geometrical model was developed in order to quantify the complexity and sinuosity of bed contours. Two different methods of measurement (polygonal approximation and two-dimensional grid cell counting) were used for 5–8 different scales. Factual similarity dimension D and other model parameters have been estimated by means of linear regression and compared for both measurement methods. From the fractal model a sinuosity coefficient of contours of the folded bed surfaces K _s and a coefficient of degree of exploration of iron ore bodies K _e were constructed. It is pointed out that parameters of the model can be used for determination of the optimal exploration length scales.     

A preliminary evaluation of the necessity of using a cumulus parameterization scheme in high-resolution simulations of Typhoon Haiyan (2013)

The Weather Research and Forecasting model was used to test the sensitivity of Typhoon Haiyan (2013) to the use of a cumulus parameterization scheme, specifically the revised Kain–Fritsch (rKF) scheme, at high horizontal resolutions with grid spacing varying from 9 to 2 km. The rKF scheme simulated the typhoon in best agreement with the observation compared with other schemes, but some fundamental drawbacks relating the rKF scheme, e.g., neglecting the momentum adjustment and being less applicable to high-resolution modeling than multi-scaled schemes, could influence the results and were discussed. Initial results showed that the typhoon track simulations benefited little from the use of the rKF scheme or a fine resolution, partially because of the similar large-scale steering flows induced by the analyzed boundary conditions used in each simulation. The influences of using the rKF scheme on typhoon intensity, size, structure, and precipitation were dependent on the grid spacing, and the most apparent changes occurred near a grid length of 4 km. At 9–4-km grid spacings, using the rKF scheme produced typhoons much stronger with more rainfall and surface latent heat flux than did using no cumulus parameterization scheme. At 3- or 2-km grid spacing, using the rKF scheme caused little changes on typhoon intensity, and the changes in precipitation and surface latent heat flux were relatively small. These results suggested that the grid spacing of 2 km for simulations using no cumulus parameterization scheme or the grid spacing of 4 km for simulations using the rKF scheme facilitated reproducing the observed Typhoon Haiyan.     

The structure of shocks in the atmospheres of pulsating stars

The structure of shocks propagating through partially ionized hydrogen gas with characteristics typical of the atmospheres of RR Lyr, W Vir, and RV Tau type variables is analyzed in terms of a self-consistent solution of the equations of gas dynamics, atomic kinetics, and radiation transfer. The solutions were obtained for shock waves with velocities 20 km/s≤U ₁≤90 km/s and unperturbed hydrogen gas with temperatures 3000 K≤T ₁≤9000 K and density ρ₁=10^?10 g/cm³. The fraction of the energy of the gas-dynamic flux converted into radiation increases with the shock amplitude, and the ratio of the radiation flux emitted by the shock to the gas kinetic energy flux is 0.4???0.8 for the velocities U ₁ considered. This ratio also increases slightly with the ambient gas temperature T ₁ due to an increase in hydrogen ionization in the radiative precursor. The flux emitted by the leading edge of the shock opposite to the gas flow is several percent higher than the flux emitted in the opposite direction by the trailing edge of the shock. Radiation is mostly concentrated in the Balmer continuum, and the region of efficient Lyman radiation transfer includes gas layers located near the viscous jump (δX=±10⁴ cm). The final gas-compression ratio in units of the limiting compression corresponding to an isothermal approximation is virtually independent of the shock amplitude, and increases with the unperturbed gas temperature from r≈0.5 at T ₁=3000 K to r≈0.9 at T ₁=9000 K.     

Compositional flow modeling using a multi-point flux mixed finite element method

We present a general compositional formulation using multi-point flux mixed finite element (MFMFE) method on general hexahedral grids. The mixed finite element framework allows for local mass conservation, accurate flux approximation, and a more general treatment of boundary conditions. The multi-point flux inherent in MFMFE scheme allows the usage of a full permeability tensor. The proposed formulation is an extension of single and two-phase flow formulations presented by Wheeler and Yotov, SIAM J. Numer. Anal. 44(5), 2082–2106 (35) with similar convergence properties. Furthermore, the formulation allows for black oil, single-phase and multi-phase incompressible, slightly and fully compressible flow models utilizing the same design for different fluid systems. An accurate treatment of diffusive/dispersive fluxes owing to additional velocity degrees of freedom is also presented. The applications areas of interest include gas flooding, CO ₂ sequestration, contaminant removal, and groundwater remediation.     

Integrated Two-Dimensional Modeling of Fluid Flow and Compaction in a Sedimentary Basin

In this paper we employ mixed finite elements and numerically study an integrated two-dimensional model of fluid flow and compaction in a sedimentary basin. This model describes a single phase incompressible flow in a two-dimensional section of a sedimentary basin with vertical compaction. At each time step, an iterative algorithm is used to solve this model. The determination of the grid movement is based on the mass conservation and movement of sediments in the basin, while the mixed method is utilized to solve the fluid flow over the moving grid. Numerical experiments are presented to verify this iterative algorithm and show representative solutions for the model under consideration.     

Spatial grid correction for short-term numerical simulation results of carbon dioxide dissolution in saline aquifers

Numerical simulation is an essential component of many studies of geological storage of carbon dioxide, but care must be taken to ensure the accuracy of the results. Unlike several other possible sources of simulation errors, which have previously been considered in detail and have well-understood techniques for mitigating their effects, comparatively little discussion of the spatial grid dependence of the dissolution rate of carbon dioxide into the formation water has appeared in the literature, despite its importance to simulation studies of geological storage of carbon dioxide in saline aquifers. In many instances, sufficient refinement of the computational grid can be a practical solution. However, this approach is not always feasible, especially for large-scale simulations in three dimensions requiring multiple realisations, which commonly feature a coarse grid due to constraints on available computational capabilities. A measure of the error in the amount of dissolved carbon dioxide introduced by the use of a finite grid is therefore of great interest. In this study, the use of finite-sized grid blocks is shown to overestimate the amount of dissolved carbon dioxide in short-term results by a factor of 1?+?V _f/V _p, where V _f is the grid block volume at the saturation front and V _p is the total grid block volume of the plume. This result can be used in a number of ways to correct the calculated short-term dissolution in coarse-scale simulations so that the amount dissolved agrees better with that obtained from fine-scale simulations.     

Performance of a flux splitting scheme when solving the single-phase pressure equation discretized by MPFA

In this article we present a series of tests to study how well suited the TPFA coefficient matrix is as a preconditioner for the MPFA discrete system of equations in an iterative solver, using a flux splitting method. These tests have been conducted for single-phase flow for a wide range of anisotropy, heterogeneity, and grid skewness (mainly parallelogram grids). We use the K-orthogonal part of the MPFA transmissibilities for a parallelogram grid to govern the TPFA transmissibilities. The convergence of the flux splitting method is for each test case measured by the spectral radius of the iteration matrix.