Integration of local–global upscaling and grid adaptivity for simulation of subsurface flow in heterogeneous formations

We propose a methodology, called multilevel local–global (MLLG) upscaling, for generating accurate upscaled models of permeabilities or transmissibilities for flow simulation on adapted grids in heterogeneous subsurface formations. The method generates an initial adapted grid based on the given fine-scale reservoir heterogeneity and potential flow paths. It then applies local–global (LG) upscaling for permeability or transmissibility [7], along with adaptivity, in an iterative manner. In each iteration of MLLG, the grid can be adapted where needed to reduce flow solver and upscaling errors. The adaptivity is controlled with a flow-based indicator. The iterative process is continued until consistency between the global solve on the adapted grid and the local solves is obtained. While each application of LG upscaling is also an iterative process, this inner iteration generally takes only one or two iterations to converge. Furthermore, the number of outer iterations is bounded above, and hence, the computational costs of this approach are low. We design a new flow-based weighting of transmissibility values in LG upscaling that significantly improves the accuracy of LG and MLLG over traditional local transmissibility calculations. For highly heterogeneous (e.g., channelized) systems, the integration of grid adaptivity and LG upscaling is shown to consistently provide more accurate coarse-scale models for global flow, relative to reference fine-scale results, than do existing upscaling techniques applied to uniform grids of similar densities. Another attractive property of the integration of upscaling and adaptivity is that process dependency is strongly reduced, that is, the approach computes accurate global flow results also for flows driven by boundary conditions different from the generic boundary conditions used to compute the upscaled parameters. The method is demonstrated on Cartesian cell-based anisotropic refinement (CCAR) grids, but it can be applied to other adaptation strategies for structured grids and extended to unstructured grids.     

Multiscale mixed/mimetic methods on corner-point grids

Multiscale simulation is a promising approach to facilitate direct simulation of large and complex grid models for highly heterogeneous petroleum reservoirs. Unlike traditional simulation, approaches based on upscaling/downscaling, multiscale methods seek to solve the full flow problem by incorporating subscale heterogeneities into local discrete approximation spaces. We consider a multiscale formulation based on a hierarchical grid approach, where basis functions with subgrid resolution are computed numerically to correctly and accurately account for subscale variations from an underlying (fine-scale) geomodel when solving the global flow equations on a coarse grid. By using multiscale basis functions to discretise the global flow equations on a (moderately sized) coarse grid, one can retain the efficiency of an upscaling method and, at the same time, produce detailed and conservative velocity fields on the underlying fine grid. For pressure equations, the multiscale mixed finite-element method (MsMFEM) has been shown to be a particularly versatile approach. In this paper, we extend the method to corner-point grids, which is the industry standard for modelling complex reservoir geology. To implement MsMFEM, one needs a discretisation method for solving local flow problems on the underlying fine grids. In principle, any stable and conservative method can be used. Here, we use a mimetic discretisation, which is a generalisation of mixed finite elements that gives a discrete inner product, allows for polyhedral elements, and can (easily) be extended to curved grid faces. The coarse grid can, in principle, be any partition of the subgrid, where each coarse block is a connected collection of subgrid cells. However, we argue that, when generating coarse grids, one should follow certain simple guidelines to achieve improved accuracy. We discuss partitioning in both index space and physical space and suggest simple processing techniques. The versatility and accuracy of the new multiscale mixed methodology is demonstrated on two corner-point models: a small Y-shaped sector model and a complex model of a layered sedimentary bed. A variety of coarse grids, both violating and obeying the above mentioned guidelines, are employed. The MsMFEM solutions are compared with a reference solution obtained by direct simulation on the subgrid.     

Numerical Calculation of Equivalent Cell Permeability Tensors for General Quadrilateral Control Volumes

A new method for upscaling fine scale permeability fields to general quadrilateral-shaped coarse cells is presented. The procedure, referred to as the conforming scale up method, applies a triangle-based finite element technique, capable of accurately resolving both the coarse cell geometry and the subgrid heterogeneity, to the solution of the local fine scale problem. An appropriate averaging of this solution provides the equivalent permeability tensor for the coarse scale quadrilateral cell. The general level of accuracy of the technique is demonstrated through application to a number of flow problems. The real strength of the conforming scale up method is demonstrated when the method is applied in conjunction with a flow-based gridding technique. In this case, the approach is shown to provide results that are significantly more accurate than those obtained using standard techniques.     

Multiscale simulations of porous media flows in flow-based coordinate system

