An Efficient Multiscale Method for the Simulation of In-situ Conversion Processes

Summary Numerical modeling of the in-situ conversion process (ICP) is a challenging endeavor involving thermal multiphase flow, compositional pressure/volume/temperature (PVT) behavior, and chemical reactions that convert solid kerogen into light hydrocarbons and are tightly coupled to temperature propagation. Our investigations of grid-resolution effects on the accuracy and performance of ICP simulations demonstrated that ICP-simulation outcomes (e.g., oil/gas production rates and cumulative volumes) may exhibit relatively large errors on coarse grids, where “coarse” means a gridblock size of more than 3 to 5 m. On the other hand, coarse-scale models are attractive because they deliver favorable computational performance, especially for optimization and uncertainty quantification workflows that demand a large number of simulations. Furthermore, field-scale models become unmanageably large if gridblock sizes of 3 to 5 m or less have to be used. Therefore, there is a clear business need to accelerate the ICP simulations with minimal compromise of accuracy. We developed a novel multiscale-modeling method for ICP that reduces numerical-modeling errors and approximates the fine-scale simulation results on relatively coarse grids. The method uses a two-scale adaptive local-global solution technique. One global coarse-scale and multiple local fine-scale near-heater models are timestepped in a sequentially coupled fashion. At a given global timestep, the global-model solution provides accurate boundary conditions to the local near-heater models. These boundary conditions account for the global characteristics of the thermal-reactive flow and transport phenomena. In turn, fine-scale information about heater responses is upscaled from the local models, and used in the global coarse-scale model. These flow-based effective properties correct the thermal-reactive flow and transport in the global model either explicitly, by updating relevant coarse-grid properties for the next timestep, or implicitly, by repeatedly updating the properties through a convergent iterative scheme. Upon convergence, global coarse-scale and local fine-scale solutions are compatible with each other. We demonstrate the much-improved accuracy and efficiency delivered by the multiscale method by use of a 2D cross-section pattern-scale ICP simulation problem. The following conclusions are reached through numerical testing: (1) The multiscale method significantly improves the accuracy of the simulation results over conventionally upscaled models. The method is particularly effective in correcting the global coarse-scale model through the use of the fine-scale information about heater temperatures to regulate the heat-injection rate into the formation more accurately. The effective coarse-grid properties computed by the multiscale method at every timestep also enhance the accuracy of the ICP simulations, as demonstrated in a dedicated test case, in which a constant heat-injection rate is enforced across models of all investigated resolutions. (2) Multiscale ICP models result in accelerated simulations with a speed-up of four to 16 times with respect to fine-scale models “out of the box” without any special optimization effort. (3) Our multiscale method delivers high-resolution solutions in the vicinity of the heaters at a reduced computational cost. These fine-scale solutions can be used to better understand the evolution of the fluids and solids (e.g., kerogen conversion and coke deposition) in the vicinity of the heaters (several-feet-long spatial scale). Simultaneously, with the fine-scale near-heater solutions, the local-global coupled multiscale model provides key commercial ICP performance indicators at the pattern scale (several-hundred-feet-long spatial scale) such as production functions.

