ROCK FORMATIONS FOR RESERVOIR SIMUUATION

BACKGROUND OF THE INVENTION

1. Field of the Invention

[0001] The present involves reservoir simulation, and in particular simulation methodologies which take into account both fluid flow and pressure distribution in the reservoir during production of fluids.

2. Description of the Related Art

[0002] In the oil and gas industries, massive amounts of data are required to be processed for computerized simulation, modeling, and analysis for exploration and production purposes. For example, the development of underground hydrocarbon reservoirs typically includes development and analysis of computer simulation models of the reservoir. A realistic simulation model of the reservoir, and the presence of its fluids, helps in forecasting the optimal future oil and gas recovery from hydrocarbon reservoirs. Oil and gas companies have come to depend on simulation models as an important tool to enhance the ability to exploit a petroleum reserve.

[0003] The underground hydrocarbon reservoirs are typically complex rock formations which contain both a petroleum fluid mixture and water. The reservoir fluid content usually exists in two or more fluid phases. The petroleum mixture in reservoir fluids is produced by wells drilled into and completed in these rock formations. Sometimes, fluids such as water or gases are also injected into these rock formations to improve the recovery of the petroleum fluids. [0004] Reservoir simulation belongs to the general domain of flow in porous media simulation. However, reservoir simulation normally involves multiple hydrocarbon components and multiple fluid phases in an underground geological formation which is under high pressure and temperature. The chemical phase behavior of these hydrocarbon fluids and the included groundwater has to be taken into account in these simulators.

[0005] The simulation models contain volumetric data which describe the specific geometry of the rock formations and the wells. The simulation models also contain reservoir properties data, such as the fluid and rock properties, as well as production and injection history pertaining to the specific reservoirs of the oil or gas field in question. The simulation models are formed by a simulator (known as a reservoir simulator) which is a suite of computer programs run on a data processing system.

[0006] The reservoir simulator which runs these models is a computer implemented numerical methodology, or coded algorithms and data constructs of an underlying mathematical model. The mathematical model which represents the physics of fluid movements in these hydrocarbon reservoirs is a system of nonlinear partial differential equations which describe the transient multiple-phase, multiple-component fluid flow. This model also represents material balance behaviors in these reservoirs induced by the production or injection of fluids, as well as the pressure- volume-temperature (PVT) relationships of the reservoir fluids.

[0007] A reservoir simulator simulates the multiphase multicomponent fluid flow and material balance in subterranean reservoirs and the included surrounding porous rock formations by subdividing the volume into contiguous cells, also known as grid blocks. In simulation models, the reservoir is thus organized into a number of individual cells. A cell or grid block is the basic finite volume where the underlying mathematical model is applied. The number of cells varies depends on the resolution needed for the simulation and the size of the reservoirs in question.

[0008] For a large reservoir, such as the type known in the industry as a giant reservoir, which may have multi-billion barrels of original oil-in-place (OOIP), the number of grid cells can be in the hundreds of millions to over a billion. This number of cells is required in order to have adequate resolution to represent flow dynamics, formation rock porosity and permeability heterogeneity, and many other geologic and depositional complexities within the reservoir. Simulation of this size reservoir can be termed giga-cell reservoir simulation. Figure 1 shows at 100 an example structured reservoir grid model of a subsurface reservoir for reservoir simulation.

[0009] The challenges in hydrocarbon reservoir simulation require the use of the latest technology to maximize recovery in a cost-effective manner. Reservoir simulators such as those provided under the trademark GigaPOWERS have been described in the literature. See, for example articles by Dogru, A.H. et al,“A Next-Generation Parallel Reservoir Simulator for Giant Reservoirs ,” SPE 119272, proceedings of the 2009 SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, Feb 2-4, 2009 and by Dogru, A. H; Fung, L. S; Middya, U.; Al-Shaalan, T. M.; Byer, T.; Hoy, H. Hahn, W.A.; Al-Zamel, N.; Pita, J.; Hemanthkumar, K.; Mezghani, M.; Al-Mana, A.; Tan, J.; Dreiman, T.; Fugl, A.; Al-Baiz, A.;“New Frontiers in Large Scale Reservoir Simulation ,” SPE 142297, Proceedings of the 2011 SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, Feb 21-23, 2011. GigaPOWERS reservoir simulation is capable of fine- scale grid simulation that exceeds a billion-cell barrier for post-processing while utilizing hundreds of GB footprint per scenario.

[0010] The total number of simulation runs for a company with a number of hydrocarbon reservoirs and appreciable reserves exceeds multiple tens of thousands per year, and one or more petabytes of high performance storage is required to host these data.

[0011] The transient solution of the multiphase multicomponent system involves the evolution of mass and energy conservation in a sequence of time steps from the initial condition of the reservoir. For each time step, the system of nonlinear discrete equations for each finite volume is linearized using what is known as the generalized Newton’s method.

[0012] In the industry, this is referred to as Newtonian iteration or nonlinear iteration. At each Newtonian iteration, a linear system of equations is constructed where the matrix, known as the Jacobian matrix, and the right-hand- side vector, known as the residuals are used to solve for the changes in the primary variables of the system.

[0013] Other factors about the reservoir further increased the complexity of giga-scale reservoir simulation. The subsurface rock in the formations is not homogenous in either permeability or porosity. For highly heterogeneous reservoirs, iterative methods may not converge in calculating reservoir pressure distribution. This introduces additional complexities to the simulation processing. [0014] Current industry practice for solving the linear system of equations is via a preconditioned iterative method. An example of preconditioning is described, for example in Saad, Yousef,“ Iterative Methods for Sparse Linear Systems,” January, 2000, Ch. 10, pp. 265 - 319. The computationally expensive preconditioner methods increase the processing time required and thus the computational cost or reservoir simulation proessing. Further, preconditioning methods do not guarantee convergence to an exact solution.

