A New Simulation Layer Optimization and Permeability Upscaling Method for Preserving Critical Reservoir Heterogeneity

Two modeling approaches are commonly utilized for simulating flow in fractured formations: the discrete fracture network (DFN) approach and the stochastic continuum (SC) approach. Although the DFN approach is the most accurate, it has computational and memory constraints. The SC approach ensures fast processing but results in system over-homogenization. The fracture continuum (FC) approach arises as an integrated technique that incorporates the merits of both approaches. The main objective of this research is to develop a computationally efficient technique based on the FC approach to simulate the radial groundwater flow towards wells through two-dimensional fractured media under both steady and transient conditions. A stochastic generation of the DFN is performed in a Monte Carlo framework taking into account wells positioning. The DFN flow system is solved by applying the mass balance equation at fracture intersections. Fracture segments are mapped onto grids of 1 × 1 m and 5 × 5 m resolution as conductivity and specific storage cells. The grid flow problem is solved via MODFLOW. Flow and head discrepancies between the proposed technique and the DFN approach (reference solution) are assessed in steady and transient conditions. A grid-conductivity correction is needed to preserve the DFN flow in the presence of wells. A porosity estimation is proposed to identify the grid-pressure transient response. Promising flow and head results are observed for fine and coarse grid models. Some of the studied cases show large discrepancies in the maximum drawdown obtained on the coarse grid model. Accordingly, a new technique is proposed to handle such discrepancies and is found efficient in transient simulations (e.g., 11% and 26.12% discrepancies are minimized to − 0.93% and − 1.03% for two studied cases). The adopted mapping technique is found efficient when the interest is to estimate the average drawdown over an aquifer as correlation coefficients of 0.89 and 0.97 are found for the coarse and fine grid models, respectively when compared to the DFN model. However, the technique has limitations in estimating the drawdown at locations of wells.

In this study, various upscaling techniques and their effects on Barik tight gas Formation simulation modelling results were investigated. The intent of this upscaling study is to recommend coarse models that provide approximately the same flow behavior or well performance as the fine grid model. The study will help to develop an awareness of the range of applicability for the upscaled coarse models. It will also allow coarse grid models to be used appropriately and with greater confidence in production optimization. These techniques comprised several alternative ways of grouping reservoir data, with respect to petrophysical rock types (herein referred to as RT). This scheme defines 5 individual rock types, with RTs 1 to 4 broadly defining pay and RT 5 defining non-pay. The layers in the simulation models are made up of single or multiple grouped RTs from the same zone. Keeping each layer in the model, ordered as it was in the well log, resulted in 85 simulation layers for the fine grid model. The upscaled or coarse models have 16 to 35 simulation layers, with the smallest being the model where RTs 1 to 4 are grouped together. Different RTs sorting were tested in each of the upscaled models. The results of this study suggest that re-ordering the log information so that all of the rock for each rock type was grouped together inside each stratigraphic unit appears to be an acceptable upscaling technique that gives reasonable efficiency and accuracy. This type of upscaling was required since the study showed that any upscaling method that averaged the properties of RT 1 with the properties of other rock types within the same simulation layer could result in optimistic EUR estimates of up to 30% relative absolute error when compared to a fine grid model. The choice of upscaling method is particularly significant in a low Kv/Kh ratio environment, but less so within environments with a high Kv/Kh ratio. The upscaling technique has only been tested for the Barik formation and should not be used elsewhere without proper testing. The upscaling technique works for the Barik formation because of several distinctive features of the reservoir and the wells. In applying this upscaling technique we assume that: The wells are all hydraulically fractured and thus the flow is generally horizontal into the fractures and from the fractures into the well.The fracture extends from the top to the bottom of each stratigraphic unit that the fracture encounters. If a fracture only partially penetrates a stratigraphic unit, this method may not work.The gas is relatively dry and thus the flow is less affected by gravity than for reservoirs containg flowing liquids.The reservoir is made up of rock with a wide range of permeability, but the flow is dominated by the well connected, higher permeability rocks (RT 1).The reservoir is very stratified, and the correlation length of the rock types has to be very much greater than the well spacing. Re-ordering of the RTs in a given upscaling method will result in acceptable accurate estimates of hydrocarbon recovery, compared to the fine grid model. It showed that the order of layering did not matter for the area studied, because a conductive fracture connects all the layers. This method probably really is only applicable for dry gas reservoirs (where gravity is not important) and in fractured wells (where horizontal flow into the fractures and into wells dominates). In oil reservoirs and rich condensate reservoirs where gravity is an important factor, the ordering of the rock types within the simulation layer may matter. It will also be shown that an approximately 90% reduction in the simulation modelling computing time could be achieved if the appropriate upscaling technique is used. To achieve this reduction in computing time, some compromises were made, including assuming RT 4 is non-pay in upscaled models where RTs 4 & 5 are grouped together in the same simulation layers, resulting in reduction of HCPV. It is important to mention that some RT 5 zones in the log have thin instances of other rock types, which are not accounted for in the upscaled models and could result in an error in average pore volume preservation.