In this paper, we propose a multiscale technique for the simulation of porous media flows in a flow-based coordinate system. A flow-based coordinate system allows us to simplify the scale interaction and derive the upscaled equations for purely hyperbolic transport equations. We discuss the applications of the method to two-phase flows in heterogeneous porous media. For two-phase flow simulations, the use of a flow-based coordinate system requires limited global information, such as the solution of single-phase flow. Numerical results show that one can achieve accurate upscaling results using a flow-based coordinate system.     

Mortar Upscaling for Multiphase Flow in Porous Media

In mortar space upscaling methods, a reservoir is decomposed into a series of subdomains (blocks) in which independently constructed numerical grids and possibly different physical models and discretization techniques can be employed in each block. Physically meaningful matching conditions are imposed on block interfaces in a numerically stable and accurate way using mortar finite element spaces. Coarse mortar grids and fine subdomain grids provide two-scale approximations. In the resulting effective solution flow is computed in subdomains on the fine scale while fluxes are matched on the coarse scale. In addition the flexibility to vary adaptively the number of interface degrees of freedom leads to more accurate multiscale approximations. This methodology has been implemented in the Center for Subsurface Modeling's multiphysics multiblock simulator IPARS (Integrated Parallel Accurate reservoir Simulator). Computational experiments demonstrate that this approach is scalable in parallel and it can be applied to non-matching grids across the interface, multinumerics and multiphysics models, and mortar adaptivity. Moreover unlike most upscaling approaches the underlying systems can be treated fully implicitly.     

Near-well upscaling for three-phase flows

Large-scale flow models constructed using standard coarsening procedures may not accurately resolve detailed near-well effects. Such effects are often important to capture, however, as the interaction of the well with the formation can have a dominant impact on process performance. In this work, a near-well upscaling procedure, which provides three-phase well-block properties, is developed and tested. The overall approach represents an extension of a recently developed oil–gas upscaling procedure and entails the use of local well computations (over a region referred to as the local well model (LWM)) along with a gradient-based optimization procedure to minimize the mismatch between fine and coarse-scale well rates, for oil, gas, and water, over the LWM. The gradients required for the minimization are computed efficiently through solution of adjoint equations. The LWM boundary conditions are determined using an iterative local-global procedure. With this approach, pressures and saturations computed during a global coarse-scale simulation are interpolated onto LWM boundaries and then used as boundary conditions for the fine-scale LWM computations. In addition to extending the overall approach to the three-phase case, this work also introduces new treatments that provide improved accuracy in cases with significant flux from the gas cap into the well block. The near-well multiphase upscaling method is applied to heterogeneous reservoir models, with production from vertical and horizontal wells. Simulation results illustrate that the method is able to accurately capture key near-well effects and to provide predictions for component production rates that are in close agreement with reference fine-scale results. The level of accuracy of the procedure is shown to be significantly higher than that of a standard approach which uses only upscaled single-phase flow parameters.     

Upscaling Hydraulic Conductivities in Cross-Bedded Formations

Modern geostatistical techniques allow the generation of high-resolution heterogeneous models of hydraulic conductivity containing millions to billions of cells. Selective upscaling is a numerical approach for the change of scale of fine-scale hydraulic conductivity models into coarser scale models that are suitable for numerical simulations of groundwater flow and mass transport. Selective upscaling uses an elastic gridding technique to selectively determine the geometry of the coarse grid by an iterative procedure. The geometry of the coarse grid is built so that the variances of flow velocities within the coarse blocks are minimum. Selective upscaling is able to handle complex geological formations and flow patterns, and provides full hydraulic conductivity tensor for each block. Selective upscaling is applied to a cross-bedded formation in which the fine-scale hydraulic conductivities are full tensors with principal directions not parallel to the statistical anisotropy of their spatial distribution. Mass transport results from three coarse-scale models constructed by different upscaling techniques are compared to the fine-scale results for different flow conditions. Selective upscaling provides coarse grids in which mass transport simulation is in good agreement with the fine-scale simulations, and consistently superior to simulations on traditional regular (equal-sized) grids or elastic grids built without accounting for flow velocities.     

Successful application of multiscale methods in a real reservoir simulator environment