[0015] A preconditioner can generally provide a reasonable initial for the unknown pressure and saturation distribution within the reservoir. This initial estimate is then improved by iterations to the correct pressure and saturation distribution for the time step. The Gauss Seidel method is one method used as a preconditioner because of its simplicity and lower comparative computational cost per iteration. One such Gauss Seidel methodology is described, for example, in Olsen-Kettle, “ Numerical Solution of Partial Differential Equations” , Chapter 4.2, The University of Queensland. It has been found, however, that Gauss Seidel preconditioning often shows slow convergence and thus extended computation time, especially for reservoirs with highly heterogeneous permeabilities and porosities in the formations.

SUMMARY OF THE INVENTION

[0016] Briefly, the present invention provides a new and improved computer implemented method of reservoir simulation of fluid flow in a producing reservoir with improved convergence of determining pressure distribution within the reservoir to determine fluid flow in the reservoir. The reservoir has rock having heterogeneous properties of permeability and porosity. The reservoir further is organized into a three dimensional (3D) grid of reservoir cells. The reservoir simulation is performed for a sequence of time steps based on fluid pressures from wells in the reservoir. The formation layers contain rock having heterogeneous permeability and porosity properties.

[0017] The computer implemented method determines fluid flow from the three dimensional grid of reservoir cells to the wells and well production rates of the wells. The computer implemented method forms an initial measure of pressure distribution in individual ones of the grid cells of the reservoir for time steps of the reservoir simulation based on at least one of the heterogeneous rock properties and distances of the individual ones of the grid cells from wells in the reservoir. A pressure distribution in the individual ones of grid cells of the reservoir for a time step of the reservoir simulation is then based on the formed initial measure of pressure distribution in the individual ones of the grid cells of the reservoir.

[0018] An estimated fluid flow rate in the grid cells of the reservoir is then determined for a time step of the reservoir simulation based on the determined an estimated pressure distribution in the grid cells of the reservoir. Next, a determination is made whether the determined fluid flow rate in the grid cells of the reservoir has converged. If so, the reservoir simulation of fluid flow in the grid cells of the reservoir is ended. If not, processing returns and repeats the steps of determining an estimated fluid flow rate in the grid cells of the reservoir for a time step of the reservoir simulation based on the determined an estimated pressure distribution in the grid cells of the reservoir, and determining whether the determined fluid flow rate in the grid cells of the reservoir has converged. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 is an isometric view of a computerized model of a subsurface reservoir structured grid for reservoir simulation.

[0020] Figure 2 is an isometric view of a one-dimensional reservoir model.

[0021] Figures 3 is a schematic diagram of a computer network for determining pressure distribution in heterogeneous rock formations for reservoir simulation according to the present invention.

[0022] Figure 4 is a schematic diagram of an application server or computer node of the computer network of Figure 3.

[0023] Figure 5 is a functional block diagram or flow chart of workflow for determining pressure distribution in heterogeneous rock formations for reservoir simulation of the present invention.

[0024] Figure 6 is a display of comparative results of pressure distributions calculated by analytical solution and by numerical solution for the reservoir model of Figure 2 with uniform or homogenous permeability and porosity for a one dimensional 11 grid block system.

[0025] Figure 7 is a display of comparative results of pressure distributions calculated by analytical solution and by numerical solution for the reservoir model of Figure 2 with heterogeneous or varying permeability and porosity for a one dimensional 11 grid block system.

[0026] Figure 8 is a display of comparative results of convergence error in pressure distribution determination according to the present invention against processing results according to the prior art for a two dimensional flat homogenous permeability reservoir model for a 5 x 5 x 1 areal grid system.

[0027] Figure 9 is a display of comparative results of convergence error in pressure distribution determination according to the present invention against processing results according to the prior art for a two dimensional flat heterogeneous permeability reservoir model for a 5 x 5 x 1 grid system.

[0028] Figure 10 is another display of comparative results of convergence error in pressure distribution determination according to the present invention against processing results according to the prior art for a two dimensional flat homogeneous permeability and porosity reservoir model for a million grid block 1,000 x 1,000 x 1 grid system.

[0029] Figure 11 is a display of comparative results of convergence error in pressure distribution determination according to the present invention against processing results according to the prior art for a large three dimensional heterogeneous permeability and porosity reservoir model a million grid block 1,000 x 1,000 x 1 grid system.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0030] In numerical reservoir simulation, based on the formulation, a pressure equation is formed either explicitly or included in the main body of the flow equations. Fluid velocities, and hence the flow of oil, water and gas is computed from the pressure distribution, hence the pressure fields must be computed properly. Often simulators face convergence problems due to difficulties in finding a proper pressure solution at a given time step. Therefore, correctly computed pressure solutions are critical for simulators to proceed in a time domain and determine reservoir performance with reasonable computational cost.

[0031] The present invention provides an improved methodology for determining pressure distribution within the reservoir to determine fluid flow in the reservoir, with improved convergence during the computer processing. This saves computational costs and solves a technological problem in the technological field of reservoir simulation. The present invention determines an accurate initial estimated measure for the unknown pressure distribution among the individual cells of a reservoir where it is known that the rock structure in the reservoir formations is heterogeneous. The present invention avoids the need to rely upon the time-consuming preconditioner condition of non-convergence likely to occur in the pressure distribution determinations using techniques such as the Gauss Seidel preconditioning for that purpose.

[0032] Underground fluid reservoirs are highly heterogeneous in rock properties. These heterogeneities are reflected in the coefficients of the partial differential equation describing the fluid flow and pressure distribution in the reservoir.

Nomenclature