Thermo-mechanical modelling of solar central receivers: Effect of incident solar flux resolution

Hybrid steam and in-situ combustion recovery processes have shown advantages over pure steam injection for recovery of oil sands resources, particularly with respect to reducing costs and lowering requirement for water and natural gas. However, it has been very challenging to predict field performance of hybrid steam and combustion processes with a reasonable degree of confidence. Usually, a combustion front has a thickness of only a few inches and high resolution grids are required to capture the steep temperature, saturation and fluid composition gradients in the vicinity of the combustion front. Using high resolution, fine grids to improve accuracy of simulation requires excessive computation time and, therefore, may be impractical for field scale modelling. It is important to have a robust simulation tool to accurately predict reservoir performance without compromising the computational efficiency. In this work, numerical modeling of a hybrid steam and combustion recovery process was performed in a typical Athabasca oil sands reservoir. A comprehensive new reaction kinetics model derived from laboratory results was incorporated to represent the complex chemical reactions in the combustion process. The hybrid recovery process utilized oxygen enriched air co-injection after several years of SAGD operation. In the numerical model, safe limits were set on producing well temperature and oxygen content of produced fluids. The initial grid size in the numerical model was at the centimeter scale resulting in large run time, and thus, in order to improve the computational efficiency, a dynamic gridding feature was applied. Parameters for controlling the creation of a dynamic grid and subsequently reverting back to the coarse grid have been examined in order to properly trigger the dynamic gridding feature in the model. Once the optimized dynamic gridding parameters were determined, several different well configurations were investigated. Comparisons were made between SAGD and hybrid steam/combustion processes in terms of cumulative water (steam) injection, cumulative oil production, and a steam-oil ratio. By comparing the simulation results from the fine grid model and the dynamic gridding model, it has been found that the temperature gradient is the best criterion to use for controlling dynamic gridding compared to fluid saturation and/or composition criteria. The threshold value for the temperature criterion was determined to be 35°C. The model locates the fine grids in close proximity to the combustion front where the temperature and fluid saturation gradients are the steepest and it places the coarse grid blocks elsewhere in the model. Comparisons are made between the computation time and the accuracy of the simulation and these demonstrate that dynamic grid amalgamation reduces the computation time significantly while maintaining reasonable computation accuracy of simulation. Compared with SAGD, the hybrid steam/in-situ combustion process reduced cumulative water usage (steam injection) by 20% to 27%, while the cumulative oil production remained the same. This paper provides a workflow for modelling of hybrid steam and combustion processes. Also, it is expected that this work will provide insights for field design of these hybrid thermal recovery processes.