For the past 10 years or so, a number of so-called multiscale methods have been developed as an alternative approach to upscaling and to accelerate reservoir simulation. The key idea of all these methods is to construct a set of prolongation operators that map between unknowns associated with cells in a fine grid holding the petrophysical properties of the geological reservoir model and unknowns on a coarser grid used for dynamic simulation. The prolongation operators are computed numerically by solving localized flow problems, much in the same way as for flow-based upscaling methods, and can be used to construct a reduced coarse-scale system of flow equations that describe the macro-scale displacement driven by global forces. Unlike effective parameters, the multiscale basis functions have subscale resolution, which ensures that fine-scale heterogeneity is correctly accounted for in a systematic manner. Among all multiscale formulations discussed in the literature, the multiscale restriction-smoothed basis (MsRSB) method has proved to be particularly promising. This method has been implemented in a commercially available simulator and has three main advantages. First, the input grid and its coarse partition can have general polyhedral geometry and unstructured topology. Secondly, MsRSB is accurate and robust when used as an approximate solver and converges relatively fast when used as an iterative fine-scale solver. Finally, the method is formulated on top of a cell-centered, conservative, finite-volume method and is applicable to any flow model for which one can isolate a pressure equation. We discuss numerical challenges posed by contemporary geomodels and report a number of validation cases showing that the MsRSB method is an efficient, robust, and versatile method for simulating complex models of real reservoirs.     

Upscaling facies models to preserve connectivity of designated facies

An upscaling algorithm has been developed that generates an irregular coarse grid that preserves flow connectivity by applying a rule-based upscaling algorithm to a fine-scale facies distribution. The algorithm is demonstrated using stochastically generated paleo-fluvial facies distributions. First, an irregular grid honoring the channel facies is created, followed by computation of effective anisotropic parameters for all coarse-grid cells. For the apparent layer-cake geometry of overbank deposits seen in outcrop, two local upscaling methods are compared: (1) the layered system approximation and (2) the mode. To assess upscaling performance, flow simulations for the original and upscaled grids are compared. The horizontal layered approximation (arithmetic mean) performs poorly, over-predicting lateral connectivity where even infrequent disconnection becomes important. Performance of the mode as an upscaling algorithm depends on the probability that a coarse-grid cell will be dominated by a single facies, and it performs surprisingly well because the upscaled grid-generation algorithm honors the channels, informing the upscaling process. Lastly, the irregular coarse grid was compared to a uniform coarse grid, showing superior performance with the irregular grid. The reduction in grid size achieved by irregular-grid generation will be a function of the geometrical complexity of the geologic objects to be honored.     

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.     

Simulation of anisotropic heterogeneous near-well flow using MPFA methods on flexible grids

Control-volume multipoint flux approximations (MPFA) are discussed for the simulation of complex near-well flow using geometrically flexible grids. Due to the strong non-linearity of the near-well flow, a linear model will, in general, be inefficient. Instead, a model accounting for the logarithmic pressure behavior in the well vicinity is advocated. This involves a non-uniform refinement of the grid in the radial direction. The model accounts for both near-well anisotropies and heterogeneities. For a full simulation involving multiple wells, this single-well approach can easily be coupled with the reservoir model. Numerical simulations demonstrate the convergence behavior of this model using various MPFA schemes under different near-well conditions for single-phase flow regimes. Two-phase simulations support the results of the single-phase simulations.     

Ensemble level upscaling for compositional flow simulation

Uncertainty quantification is typically accomplished by simulating multiple geological realizations, which can be very expensive computationally if the flow process is complicated and the models are highly resolved. Upscaling procedures can be applied to reduce computational demands, though it is essential that the resulting coarse-model predictions correspond to reference fine-scale solutions. In this work, we develop an ensemble level upscaling (EnLU) procedure for compositional systems, which enables the efficient generation of multiple coarse models for use in uncertainty quantification. We apply a newly developed global compositional upscaling method to provide coarse-scale parameters and functions for selected realizations. This global upscaling entails transmissibility and relative permeability upscaling, along with the computation of a-factors to capture component fluxes. Additional features include near-well upscaling for all coarse parameters and functions, and iteration on the a-factors, which is shown to improve accuracy. In the EnLU framework, this global upscaling is applied for only a few selected realizations. For 90 % or more of the realizations, upscaled functions are assigned statistically based on quickly computed flow and permeability attributes. A sequential Gaussian co-simulation procedure is incorporated to provide coarse models that honor the spatial correlation structure of the upscaled properties. The resulting EnLU procedure is applied for multiple realizations of two-dimensional models, for both Gaussian and channelized permeability fields. Results demonstrate that EnLU provides P10, P50, and P90 results for phase and component production rates that are in close agreement with reference fine-scale results. Less accuracy is observed in realization-by-realization comparisons, though the models are still much more accurate than those generated using standard coarsening procedures.     

Implementation of a Locally Conservative Numerical Subgrid Upscaling Scheme for Two-Phase Darcy Flow

