DEPARTMENT OF PETROLEUM ENGINEERING DWIGHT LOOK COLLEGE OF ENGINEERING TEXAS A&M UNIVERSITY

Model Calibration and Efficient Reservoir Imaging (MCERI)

A JOINT INDUSTRY PROJECT

Home
I.             Past Work

3-D streamline simulation is widely considered as one of the major developments in petroleum reservoir simulation and performance forecasting in the last decade. The technology has been rapidly assimilated by the industry for highly detailed flow simulation, reservoir management, history matching and uncertainty assessment. Several commercial reservoir simulators have been developed based on the streamline ‘time of flight’ concept introduced by Dr. Datta-Gupta and King (1995). With the advancement in high resolution data acquisition and seismic technolgies, geologic models now routinely consist of multimillion cells. This resulted in a widening gap between geologic modeling, flow simulation and uncertainty assessement. Streamline simulation has effectively bridged this gap. Our work developing streamline flow simulation and history matching of high resolution geologic models has been pivotal to these developments.

Links to Previous Research Works



Π.  Ongoing work

 ·        Fast Marching Methods: Application to Unconventional Reservoirs

The concept of depth of investigation is fundamental to well test analysis and is routinely used to design well tests and to understand the reservoir volume investigated. The depth of investigation can also be useful in identifying new well locations and planning, designing and optimizing hydraulic fractures in unconventional reservoirs. One way to calculate the depth of investigation is through numerical simulation. The numerical approach is very general and can be applied to arbitrary reservoir and well conditions. The computation time and expenses, however, make the numerical simulation approach unfeasible for routine applications.

As an alternative to numerical simulation, we can directly solve for the propagation equation for the pressure ‘front’ defined as the maximum of the pressure response for an impulse source.  Such propagation equation can be derived using asymptotic ray theory which has been used extensively in electromagnetic and seismic wave propagation (Virieux et al., 1994). The asymptotic method draws upon an analogy between the propagating pressure ‘front’ and a propagating wave front, and many of the concepts such as rays and propagating fronts have their counterparts in Petroleum Engineering (Datta-Gupta and King, 2007). Specifically, a high frequency asymptotic solution of the diffusivity equation leads to the following equation for a propagating pressure ‘front’ for an impulse source or sink (Vasco et al., 2000; Datta-Gupta and King, 2007):

Overview
News & Events
JIP Meeting
People
Research
Recent Project
Publications
Software
Members Only*
Contact

In the above equation,  is the diffusivity (consisting of  = permeability, = porosity, µ = fluid viscosity, ct = total compressibility) which can be function of position and  is the propagation time of the pressure ‘front’ (also called a ‘diffusive time of flight’). It simply tells us that the pressure ‘front’ propagates in the reservoir with a velocity given by the square root of diffusivity. This propagation depends on reservoir and fluid properties and is independent of flow rate. Also,  has unit of square root of time which is consistent with the scaling behavior of pressure diffusion. In fact,  is related to physical time through a simple expression of the form  where the pre-factor c depends on the specific flow geometry (Kim et al., 2009). For example, for radial flow, c = 1/4, for linear flow, c = 1/2 and for spherical flow, c = 1/6.

The pressure ‘front’ equation is a form of the Eikonal equation which appears in many contexts such as elastic and electromagnetic wave propagation; its properties are well developed in the literature (Sethian, 1999). Most importantly, the Eikonal equation can be solved very efficiently (in seconds opposed to hours of simulation time for comparable problems) by a class of solution called Fast Marching Methods (FMM) (Sethian, 1999). Following are three examples of depth of investigation: 1) generalization to heterogeneous medium; 2) generalization to anisotropic permeability and unstructured grid; and 3) a horizontal well with multistage hydraulic fractures.


Figure 1: Generalization of depth of investigation to heterogeneous medium (a) log permeability field (md) (b) depth of investigation in days (log scale) from FMM for a vertical well and (c) for a vertical well with an infinite conductivity fracture.

 