[0033] Set forth below for ease of reference and understanding is a listing of the nomenclature used in the Equations which represent the various parameters, measurements and units of measure according to the present invention. m = fluid viscosity, centipoise (cP) p= fluid density, grams/cubic centimeter (g/cm ³⁾ p=pressure, lbs/in ² , (psi)

F= fluid potential, barrels/day/psi (STB/d/psi) q= well production or injection rate, barrels/day (b/d)

t= time

K = permeability tensor

k = absolute permeability of the reservoir rock, millidarcy (mD) k _xx = Value of x directional rock permeability at any x location k _xy = Value of x directional rock permeability at x and y location k _xz = Value of x directional rock permeability at x and z location Elements of permeability tensor K at y and z locations interpreted in a similar manner x, y, z = space coordinate system

i, j, k = grid block index in x, y and z directions

cp= porosity

f _re = reference porosity

C= compressibility, per psi (l/psi)

Cp= Pore Compressibility

CL= Liquid Compressibility

Co = oil compressibility

Pref = reference pressure

T= Transmissibility between grid blocks

A= Cross Sectional Area perpendicular to flow Ax, Ay, Dz = Grid Block Size in x, y, z directions

V _p = Pore Volume

G = Green’s Function

_< 5(a— b) = Dirac'sDelta function (= 1 if a = b, otherwise = 0) [0034] The known pressure equation defining pressure conditions in a body such as a subsurface reservoir for a single phase or multiple fluid phases with pressure dependent fluid density can be written as below:

[0035] In Equation ( 1 ), K is a tensor of the reservoir parameter or attribute of rock permeability, defined as below:

[0036] For the rock permeability tensor K of Equation (1), heterogeneous rock permeability is taken into account by a three-dimensional matrix of measures of the permeability according to Equation (2). In Equation (2), k is the rock permeability, subscript x indicates the x direction, subscript y indicates y direction and subscript z indicates the z direction (vertical) in a three dimensional space.

[0037] In Equation (1), //describes fluid viscosity; //is the fluid density representing density of a single fluid such as oil, water or gas or an average density of the fluid mixture in the pore; F represents rock porosity; q is the volumetric mass injection or production rate (well); and t represents the time. Fluid potential F in Equation (1) is represented by Equation (3): f = P - (

yc P ^Z 3)

[0038] In Equation (3), p is the fluid pressure, g is the gravitational force and g _c is the gravitational conversion factor, p is the average density between two spatial locations and is the vertical distance from a reference point.

[0039] In Equation (1), Eis the divergence operator and is defined by Equation (4) as follows:

[0040] In Equation (4), i, j and k are the unit vectors in x, y, and z directions. [0041] A number of physical properties or measures relating to a subsurface reservoir which are expressed in the discussion above are not constants, but rather dependent on pressure, temperature or other reservoir conditions.

[0042] For example, in Equation (1), F is the rock porosity and is defined by Equation (5) as follows:

[0043] In Equation (5), F _tb† is the porosity measured at a reference pressure, p _ref and c _P is the pore compressibility measured in laboratory conditions, at constant temperature.

[0044] Similarly, fluid density p is defined for fluids oil and water above the fluid mixtures bubble point pressure by Equation (6) as follows:

[0045] In Equation (6), c _L is the liquid (oil or water) compressibility, and subscript L stands for liquid at constant temperature.

[0046] For gas reservoirs fluid density is defined by Equation (7) as follows:

[0047] In Equation (7), M is the molecular weight of the gas, p is the gas pressure, R stands for universal gas constant, T is the absolute temperature and Z(p,T ) is the gas deviation factor which can be obtained from laboratory or reference charts or by using equation of state calculations for a given gas composition.

[0048] The fluid potential F of Equation (1) defining pressure conditions in a subsurface reservoir is thus pressure dependent. To take this into account, Equation (3) is substituted into Equation (1), resulting in the following Equations (8a) and (8b):

[0049] In Equations (8a) and (8b), p in the second term on the left hand side representing the flow due to gravity forces is defined as:

P = & i/c P (8c)

This term is important for reservoirs, taking into account elevation differences and three dimensional effects in general.

Linearization the Pressure Equation (Equation (1))

[0050] Using the pressure depending property formulas of Equations (5) through (7), Equation (1) becomes Equation (9) as follows:

[0051] For slightly compressible oil and water, Equation (9) becomes Equation (10) as follows:

[0052] In Equation (10), ct is the total pore compressibility including the rock (pore) and the fluids, and q is the volumetric production or injection rate. Equation (10) represents a“volume balance” relationship of pressure conditions in a subsurface reservoir as compared to Equation (1) which represents a“mass balance” relationship of subsurface reservoir pressure conditions.

[0053] For flat reservoirs, Equation (1) becomes Equation (11) as follows:

which is a linear partial differential equation. The interrelation between the numerous variables representing pressure conditions in a subsurface reservoir with pressure dependent fluid density are thus in a form in which it is possible to determine the pressure distribution in the reservoir by analytical solution for a given set of initial and boundary conditions.

[0054] For a diagonal tensor of the type expressed in Equation (2) , and assuming only property changes in x, y, and z directions, a partial differential Equation (12) can be written representing pressure conditions in a body such as a subsurface reservoir according to Equation (11) as follows:

9 , 9 p 9 , dp 9 , 9 p , dp

— k _Y — +— k _v — +— k ₇ - \- q = c _† — (12) dx ^x dx dy dy dz ² dz “ ^c dt C _t — Cp + c ₀

Discretization

[0055] An initial condition and boundary conditions are needed to determine the pressure distribution in the reservoir at any time t and location (x, y, and z) according to Equation (12). In order to solve Equation (12) for a given reservoir, the reservoir is divided into grid blocks, which can be rectangular boxes, and into time steps to cover the past reservoir history and future performance in time.

[0056] The following example and explanation of discretization assume a two dimensional flat reservoir. However, the methodology and principles applied are of general applicability and can be easily extended into a three dimensional system.

[0057] Assuming that reservoir such as that shown at 100 in Figure 1 is divided into N _x number of blocks in x direction and N _y number of block in y direction and i represents the grid block number in x direction, z=l, 2, ... N _x and j =1, 2, ... N _y in y direction and 1 grid block in the z direction, Nz=l.