Fractured reservoirs are characterized by complex fractures on multiple scales and are quite difficult to model. Numerical simulation of fractured media is usually done based on dual-porosity model and equivalent continuum model. Dual-porosity model treats matrix system and fracture system as two parallel continuous systems coupled by crossflow function. This model is only valid for reservoirs with highly developed fractures. Equivalent continuum model treats fractured media as a continuum media and equivalent absolute permeability tensors for each grid block are calculated to describe the heterogeneity of the reservoir. This model is efficient only when there exist representative element volume and the equivalent permeability are difficult to decide. Both models treat the fractured media as a simplified model and cannot describe the multiscale flows exactly, because they cannot precisely consider the diversion effect of the fractures. Although the discrete fracture network (DFN) model can provide a detailed representation of flow characteristic, traditional numerical method does not suitable for DFN. The major difficulty is the size of the computation. A tremendous amount of computer memory and CPU time are required, and this can exceed limit of today’s computer resources. Upscaling methods are generally used to reduce the computational cost. However, it is not possible to have a priori estimates of the errors that are present when complex flow processes are investigated using coarse models constructed via simplified settings. In this paper, multiscale mixed finite element method (MsMFEM) is proposed to simulate water/oil two phase flow in discrete fracture media. By combining MsMFEM with the discrete fracture model, we aim towards a numerical scheme that facilitates fractured reservoir simulation without upscaling. The MsMFEM uses a standard Darcy model to approximate pressure and fluxes on a coarse grid. The multiscale basis functions are constructed numerically by solving local differential equations on the fine-scale grid. The advantage of MsMFEM is that the basis functions capable of reflecting information about fractures within elements. Therefore, this method can capture the fine-scale effects on the coarse grid, that is, multiscale method can reduce the computational cost and keep high calculation accuracy at the same time. Traditional numerical methods generally difficult to deal with complex grid element, in this paper, mimetic finite difference (MFD) method is used to construct the multiscale basis functions due to its local conservativeness and applicability of complex grids. Compared with traditional multiscale mixed finite element methods, this method is suitable for arbitrary complex grid system. This paper introduced fundamental principles of the multiscale mixed finite element method and described the numerical scheme of discrete fracture model based on MsMFEM in detail. Then we deduced discrete fracture model computing formulate for the multiscale basis function by using mimetic finite difference method. Oversampling technique is applied to get more accurate small-scale details. IMPES scheme is used in the two-phase flow simulation. Physical experiment is used to prove the validity the multiscale method. The numerical results show that compared with traditional numerical method, the MsMFEM can represent the fine-scale flow in fracture networks exactly and meanwhile has a higher computational efficiency.

Summary A robust and efficient simulation technique is developed on the basis of the extension of the mimetic finite-difference method (MFDM) to multiscale hierarchical-hexahedral (corner-point) grids by use of the multiscale mixed finite-element method (MsMFEM). The implementation of the mimetic subgrid-discretization method is compact and generic for a large class of grids and, thereby, is suitable for discretizations of reservoir models with complex geologic architecture. Flow equations are solved on a coarse grid where basis functions with subgrid resolution account accurately for subscale variations from an underlying fine-scale geomodel. The method relies on the construction of approximate velocity spaces that are adaptive to the local properties of the differential operator. A variant of the method for computing velocity basis functions is developed that uses an adaptive local-global (ALG) algorithm to compute multiscale velocity basis functions by capturing the principal characteristics of global flow. Both local and local-global methods generate subgrid-scale velocity fields that reproduce the impact of fine-scale stratigraphic architecture. By using multiscale basis functions to discretize the flow equations on a coarse grid, one can retain the efficiency of an upscaling method, while at the same time produce detailed and conservative velocity fields on the underlying fine grid. The accuracy and efficacy of the multiscale method is compared with those of fine-scale models and of coarse-scale models with no subgrid treatment for several two-phase-flow scenarios. Numerical experiments involving two-phase incompressible flow and transport phenomena are carried out on high-resolution corner-point grids that represent explicitly example stratigraphic architectures found in real-life shallow-marine and turbidite reservoirs. The multiscale method is several times faster than the direct solution of the fine-scale problem and yields more accurate solutions than coarse-scale modeling techniques that resort to explicit effective properties. The accuracy of the multiscale simulation method with ALG-velocity basis functions is compared with that of the local velocity basis functions. The multiscale simulation results are consistently more accurate when the local-global method is employed for computing the velocity basis functions.

A robust and efficient simulation technique is developed based on the extension of the mimetic finite volume method to multiscale hierarchical hexahedral (corner-point) grids via use of the multiscale mixed finite element method. The implementation of the mimetic subgrid discretization method is compact and generic for a large class of grids, and thereby, suitable for discretizations of reservoir models with complex geologic architecture. Flow equations are solved on a coarse grid where basis functions with subgrid resolution accurately account for subscale variations from an underlying fine-scale geomodel. The method relies on the construction of approximate velocity spaces that are adaptive to the local properties of the differential operator. A variant of the method for computing velocity basis functions is developed that utilizes an adaptive local-global algorithm to compute multiscale velocity basis functions by capturing the principal characteristics of global flow. Both local and local-global methods generate subgrid-scale velocity fields that reproduce the impact of fine-scale stratigraphic architecture. By using multiscale basis functions to discretize the flow equations on a coarse grid, one can retain the efficiency of an upscaling method, while at the same time produce detailed and conservative velocity fields on the underlying fine grid.The accuracy and efficacy of the multiscale method is compared to that of fine-scale models and of coarse-scale models with no subgrid treatment for several two-phase flow scenarios. Numerical experiments involving two-phase incompressible flow and transport phenomena are carried out on high-resolution corner-point grids that explicitly represent example stratigraphic architectures found in real-life shallow marine and turbidite reservoirs. The multiscale method is several times faster than the direct solution of the fine-scale problem and yields more accurate solutions than coarse-scale modeling techniques that resort to explicit effective properties. The accuracy of the multiscale simulation method with adaptive local-global velocity basis functions are compared to that of the local velocity basis functions. The multiscale simulation results are consistently more accurate when the local-global method is employed for computing the velocity basis functions.