Figure 2: Generalization to anisotropic permeability and unstructured grid
 

                                                                          

Figure 3: Depth of investigation at various times for a horizontal well with multistage hydraulic fractures (a) well with hydraulic fractures (b) depth of investigation at 3 months (c) depth of investigation at 6 months (d) depth of investigation at 1 year.

·         Multiscale Parameterization and Model Calibration with Grid Connectivity Based Transforms (GCT)

The motivation of this topic is to expand upon the current state of the art in parameterization by linear transformation for history matching. In addition to theoretical developments, an emphasis is placed on the development of parameterization methods that are applicable to large reservoir models, on par with current industry standards and computational modeling capabilities, and also varying levels of prior geologic complexity.

                                             

Figure 4: Parameterization of high resolution permeability field with GCT

Methods of linear transform dependent on prior model assumptions have commonly utilized the parameter covariance matrix. The Karhunen-Loeve transform (KLT) or principle component analysis (PCA) has been widely used [Gavalas et al., 1976; Karhunen, 1947; Loéve, 1978; Reynolds et al., 1996; Li and Cirpka, 2006; Ma et al., 2008; Jafarpour and McLaughlin, 2009]. In this approach the estimable property defined at any grid cell is represented as the linear expansion of the weighted eigenvectors of the property covariance matrix. The eigenvectors form the transform basis and are typically ranked by their corresponding eigenvalues, from largest to smallest, that are related to the variance contribution of each eigenvector to the total parameter variance. The expansion is optimal in the mean-squared-error sense among the linear class of transformations; therefore, the KLT coefficients present the fewest number of parameters that capture the maximum amount of variation for any low-rank approximation. Accordingly, the sorted eigenvectors convey larger to smaller scales of spatial parameter variation. A limitation of this approach is that in realistic problems with high-resolution models, the covariance can be unknown (or uncertain), resulting in misleading basis functions when prior model assumptions are incorrect [Jafarpour and McLaughlin, 2009]. Eigendecomposition of the covariance matrix, required for each parameter update, is also prohibitively expensive for high-resolution models relative to current computational capability. Last, the preservation of moments beyond second-order (i.e., covariance), as well as complex continuous geologic structures (e.g., channel systems), are not guaranteed in the KL expansion during the parameter updates.

A new model-independent basis constructed from grid-connectivity information is developed here. The parameterization is applicable to any grid structure and domain geometry. The development of the grid-connectivity-based transform (GCT) basis begins from first principles, merging discrete Fourier analysis and spectral graph theory. The basis functions represent the modal shapes or harmonics of the grid, are defined by a modal frequency, and converge to the discrete cosine transforms (DCT) for certain grid geometries and boundary assumptions; therefore, reservoir model calibration is performed in the spectral domain. Using an adaptive multiscale workflow, the GCT parameterization is successfully applied for history matching of several synthetic reservoir models of varying geometry and geologic complexity and also of a field case.

                                                           

Figure 5: Multiscale Model Calibration Workflow Using GCT

 

·        A Hierarchical Multi-scale Approach to History Matching

We develop a hierarchical multiscale history matching approach with two levels of hierarchy for inversion in parameter space and geometric space. The approach follows a sequence of calibrations from global to local parameters in coarsened and fine scales. At first, we identify the ‘heavy hitters’ in the large scale static and dynamic parameters and then calibrate them using an evolutionary algorithm. This global parameter calibration, matching and balancing field level energy, is followed by streamline-assisted multiscale inversion to match well by well production history with local parameter updates. Starting with a relatively coarse grid, we match the production data at wells by dynamically refining the reservoir grid. This multiscale data integration offers not only benefits of enhancement in computation efficiency and but also effective iterative minimization. The approach has been successfully applied to a variety of field cases in collaboration with operating companies.

 

Figure 6: Overview of hierarchical multiscale history matching

 