[0058] For any grid block located at grid coordinates (z, j ) in the reservoir, finite volume discretization yields the following volume balance relationship expressed in Equation (13):

where T represents the interblock transmissibility between two neighboring grid blocks and defined as below.

[0059] In Equation (13), the superscript n represents the preceding or old time step, while the new or current time step is represented by the superscript n+1.

[0060] For example, between grid blocks (z-l, j) and (i ,j) interblock transmissibility T can be expressed according to Equation (14) as follows: where A is the surface area perpendicular to flow between the grid blocks and Ax is the length of the grid block in x direction. The term ki- _± / _2,j is ^a harmonically averaged interface permeability between grid blocks (z-l, j) and (z, j).

[0061] The pore volume V _P of the grid block is defined by Equation 15 as follows:

V _{Pi .} = Ax _u j Ay ₀ j F _{ί }} (15)

[0062] In Equation (13) pressures for each grid block p ⁿ i,j are known (starting from an initial reservoir condition) at the start of production. Reservoir simulation is performed to compute the unknown pressures for the new time step (n+l), namely p ⁿ⁺¹ i,j for all grid blocks z=l, 2, . . N _x j=l, 2, ..N _y .

[0063] For this purpose, Equation (13) is in matrix form for N _x x N _y grid blocks, using the proper boundary conditions (usually no flow from regions around the reservoir) and with given well rates qi . Equation (13) yields an N _x x N _y linear system with a sparse coefficient matrix for the unknown pressures p ⁿ⁺¹ i, j.

[0064] The linear system of Equation (13) is solved for the unknowns with several known types of iterative or direct methods (for smaller systems). Current industry practice for solving the linear system of equations is via a preconditioned iterative method.

[0065] However, as has been noted, for highly heterogeneous media in a reservoir, iterative methods may have difficulty in obtaining convergence in calculating reservoir pressure distribution. To prevent this, computationally expensive preconditioner methods have previously been used. So far as is known, however, these methods increase computational costs, and do not guarantee convergence of the pressure distribution determination for the reservoir cells to exact solution.

The Present Invention

[0066] The present invention improves functioning of the data processing system S in to solve this linear system of equations faster and determine pressure distribution for use in reservoir simulation. Many iterative methods require large cost to precondition the system which becomes computationally expensive. Several iterative methods will not converge to the solution for a given time step At, therefore, ending up reducing the time step number to a small number. Small time steps sizes can increase the total number of time steps to complete the simulation. However, this can be very expensive in terms of unnecessary computer usage and processing time and costs.

[0067] On the other hand other iterative methods require a close initial estimate of reservoir pressure for the determination of the unknown fluid flow. The present invention provides a pressure solution closer to the unknown solution, so that simple iterative methods can determine unknowns with less computational cost.

[0068] A single phase pressure equation is selected for demonstration purposes. The present invention provides a new approximate initial pressure estimate for the linearized pressure equation with highly heterogeneous coefficients. Numerical examples are presented below to demonstrate that the present invention can solve the pressure equation effectively by reducing the numerical error per iteration an order to several orders of magnitude compared to a basic iterative solver.

Solution Structure for Initial Estimate of Reservoir Pressure

[0069] Reservoir simulation models usually implement no flow boundary conditions at the boundaries of the reservoir. Wells provide an inner boundary condition. The present invention provides an improved methodology for an initial estimate of reservoir pressure which covers such types of boundary conditions.

[0070] An illustrative example for fluid production and related rock properties of Equation (12) in conjunction with Equation (13) is presented. Consider for a simplified example (Figure 2) a one dimensional linear porous media 110 (a simplified reservoir) with constant permeability, constant porosity, flow cross sectional area, and constant fluid viscosity. In Figure 2, the porous media 110 has a length L along a linear or x dimension, and a well 112 at location at X* located between x= 0 and x = L. Both an east boundary 114 at x = 0 and a west boundary 116 at x = L are considered sealed (no flow).

[0071] By discretizing the time derivative and leaving the space derivatives at the new time level (n+l), Equation (12) becomes:

[0072] Removing the source term results in:

assuming p ⁿ = 0, Equation (16) becomes:

[0073] Since no time derivative is involved, Equation (17) can be written in an ordinary differential equation format as Equation (18):

[0074] Defining a =

, Equation (18) results in Equation (19) as follows:

[0075] Equation (18) needs two boundary conditions to have a unique analytical solution. Considering the solution to the left of the source x ^* ,

Left and right Boundary Conditions

[0076] At the left and right of the source x ^* the well production or injection flow q is thus:

[0077] Defining normalized variables results in Equations (22a) and (22b) as follows:

x = ^ , x < x ^* (22a)

[0078] Approximate analytical solutions for Equation (18) for the porous media 110 of Figure 2 and the above stated boundary conditions are as follows:

Left solution, x < x ^* : Right solution, x > x ^* :

[0079] In Equations (23) and (24):

and

and A is the cross sectional area to flow in Equations (14), (25) and (26).

Multiple Dimensions

[0080] A more complex reservoir model than that of Figure 2 can now be considered. For a given time step, the unknown pressure p(x, y, z) for a porous medium of multiple dimensions (x, y, z) can be assumed to be the product of the solution in x, y and z directions. Since the pressure equation is linear, a product assumption for the solution can be used:

[0081] For an x direction pressure solution p _x (x ) is presented in Equations (23) and (24). A solution in the y direction, p _y (x) is constructed similarly to Equations (23) and (24) except variables y, y ^* and y are used in place of the x variables. A pressure solution p constructed similarly.

[0082] For purposes of simplification, the following discussion considers only x, y directions with no gravity effects. A complete three dimensional solution can be obtained by using Equation (28).