Multiple methodologies and grid geometries were applied to model the production performance and uncertainties in the planning and execution of a pilot horizontal well with multiple hydraulic fractures. The well is located in a tight carbonate formation and the fracture will be induced hydraulically with acid. A methodology was established to scale the observations of the single well model to a full field scale. This, when calibrated against actual performance from the pilot well, will enable improvements in the planning the development of the target reservoir at full field scale. The methodologies described in this paper consist of fine grid models to explicitly represent the fractures, a method to calibrate and represent the production performance on a coarse grid full field model. In the fine grid model, the grid was geometric refined to fracture aperture grid size to capture the physics and flow behavior induced by the acid-fracs. Fracture plane conductivity is based on laboratory measurement of acid-frac experiments. The fine scale model was also used in sensitivity analysis to optimize the pilot well design such as number of fractures and desired fracture height growth, identify risk of water production and help formulate of mitigation strategies. To enable modeling at full field scale – which is impractical with fine grid models – a novel methodology was established to calibrate a coarse grid field scale model to consistently reproduce the fine grid model behavior. This enables early planning of full field development should the performance of the pilot well support the application of multi acid-frac horizontal well as the preferred technology.

The process of evaluating and assessing a reservoir early in its development is rather difficult and uncertain although the results of the evaluation may have far-reaching consequences for the reservoir. This paper proposes and describes a methodology for enhancing this assessment. The method involves several basic steps including: identification of key heterogeneities, their quantification, the development of models representing the full breadth of variations in heterogeneity, the detailed modeling of limited areas of the reservoir to assess the effects of the model and geologic uncertainties, and the use of these results for calibrating up-scaled full-field models. The methodology is applied for a gas condensate reservoir appraised by only two wells in which the choice of recovery process (primary depletion or gas-cycling) has not yet been made. Given that the key heterogeneity is unknown at this early stage, stochastic modeling is then undertaken on a grid emphasizing vertical discretization. In order to provide contrasted models clearly showing the uncertainty associated with the geological model, two "end-member hypotheses" are used, based on the uncertainty as to the frequency and range of the heterogeneities. Simulations carried out on only part of the reservoir show how sensitive the results are to the range of the heterogeneities. This is illustrated particularly clearly in the gas-cycling option in which gas injection is much more effective in the case of low range heterogeneities than for long range heterogeneities. After optimizing the clustering of fine grid cells into coarse grid cells, it is demonstrated that simulation results obtained on a fine grid model can be reproduced with a significantly smaller number of layers.

A flood inundation modelling using v-support vector machine regression model

Baram field is one of the nine major oil fields in Baram Delta Operation (BDO) area which is located 28 km northwest of Lutong, in the offshore area of Miri, Sarawak in East Malaysia. The field has been producing since 1967 with multiple stacked sandstone reservoirs in a shallow offshore environment. Based on the BDO EOR feasibility study conducted in 2004 the field has been selected as one of the major EOR candidates in the area. This paper describes a successful application of the PETRONAS upgridding technique in Baram field. The so called "Optimum Grid Coarsening Scheme" uses fine grid static and dynamic reservoir variability to identify the non-uniform grid cell configuration at the coarse grid model. Averaged saturation equation is the key driver in calculating the sub-grid properties in the fine level that subsequently will be the basis for determining the coarse grid scheme. It considers the balance of all forces including gravity, viscous and capillary that contribute to fluid flow in the reservoir. The upscaled model has been successfully history matched with the past 40 years of production and pressure history in a timely manner with limited adjustment of the static parameters. The history matched model was then used in our study to investigate many possible development scenarios in order to improve the ultimate recovery factor of the field. In our study, Water Alternative Gas (WAG) injection has shown superiority over other studied methods as the most promising EOR scenario for the field. The study shows that the recovery factor of the S sand in Baram can be increased by an additional of 16% above infill drilling and 7% above water flooding options using WAG application.

Flow fields in a centrifugal fan for an indoor unit of an air-conditioner were calculated with finite element method-based large eddy simulation (LES) with the aim of predicting fan performance and aerodynamic noise in this study. The numerical simulation code employed throughout the LES was called FrontFlow/blue (FFB). We compared 10M grid [coarse grid] and 60M grid [fine grid] calculation results for investigation of influence of grid resolution. In the fine grid, the number of grid elements in blade-to-blade direction, and of region between the shroud and the bell mouth increased in particular. By calculating with the fine grid, calculated distributions of absolute velocities at blade exit reasonably agreed with experimental results. Because of this, maximum absolute velocity by fine grid near hub decreased as compared to those by coarse grid. Calculated sound pressure level by fine grid was therefore smaller than that by coarse grid, and the overestimation of sound pressure was suppressed by calculating with fine grid. This decrease of the absolute velocity was a first factor for the improvement of calculation accuracy. Moreover, number of captured streaks on the blade, hub, and shroud surfaces by fine grid increased as compared to those by coarse grid. As a result, size of streak by fine grid became smaller than that by coarse grid. Static pressure fluctuations by fine grid on the blade, hub, and shroud surfaces therefore reduced as compared to those by coarse grid. Aerodynamic noise was related to static pressure fluctuations according to Curle’s equation. This reduction of static pressure fluctuations was therefore a second factor for improvement of calculation accuracy.

