For many field-scale reservoir simulation models, aquifer regions (both peripheral and bottom) can represent a large percentage of total active grid blocks. The objective of this work is to build an accurate method respecting the fine-scale heterogeneities, but significantly reduce the computational cost of keeping the aquifer grids in the model, where they are deemed necessary to better describe the variable strength of water influx.
The present method can be considered a specialized two-level grid coarsening or upscaling method. Within the aquifer regions, each column of fine cells are coarsened vertically based on fine-scale z-transmissibility. A coarsened column may consist of a single amalgamated aquifer cell or multiple vertically disconnected aquifer cells separated by flow barriers. The pore volume (PV), compressibility, and lateral flow terms of the coarse cell are restricted from the fine-grid cell. The lateral connectivities within the aquifer regions and between the aquifer and the reservoir are honored, inclusive of the fine-scale description of faults, pinchouts, and null cells. Reservoir regions are not coarsened. Two alternatives exist for the fine-scale pressure reconstruction from the coarse-grid solution. The first method uses the vertical equilibrium concept. The second method performs a 1D solution for the fine-scale pressure. A spatially variable multi-porosity permeability model can also be accommodated. In parallel implementation, besides the load balancing method, the memory balancing method is also necessary. All these aspects will be explained in the paper.
The method has been applied successfully to several complex full-field simulation models where the transient aquifer water influx has been identified as a key factor. These models include dual porosity, dual permeability models with a complex geologic description. We compare simulation results using the accelerated aquifer calculation method against the original fine-grid simulation model without acceleration and demonstrate the efficacy of the new method. For the models tested, the actual speedup factors achieved are also compared against the theoretically computed improvement factors.
The dual grid method strives to significantly reduce computational cost but retain the fine-scale heterogeneity data to accurately represent the water movement within the water zones of the simulation model. The method differs from the standard upscaling and grid coarsening method where the coarse-grid properties are computed a priori. Instead, the fine-scale information is restricted to the coarse grid during Newton's iteration to represent the fine-scale flow behavior. This is analogous to the multiscale method, but we apply this only vertically within the aquifer regions.
Ma, Xiang (ExxonMobil Upstream Research Company) | Hetz, Gill (Texas A&M University) | Wang, Xiaochen (ExxonMobil Upstream Research Company) | Bi, Linfeng (ExxonMobil Upstream Research Company) | Stern, Dave (ExxonMobil Upstream Research Company) | Hoda, Nazish (ExxonMobil Upstream Research Company)
Many recent developments in generating history matched reservoir models that approximately characterize subsurface uncertainty are associated with the ensemble smoother (ES) method. It is much better suited for practical history matching applications because it does not require updating of the dynamical variables and thus the frequent simulation restarts required by ensemble kalman filter (EnKF) are avoided. However, the performance of original single update scheme of ES is poor for strongly nonlinear problems and therefore iterations may be needed. Several iterative forms of ES were proposed in the past few years, most of which combine ideas from random maximum likelihood (RML) and ensemble-based techniques. Unlike previous implementations, we pose the history matching problem as a full nonlinear least squares optimization problem and classical Levenberg-Marquardt (LM) algorithm is used as the optimization solver. By showing the that solution of the linearized least squares subproblems arising from each iteration has similar structure to that of standard ES update equation, we propose to use ES as the linear least squares solver to avoid the expensive adjoint calculation. In this way, the proposed algorithm can be considered as an iterative ES and the regularization parameter can be updated following the standard LM rule. Furthermore, because it is casted as an optimization problem, it is straightforward to extend it to robust nonlinear least squares method that can automatically estimate the measurement noise level and reduce the effect of outliers in the data that is essential for field applications. Two synthetic reservoir models are used to showcase the effectiveness and robustness of the newly developed algorithm.
Wang, Shihao (Colorado School of Mines) | Zhang, Juncheng (CNPC Chuanqing Driling Company) | Yang, Zhenzhou (CNPC USA) | Yin, Congbin (CNPC Chuanqing Driling Company) | Wang, Yonghong (CNPC USA) | Zhang, Ronglei (Colorado School of Mines) | Winterfeld, Philip (Colorado School of Mines) | Wu, Yu-Shu (Colorado School of Mines)
We present the development and application of a multi-physical simulator for evaluating the combined thermal-hydraulic-mechanical behaviors of petroleum reservoirs. The simulator combines non-isothermal multiphase compositional modeling with coupled geomechanical simulation module. The simulator consists of two major modules, namely, the fluid and heat flow module and the geomechanical module. An isenthalpic flash calculation approach is implemented in the fluid and heat flow module. In the flash calculation module, a nested approach is adopted, in which PT flash calculations are conducted in the inner loop and temperature is updated in the outer loop. The iteration is continued until both the fugacity and energy stopping criteria are satisfied. An improved version of the Beltrami-Michell equation, called extended Beltrami-Michell equation, has been derived and implemented in the geomechanical simulation module to simulate heterogeneous and plastic behavior of formation rocks. The three normal stress components inside the stress tensor are solved simultaneously with the pressure and enthalpy in the fluid/heat module, ensuring the mass/energy conservation. The newly-derived extension of the Beltrami-Michell equation is capable of handling materials with changing mechanical properties. This way, the simulator is able to capture the phase change as well as the poro-mechanical effects on rock deformation induced by fluid injection/extraction. The multi-physics simulator is built on an object-oriented parallel simulation framework, with a speedup factor up to hundreds.
Reservoir simulation plays an important role in the petroleum industry. Today, there is a specific demand to run ensembles of megaand even giga-cell models. The iterative solution of the large systems of nonlinear governing equations, which describe the multiphase mass transfer in the subsurface, takes the most of the simulation time. The linearization part of the solution process occupies a significant fraction of that time, especially in compositional models. Moreover, the implementation of the linearization step usually embodies the most substantial, complex, and specific part of the computational loop in modern simulators, defining which physical mechanisms and assumptions are employed. This significantly complicates the implementation of simulation codes for heterogeneous computing hardware, which promises significant improvements in simulation time. In this work, we use the recently proposed Operator-Based Linearization (OBL) approach to develop a general purpose reservoir simulation code aiming to substantially decrease the simulation time. OBL offers a simplified linearization method, enhancing the computational performance of simulation and providing an opportunity of a painless porting to heterogeneous computing architectures. To distinct the contribution of both factors, we developed two versions of the compositional simulation prototype code: for traditional CPU and GPU-accelerated hardware architectures. While the former allowed us to speed the linearization stage up by an order of magnitude in comparison with the conventional approach, based on Automatic Differentiation (AD), the latter improved it further by another order of magnitude. The developed prototype realizes the potential of the OBL approach and GPU computing architecture, proving significant improvement in general purpose simulation performance.
This paper presents numerical simulation results of pulse testing experiments carried out at a test site of a carbon capture and geological storage project in Mississippi, USA. The primary objective of this study is to validate the effectiveness of pulse testing as a monitoring tool for detecting potential CO2 leakage pathways with application to the test site. Detrending followed by Fourier transform is adopted to decompose sinusoidal pressure anomalies induced by a periodic injection of CO2 into frequencies used as target parameters of history matching. The secondary objective is to calibrate the geologic model of the test site by reducing the discrepancy between observed and simulated Fourier parameters and assess uncertainties associated with the compositional brine-CO2 flow. An assisted history matching tool that mounts global- and multi-objective evolutionary algorithms is developed, integrated with an in-house flow-geomechanics simulator, and employed to manage pulse testing simulations with a low computational cost in high-performance parallel computing environments. Grid cells in the test site are locally refined using enhanced-velocity that allows nonmatching grids on interfaces between subdomains. Experiments performed with one pulser well and two monitoring wells under steady-state conditions are considered baselines for subsequent experiments that convert one monitoring well into a production well as an artificial CO2 leakage pathway. The difference between the pressure anomalies obtained from the baseline and leak experiments are captured as a signal of CO2 leakage detection with reliability in the simulation results. A trade-off relationship between the matching qualities at the two monitoring wells is revealed more clearly by invoking multi-objective history matching than conventional global-objective history matching. This comparative study to investigate the significance of multi-objective optimization in subsurface characterization represents that diversity-preservation in the ensemble composed of qualified geologic models has the advantage of reducing bias for uncertainty quantification.
Given the high viscosity of the oil, bitumen from oil sands reservoirs in western Canada is recovered by using steam which, due to its temperature, lowers its viscosity. One of the key issues faced by the operators is the steam conformance of the depletion chamber around wells. The greater the fingering phenomena of steam at the edge of chamber, the worse is the chamber uniformity and utilization of the well, and the greater are the green house gas emissions and water use per unit oil recovered. Fingering has long been explained as the penetration of steam phase into the oil phase which arises from an unfavourable mobility ratio. In this paper, we introduce linear instability analyses (Orr-Sommerfeld and Rayleigh-Taylor/Saffman-Taylor instability) of the interface between steam and oil layers and conduct a series of numerical simulations to reveal that fingering in the steam-assisted heavy oil recovery at the top of the steam chamber is created due to solution gas exsolution whereas fingering at the chamber edge is due to viscous shear instability. The results show that non-ideal steam conformance is inevitable even in homogeneous reservoirs.
Earth models are important tools for support of decision making processes for optimal exploitation of subsurface resources. In earth modelling applications, massive amounts of information result from collecting and interpreting measurements and simulating physical scenarios at various scales and resolutions. The size of an earth model grid is a consequence of its scale and resolution. The grid size is a major factor for the computational efficiency when managing the model and performing various model based simulations. The grid is controlled by the geological structure, and the resolution of the grid is therefore controlled by the resolution of the geological structure (the number of structural elements).
We discuss principles for a multi-resolution framework for earth model gridding where the aim is to allow automatic, local control of the resolution of a populated grid by locally controlling the resolution of the geological structure. The structure is ordered in a hierarchical manner, and splits the subsurface into a nested set of regions that are ordered in a corresponding hierarchical fashion. Each region can be separately gridded. The method is based on a recent method for local updates of the connectivity (topology) of the geological structure and of the resolution of a populated grid. Similar to existing strategies for multi-resolution modelling in other sciences and applications, the proposed method allows local control with the trade-off between numerical accuracy and grid size. The method should allow to locally change the model resolution even during a modelling exercise; subsurface volumes of high interest can be represented at fine resolution, whereas in volumes of less interest, details can be omitted. The method also enables local updates and local scale uncertainty management of the geological structure.
The aim of the proposed method is to develop an effective methodology that supports real-time earth model based workflows such as geosteering where multiple grid realizations are never fixed and always updated with the most recent measurements and interpretations, and where each realization is always kept at an optimal resolution.
This paper presents new coupled three-phase relative permeability and capillary pressures models for implementation in a four-phase flow compositional EOS simulator. Identification of hydrocarbon and aqueous phases based on their molar Gibbs free energy is a key feature of the new model. Instead of using phase labels (gas/oil/solvent/aqueous) to define parameters such as relative permeability endpoints, residual saturations and exponents, the parameters are continuously interpolated between reference values using the Gibbs free energy of each phase at each time step. Models independent of phase label have many advantages in terms of both numerical stability and physical consistency. The models integrate and unify relevant physical parameters including trapping number and hysteresis into one rigorous algorithm for application in numerical reservoir simulators. The robustness of the new models is demonstrated with simulations of miscible WAG and solvent stimulation to remove condensate and/or water blocks in both conventional and unconventional formations.
In this paper, we formulate an adaptive time-stepping strategy in a control system framework for efficient reservoir simulation. Current simulators achieve a trade-off between stability and computational effort by means of heuristically-based time-stepping strategies during the simulation runs of complex reservoirs. This can lead to larger computational time and unrealistic time-steps. To this end, we assess the effectiveness of a PID (proportional-integral-derivative) controller at regulating the time steps used during the time discretization of a simulator. We also evaluate the robustness of a PID controller using various metrics - computing time, total number of time steps, and number of rejected steps - and compare it to adaptive time-stepping techniques currently used in the industry. Feedback control strategies are extensively applied in the upstream and downstream of the oil and gas industry in the sense of controlling fluid flow displacement. In particular, PID controllers are preferred for their ease of implementation and robustness. We test our methodology using a two-phase flow reservoir simulator. The goal is to determine if the use of the PID controller can significantly reduce the computational cost while maintaining the same level of accuracy. The results show that the PID controller can improve the efficiency of the reservoir simulator as long as it is tuned properly. We recorded a lower rise time in the PID time-step response, a lower number of total time-step used for the time discretization and a near identical accuracy in the solution obtained. In most of the scenarios run, the PID controller consistently outperformed the basic controller while achieving the same accuracy in the pressure and saturation solutions. This trend was observed in a wide range of permeability and well locations.
Standard reservoir simulation schemes use single-point upstream weighting for computing the convective fluxes when multi-phase/component fluids are present. These schemes are only first-order and may cause high viscosity effect. Second-order schemes provide a better resolution and reduce the smoothing near the shocks. In reservoir simulation practice, implicit discretisations capable of taking large time steps are preferred over explicit schemes. However, even for a simple first-order method, solving a large nonlinear system is often very expensive; Extra coupling and nonlinearity of the discretised equations can be introduced from a higher-order spatial discretisation. It has been shown that strong nonlinearity as well as lack of continuous differentiability in numerical flux function and flux limiter can cause serious nonlinear convergence problem.
Cell-centered finite-volume (CCFV) discretizations may offer several attractive features, especially for fluid flow in heterogeneous or fractured porous media. The objectives of this work are to develop a fully-implicit CCFV framework that could achieve second-order spatial accuracy on smooth solutions, while at the same time maintain robustness and nonlinear convergence performance. We develop a novel multislope MUSCL method which interpolates the required values at the edge centroids in a simpler way by taking advantage of some geometric features of the triangular mesh. Through the edge centroids, the numerical diffusion caused by mesh skewness is significantly reduced, and optimal second-order accuracy can be achieved. An improved gradually-switching piecewise-linear flux-limiter is introduced according to mesh non-uniformity in order to prevent spurious oscillations. The smooth flux-limiter can achieve high accuracy without degrading nonlinear convergence behavior.
For the discretization of pressure and Darcy velocities, a mimetic finite difference method that provides flux- continuity and an accurate total velocity field is used. The fully-coupled MFD-MUSCL framework is adapted to accommodate a lower-dimensional discrete fracture-matrix (DFM) model. Several numerical tests with discrete fractured system are carried out to demonstrate the efficiency and robustness of the numerical model. The results show that the high-order MUSCL method effectively reduces numerical diffusion, leading to improved resolution of saturation fronts compared with the first-order method. In addition, it is shown that the developed multislope method and adaptive flux limiter exhibit superior nonlinear convergence compared with other alternatives.