Multiple wells and Superposition [0083] The analytical solutions given by Equations (23) and (24) for a single well located at grid block k with a production or injection rate q by Equation (27a) can be abbreviated as follows: O, <7 _f c) = ¾¾( ) (27a)

[0084] For multiple wells, the solutions obtained due to a single source located at grid block k can simply be superposed as expressed in Equation (27b) below:

Calibration

[0085] In order to match an analytical solution to a numerical solution, a calibration has been found helpful. First, an average pressure _m¾e is determined according to Equation (28a) based on material balance conditions, as follows:

[0086] In Equation (28a), Vp,t is the total pore volume of the reservoir.

[0087] Average reservoir pressure p _AS for the analytical solution is determined according to Equation (28b), as follows:

[0088] The corrected analytical solution for pressure for numerical solution in successive ones of n iterations is set forth in Equation (28c) as follows:

n+l,corr n+1 Pmbe

Pi = Pi (28c)

PAS

i = 1, 2, ...n

Verification

[0089] For verification, a one simplified one dimensional reservoir with uniform rock and fluid properties such as the simplified model 110 (Figure 1) is used. Equation (18) is discretized in space and time in a manner similar to Equation (19). This results in a linear system of equations with a tridiagonal coefficient matrix for the unknown pressures at each grid block for the new time step (n+l) as shown below in Equation (29):

In Equation (29) for any grid block i, Equations (30) and (31) below define the west neighbor transmissibility.

[0090] The east transmissibility T _{e i} is defined similarly for the i + 1/2 cell face.

Boundary Conditions

[0091] West and east end boundaries are assumed to exhibit no fluid flow. Therefore,

T _w,i = 0, T _e,i = o (32)

[0092] The well production or injection rates q _l q ₂ —q _n f° ^{r a} grid cell i are set forth in the right side of matrix Equation (29) above.

[0093] Equation (30) is solved for the first time step pressure distribution for the one dimensional reservoir by using Gaussian elimination or another suitable linear solver method.

[0094] Set forth below are an example set of parameter values for such a pressure distribution in a reservoir model of comparable dimensions to the model 110 of Figure 2 organized into eleven grid blocks. Also included below are units of measure for these parameters which are applicable unless otherwise indicated for the models and examples given in the present application. The pressure distribution in the reservoir model is determined by both analytical solution and numerical solution to verify conformity between the analytical and the reservation simulation runs:

[0095] Table 1 shows the results of pressure distribution calculated by the analytical solution and numerical solution described above. The numerical solution was obtained by solving Equation (30) in the data processing system S.

Table 1

[0096] Table 1 verifies that the proposed analytical solution yields a very close solution to the numerical solution. Figure 6 is a display of comparative results of pressure distributions for a succession of simulation runs for time steps Dt of one day for a homogeneous permeability reservoir as calculated by analytical solution and by numerical solution for the reservoir model of Figure 6 having the example set of parameter values identified above. The comparative results of the pressure distribution determination runs for the analytical solution and the numerical solution are so close as to overlap in the plot of Figure 6.

Initial Pressure Determination for Heterogeneous Reservoirs

[0097] Reservoirs in nature are heterogeneous in rock properties such as permeability and porosity in the nature. For actual heterogeneous reservoirs, there is no analytical solution. A solution can only be obtained by numerical methods such as the one described by Equation (30). The present invention provides an initial pressure determination for heterogeneous reservoir rock properties of permeability and porosity using the same structure offered by Equations (23) and (24). [0098] In determining the pressure at a location x with according to Equations (23) and (24), due to a source at x*, a harmonically averaged permeability k _{i k} between x and x* is used, and a pore volume averaged porosity <p _{i k} between the same two points is also used.

[0099] Assuming that well location x* corresponds to grid block k and location x has a grid block index i, the harmonically averaged permeability k _{i k} can be expressed for each individual grid block i as: åU* ^x i

ki,k — . DC :

y 7 (33)

/=l fey

[00100] Similarly, for porosity, a pore volume averaged porosity fi _k can be expressed for each individual grid block i with constant cross sectional area A is as follows:

Example for Heterogeneous Reservoir

[00101] The same example grid blocks used for the verification described above were adjusted to have heterogeneous permeability and porosity. In the model, it is assumed that the reservoir has one injector at the left end and one producer at the grid block 11 on the right, and with constant rates of 500 b/d injection and 300 b/d production. The reservoir permeability and porosity distributions among the eleven grid blocks are given in Table 2.

Table 2

[00102] Pressure distributions were calculated by both the analytical solution and the numerical solution. Results of the pressure distribution runs for both the analytical solution and the numerical solution are plotted in Figure 7. As seen, an analytical solution captures the numerical solution trend and yields a highly reasonable approximation.

Data Processing System

[00103] The processing and workflow described above are suitable for deployment on a variety of high performance computer (or HPC) cluster computer/processor node computer systems. These are typically rack mounted hardware with several compute nodes which contains multiple CPUs with multi-core architecture. Nodes are typically interconnected with conventional low latency high bandwidth networks, switches, and routers.

[00104] It should be understood that the present invention may be performed on a mainframe computer of any conventional type of suitable processing capacity such as those available from International Business Machines (IBM) of Armonk, N.Y. or other source.

[00105] The typical HPC environment for use with this simulation system is today’s multi-node, multi-CPU, multi-core compute clusters. An example such cluster is illustrated at C in the data processing system S of Figure 3. The cluster C is formed of a plurality of computer nodes 150 (Figures 3 and 4) which are provided with data in parallel as indicated by and arrow 152 by a router server or servers 154. If desired, several such router servers may be used for this purpose. Original simulation or input data of the types described above is stored in a suitable number of data storage/file servers 156. The data storage/file servers 156 store operating instructions, control information and database records for use by the data processing system S. [00106] The router servers 154 operate under control of computer code 155 stored in memory and transfer input simulation data in parallel from the storage servers 156, as well as simulation processing results as indicated by an arrow 158, to and from the computer nodes 150 of the cluster C.