When reliable a priori information is not available, it is difficult to correctly predict near-surface S-wave velocity models from Rayleigh waves through existing techniques, especially in the case of complex geology. To tackle this issue, we have developed a new method: two-grid genetic-algorithm Rayleigh-wave full-waveform inversion (FWI). Adopting a two-grid parameterization of the model, the genetic algorithm inverts for unknown velocities and densities at the nodes of a coarse grid, whereas the forward modeling is performed on a fine grid to avoid numerical dispersion. A bilinear interpolation brings the coarse-grid results into the fine-grid models. The coarse inversion grid allows for a significant reduction in the computing time required by the genetic algorithm to converge. With a coarser grid, there are fewer unknowns and less required computing time, at the expense of the model resolution. To further increase efficiency, our inversion code can perform the optimization using an offset-marching strategy and/or a frequency-marching strategy that can make use of different kinds of objective functions and allows for parallel computing. We illustrate the effect of our inversion method using three synthetic examples with rather complex near-surface models. Although no a priori information was used in all three tests, the long-wavelength structures of the reference models were fairly predicted, and satisfactory matches between “observed” and predicted data were achieved. The fair predictions of the reference models suggest that the final models estimated by our genetic-algorithm FWI, which we call macromodels, would be suitable inputs to gradient-based Rayleigh-wave FWI for further refinement. We also explored other issues related to the practical use of the method in different work and explored applications of the method to field data.

Compositional upscaling in fractured reservoirs during gas recycling

This paper presents a multidomain (MD) computational approach to multiscale (MS) compositional reservoir simulation amenable to parallel processing. The main objective is to minimize the upscaling process and preserve the pertinent flow characteristics of the fine-scale geological features of the reservoir. The preservation of the fine-scale reservoir characteristics yields higher-accuracy solution compared to an upscaled model. The computing algorithm involves breaking down the fine-grid (FG) computing mesh into an underlying coarse-grid (CG) subdomain. Unlike other MS techniques, the CG subdomain only serves as an assemblage of the FG cells and not as a computational CG overlaying the FG cells. The flow equations are then restructured into a global pressure equation which is solved on the FG using a large time step without solving for the compositions, saturations, etc. The convergence tolerance for the global pressure solution is soft (no iterations) to reduce computing time. The pressure solution is then used to calculate the molar flux of components only at the FG boundaries of the CG subdomain cells. Finally, each CG subdomain subject to its boundary conditions is solved independently using a smaller time step to compute cell pressures, compositions, and phase saturations. This procedure is amenable to parallel processing. This MS scheme was applied to time domain, in which the global pressure solution is obtained with a larger time step than in the local FG solution. This means that the molar flux boundary condition at the CG subdomain boundaries is not updated every time step, thus, buying additional computational advantage. This new technique was tested on several volatile oil and gas-condensate systems and excellent agreement was obtained with FG simulation with about 60% reduction in computing time using a single processor. Much greater time reduction is expected in a multi processing environment.