·        Flood Rate Optimization Using Streamlines

It is well recognized that field-scale rate optimization problems often involves highly complex reservoir model, production and facility related constraints and geological uncertainty. Deployment of smart well completions with inflow control valves (ICV) to control production/injection rates for various segments along the wellbore further compounds the complexity of the problem. All these make optimal reservoir management via rate control difficult without efficient optimization algorithms.

Alhuthali et al. (2007, 2008, 2010) proposed a streamline-based approach to maximize waterflood sweep efficiency. The approach focuses on equalizing arrival times of the waterfront at all producers within selected sub-regions of a waterflood project. This results in delayed water breakthrough and reduces field water cut after water breakthrough. The optimization was performed under production and facility constraints using a sequential quadratic programming approach. A major advantage of the streamline-based approach is the analytical computation of the sensitivities of the waterfront arrival times at the producers with respect to well production/injection rates and also, analytic calculations of the gradient and Hessian of the objective function. This makes it computationally efficient and suitable for large field cases. It also accounts for geological uncertainty via a stochastic optimization framework based on a combination of the expected value and variance of a performance measure from multiple realizations. We have also studied the tradeoff between maximizing sweep efficiency and accelerating production (maximizing NPV) using an augmented objective function. In addition to waterflooding, we have applied this method to EOR processes like polymer and CO2 flooding.

Figure 7: We are showing application of streamline based optimization to a giant middle eastern field. On left base case injection (blue) and production rates (green) are shown. On right these rates are shown after optimization. Optimization results in higher rate allocation to producers in the center.

 

·        Robust Streamline Tracing in Unstructured and Locally Refined Grids

We have developed robust, fast, and easy-to-implement streamline tracing algorithms to be used in unstructured PEBI grids and grids with locally refinement (LGR) and faults. The proposed algorithms work as a post processing procedure of regular finite difference/finite volume simulation. The calculated streamline trajectories and time of flights can be used for visualization, streamline-based history matching, streamline-based reservoir management, and many other applications. The PEBI grid tracing algorithm proposed assumes no constitutive properties and maximizes the velocity field continuity. The LGR and faulted cell tracing algorithms proposed are based on “local boundary layer approximation” which conserves fluxes, computationally inexpensive, and robust for various geometric configurations.

Figure 8:  Tracing in 2.5D PEBI grid

 


Figure 9:  Nested LGR grid with multilateral well completed in the LGR region 
 

Figure 10: Streamlines near a branch of a multilateral well

 

·        Improved Upscaling for Flow Simulation of Tight Gas Reservoir Models

Tight gas and shale gas reservoirs provide almost half of the current U.S. domestic gas production, with significant projected increases in the next several decades. Reservoir simulation, together with detailed 3D geologic models, may be used for improved management, but only if the flow simulation through these models can capture the essential heterogeneity and 3D continuity of these intermittently connected rock packages. The current work will examine tight gas reservoirs, but we recognize that similar issues will arise in source rock plays.

The current approach has three key technical elements: (1) We explicitly preserve the local continuity of the reservoir sands in the geologic model through the design of the simulation grid (Fig.1). (2) We have developed a “Well Index” based upscaling, to determine the upscaled permeability and which preserves the local reservoir quality. (3) The heterogeneity is preserved through the use of transmissibility upscaling, which we show performs systematically better than the usual permeability upscaling. In combination, these three elements provide simulation results which are almost indistinguishable from the fine scale model (Fig.2).

As a posteriori analysis, well test depth of investigation predicts local horizontal & vertical connectivity within the model, which will lead our future areal adaptive work.



Figure 11:  Vertical Adaptive Grid Design: Color shows final pressure after 24 years depletion. Cells follow the pay/non-pay contrast in each column of the model.

                         Figure 12:  Recovery Curves Comparison: The recovery curve of 1x1xN adaptive coarsening (with transmissibility based upscaling) overlaps with the fine scale model