[00107] The program codes 155 according to the present invention is in the form of non- transitory computer operable instructions causing the server or servers 154 to index, order, process and transfer the reservoir simulation results. Typically, the data processing system S includes a set of work stations 157 of suitable, conventional type which are connected to the system by a network 159.

[00108] The computer nodes 150 of the cluster C include a plurality of processors or cores 160 of the type illustrated in Figure 6 operating in parallel under instructions of a computer code or program product 162 stored in stored in memory 164 of the computer nodes 150. The program code 162 according to the present invention is in the form of non-transitory computer operable instructions causing the data processors 160 to perform reservoir simulation in reservoir simulation modules 166.

[00109] It should be noted that program codes 155 and 162 may be in the form of microcode, programs, routines, or symbolic computer operable languages that provide a specific set of ordered operations that control the functioning of the data processing system D and direct its operation. The instructions of program codes 155 and 162 may be stored in memory of the servers 154 or processor nodes 150, or on computer diskette, magnetic tape, conventional hard disk drive, electronic read-only memory, optical storage device, or other appropriate data storage device having a non-transitory computer usable medium stored thereon. Program code 160 may also be contained on a data storage device such as server 156 as a computer readable medium, as shown.

The Present Invention

[00110] The present invention provides improved determination of pressure distribution and pressure values over the three dimensions of heterogeneous rock formations of a reservoir such as the reservoir 100 (Figure 1). The present invention further permits prompt initial pressure measures for pressure determination through the reservoir model during reservoir simulation for the subsurface reservoir 100. The present invention reduces convergence problems with existing reservoir simulation. The initial pressure distribution which is determined represents a correct pressure profile within the reservoir. The determined initial pressure in the reservoir is provided as input pressure measures for reservoir simulation by iterative methods.

[00111] The estimate of initial pressure is proven close to the actual solution based on analytical solutions in simplified reservoir models. The determined initial pressure estimate enables the reservoir simulator to converge to a desired solution within specified accuracy limits. Determination of the initial pressure distribution according to the present invention also converges with further computer processing iterations then other iterative determination methods. An example iterative determination method is Gauss Seidel methodology, which is described at: Olsen-Kettle, ‘ Numerical Solution of Partial Differential Equations” , previously mentioned. Gauss Seidel methodology was used as iterative method of determination of the initial pressure distribution for convergence time benchmark or comparison purposes.

[00112] Figure 5 is a schematic diagram of a workflow F outlining a preferred embodiment of a sequence according to the present invention for reservoir simulation based on improved determination of pressure distribution over the three dimensions of heterogeneous rock formations of a reservoir.

[00113] During a step 200, simulation begins by reading the reservoir and production data. Reservoir and production data read in during step 200 are of the types discussed above. The reservoir simulator is also initialized during step 200, setting the simulation day and the time step to zero.

[00114] Reservoir data includes reservoir geometric information-its size, extent (length) in x, y, and z directions, reservoir properties such as permeability distribution, porosity distribution, thickness of layers, relative permeability data, capillary pressure data, fluid property data such as fluid density tables, formation volume factor tables, viscosity tables, location of wells, and location of well perforations within the reservoir.

[00115] Production data includes oil, water, and gas production rate measured or specified for the wells defined in the previous step. Production data also includes minimum bottomhole pressure for each well.

[00116] During step 202, pressure distribution for the cells of the reservoir model are determined solution using the same reservoir cell matrix structure offered by Equations (23) and (24). For the initial iteration of step 202 or time step n = 1 of the reservoir simulation, the pressure change is set to zero. Improved determination of pressures according to the present invention in heterogeneous rock formations permits acceptable initial pressure measures for pressure determination during reservoir simulation and also reduces convergence problems with existing reservoir simulation.

[00117] During step 204, the determined initial pressures for the cells of the reservoir model resulting from step 202 are used to determine final pressure distribution in the reservoir cells organized in multiple dimensions as shown in Equation (28) to calculate the estimate for the pressure distribution for the next fluid flow reservoir simulation time step (n+l). The iterative solver methodology for is preferably of the Gauss Seidel type described above, although it should be understood that other iterative solvers may be used. Processing by the iterative solver during continues until convergence is achieved. The pressures so determined then are provided as part of the reservoir simulation during step 206 to determine fluid flow in the grid cells of the reservoir. The pressure distribution in the reservoir cells in the multiple dimensions determined during step 204 is then available for performance of step 206.

[00118] During step 206, the pressure estimate determined during step 204 is used in fluid flow reservoir simulation based on improved determination of pressures over the three dimensions of heterogeneous rock formations of a reservoir. The fluid flow reservoir simulation is determined in connection with an iterative solver of the conventional type to perform multiple iterations and determine pressures by numerical solution to be used as pressures at time step (n+l).

[00119] As will be set forth, the present invention provides improved determination of pressures in heterogeneous rock formations. The present invention permits acceptable initial pressure measures for pressure determination during reservoir simulation, and reduces convergence problems with existing reservoir simulation.

Simplified Reservoir Models Initial Pressure Distribution Determination Examples

[00120] Comparative example cases of initial pressure determination for simplified reservoir models are presented: 1) a homogeneous reservoir; and 2) a heterogeneous reservoir.

[00121] For the two simplified models, a two dimensional flat square reservoir with the following properties was used as the reservoir model. The following chart sets forth reservoir parameters common to both examples.

1. Homogenous Reservoir Example

[00122] For the homogenous reservoir example, the formation rock was assigned uniform heterogeneous properties of permeability = 100 mD, and porosity = 0.20. The time step for the simulation, Dt = 15 days.