Seismic data are routinely used to infer in situ properties of earth materials on many scales, ranging from global studies to investigations of surficial geological formations. While inversion and imaging algorithms utilizing these data have improved steadily, there are remaining challenges that make detailed measurements of the properties of some geologic materials very difficult. For example, the determination of the concentration and orientation of fracture systems is prohibitively expensive to simulate on the fine grid and, thus, some type of coarse-grid simulations are needed. In this paper, we describe a new multiscale finite element algorithm for simulating seismic wave propagation in heterogeneous media. This method solves the wave equation on a coarse grid using multiscale basis functions and a global coupling mechanism to relate information between fine and coarse grids. Using a mixed formulation of the wave equation and staggered discontinuous basis functions, the proposed multiscale methods have the following properties. • The total wave energy is conserved. • Mass matrix is diagonal on a coarse grid and explicit energy-preserving time discretization does not require solving a linear system at each time step. • Multiscale basis functions can accurately capture the subgrid variations of the solution and the time stepping is performed on a coarse grid. We discuss various subgrid capturing mechanisms and present some preliminary numerical results.

Reservoir simulation models often require grids sufficiently resolved to capture the complexities of the geological structures present so that pressure and saturation profiles are deemed reasonable enough for decisions to be made with confidence. However, high resolution often implies long simulation run times, which is compounded in optimization problems. Reservoir optimization problems usually involve localized decisions, e.g., production control for a specific time-dependent objective function. With this tighter focus on optimization objectives, the reservoir simulator can often be replaced by an approximate proxy function, e.g., decline curves, neural networks and polynomial response surfaces. This accelerates the run time, allowing the optimizer to make more calls to the simulator, but at the cost of simplified physics and some consequent loss of accuracy. Recent work in numerical solutions of partial differential equations (PDE) by Druskin and Knizhnerman (2000) introduced a finite-volume optimal grid approach to improve solution accuracy of linear PDEs at targeted points of a computational domain. This approach was applied to the solution of electromagnetic and acoustic forward and inverse problems arising in oil exploration and medical imaging. We extend this work to non-linear PDEs by introducing a reservoir simulator that is used as its own proxy. We distinguish our optimal grid from an optimal upscaling in that it may not preserve local features of the model, but will preserve some aggregate objective value used in production optimization. This narrow focus allows greater coarseness than a properly upscaled grid by replacing the fine-grid model with a coarser grid whose coarsening is determined by an optimizer that preserves the accuracy of a predefined aggregate objective value. In this work we optimize injection and production rates to maximize the field oil production total (FOPT) over a fixed period of time. The proxy optimizer will seek the best coarsening that preserves this objective value over a range of injection/production rates for each well. The resulting coarse-grid model has a similar physical behavior in these output values to that of the fine-grid model. The grid coarsening problem is a non-linear discrete optimization problem, where the optimizer adjusts the grid coarsening while preserving the fine-grid response of interest. Although we focus on FOPT here, other reservoir objectives can also be handled.

Flow simulation of multimillion grid cell models with hundreds of wells and decades long production history can be extremely time-consuming. This often limits the applicability of assisted history matching techniques. A pragmatic solution to this problem is grid coarsening which is now embedded in many commercial reservoir simulators. Instead of up-scaling geological models in external packages, grid cells are automatically amalgamated within the simulator while preserving flux distribution and reducing the total number of active cells. The resulting speedup can be significant, often only at small loss of accuracy. Both characteristics are essential elements of any inversion technique in a multimillion grid cell environment. For water-flood history matching, we have utilized a commercial finite-volume simulator and the streamline-based generalized travel time inversion whereby water-cut behavior is matched by adjusting inter-well permeabilities. To apply the assisted history matching technique on high resolution models, a flux reconstruction method is devised which makes full use of the benefits of automatic grid coarsening. It approximates fluxes in the original geological grid by recalculating transmissibility at the fine scale and redistributing coarse scale fluxes accordingly. The inversion is thus conducted on the fine scale grid, while the forward simulation model uses the coarse grid. While solving for an inexpensive coarse model, the calibration takes into account fine grid resolution. Upscaling is not required and the degree of coarsening can be adjusted during the matching process. The proposed method was successfully tested on a supergiant carbonate oilfield with about hundred wells and large-scale water injection. The history match improved dramatically at relatively low numerical cost, which also allowed for investigating multiple sensitivities. The results were verified against non-coarsened model. The significant increase in efficiency makes this a potential method of choice for cases, where previously assisted history matching techniques could not be deployed due to excessive run times.