Rapid and accurate simulation of the In-situ Conversion Process using upscaled dynamic models

Spectral multiscale finite element for nonlinear flows in highly heterogeneous media: A reduced basis approach

Summary A large number of multiscale methods have been developed based on the same basic concept: Solve localized flow problems to estimate the local effects of fine-scale petrophysical properties. Use the resulting multiscale basis functions to pose a global flow problem a coarser grid. Reconstruct conservative fine-scale approximations from the coarse-scale solution. By extending the basic concept with iteration cycles and additional local stages, one can systematically drive the fine-scale residual towards machine precision. Posed algebraically, this can be seen as a set of restriction operators for computing a reduced global problem and a set of prolongation operators for constructing conservative fine-scale approximations. Such multiscale finite-volume methods have been extensively developed for Cartesian grids in the literature. The industry, however, uses very complex with unstructured connections and degenerate cell geometries to represent realistic structural frameworks and stratigraphic architectures. A successful multiscale method should therefore be able to handle unstructured polyhedral grids, both on the fine and coarse scale, and be as flexible as possible to enable automatic coarse partitionings that adapt to wells and geological features in a way that ensures optimal accuracy for a chosen level of coarsening. Herein, we will discuss a compare a set of prolongation operators that can be combined with finite-element or finite-volume restriction operators to form different multiscale finite-volume methods. We consider the MsFV prolongation operator (developed on a dual coarse grid with unitary at coarse block vertices), the more recent MsTPFA operator (developed on primal grid with unitary flux across coarse block faces), as well as a simplified constant prolongation operator. The methods will be compared on a variety of test cases ranging from simple synthetic grids to highly complex, real-world, field models. Our discussion will focus on flexibility wrt (coarse) grids and tendency of creating oscillatory approximations. In addition, we will look at various methods for improving the methods’ convergence properties when used as preconditioners, as well as for generating novel prolongation operators. This is relevant for oil recovery because: Multiscale methods may provide a way to significantly speed up reservoir simulation and make previously intractable problems possible to solve. The extension of such methods to industry standard grids used for reservoir modelling enables the evaluation of the methods on real world models The construction of basis functions for multiscale methods may have direct connections to the process of upscaling rock derived properties such as transmissibility

We introduce a novel multiscale approach for reservoir simulation as an alternative to industry-standard upscaling methods. In our approach, reservoir pressure and total velocity is computed separately from the fluid transport. Pressure is computed on a coarse grid using a multiscale mixed-finite element method that gives a mass-conserving velocities on a fine subgrid. The fluid transport is computed using streamlines on the underlying fine geogrid.

Previous research has shown that multiscale methods are robust and capable of providing more accurate solutions than traditional upscaling methods. Multiscale methods solve the pressure equation on a coarse grid, but capture the effects from fine-scale heterogeneities through basis functions computed numerically from local single-phase problems on the underlying geocellular grid. Published results have so far been limited to simple Cartesian grids and/or incompressible flow. Here, we present a multiscale mixed finite-element method for three-phase black-oil flow on geomodels with industry-standard complexity. In particular, we discuss which effects can be incorporated in the multiscale basis functions and which effects should be modeled only on the coarsened simulation grid. Moreover, we describe how to handle degenerate hexahedral cells and non-matching interfaces that occur across faults. Finally, we present results of flow simulations on models of industry-standard complexity and demonstrate how multiscale methods can be used to simulate three-phase black-oil flow directly on high-resolution geomodels. The multiscale methods presented herein enable varying resolution and provide a systematic procedure for coarsening or refining the simulation model.