[00123] Figure 8 is a display of comparative results of convergence error in pressure distribution determination according to the present invention against processing results according to the prior art for the selected two dimensional flat homogenous permeability reservoir model as defined by properties stated above. Figure 8 shows an error in barrels per day averaged over the 25 grid blocks of the reservoir model versus number of fixed point Gauss Seidel iterations. Figure 8 shows at 220 results of convergence as a function of iteration count a standard Gauss Seidel solution with zero pressure change as an initial estimate of the initial pressure gradient for the reservoir model. Figure 8 shows at 222 results of convergence as a function of iteration count a standard Gauss Seidel solution with a pressure distribution initial estimate formed according to the present invention.

[00124] As can be seen from Figure 8, the results shown at 222, pressure determination according to the present invention is superior to standard methods with a zero initial pressure estimate. Residual error is reduced to 0.01 b/d in 10 iterations as shown at 222a and 0.00001 b/d in 20 iterations as indicated at 222b. In contrast, standard processing method a shown at 220a requires 210 iterations to reduce the convergence error to 0.1 b/d, which is ten times the residual error from the present invention, and requiring more than twenty times the number of iterations than the present invention.

[00125] The convergence rate is thus much more rapid and the number of computer processing iterations is greatly reduced with the new present invention. It is also to be noted that conventional methods with a pressure estimate of zero yields much flatter curve as shown at 220 for relatively large time step sizes, indicating the residual is changing relatively little for a considerable number of iterations.

2. Heterogeneous Reservoir Pressure Example

[00126] A highly heterogeneous reservoir model was set up using the same grid dimensions for the same reservoir presented in the previous example. Permeability was varied randomly from 1 mD to 4,000 mD. Porosity was varied from ½ percent to 20 percent. The same time step size of 15 days was used.

[00127] Results are displayed in Figure 9. As shown at 230, the standard Gauss Seidel with zero as an initial estimate for the reservoir model yields during a continuing series of iterations an almost flat curve, indicating extreme difficulty in reducing the residual error. This is a typical performance of a non-preconditioned iterative method. On the other hand, as indicated 232, the present invention reduces the residual error systematically. The present invention as indicated at 232a reduces residual error for the heterogeneous model to less than ½ b/d in 20 iterations while as indicated at 230a the standard Gauss method has not reduced the error to 1 b/d in 100 iterations. Further iterations with the present invention converge to an exact solution at a given accuracy.

[00128] It is to be noted that as indicated in the legend of Figure 9, the present invention accommodated an increase in pressure for the heterogeneous model of as much as 764 psi for a time step. This occurred with a time step pressure set constant of zero while determining a pressure distribution within the reservoir at 15 days (Dt=l5 days) based on injection at the center well as shown in Figure 9. The resulting pressure distribution within the reservoir-in these conditions was determined that to have increased from 0 to 764 psi. This is considered a significant and unusually large change in pressure in the context of reservoir simulation. Typical reservoir simulators, due to convergence problems, cannot take large time steps of sizes which would result from a pressure change of this magnitude. Rather, to avoid convergence problems, the reservoir simulators limit the time step size so that the maximum pressure change can remain within a much smaller pressure change, typically in the range of 100 psi. However, the reduced time step size increased the number of time steps and thus the computation time and cost of reservoir simulation.

Highly Heterogenous Reservoir Example: Estimation of Pressure Distribution

[00129] For a highly heterogeneous reservoir model, a two dimensional flat square reservoir with the following properties was used. The reservoir model lateral dimensions were 3,000 ft. by 3,000 ft., and thickness of 1 ft. The number of grid cells were selected to be N _x = 100, and N _y =l00, and grid cell dimensions were: D _x = D _y = 30 ft, and D _z = 1 ft. A well was located the center, injecting 50 b/d. Permeability was varied randomly among grid cells, varying between 5.e-4 mD to 100 mD, and porosity varied randomly between 2e-5 to 0.20. The time step Dt = 1 day.

[00130] Table 3 below shows the performance of Gauss Seidel method with different Initial pressure estimates.

Table 3

[00131] An analytical solution was used to generate the initial pressure estimate and the results are indicated in Table 4.

Table 4

[00132] As seen from the comparison of Table 3 and 4, Gauss Seidel iterations are sensitive to the pressure estimate assigned at the start of pressure determination, while the results obtained according to the present invention are closer to the actual solution. If the starting estimate is not properly chosen, the Gauss Seidel pressure determination persists closely to the initially selected pressure over continuing iterations and increased computer processing times.

[00133] As will be set forth, the present invention requires a much smaller number of iterations to find a solution. Gauss Seidel method can require excessive number of iterations if the initial estimate (usually unknown) is selected to be widely different from the actual solution, as Tables 3 and 4 illustrate. For example, as shown in Table 3 if the initial estimate is 2,000 psi it would take 6,600 iterations to reach the solution, while the present invention requires 8 times less the number of iterations.

Large Reservoir Example (Million Cell) Initial Pressure Determination

[00134] For a much larger reservoir model, a two dimensional flat square reservoir with the following dimensional properties was used as the reservoir model was used with lateral dimensions of 30,000 ft. by 30,000 ft., and with a thickness of 10 ft. The number of grid cells were selected = 1,000; N _y = 1,000, and Nz = 1. Grid cell dimensions were: D _x = D _y = 30 ft, and D _z = 1 ft. The total number of grid cells was thus one million.

[00135] The following formation data were assigned: pore compressibility Ct= 5e-6 l/psi; oil viscosity = 1 cP. Flow rate was set at 1,000 b/d injection at the center. Initial reservoir was assigned pressure = 0 psi. The time step Dt was established at 10 days.

Homogeneous Large Reservoir [00136] As an example the larger reservoir model was assigned homogeneous permeability, K _x = K _y = 100 mD for each grid cell, and homogeneous porosity = 0.20. The other parameters were the same as in the previous examples.

[00137] A reservoir simulation was run by using the Gauss Seidel method again as a benchmark to determine grid cell pressure distribution at the end of the first time step (10 days). The processing of the present invention according to Figure 10 was run.