We present a locally mass conservative scheme for the approximation of two-phase flow in a porous medium that allows us to obtain detailed fine scale solutions on relatively coarse meshes. The permeability is assumed to be resolvable on a fine numerical grid, but limits on computational power require that computations be performed on a coarse grid. We define a two-scale mixed finite element space and resulting method, and describe in detail the solution algorithm. It involves a coarse scale operator coupled to a subgrid scale operator localized in space to each coarse grid element. An influence function (numerical Greens function) technique allows us to solve these subgrid scale problems independently of the coarse grid approximation. The coarse grid problem is modified to take into account the subgrid scale solution and solved as a large linear system of equations posed over a coarse grid. Finally, the coarse scale solution is corrected on the subgrid scale, providing a fine grid representation of the solution. Numerical examples are presented, which show that near-well behavior and even extremely heterogeneous permeability barriers and streaks are upscaled well by the technique.     

Modeling semi-steady state near-well flow performance for horizontal wells in anisotropic reservoirs

This paper presents a novel methodology to model semi-steady state horizontal well flow performance in an anisotropic reservoir taking into account flow in the near-well region for an arbitrary well trajectory. It is based on an analytical productivity model describing coupled axial reservoir flow and radial well inflow. In order to apply this model in an anisotropic reservoir, the permeability field relative to the radial direction perpendicular to the well trajectory and the axial direction along the well trajectory must first be determined. A classical space transformation is used in concert with rotational transforms to obtain a virtual isotropic model. The transformation preserves the volumes and pressures. It is not a novel concept, but different from previous approaches in the sense that it is only applied in the near-well domain to formulate an equally isotropic media. As a result, the use of this virtual isotropic model requires the Dietz shape factor for an ellipse, transformed from the original cylindrical near-well domain. The Dietz shape factors are determined numerically in this research. The semi-steady state well/near-well model is implemented in a numerical simulator incorporating formation anisotropy and wellbore hydraulics. The specific productivity index along the well trajectory is generated using the virtual configuration. Numerical results for different anisotropy ratios and also incorporating frictional losses in the well are presented. Furthermore, the well/near-well model is applied in coupling with streamline reservoir model for a water flooding case. This appears to be the first coupling of a well hydraulics model and a streamline simulator. It presents the application of the well/near-well model in integrated reservoir simulation in an efficient and accurate manner. The results demonstrate that the coupling approach with a streamline reservoir model and the well/near-well is of great potential for advanced well simulation efficiently.     

A finite volume scheme with improved well modeling in subsurface flow simulation

We present the latest enhancement of the nonlinear monotone finite volume method for the near-well regions. The original nonlinear method is applicable for diffusion, advection-diffusion, and multiphase flow model equations with full anisotropic discontinuous permeability tensors on conformal polyhedral meshes. The approximation of the diffusive flux uses the nonlinear two-point stencil which reduces to the conventional two-point flux approximation (TPFA) on cubic meshes but has much better accuracy for the general case of non-orthogonal grids and anisotropic media. The latest modification of the nonlinear method takes into account the nonlinear (e.g., logarithmic) singularity of the pressure in the near-well region and introduces a correction to improve accuracy of the pressure and the flux calculation. In this paper, we consider a linear version of the nonlinear method waiving its monotonicity for sake of better accuracy. The new method is generalized for anisotropic media, polyhedral grids and nontrivial cases such as slanted, partially perforated wells or wells shifted from the cell center. Numerical experiments show noticeable reduction of numerical errors compared to the original monotone nonlinear FV scheme with the conventional Peaceman well model or with the given analytical well rate.     

Grid adaptation for the Dirichlet–Neumann representation method and the multiscale mixed finite-element method

A Dirichlet–Neumann representation method was recently proposed for upscaling and simulating flow in reservoirs. The DNR method expresses coarse fluxes as linear functions of multiple pressure values along the boundary and at the center of each coarse block. The number of flux and pressure values at the boundary can be adjusted to improve the accuracy of simulation results and, in particular, to resolve important fine-scale details. Improvement over existing approaches is substantial especially for reservoirs that contain high-permeability streaks or channels. As an alternative, the multiscale mixed finite-element (MsMFE) method was designed to obtain fine-scale fluxes at the cost of solving a coarsened problem, but can also be used as upscaling methods that are flexible with respect to geometry and topology of the coarsened grid. Both methods can be expressed in mixed-hybrid form, with local stiffness matrices obtained as “inner products” of numerically computed basis functions with fine-scale sub-resolution. These basis functions are determined by solving local flow problems with piecewise linear Dirichlet boundary conditions for the DNR method and piecewise constant Neumann conditions for MsMFE. Adding discrete pressure points in the DNR method corresponds to subdividing faces in the coarse grid and hence increasing the number of basis functions in the MsMFE method. The methods show similar accuracy for 2D Cartesian cases, but the MsMFE method is more straightforward to formulate in 3D and implement for general grids.