Compositional kinetics for hydrocarbon evolution in the pyrolysis of the Chang 7 organic matter: Implications for in-situ conversion of oil shale

Accurate multiscale finite element method for numerical simulation of two-phase flow in fractured media using discrete-fracture model

Summary Multiscale methods have been developed as a robust alternative to upscaling and to accelerate reservoir simulation. In their basic setup, multiscale methods use both a restriction operator to construct a reduced system of flow equations that can be solved on a coarser grid and a prolongation operator to map pressure unknowns from the coarse grid back to the original simulation grid. When combined with a local smoother, this gives an iterative solver that can efficiently compute approximate pressures to within a prescribed accuracy and still provide mass-conservative fluxes. We present an adaptive and flexible framework for combining multiple sets of such multiscale approximations. Each multiscale approximation can target a certain scale; geological features such as faults, fractures, facies, or other geobodies; or a particular computational challenge such as propagating displacement and chemical fronts, wells being turned on or off, and others. Multiscale methods that fit the framework are characterized by three features. First, the prolongation and restriction operators are constructed by use of a nonoverlapping partition of the fine grid. Second, the prolongation operator is composed of a set of basis functions, each of which has compact support within a support region that contains a coarse gridblock. Finally, the basis functions form a partition of unity. These assumptions are quite general and encompass almost all existing multiscale (finite-volume) methods that rely on localized basis functions. The novelty of our framework is that it enables multiple pairs of prolongation and restriction operators—computed on different coarse grids and possibly also by different basis-function formulations—to be combined into one iterative procedure. Through a series of numerical examples consisting of both idealized geology and flow physics as well as a geological model of a real asset, we demonstrate that the new iterative framework increases the accuracy and efficiency of the multiscale technology by improving the rate at which one converges the fine-scale residuals toward machine precision. In particular, we demonstrate how it is possible to combine multiscale prolongation operators that have different spatial resolution and that each individual operator can be designed to target, among others, challenging grids, including faults, pinchouts, and inactive cells; high-contrast fluvial sands; fractured carbonate reservoirs; and complex wells.

We propose a multiscale finite element method for solving second order elliptic equations with rapidly oscillating coefficients. The main purpose is to design a numerical method which is capable of correctly capturing the large scale components of the solution on a coarse grid without accurately resolving all the small scale features in the solution. This is accomplished by incorporating the local microstructures of the differential operator into the finite element base functions. As a consequence, the base functions are adapted to the local properties of the differential operator. In this paper, we provide a detailed convergence analysis of our method under the assumption that the oscillating coefficient is of two scales and is periodic in the fast scale. While such a simplifying assumption is not required by our method, it allows us to use homogenization theory to obtain a useful asymptotic solution structure. The issue of boundary conditions for the base functions will be discussed. Our numerical experiments demonstrate convincingly that our multiscale method indeed converges to the correct solution, independently of the small scale in the homogenization limit. Application of our method to problems with continuous scales is also considered.

Characteristics and quantitative models for hydrocarbon generation-retention-production of shale under ICP conditions: Example from the Chang 7 member in the Ordos Basin

A new multiscale method for wave simulation in 3D heterogeneous poroelastic media is developed. Wave propagation in inhomogeneous media involves many different scales of media. The physical parameters in the real media usually vary greatly within a very small scale. For the direct numerical methods for wave simulation, a refined grid is required in mesh generation to maintain the match between the mesh size and the material variations in the spatial scale. This greatly increases the computational cost and computer memory requirements. The multiscale method can overcome this difficulty due to the scale difference. The basic idea of our multiscale method is to construct computational schemes on two sets of meshes, i.e., coarse grids and fine grids. The finite-volume method is applied on the coarse grids, whereas the multiscale basis functions are computed with the finite-element method by solving a local problem on the fine grids. Moreover, the local problem only needs to be solved once before time stepping. This allows us to use a coarse grid while still capturing the information of the physical property variations in the small scale. Therefore, it has better accuracy than the single-scale method when they use the same coarse grids. The theoretical method and the dispersion analysis are investigated. Numerical computations with the perfectly matched layer boundary conditions are completed for 3D inhomogeneous poroelastic models with randomly distributed small scatterers. The results indicate that our multiscale method can effectively simulate wave propagation in 3D heterogeneous poroelastic media with a significant reduction in computation cost.