[00138] Figure 10 is a plot of comparative results of convergence error in pressure distribution determination according to the present invention against processing results according to the Gauss Seidel prior art benchmark for a two dimensional million cell homogenous permeability reservoir model. Figure 10 illustrates the error reduction by the iterations. Both methods converged to the same pressure solution with a determined well bottomhole pressure of 707 psi and a residual error of l.e-7.

[00139] Processing according to the present invention formed an initial estimate of the well bottomhole pressure as 221 psi and took 1,993 iterations to converge on the 707 psi well bottomhole pressure. The conventional Gauss Seidel benchmark with zero change as an initial estimate took 6,030 iterations to achieve the 707 psi well bottomhole pressure solution, or more than three times the number of processing iterations for one time step.

Heterogeneous Large Reservoir

[00140] The same reservoir dimensions and parameters as for the homogenous larger reservoir model described above were used for this example. However, random permeability and porosity distribution were assigned. Minimum permeability for the random distribution was 0.0001 mD, maximum permeability was 100 mD; and average permeability was 50 mD. Minimum porosity of = le-7 was established; maximum porosity was = 0.20, and average porosity was = 0.10.

[00141] The present invention and the Gauss Seidel benchmark method were again run to calculate reservoir pressure distribution at the end of first time step and then compared. Performance of both methods is presented in Figure 11. For the Gauss Seidel benchmark an initial estimate of well bottomhole pressure was 0 psi. Initial pressure determined according to the present invention was estimated at 1863 psi. As shown at 260, the Gauss Seidel benchmark converged to a well bottomhole pressure of 1764 psi with a residual error of l.e-7. As shown at 262, the present invention converged to the same well bottomhole pressure of 1764 psi and residual error of l.e-7 after significantly fewer iterations.

[00142] However, the Gauss Seidel benchmark required 5,185 to achieve convergence for the time step, while the present invention obtained the same results based on the determined initial pressure in 3,125 iterations, about 40% less processing time.

Reservoir Fluid Flow Determination

[00143] Matrix Equation (35) below can be used according to the present invention for consecutive time steps to obtain pressure values during step 204. As noted a constant pressure of zero is assumed for an initial step, which is identical to defining existing pressure as the difference from a constant initial pressure.

[00144] Assuming the right hand side vector of Equation (35) to be new source (well) terms, an analytical solution can be constructed by the superposition of the analytical solutions. Grouping the right hand side vector elements as:

[00145] The superposition can be expressed in the following manner: [00146] Without improved determination of pressure distribution over the three dimensions of heterogeneous rock formations of a reservoir according to the present invention, the superposition presented in Equation (37) can be expensive due to the large number n of terms involved in the processing. This is especially for large scale simulation operations which almost certainly prove to be computationally expensive.

Reducing Computational Cost for Large Systems

[00147] Since the computational cost involved is in the evaluation of the analytical solution (G/c(x)) involving exponential terms contained in the pressure relationships expressed in Equations (23) and (24), an approximate radial solution to the pressure equation is presented.

[00148] Let x*, y* be the x, y coordinate of a source (well) in the (x, y) grid. The radial distance between a point (x, y) and (x*,y*) can be defined by:

[00149] Finding an average area of an equivalent radial system to L _x L _y Cartesian system can be expressed as Equation (39):

L _X Ly

r _e (39)

[00150] The basic solution becomes one expressed in Equation (4):

where

and r ^* is the radial coordinate of the source, and k _r is the radial permeability calculated from k _x and k _y .

Approximation to Exponential Function

[00151] The exponential expressions in Equations (23), (24) and (40) can be approximated by Equation (42):

e ^x = l + x (42) [00152] Depending on the value of x, Equation (42) should further reduce the computational cost for pressure estimate for large scale applications.

Parallel Computing

[00153] Since the analytical solution is expressed in closed form expressions, no dependency in the solution exists. This provides a significant advantage over the conventional solution methods which deals with dependency, or reclusiveness in the mathematical expressions. Therefore, the method presented here is highly parallelizable.

[00154] The present invention provides reservoir simulation based on improved determination of pressure distribution over the three dimensions of heterogeneous rock formations of a reservoir. It is adapted for and implemented in multiple dimensions. The present invention provides a good initial approximate solution to pressure distribution within the reservoir. Based on the determination of an initial approximate pressure, a simple iterative method then is performed to obtain the final pressure (exact) distribution.

[00155] It has been found that the that conventional Gauss Seidel processing methodology with no realistic initial pressure for heterogeneous reservoirs converges to different solutions than the actual depending on the starting pressure solution estimate. In contrast, the methodology of the present invention analytical solution converges to a correct and demonstrably accurate answer. Numerical experiments described above have shown that the present invention works well for two dimensional reservoirs, whether homogenous and heterogeneous.

[00156] Advantages of the present invention are twofold: the present invention develops a close estimate of pressure distribution for reservoir simulation. The iterative method used for reservoir simulation runs after determining an initial estimate of pressure distribution can be very simple. Therefore, the total process is attractive and cost effective. The present invention provides improved reservoir simulation based on improved determination of pressure distribution over the three dimensions of heterogeneous rock formations of a reservoir, and avoids prior problems with existing simulation techniques.

[00157] The invention has been sufficiently described so that a person with average knowledge in the field of reservoir modeling and simulation may reproduce and obtain the results mentioned in the invention herein. Nonetheless, any skilled person in the field of technique, subject of the invention herein, may carry out modifications not described in the request herein, to apply these modifications to a determined structure and methodology, or in the use and practice thereof, requires the claimed matter in the following claims; such structures and processes shall be covered within the scope of the invention.

[00158] It should be noted and understood that there can be improvements and modifications made of the present invention described in detail above without departing from the spirit or scope of the invention as set forth in the accompanying claims.