PEH:Reservoir Simulation

Publication Information

Petroleum Engineering Handbook

Larry W. Lake, Editor-in-Chief

Volume V – Reservoir Engineering and Petrophysics

Edward D. Holstein, Editor

Chapter 17 – Reservoir Simulation

Rod P. Batycky ,Marco R. Thiele, StreamSim Technologies, Inc.; K.H. Coats, Coats Engineering, Inc.; Alan Grindheim, Dave Ponting, Roxar Software Solutions; John E. Killough, Landmark Graphics; Tony Settari, U. of Calgary and Taurus Reservoir Solutions Ltd.; L. Kent Thomas, ConocoPhillips; John Wallis, Wallis Consulting Inc.; J.W. Watts, Consultant; and Curtis H. Whitson, Norwegian U. of Science and Technology and Pera

Pgs. 1367-1398

ISBN 978-1-55563-120-8
Get permission for reuse

The Merriam-Webster Dictionary defines simulate as assuming the appearance of without the reality. Simulation of petroleum reservoir performance refers to the construction and operation of a model whose behavior assumes the appearance of actual reservoir behavior. The model itself is either physical (for example, a laboratory sandpack) or mathematical. A mathematical model is a set of equations that, subject to certain assumptions, describes the physical processes active in the reservoir. Although the model itself obviously lacks the reality of the reservoir, the behavior of a valid model simulates—assumes the appearance of—the actual reservoir.

The purpose of simulation is estimation of field performance (e.g., oil recovery) under one or more producing schemes. Whereas the field can be produced only once, at considerable expense, a model can be produced or run many times at low expense over a short period of time. Observation of model results that represent different producing conditions aids selection of an optimal set of producing conditions for the reservoir.

The tools of reservoir simulation range from the intuition and judgment of the engineer to complex mathematical models requiring use of digital computers. The question is not whether to simulate, but rather which tool or method to use. This chapter concerns the numerical mathematical model requiring a digital computer. The Reservoir Simulation chapter in the 1987 edition of the Petroleum Engineering Handbook^[1] included a general description of reservoir simulation models, a discussion related to how and why they are used, choice of different types of models for different-reservoir problems, and reliability of simulation results in the face of model assumptions and uncertainty in reservoir-fluid and rock-description parameters. That material is largely omitted here. Instead, this chapter attempts to summarize current practices and trends related to development and application of reservoir simulation models.

Models have been referred to by type, such as black-oil, compositional, thermal, generalized, or IMPES, Implicit, Sequential, Adaptive Implicit, or single-porosity, dual-porosity, and more. These types provide a confusing basis for discussing models; some refer to the application (e.g., thermal), others to the model formulation (e.g., implicit), and yet others to an attribute of the reservoir formation (e.g., dual-porosity). The historical trend, though irregular, has been and is toward the generalized model, which incorporates all the previously mentioned types and more. The generalized model, which represents most models in use and under development today, will be discussed in this chapter. Current model capabilities, recent developments, and trends will then be discussed in relation to this generalized model.

The 10 SPE Comparative Solution Project problems, SPE1 through SPE10,^[2]^[3]^[4]^[5]^[6]^[7]^[8]^[9]^[10]^[11] are used for some examples below. Table 17.1 gives a brief description of those problems.

Table 17.1

Overview

The Generalized Model

Any reservoir simulator consists of n + m equations for each of N active gridblocks comprising the reservoir. These equations represent conservation of mass of each of n components in each gridblock over a timestep Δt from tⁿ to tⁿ⁺¹. The first n (primary) equations simply express conservation of mass for each of n components such as oil, gas, methane, CO₂, and water, denoted by subscript I = 1,2,…, n. In the thermal case, one of the "components" is energy and its equation expresses conservation of energy. An additional m (secondary or constraint) equations express constraints such as equal fugacities of each component in all phases where it is present, and the volume balance S_w + S_o + S_g + S_solid = 1.0, where S solid represents any immobile phase such as precipitated solid salt or coke.

There must be n + m variables (unknowns) corresponding to these n + m equations. For example, consider the isothermal, three-phase, compositional case with all components present in all three phases. There are m = 2n + 1 constraint equations consisting of the volume balance and the 2n equations expressing equal fugacities of each component in all three phases, for a total of n + m = 3n + 1 equations. There are 3n + 1 unknowns: p, S_w, S_o, S_g, and the 3(n – 1) independent mol fractions x_ij , where i = 1,2,..., n – 1; j = 1,2,3 denotes the three phases oil, gas, and water. For other cases, such as thermal, dual-porosity, and so on, the m constraint equations, the n + m variables, and equal numbers of equations and unknowns can be defined for each gridblock.

Because the m constraint equations for a block involve unknowns only in the given block, they can be used to eliminate the m secondary variables from the block’s n primary or conservation equations. Thus, in each block, only n primary equations in n unknowns need be considered in discussions of model formulation and the linear solver. The n unknowns are denoted by P_i1, P_i2 ,…, P_in, where P_in is chosen as pressure p_i with no loss of generality. These primary variables may be chosen as any n independent variables from the many available variables: phase and overall mol fractions, mol numbers, saturations, p, and so on. Different authors choose different variables.^[12]^[13]^[14]^[15] Any sensible choice of variables and ordering of the primary equations gives for each gridblock a set of n equations in n unknowns which is susceptible to normal Gaussian elimination without pivoting. The (Newton-Raphson) convergence rate for the model’s timestep calculation is independent of the variable choice; the model speed (CPU time) is essentially independent of variable choice.

The Ith primary or conservation equation for block i is

....................(17.1)

where M_iI is mass of component I in gridblock i, q_ijI is the interblock flow rate of component I from neighbor block j to block i, and q_iI is a well term. With transposition, this equation is represented by f_iI = 0, the Ith equation of gridblock i. All n equations f_iI = 0 for the block can be expressed as the vector equation F_i = 0 where f_iI is the Ith element of the vector F_i. Finally, the vector equation

....................(17.2)

represents the entire model, where the ith element of the vector F is F_i. F is a function of the N vector unknowns P_i, where the Ith scalar element of P_i is P_iI. Application of the Newton-Raphson method gives

....................(17.3)

where δP is P^{l +1}–P^l and the N × N matrix A represents the Jacobian ∂F/∂P. The element A_ij of A is itself an n × n matrix ∂F_i/∂P_j with scalar elements a_rs = ∂f_ir /∂P_js, r and s each = 1,2,..., n. Eq. 17.3 is solved by the model’s linear solver. The matrix A is very sparse because A_ij is 0 unless block j is a neighbor of block i.

The calculations for a timestep consist of a number of Newton (nonlinear or outer) iterations terminated by satisfaction of specified convergence criteria. Each Newton iteration requires:

(a) Linearization of the constraint equations and conservation Eq. 17.1.
(b) Linear algebra to generate the A matrix coefficients.
(c) Iterative solution of Eq. 17.3 (inner or linear iterations).
(d) Use of the new iterate P^l+1 to obtain from Eq. 17.1 the moles of each component in the gridblock.
(e) A flash to give phase compositions, densities, and saturations which allow generation of the A matrix coefficients for the next Newton iteration.

Model Formulations

A major portion of the model’s total CPU time is often spent in the linear solver solution of Eq. 17.3. This CPU time in turn reflects the many multiply operations required. The model formulation has a large effect on the nature and expense of those multiplies.

Implicit vs. Explicit. The interblock flow term in Eq. 17.1,

....................(17.4)

uses phase mobilities, densities, and mol fractions evaluated at the upstream blocks. A gridblock is implicit in, say, the variable S_g if the new time level value S_gⁿ⁺¹ is used to evaluate interblock flow terms dependent upon it. The block is explicit in S_g if the old time level value S_gⁿ is used.

The Implicit Forumulation. The implicit formulation^[16] expresses interblock flow terms using implicit (new time level) values of all variables in all gridblocks. As a consequence, all nonzero A_ij elements of the A matrix of Eq. 17.3 are full n × n matrices. The resulting multiplies in the linear solver are then either matrix-matrix or matrix-vector multiplies, requiring work (number of scalar multiplies) of order n³ or n², respectively.

The IMPES Formulation. Early paper^[17]^[18]^[19] presented the basis of the IMPES (implicit pressure, explicit saturations) formulation for the black-oil case: take all variables in the interblock flow terms explicit, except for pressure, and eliminate all nonpressure variables from the linearized expressions for M_iIⁿ⁺¹ in Eq. 17.1. The obvious extension to any type model with any number of components was presented later,^[20] and numerous IMPES-type compositional models have been published.^[13]^[14]^[15]^[21]

The model Eq. 17.3 can be written as:

....................(17.5)

If all variables but pressure are explicit in the interblock flow terms, then all entries but those in the last column of the n × n A_ij (j ≠ i) matrix are zero (recall, the n th variable in each gridblock, P_in, is pressure p_i). This allows elimination of all nonpressure variables and reduction of the vector Eq. 17.5 to the scalar equation in pressure only^[22]:

....................(17.6)

or

....................(17.7)

where A is now a scalar N × N matrix and the P and F vectors have N scalar elements p_i and f_i , respectively. The multiplications required in solution of the IMPES pressure Eq. 17.7 are scalar multiplications, requiring a small fraction of the work of the matrix-matrix and matrix-vector multiplications of the implicit formulation. Thus, the model CPU time per gridblock per Newton iteration for moderate or large n is much less for the IMPES formulation than for the implicit formulation.

The Sequential Formulation. The stability of the IMPES formulation for the two-phase water/oil case was improved by following the IMPES pressure equation solution with solution of a water saturation equation using implicit saturations (mobilities).^[23] This concept was extended to the three-phase case and called the sequential formulation.^[24] For each Newton iteration, this method requires solution of the IMPES pressure Eq. 17.7, followed by solution for two saturations from a similar equation where the A_ij elements of A are 2 × 2 matrices.

A sequential compositional model was described^[15] and mentioned the desirability of a sequential implicit treatment of mol fractions in addition to saturations.

The Adaptive Implicit Forumlation. The Adaptive Implicit Method (AIM) ^[25] uses different levels of implicitness in different blocks. In each gridblock, each of the n variables may be chosen explicit or implicit, independent of the choices in other gridblocks. The choices may change from one timestep to the next. This results in the same equation AδP = – F^l as the Implicit formulation except that the elements A_ij of the A matrix are rectangular matrices of variable size. The numbers of rows and columns in A_ij equal the numbers of implicit variables in blocks i and j, respectively; all A_ii are square matrices. The CPU expense per Newton iteration of an AIM model lies between those of IMPES and Implicit models, tending toward the former as more blocks are taken implicit in pressure only.

Choice of Formulation. For a given problem, the previous four formulations generally give widely different CPU times. Generalizations regarding the best formulation have many exceptions. Arguably, the trend is or should be toward sole use of the AIM formulation. This is discussed in the Stable Step and Switching Criteria sections to follow. Current simulation studies use all of these formulations. The Implicit formulation is generally faster than IMPES for single-well coning studies, and for thermal and naturally fractured reservoir problems. For other problems, IMPES is generally faster than Implicit for moderate or large n (say, n > 4). Most participants used IMPES for SPE Comparative Solution Project problems SPE1, SPE3, SPE5, and SPE10. All participants used the Implicit formulation for SPE2, SPE4, SPE6, and SPE9. No participants in SPE1 through SPE10 used a Sequential model, and, with few exceptions, none used AIM.

A frequently stated generalization is that numerical dispersion error is significantly larger for Implicit than for IMPES formulations. Truncation error analysis^[26] shows this error to be proportional to Δx + uΔt for Implicit and Δx – uΔt for IMPES. Real problem nonlinearities and heterogeneity render the analysis approximate and the generalization of limited merit. For example, Figs. 17.1 through 17.3 show virtually identical Implicit and IMPES results for the black-oil 9,000-block SPE9 and 2,000-block gas/oil SPE10 problems. For SPE9 (SPE10), the average timestep was 67 (9.7) times larger for Implicit than for IMPES. The percentage of total CPU time spent in the linear solver for IMPES (Implicit) was 23.7 (57.3) for SPE9 and 35.4 (73.4) for SPE10.

Fig. 17.1 – IMPES vs. implicit oil rate and cumulative oil production, SPE9
Fig. 17.2 – IMPES vs. implicit GOR and water cut, SPE9.
Fig. 17.3 – IMPES vs. implicit oil rate, cumulative oil production, and GOR, SPE10 Model 1.

Implementation is an important factor in the relative efficiencies of different formulations. For a given problem, different models using the same formulation can give widely different CPU times. For example, the IMPES CPU times reported by different participants in the six-component compositional SPE5 problem^[6] differed by a factor of over 50.

Advances in Model Forumlations

The IMPES formulation was improved by concepts of relaxed volume,^[13]^[14]^[15] better choice of variables, 13 and "adaptive" flash calculations.^[13]

Relaxed Volume. The relaxed volume concept relates to the timestep calculation Steps (d) and (e) given previously. Step (d) gives the mass of each component in the gridblock, , which in turn gives overall composition {z_I}^l+1. The Step (e) flash then gives phase amounts and densities which in turn give new iterate S_w, S_o, and S_g values. These saturations do not sum to 1.0 because of the nonlinear nature of the conservation Eq. 17.1. If the saturations are altered (e.g., divide each by their sum) to exactly satisfy the volume balance Σ_J S_J = 1, then an incremental (timestep) mass-balance error occurs. If the saturations are not altered, then mass is conserved but there is a volume-balance error Σ_J S_J – 1. The authors^[13]^[14]^[15] chose to preserve mass and carry forward the volume balance error from iterate to iterate and step to step. The volume balance going into iterate . This in effect conserves both mass and volume because there is no permanent or accumulating volume error—only that of the given timestep. Equally important, there is no need to iterate out the volume error to a "tight" tolerance, and Newton iterations and model CPU are reduced. In contrast, the previous or historical IMPES procedure reset saturations to preserve volume and iterated out the mass-balance error. Because the latter error was not carried forward, more Newton iteration (and CPU time) was required to keep the permanent, accumulating mass balance error tolerably low. This use of relaxed volume with carryover also reduces Newton iterations and CPU time in the Implicit formulation.^[21]

This discussion implies some fundamental advantage of preserving mass and iterating out volume error as opposed to preserving volume and iterating out mass error. In the writer’s opinion, that is not true provided the error is carried forward in both cases. The Newton iteration requirement and CPU time should be similar if "equivalent" mass and volume error tolerances are used as convergence criteria.

Variable Choice. The linear algebra required to reduce the gridblock’s n conservation equations to the IMPES pressure equation is influenced by the choice of variables. The influence is absent for black oil, moderate for "moderate" n and up to a factor of three for large n (say, > 15).^[22] The choices of p and mol fractions {z_I}^[13] or mol numbers^[14]^[15] are better than the choice of p, saturations, and phase mol fractions^[12] for large n. The effect of this variable choice on total CPU time is often small because the affected work is often a small part of total CPU time. This IMPES reduction is absent in the Implicit formulation and the last of the above variable choices is arguably preferable.^[22]

Adaptive Flash Calculations.^[13] The work of EOS flash calculations, including the generation of fugacities and their derivatives, can significantly affect model efficiency when the linear solver does not dominate total CPU time. There may be little need to perform (most of) that work in a gridblock when p and composition are changing slowly. Use of internal, intelligent criteria dictating when that work is needed can significantly reduce the total-run flash calculation CPU time.^[13] This is similar in principle to the AIM selection of explicit variables for gridblocks which are quiescent in respect to throughput ratio.

Stable Timestep and Switching Criteria

This topic relates to the observation that lower run turnaround time can increase benefits from a reservoir study allotted a budgeted time period. As a corollary, time spent in repeated runs fighting model instabilities or time-stepping is counterproductive. While many factors affect this run time, it always equals the product (CPU time/step) × (number of timesteps). The first factor is "large" and the second "small" for the Implicit formulation, and conversely for the IMPES formulation. IMPES is a conditionally stable formulation requiring that Δt < Δt* to prevent oscillations and error growth, where Δt* is maximum stable timestep. The conditional stability stems from the explicit treatment of nonpressure variables in the interblock flow terms. Mathematicians performed stability analyses for constant-coefficient difference equations bearing some resemblance to IMPES. Authors in our industry extended and applied their results to derive expressions for Δt*, in particular,^[27]

....................(17.8)

for the black-oil 3D case of gas/oil flow. This shows that stable step Δt* is dependent upon flow rates, phase mobility, and capillary pressure derivatives, which of course vary with time and from one gridblock to another. Thus, at a given timestep, there are block-dependent stable step values Δt*_i , where 1 < i < N, and the IMPES stable step is Min(i) Δt*_i . An IMPES model using this internally determined stable step will run stably but may suffer from the weakest-link principle. As an extreme example, consider a 500,000-gridblock problem where, over a 100-day period, the Δt*_i value is 0.01 day for one block and > 30 days for all other blocks. The IMPES model will require 10,000 timesteps over the 100-day period.

In the AIM formulation, the stable step Δt*_i depends upon the number and identities of variables chosen explicit in block i; theoretically, Δt*_i = ∞ if all block i variables are chosen implicit. In the previous example, all nonpressure variables could be chosen implicit in the block where Δt*_i = 0.01 and explicit in all other blocks. The AIM model would then require CPU time/step essentially no greater than the IMPES model but would require only three timesteps for the 100-day period.

Numerous papers^[28]^[29]^[30]^[31]^[32]^[33] address the problem of determining expressions for the Δt*_i for use internally as switching criteria to select block variables as explicit or implicit in the AIM model. The stability analyses involved are complex and may be impractically complex when allowing the implicit vs. explicit variable choice to include all permutations (in number and identity) of the n variables. The most reliable and efficient AIM models in the future will stem from continuing research leading to the following: (a) Δt*_i estimates which are "accurate," and (b) implicit vs. explicit variable choices, block by block, which are near-optimal^[34] and minimize total CPU time, (CPU time/step) × (number of steps).

The Linear Solver

Preconditioned Orthomin^[35] is the most widely used method for iterative solution of Eqs. 17.3 or17.7. Nested Factorization (NF) ^[36] and incomplete LU factorization [ILU(n)] ^[37] are the two most widely used preconditioners. The term "LU factorization" refers to the factoring of the matrix A into the product of a lower triangular matrix L and an upper triangular matrix U. That is an expensive operation but is straightforward, involving only Gaussian elimination. The term "ILU(n)" denotes incomplete LU factorization, where only limited fill-in is allowed and n is the "order of fill."^[37] NF performs exceptionally well when transmissibilities associated with a particular direction (in a structured grid) dominate those in other directions uniformly throughout the grid. In general, ILU(n) or red-black ILU(n)^[38] [RBILU(n)] is less sensitive than NF to ordering of the blocks and spatial variation of the direction of dominant transmissibilities. In addition, RBILU(n) or ILU(n) have the parameter n (order of allowed infill) which can be increased as needed to solve problems of any difficulty.

A literature search and discussions with numerous developers and users have failed to establish consensus on whether NF or ILU preconditioning is better. Some are strong advocates of one method and others are just as adamantly supportive of the other. But many find, like this writer, that the better method is problem-dependent and it is difficult to find a reliable a priori indicator for making an up-front choice. In the writer’s experience, (a) when NF works well, it is faster than ILU methods, (b) RBILU(0) with no residual constraint is frequently the best of the ILU variants and a good default choice, and (c) in some cases, global residual constraint with the ILU or RBILU method is beneficial.

Cartesian Grids and Reservoir Definition

For many years, simulation used orthogonal Cartesian grids. In the past 15 years, numerous papers have described local grid refinement and various non-Cartesian grids, as discussed in the Gridding section. These papers show that non-Cartesian grids can reduce grid-orientation effects and provide definition and accuracy near wells, faults, highly heterogeneous areas, and so on more efficiently than Cartesian grids. The premise that Cartesian grids cannot provide required accuracy efficiently in these respects has come to be accepted as a fact. In addition, advances in geophysics have led to geostatistical description of permeability and porosity on a fine scale once unimaginable. Increasingly, our papers include examples using thousands of gridblocks for two- or few-well "patterns," in part to reflect these geostatistical descriptions. The purpose of this section is to show, using a few examples, that Cartesian grids can provide adequate accuracy and reservoir and near-well definition efficiently in some cases, even without local grid refinement. No generalizations from the examples used are intended. For the most part, the examples are taken from the literature.

SPE7 is an x-direction horizontal well problem with a 9 × 9 × 6 Cartesian grid representing a 2,700 × 2,700 × 160-ft reservoir section. The specified block Δy values decrease from 620 to 60 ft at the well, presumably to increase near-well definition and accuracy of results. The Δx are uniformly 300 ft. Fig. 17.4 compares Case 1A results for the SPE7 grid with results using uniform areal spacing Δx = Δy = 300 ft. The near-well y-direction refinement of the specified grid has no effect and is not necessary in this problem.

Fig. 17.4 – Effect of near-well grid refinement, SPE7 Case 1A.

SPE8 is a gas/oil problem with one gas injector and two producers on the corners of a 5,000 × 5,000 × 325-ft square reservoir. A 10 × 10 × 4 Cartesian grid with uniform Δx = Δy = 500 ft is specified. Five participants compared their results for that grid with results from their (areally) locally refined or unstructured grids. They showed good agreement for grids having approximately four times fewer blocks than the 10 × 10 × 4 grid. Fig. 17.5 shows equally good agreement for a 5 × 5 × 4 (Δx = Δy = 1,000 ft) Cartesian grid with no local refinement.

Fig. 17.5 – Effects of coarser grids on producing GOR.

SPE10 (Model 1) is a 2D cross-sectional gas/oil problem with a geostatistical permeability distribution given on a 100 × 1 × 20 Cartesian fine grid. Coarse-grid submittals included results using upscaling and local grid refinement. A homogeneous 5 × 1 × 5 Cartesian grid with no alteration of relative permeability matched the 100 × 1 × 20 results nearly exactly.^[11]

SPE10 (Model 2) is a 3D water/oil problem with a 1.122 million-cell geostatistical grid. Some coarse grid submittals included sophisticated upscaling and gridding techniques with no pseudoization of relative permeability and grids from 4,810 to 70,224 blocks. Others used simple flow-based upscaling to 75- to 2,000-block Cartesian grids with moderate k_r changes. In general, the latter submittals showed the best agreement with the fine-grid solution.^[11]

Numerous papers show that non-Cartesian grids can significantly reduce the grid-orientation effects of Cartesian grids. However, most of the examples used to study those effects are highly adverse mobility ratio displacements in homogeneous, horizontal reservoirs. In reservoirs with more normal fluid mobilities, areal fluid movement is more strongly affected by heterogeneity and/or gravity forces associated with reservoir structure (variable dip), and grid-orientation effects tend toward a second-order effect. As an example, the SPE10 (Model 2) water/oil problem reservoir is highly heterogeneous. Fig. 17.6 compares five-point and nine-point field results for an upscaled 28 × 55 × 85 Cartesian grid. The close agreement indicates an absence of grid-orientation effects even though the unfavorable oil/water viscosity ratio is 10 and there is no dip.

Fig. 17.6 – 9-point vs. 5-point, SPE10 Model 2.

Example 17.1

Table 17.2 gives data for Example 17.1, a ¼ five-spot, vertical-well problem. Fig. 17.7 shows two block-centered grids (a) and (b) used for this type of problem. The four-fold smaller well blocks of grid (b) provide finer well definition and presumably increase the accuracy of results. Fig. 17.8 shows the identical results for 10 × 10 grid (a) and 11 × 11 grid (b). Results are nearly identical for the 5 × 5 grid (a) and 6 × 6 grid (b), and Fig 17.9 shows insignificant difference between 3 × 3 grid (a) and 4 × 4 grid (b) results. The grid (b) doubles the grid (a) IMPES run CPU time but contributes no greater accuracy. Well-index effects are not important here. When they are, a single one-layer single-phase run can be made to determine the index correction factor for grid (a) wells located in the corners of their gridblocks.

Table 17.2
Fig. 17.7 – Block-centered 3 x 3 and 4 x 4 Cartesian grids used in Example 17.1
Fig. 17.8 – Effect of (a) 10 x 10 vs. (b) 11 x 11 grid on Example 17.1 results.
Fig. 17.9 – Effect of (a) 3 x 3 vs. (b) 4 x 4 grid on Example 17.1 results.

Fig. 17.10 shows a small effect of grid refinement on Example 17.1 results for grids from 20 × 20 to 3 × 3. The results indicate little need to enhance near-well definition by unstructured grids or by grid refinement (global or local) for grids finer than 3 × 3 for this problem.

Fig. 17.10 – Effects of Cartesian grid coarsening, Example 17.1 results.

Example 17.2

Flexible non-Cartesian grids are shown to significantly reduce the required number of gridblocks.^[39] An example provided^[40] was water/oil coning in a horizontal well in a 600 × 300 × 230 m homogeneous reservoir. Results were: (a) a 25,823-block 31 × 17 × 49 Cartesian grid was required to obtain a converged solution, and (b) a 3D 2,066-block CVFE unstructured grid gave correct results. Table 17.3 gives data for Example 17.2, a similar problem. Fig 17.11 compares Example 17.2 results for 60 × 31 × 48 and 10 × 7 × 9 Cartesian grids. The 630-block coarse Cartesian grid results here agree as well with the Cartesian 60 × 31 × 48 fine-grid results as the reported 2,066-block CVFE results agree with the 31 × 17 × 49 Cartesian fine-grid results.

Table 17.3
Fig. 17.11 – Fine vs. coarse Cartesian grids, Example 17.2

Non-Cartesian grids are argued to define irregular reservoir boundaries more efficiently than Cartesian grids. This is not necessarily true. For over 30 years, many models have used active-block coding. While the Cartesian grid extends past boundaries to numerous inactive blocks, those inactive blocks are dropped by the model and require no computer storage or CPU. These numerous inactive blocks pose a problem only for models, if any, that do not use active-block coding.

Linear Solver—John Wallis and J.W. Watts

The linear equation solver is an important component in a reservoir simulator. It is used in the Newton step to solve the discretized nonlinear partial differential equations. These equations describe mass balances on the individual components treated in the model. For nonisothermal problems, an energy balance is added to the system. The matrix problem involves solving Ax=b, where A is typically a large sparse matrix, b is the right-side vector, and x is the vector of unknowns. In the IMPES formulation, there is a single unknown per cell pressure. In the fully implicit formulation, there is a fixed number n of unknowns per cell where n ≥ 2. In the adaptive implicit formulation, there is a variable number of unknowns per cell. In most formulations, pressure is an unknown for each cell. The matrix A typically has associated well constraint equations and well variables and may be partitioned in block 2 × 2 form as

....................(17.9)

where x_w is the well variable-solution vector and x_R is the reservoir variable-solution vector. The matrix A_ww is often diagonal. In this case, the well variables may be directly eliminated, and the iterative solution is on the implicitly defined matrix system

....................(17.10)

The well variables are then obtained by back substitution as

....................(17.11)

If A is large, solution of the matrix equations is impractical using direct methods such as Gaussian elimination because of computer storage or CPU time requirements. Iterative solution based on projection onto Krylov subspaces is typically used. These Krylov subspaces are spaces spanned by vectors of the form p(A)v, where p is a polynomial. Basically, these techniques approximate A^-1b by p(A)b . The commonly used methods for constructing p(A)b are ORTHOMIN^[36] and GMRES.^[41] Both methods minimize the residual norm over all vectors in span{b, Ab, A2b, …, A^m-1b} at iteration m. They should yield identical results. From a practical standpoint, it does not matter which is used.

A technique known as preconditioning can improve both the efficiency (speed in a typical problem) and robustness (ability to solve a wide range of problems at least reasonably well) of Orthomin or GMRES. Preconditioning involves transforming the original matrix system into one with the same solution that is easier to solve. As a rule, the robustness of the iterative scheme is far more dependent on the preconditioning than on the specific Krylov subspace accelerator used. The preconditioner M is a matrix that approximates A, and has the property that the linear systems of the form Mx = b are easily and inexpensively solved. For most linear solvers the following preconditioned system is solved:

The preconditioned ORTHOMIN(k) algorithm, which retains the last k A-orthogonal direction vectors, is given by

Compute r₀ = b- Ax₀. Set p₀ = M^-1r₀.
For j = 0,1,2, …, until convergence Do :
x_j+1 = x_j + α_jp_j
r_j+1 = r_j - α_jAp_j
End Do

Because of the nature of the reservoir simulation equations, only certain preconditioners are effective in solving them. Reservoirs are typically shaped like pancakes, being much broader than they are thick. This geometry leads to strong vertical connectivity. Some preconditioners exploit this property. The most commonly used such preconditioner is NF.^[37] The convergence of NF is sensitive to the cell ordering. Best results are usually obtained by ordering the cells first along the direction of highest transmissibility and then successively along directions of decreasing transmissibility. This nearly always means that NF should be ordered first in the vertical direction.

The other commonly used preconditioners are incomplete lower triangular/upper triangular factorizations of the matrix, or ILU(n), where n is the level of infill that is retained during the elimination process. Performance of these can be improved by using a red-black checkerboard ordering (also called D4 ordering) of the nodes.^[42] Red-black ordering on a five-point (in 2D) or seven-point (in 3D) grid leads to direct elimination of the unknowns at the red cells, leaving a system containing only the unknowns at the black cells. The result is a halving of the number of unknowns. An ILU preconditioner using red-black ordering with zero infill on the reduced system is referred to as RBILU(0) and is the most frequently-used form of ILU.

In IMPES models with either no faults or vertical faults only, RBILU(0) or ILU(1) combined with z-line additive corrections,^[43]^[44] typically converges very rapidly. NF can also be used effectively in cases involving vertical faults and pinchouts because the matrix retains the structure required by NF. Nonvertical faults interrupt the matrix structure that makes red-black orderings attractive. In models containing them, ILU is the method of choice. Nine-point discretizations also cause problems for red-black orderings, but cause no difficulty for NF.

NF and RBILU(0) are commonly used in implicit models. Another very effective approach exploits the fact that pressure is the "stiff" variable. The Combinative^[45] or CPR^[44] method is a two-step preconditioning that extracts a pressure equation from the implicit matrix. It iteratively solves for a pressure correction at each iteration, uses the pressure correction to form a new residual, applies an inexpensive implicit preconditioning such as diagonal scaling or line Gauss-Seidel to the new residual, and then uses the sum of the two steps as the approximate solution. In compositional models, this two-step method can be much faster than one-step methods.

Many models include well-constraint equations that add well pressures to the set of unknowns. A simple but effective way of dealing with these equations is to order them first in the global matrix and then eliminate the well pressures from the set of unknowns.^[46] In this approximate elimination, any infill terms are either column-summed or row-summed into the main diagonal of the reservoir matrix, which is then factored using NF. The reservoir matrix problem is then solved iteratively and the well variables are obtained by back substitution.

Parallel iterative solution typically uses a domain decomposition approach in which the grid is partitioned into domains that contain approximately the same number of cells. The partitioning should be done such that coupling is strongest between cells within a domain. As a result, domains normally are groups of columns of cells. One way to accelerate the iteration is to color the domains in red-black fashion and apply an NF-type procedure in which the outer level of nesting is the coupling between domains and either NF or an ILU variation is used to factor the individual domains.^[47]^[48]

The solution of matrices arising from unstructured grids typically involves some variant of ILU with cell orderings such as Reverse Cuthill-McKee (RCMK)^[49] or Minimum Degree Fill (MDF).^[50]^[51]

Gridding in Reservoir Simulation

Introduction

The aim of gridding in reservoir simulation is to turn the geological model of the field into a discrete system on which the fluid-flow equations can be solved.

The basic structure of an oil reservoir is a set of geological horizons representing bedding planes. The reservoir may contain faults, at which the strata are displaced. It is usually possible to identify many more layers in the geological model than it is practical to include in reservoir flow simulation, so some upscaling of rock properties will normally be carried out. Even after this process, the geology to be represented is rarely homogeneous at the scale of the simulation grid.

Two related issues are involved in choosing a grid for reservoir simulation: the accuracy with which the geological description of the reservoir is matched, and the discretization of the flow equations. In a classical finite-difference scheme, the point values of pressures and saturations are used as solution variables, and the differential operators that appear in the fluid-flow equations may be expanded as difference expressions of these point values to some order. An alternative approach is to use an integral finite-difference^[52] or finite-volume^[53] method in which the fluid-flow equations are integrated over a set of cell volumes. This yields a set of equations in which the mass conservation conditions for the fluid in the simulation cell volumes are related to the flows through the interfaces between those cell volumes. Rock properties such as porosity are assumed constant over the cell or controlled volume. This yields a discretization scheme which is conservative (each outflow from one cell is an inflow to another) and for which the fluid in place may be obtained straightforwardly. The mass conservation equations for a timestep from T to T + ΔT then become:

....................(17.12)

where V_pa is the pore volume of cell a, m_ca is the density of conserved component c in cell a, Q_ca is the injection or production rate of component c because of wells, and F_cpab is the flow rate of component c in phase p from cell a to its neighbor b. In general, the flows F_cpab may involve the solution values of a number of cells, the number of cells involved defining the stencil of the numerical scheme. The linear pressure dependence of flows given by Darcy’s law leads to an expression of the type:

....................(17.13)

M_cpax is the mobility of component c in phase p for the contribution to the flow between a and x, given by x_cp.K_rp/μ_p, where x_cp is the concentration of component c in phase p, K_rp is the relative permeability of phase p, and μ_p is the viscosity of phase p. This is often set to an upstream value of the mobility, depending upon the sign of the potential difference.

ΔФ_pax is the potential difference of phase p between cell a and cell x, which includes pressure, gravity and capillary pressure contributions:

....................(17.14)

The constant coefficients of mobility and potential difference products, T_ax , are commonly termed the transmissibilities.

When the flows between two cells a and b can be expressed as a function of the solution values in just those two cells, so that the summation over cells includes just x = b, the flow expression takes a two-point form. The flow expression then takes a simple form:

....................(17.15)

When solution values from other cells are required, the flow takes a multipoint form.^[54]

Other options for discretization are available, such as Galerkin finite elements^[55]^[56]^[57] and mixed finite-element.^[58] It is sometimes possible to cast a finite-element Galerkin discretization into the upstreamed transmissibility-based form.^[57]

Regular Cartesian Grids

A simple 3D grid is the regular Cartesian grid (Fig. 17.12). Cells in such a grid may be simply identified using their (i,j,k) index values.

Fig. 17.12 – Regular Cartesian grids.

Each of the grid elements will be assigned a single permeability or porosity value. In this case, it is possible to obtain the transmissibility value as a harmonic average:

....................(17.16)

where cell b is the neighbor to cell a in some direction and K is the cell permeability in that direction. A is the area of the cell orthogonal to the direction of flow, and d the dimension of the cell in that direction. Such a two-point transmissibility assumes a permeability tensor with primary axes aligned along the grid axes.

Although regular grids are normally defined in normal Cartesian coordinates, it is also possible to use an (r, Ф, z) radial system.^[52] The resulting grid is cylindrical and is important for the special case of near-well studies dominated by radial inflow. For a 3D system, regular grids yield seven-point schemes, in which the flow equations for a cell involve solution values for just the cell and its six neighbors. Not all the elements in the grid need represent active solution variables in the simulation. Some cells may be inactive, representing volumes of the reservoir with zero porosity. Such inactive cells are usually compressed out of reservoir simulation solution arrays prior to the memory and time-intensive flow solution stage, and enable reservoirs with irregular boundaries to be represented within extended simulation grids.

The horizons that delimit rock strata are generally not horizontal, but are dipped, curved, or faulted. Unless extremely fine, a true regular grid that is orthogonal in all three axes will be unable to assign rock properties accurately to cell volumes. Such a layer-cake structure can be used, but will generally misalign property values (Fig. 17.13) in which the orthogonal grid provides a rather poor match to the dipping strata represented by the shaded layers. However, it is possible that improving computer power will bring such rasterized grids to a level of refinement at which a sufficiently good representation may be obtained.

Fig. 17.13 – Orthogonal grid used to represent dip.

Dip-Normal Geometry. A simple variation of a regular grid, in which the regular grid is rotated to bring the layers of cells into alignment with the bedding planes. Such a description would only suit a reservoir with a single, constant angle of dip. As geological descriptions have improved, fewer and fewer model reservoirs are found to fit this simple pattern, and something more flexible is required.

Block-Center Geometry. A simple model in which transmissibility between blocks is calculated on the basis of linear interpolation between the center values of the cells. This is a simple way of representing variable dip, but is difficult to represent graphically in a consistent way. Pore volumes are calculated on the basis of a series of flat regular cells with variable depths (Fig. 17.14a), but transmissibilities are calculated on the basis of interpolated values (Fig. 17.14b). The areal grid is rectangular.

Fig. 17.14 – (a) Obtaining pore volumes; (b) obtaining transmissibilities.

Thus, for the pair of cells illustrated,

....................(17.17)

where A is the average area over which flow occurs and c is a dip correction given by cos²θ, where θ is the angle of dip of a line joining the cell centers to the horizontal. Such a block-center option is suitable for unfaulted reservoirs and is commonly supplied as a simulator option.

Hexahedral Grids

Further improvements in geological modeling threw an emphasis on describing faults, and made it important to distinguish depth displacements due to dip and faulting. This is difficult in block centre geometry in which the cell is positioned by its centre depth and Δx, Δy, Δz dimensions. To define faulting more precisely it is useful to define the position of grid cell by its corner point locations. A hexahedral shape with eight corners and bilinear planes as surfaces then describes the cell geometry. Faults, both vertical and inclined, may be described precisely (Fig. 17.15). Such grids are often called corner-point grids.

Fig. 17.15 – Hexahedral grid system.

In both the dipped and general hexahedral grids, the orthogonality of a completely rectangular grid no longer exists, and the result is that the two-point property of the flows between the cells is lost—the flow between cell a and cell b is not just a function of the solution values of cells and a and b.^[53]^[59]^[60]^[61]^[62] Typically, the result is a 27-point scheme in three dimensions. However, if the grid distortion is mild, it may be possible to ignore some additional couplings and use a low-order transmissibility scheme. This is normally done for extra couplings introduced by dip angles, which are often small.

Although this corner-point description handles the fault issue, the basic coordinate system remains a regular grid (i.e., the grid is structured). Fitting such a basically regular system to the irregular shapes of a reservoir remains a difficulty that may be solved in two basic ways—either by distorting the grid and fitting the cells into the geometry, or by truncating the grid to the reservoir position.

Multiple-Domain Hexahedral Grids

In some cases, a single structured grid system cannot match the overall structure of a reservoir, so a block grid or domain-based grid is used.^[63] This consists of a number of subgrids, each with a local regular (i,j,k) structure, but linked together to model the entire reservoir. The block hexahedral system gives rise to multiple (i,j,k) indexing systems—(i,j,k,l), where the l index specifies the grid system. These comprise a series of regular grids. Such regular gridding systems have advantages for upscaling and downscaling—for example, a natural coarsening of a regular grid may be simply defined by grouping sets of coordinates in each direction.

Grid Refinement

A common requirement in reservoir simulation is an increased level of detail around an item of interest such as a well. This is frequently obtained in structured grids by local grid refinement, replacing a set of cells in the original grid by a finer grid (Fig. 17.16). The inserted grid may be Cartesian (center) or radial (upper left). Local refinement may be regarded as a form of multiple domain structured grid, in that it consists of a number of linked structured grids. Flows at the edges of local refinements generally take a multipoint form.^[54]

Fig. 17.16 – Local grid refinement.

Unstructured Grids

The problems involved in using a regular structured grid to represent reservoir geometry can be avoided by using an unstructured grid. This is constructed around a set of solution points that need have no particular indexing scheme. These points may be triangulated into a mesh of triangles or tetrahedrons. A control volume is constructed around the nodes of the resulting mesh to define the simulator cell volumes. The perpendicular bisector (PEBI) method introduced into reservoir simulation by Heinemann^[64]^[65] used a technique also known as a Voronoi grid.^[58]^[66] Starting from any set of solution points, the PEBI cell volumes are defined by the perpendicular bisection planes between these points. The resulting control volume is defined by the perpendicular planes—it is the set of points closer to the node than any other. This is shown in Fig 17.17, in which the bisectors to the heavy lines joining the solution points enclose the control volume, represented by the shaded area. The grid is locally orthogonal, and the desirable property of two point flows is obtained. The actual cell volumes may have a variety of shapes, depending on the exact placement of the solution points, but are typically hexagonal in two dimensions. Grid refinement occurs naturally in areas where solution points are closely spaced.

Fig. 17.17 – Unstructured grids may result in hexagonal cells.

The two-point property is not naturally preserved in anisotropic reservoirs, although it can be regained by transforming to a K-orthogonal grid in which the geometry is transformed so that K_n is parallel to the vector joining the solution nodes, where K is the permeability tensor and n is the normal to the cell volume surface.^[67] For nonisotropic cases in which the grid is not K-orthogonal, the flows will be functions of the solution values in more than two cells, as in the general hexahedral case.

An unstructured grid may be defined in two dimensions, and then applied to each layer of a reservoir model, so that a typical cell is a hexagonal prism. This is sometimes termed a 2½D unstructured grid.^[68] Alternatively, a full 3D unstructured grid may be defined. The 3D approach is most effective when applied to model a local structure such as a branched well.

Unstructured grids yield an elegant and flexible grid description. However, the ability to identify cells by a simple set of indices is lost, and items such as wells need to be positioned in true space terms. The systems of linear equations generated by unstructured grids are also commonly regarded as more difficult to solve than those produced by structured grids. However, it may be more true to say that optimal solution schemes are simpler to find for structured grids, where the row and plane order provides a natural ab initio solution variable-ordering scheme.

Truncated Regular Grids^[69]

The truncated approach fits in well with the rectangular grids used in geological modeling. A simple rectangular grid is always used in the areal direction, but faults may subdivide the rock volume in a given column. The areal grid is not modified to match the faults. Thus the two marked volumes in Fig. 17.18 represent different cells, but may have the same i, j indices, so this creates a multiple-domain grid. A disadvantage is the more complex shape of cells at the edge of the grid. Transmissibilities for such cells may need to be calculated numerically. Apart from the truncated cells, all the grid cells are hexahedra that are rectangular in plan.

Fig. 17.18 – A fault creating two cell volumes in a truncated grid.

Other Gridding Systems

Triangular or Tetrahedral Grids. The underlying solution points of a PEBI mesh can be linked together into a Delaunay triangulation. In 2D, this creates triangles, and in 3D it creates tetrahedra. One option would clearly be to use triangular or tetrahedral cells directly and associate cell volumes with these. This technique is rather rarely used in reservoir simulation. Partly this may be historical, but the Delaunay triangulation is rather less stable under grid changes than a Voronoi grid, and triangulation can more often lead to "sliver" cells with a high surface area but a small volume.

Curvilinear Grid Systems. In some special cases, a transformed coordinate system may be used, based around an expected flow pattern. Such grids are not well adapted to represent geological data, and have been used less frequently as more detailed reservoir descriptions have become available.

Future Directions

Two themes emerge from current trends in reservoir simulation gridding. The increasing sophistication of data preparation and solver technology indicates a move towards unstructured grids as a general method of solving the flow equations for a given reservoir simulation problem. On the other hand, reservoir simulation is increasingly seen as part of a decision-making process rather than an isolated activity, so the ability to map easily onto the generally regular data structures used in seismic and geological modeling becomes an important issue. In this role, structured grids may have advantages of simplicity and scalability.

An ideal is to separate the construction of the flow-simulation grid from the description of the reservoir geometry. This ties in with a further ideal, inherent in many discretization schemes, that the scale of the simulation grid should be below the scale of the problem structure.

For more complex shape-dominated problems, the unstructured approach looks general and flexible, providing that the data-handling and cell-identification methods can be moved to true x,y,z space preprocessing software.

Upscaling of Grid Properties—Alan Grindheim

Definition

Upscaling, or homogenization, is substituting a heterogeneous property region consisting of fine grid cells with an equivalent* homogeneous region made up of a single coarse-grid cell with an effective property value. It is performed for each of the cells in the coarse grid and for each of the grid properties needed in the flow-simulation model. Therefore, the upscaling process is essentially an averaging procedure in which the static and dynamic characteristics of a fine-scale model are to be approximated by that of a coarse-scale model. A conceptual illustration of the upscaling process is shown in Fig. 17.19.

Fig. 17.19 – Concept of upscaling.

* Either volume or flux vice, depending on the type of property that is to be upscaled.

Can Upscaling Be Avoided?

Typically, 3D geological models contain detailed descriptions of the reservoir that can be hard to capture properly with a significantly coarser model. Therefore, it would be preferable if upscaling could be avoided. Currently, an average-sized flow simulation model consists of approximately 100,000 active grid cells. This is to ensure that the CPU consumption of a simulation run will be reasonable (i.e., within practical limits).

Because a typical 3D geological model may consist of approximately 10 million active grid cells, it is obviously infeasible to run fluid-flow simulations directly on the geological model. Hence, upscaling is a required part of current reservoir modeling workflows.

Seen through the eyes of the geologist, the upscaling task may be a painful experience because all the geological details that were put into the model seem to be lost in the process. For a reservoir engineer, on the other hand, effective properties might be all that matter.

For volumetric (additive) properties such as porosity and saturation, the effective flow-cell value is simply given by the bulk and pore volume weighted arithmetic average, respectively, of the geo cells inside it. For the permeability, which is intrinsic (nonadditive) by nature, no such simple averaging method exists. The complexity one needs to take into account when upscaling permeability is considerable; therefore, all current techniques provide only an approximation of the true effective cell permeability. This approximation may range from very good to very poor, depending on the complexity of the fine-scale permeability distribution as well as the upscaling method used.

Upscaling Techniques for Absolute Permeability

Homogenization of absolute permeability does not have an exact analytical solution, except for in a few idealized cases. The challenge of computing an accurate effective permeability has resulted in a large number of upscaling techniques. These techniques range from simple statistical averages to advanced numerical methods.

Tensor methods are the most accurate techniques available for computing the effective cell permeability. These are based on solving a second-order elliptic differential equation describing single-phase, incompressible, steady-state fluid flow in a region without sources and sinks (i.e., wells). Some flow-based methods may provide a full permeability tensor. However, because most multiphase flow simulators can only handle a diagonal permeability tensor because of the use of a seven-point stencil in 3D, diagonal tensor methods are most frequently used whether directly or indirectly (through a diagonalization of a full tensor). For a diagonal tensor, only the effective permeability in the principal directions of flow (x, y, and z) will be nonzero.

The flow equation is usually discretized with a finite-difference scheme, although finite-element methods are also applied occasionally. To compute all the directional components of the permeability tensor, the discretization and solution of the flow equation must be performed for each of the principal flow directions (i.e., three separate single-phase simulations need to be performed). Each simulation involves the iterative solution of a linear equation system (typically, the linear solver is a conjugate gradient method, preconditioned by incomplete Cholesky or LU factorization). The unknowns in this equation system are the geo-cell pressures inside the flow cell, whereas known quantities are the geo-cell dimensions and permeabilities, as well as the pressure conditions along the faces of the flow cell. When the numerical solution of the fine-scale pressure distribution has converged, an effective permeability is computed by equating expressions for the flux through the heterogeneous geo cells with the flux through the equivalent homogeneous flow cell using some form of Darcy’s law.

The pressure field is usually solved locally—that is, for one flow cell at a time. However, as discussed in the next subsection, the size of the computational region may not necessarily be limited to that of the upscaling region (i.e., the flow cell).

Upscaling Schemes for Absolute Permeability

Based on the size of the computational region, the single-phase upscaling process may either be described as local, regional, or global. With local upscaling techniques, the computational region is identical to the upscaling region (i.e., only geo cells inside the flow cell are considered in the upscaling computations). For regional upscaling, the computational region is expanded beyond that of the flow cell to include a buffer region of neighboring geo cells. In the case of global upscaling, the computational region is that of the entire geo model. Fig. 17.20 provides a schematic drawing of how the computational region varies with the different upscaling schemes. These are further discussed in the subsections that follow.

Fig. 17.20 – Upscaling schemes and the size of the computational region (geo grid in white, flow grid in black).

It should be noted that the different upscaling schemes are only relevant when considering flow-based (tensor) methods. It is also important to realize that even though the computational region may vary according to the scheme used, the upscaling region remains unchanged and is of course defined by the flow cell, as in the case of the simple, analytical upscaling techniques.

Local Upscaling

Because it used to be too time-consuming to compute the fine-scale pressure field for the complete geo grid in a single operation, the flow-based methods have traditionally been restricted to solving the pressure field locally—that is, for a single flow cell at a time. Hence, the effective cell permeability is computed separately and independently of the other flow cells, which may or may not be correct depending on how representative the imposed pressure conditions along the faces of the flow cell are.

Different types of artificial boundary conditions for the flow cell have been suggested over the years, all with the objective of providing as good an approximation of the real boundary conditions as possible. An important design criterion for the artificial boundary conditions is the conservation of flux in and out of the flow cell.

The first type of boundary conditions proposed for the local solution of the pressure equation was published by Warren and Price in 1961.^[70] Their approach is to impose a constant pressure gradient in a selected direction of flow by specifying a pressure of 1 on the inflow face and a pressure of 0 on the outflow face. By allowing no flow to pass through the sides of the cell, all fluxes are forced to go in the principal direction of flow. Therefore, this type of boundary conditions is often referred to as the no-flow or sealed-sides boundary conditions. The sealed-sides boundary conditions are graphically illustrated in Fig. 17.21 for flow in the vertical direction (here in the case of a flow cell containing a barrier).

Fig. 17.21 – The sealed sides boundary conditions (in z-direction).

The choice of boundary conditions emulates the way core permeability is measured in the lab. This is hardly a coincidence. As in the coreflood experiment, the local numerical flow simulation is in effect 1D because the cell faces parallel to the main flow direction are sealed. This implies that the estimated effective permeability will be scalar. Hence, the maximum number of directional permeability components that can be obtained with this type of boundary conditions is three, one for each of the principal directions of flow. In practice, the diagonal permeability tensor is derived by setting up the boundary conditions for x, y, and z directions, respectively, in three independent single-phase simulations.

As documented in two different sources^[71]^[72], a tensor technique based on the sealed-sides boundary conditions tends to bias the estimated effective permeability toward a low value. The physical implication of this is most clearly seen in the case of a bimodal permeability system of sand and shale. This is because the sealed-sides method consistently underestimates the reservoir flow characteristics by thickening shale barriers and narrowing sand channels. The latter effect also has a tendency of disconnecting stacked sand channels.

Take, for example, the flow illustrated by Fig. 17.21. Because the barrier extends across the entire length of the local upscaling region, the resulting effective permeability (in the z-direction) will be zero. For vertical flow, the result, therefore, is a thickening of the shale in the flow model equal to the thickness of the flow cell. Depending on which factors that affect fluid flow in the region of the cell, this may or may not be a representative value for that particular flow cell.

Strictly speaking, the sealed-sides boundary conditions are only valid if no wells are present and the flow cells are symmetric in each direction of the grid as illustrated in 2D by Fig. 17.22. Hence, the sealed-sides boundary conditions assume that the flow cell is surrounded by mirror images of itself.

Fig. 17.22 – The fine-scale permeabilities are symmetric on the scale of the flow cell (geo grid in white, flow grid in black).

By the end of the 1980s, 3D geological models had started to appear more regularly on the modeling scene. This resulted in a new demand for advanced upscaling. In this renewed effort, two alternative boundary conditions for solving the local pressure solution in a flow-based method were suggested more or less at the same time. One was based on linear boundary conditions, the other on periodic boundary conditions.

The use of linear boundary conditions in flow-based upscaling was suggested by Guerillot et al.^[73] in 1989 and Samier^[74] in 1990 to enable the computation of a full-permeability tensor. Instead of setting the flow through the sides of the cell to zero, the pressure along the sides is allowed to vary in a linear fashion that matches the constant pressure on the two cell faces perpendicular to the flow. Hence, the imposed pressure gradient is still constant, but the flow is allowed to enter and leave the cell at any point along the sides parallel to the main flow direction. Therefore, this type of boundary conditions is also referred to as the open-sides boundary conditions. The situation is graphically illustrated in Fig. 17.23 for flow in the vertical direction (here in the case of a flow cell containing a barrier).

Fig. 17.23 – The open sides boundary conditions (in z-direction).

As with the sealed-sides boundary conditions, three independent single-phase simulations, with the main flow direction in x, y, and z, respectively, are needed to yield all of the components of the permeability tensor. With open-sides boundary conditions, however, also the off-diagonal components will generally be nonzero. Hence, unlike the sealed-sides boundary conditions where the effective permeability is limited to that of a diagonal tensor, the open-sides boundary conditions, as previously mentioned, give a full permeability tensor. The resulting full tensor may be either symmetric or nonsymmetric depending on the properties of the method under consideration.

As documented in two different sources^[75]^[76], a tensor technique based on the open-sides boundary conditions tends to bias the estimated effective permeability toward a high value. The physical implication of this is most clearly seen in the case of a bimodal permeability system of sand and shale. This is because the open-sides method consistently overestimates the reservoir flow characteristics by narrowing shale barriers and thickening sand channels. The latter effect also has a tendency of connecting isolated sand channels.

Take, for example, the situation illustrated by Fig. 17.23. Even though the barrier extends across the entire length of the local upscaling region, the resulting effective permeability (in the z-direction) will be significantly larger than zero. For vertical flow, the result is therefore a narrowing of the shale in the flow model equal to the horizontal dimensions of the flow cell. Depending on which factors affect fluid flow in the region of the cell, this may or may not be a representative value for that particular flow cell.

The use of linear boundary conditions has its origin in the effective medium theory,^[77] which states that any region of permeability behaves as if embedded within the average medium. Strictly speaking, these boundary conditions are therefore only valid if the neighboring flow cells are of a uniform, nonzero permeability. This is illustrated in 2D by Fig. 17.24.

Fig. 17.24 – The fine-scale properties are uniform and nonzero outside the upscaling region (the geo cells outside the upscaling region are drawn in gray to reflect that the average value need not be known).

The use of periodic boundary conditions originates from the volume averaging theory, and its use in flow-based upscaling was first introduced by Durlofsky and Chung in 1990 and by Durlofsky in 1991.^[75]^[76] Durlofsky used periodic boundary conditions, together with Darcy’s law and the classic requirement of flux conservation, to derive a full permeability tensor. A somewhat different approach, which also uses a periodic pressure field around the flow cell, was proposed by Øistein Bøe et al. in 1994.^[78]^[79] This uses a weak form of Darcy’s law to prove that periodic boundary conditions result in a full permeability tensor that is both symmetric and positive definite. The Norsk Hydro tensor method is based on the conservation of dissipation (mechanical energy per unit weight of fluid), although it turns out that fluxes are conserved as well.

A simplistic illustration of the periodic boundary conditions is given in Fig. 17.25 for flow in the vertical direction (here in the case of a flow cell containing a barrier).

Fig. 17.25 – The periodic boundary conditions in z-direction (P_1k = P_nx,k, P_i1 = P_i,nz).

Although the periodic boundary conditions generally result in an effective permeability that is higher than that computed with the sealed-sides boundary conditions, the effective vertical permeability for the upscaling region illustrated in Fig. 17.25 will also be zero.

Strictly speaking, the periodic boundary conditions are only valid if no wells are present and the fine-scale medium is periodic on the scale of the flow cells (i.e., the fine-scale property distribution inside each flow cell must be identical). This is illustrated in 2D by Fig. 17.26. Please note that if a medium is symmetric on the scale of ΔL, then it will be periodic on the scale of 2ΔL.

Fig. 17.26 – The fine-scale properties are periodic on the scale of the flow cell (geo grid in white, flow grid in black).

The relative performance of the tensor methods that is caused by the various boundary conditions has proven to be of considerable interest. As it happens, the sealed-sides method provides a lower bound and the open-sides method an upper bound of the effective permeability. The periodic-based method turns out to give an effective permeability estimate that generally lies in between the two previous methods.

With regard to the outer bounds of effective permeability, it is well known that the harmonic and arithmetic means provide the absolute lower and upper limit of the effective permeability, respectively. It is less known that the uncertainty range in the effective permeability may be narrowed using the composite averages. In fact, it may be mathematically proven that the harmonic-arithmetic average provides a closer lower limit than the pure harmonic mean, whereas the arithmetic-harmonic average provides a closer upper limit than the pure arithmetic mean (truly valid only for regular grids). In this context, it is important to realize that the two flow-based methods (sealed and open sides) provide an even narrower uncertainty band for the effective permeability, but at the expense of increased CPU time.

The relative performance of the most important local upscaling techniques is shown in Fig. 17.27. Fig. 17.27 also indicates the "inner" uncertainty range of the true effective permeability for the sake of comparison.

Fig. 17.27 – The relative performance of local upscaling techniques.

Regional Upscaling

Regional upscaling is applied to reduce the influence of the artificial boundary conditions on the effective permeability estimate by moving the boundary of the computational region away from the flow cell. This implies that the influence of neighboring geo cells is taken into account in addition to the geo cells inside the flow cell. In other words, regional upscaling represents an expansion of the local computational region outside the volume of the flow cell. The size of the so-called buffer or skin around each flow cell is usually given in number of neighboring geo cells to either side of the flow cell and must be specified by the modeler for each of the three coordinate directions.

The permeability estimate of a regional upscaling method will improve as the size of the buffer region increases, and it will ultimately be equal to the "true" effective permeability when the buffer size has reached the boundaries of the geo model for all three directions. The gain in accuracy is largest in the beginning (i.e., for small buffer values). This is illustrated by Fig. 17.28, showing the behavior of the lower- and upper-bound tensor techniques in the case of increasing buffer size. Please note that in Fig. 17.28, the outer bounds are shown to be symmetric around the "true" effective permeability. Generally, this is not the case.

Fig. 17.28 – The outer bounds of effective permeability and their behavior with increasing buffer size.

Global Upscaling

Strictly, the fine-scale pressure field must be determined for the entire geo grid simultaneously to compute "exact" effective permeabilities for the flow cells. In the past, however, this has been too CPU-intensive to be performed in practice. With the introduction of new and promising solution algorithms such as the Output Least Squares (OSL) method, global upscaling schemes can now be realized. In the paper by Holden and Nielsen,^[80] the OSL method is used to minimize the difference in pressure, as well as velocity, between the geo and flow grids in an iterative process. Because the CPU consumption of the applied equation solver is proportional to the number of geo cells, a global solution will use approximately the same amount of computational time as the sum of all the local computations. Therefore, the new global upscaling scheme is just as fast as any local method.

An obvious advantage with the global upscaling approach is that one avoids using artificial boundary conditions around the upscaling region (i.e., instead of guessing what the boundary conditions for the flow cells might be, the pressure conditions surrounding the cells are explicitly known). Another important benefit is that a poor separation of scales in the upscaling will no longer occur because the size of the computational region is the same as the geo model.

Although still in its research stage, global upscaling has much potential for improving today’s permeability estimators, especially for models containing a complex facies architecture with large permeability contrasts between facies. In fact, according to Holden and Nielsen,^[80] preliminary results show an improvement factor of 10 in some cases.

Still, as discussed by Holden and Nielsen,^[80] the global upscaling approach is not enough to ensure maximum accuracy in the modeling of the effective permeability. Because the value of the effective permeability is influenced by changes in the pressure field, the flow-cell permeabilities should strictly be recomputed by the global method for every timestep taken by the multiphase flow simulator. In practice, though, it might be good enough to update the effective permeability field whenever a significant change occurs as a result of altering the well configuration or production/injection rates and so on. Hence, the ultimate upscaling scheme for the absolute permeability might be the one that is coupled with the multiphase flow simulator and automatically updates the absolute effective permeability field for each timestep. With the current computer power and the lack of proper integration between the geological model and the simulation model, this is hardly achievable yet.

Best Practice Guidelines

As may be understood from the previous sections, there exists no single upscaling method for absolute permeability that is superior to all other methods in all situations, at least not until it has been fully established that the global upscaling scheme represents the ultimate method of choice. Selecting the proper upscaling method from the many available choices can be quite a challenge. The choice of sophistication in the upscaling method generally depends on one or several of the following factors:

The complexity of the fine-scale permeability distribution (i.e., the geo model).
The degree of upscaling that needs to be performed (i.e., the coarsening factor).
The number of permeability realizations that need to be upscaled.
The time available to the project for performing upscaling.
The intended use of the flow model.

Because an exact validation of the upscaling process cannot be performed unless a multiphase flow simulation is carried out on the geo model itself, two alternative upscaling approaches for identifying the proper homogenization method are presented here.

The Absolute Upscaling Approach. This approach assumes that there exists a way to properly validate the absolute performance of an upscaling method without resorting to an extremely time-consuming (if at all feasible) finite-volume simulation of the geo model. As documented in Samier et al.^[81], streamline simulation offers a very efficient way of validating the performance of upscaling methods. The validation process is carried out by first running a streamline simulation on the geo model to compute the reference solution. Then a streamline simulation is run on the flow model for each of the upscaling methods that are to be evaluated. The simulated performance of the various upscaling methods is then compared to that of the geo model. The validation of upscaling methods is best done under single-phase flow conditions to avoid introducing other model parameters (e.g., relative permeabilities and associated rock types) that also need to be upscaled in one way or the other. A higher confidence may also be obtained for the validation process if the actual well pattern is used in the streamline simulations.

Using the previously described validation scheme, the modeler may choose to evaluate any upscaling method until one with a satisfactory performance is found. Still, a more systematic way of identifying the optimum upscaling method is desirable. With the absolute upscaling approach that is presented here, the modeler is offered a multistep procedure that is to be terminated as soon as a satisfactory upscaling method has been identified. The recommended procedure involves the following steps:

Compute the upper and lower bounds of the effective permeability using the arithmetic-harmonic (or pure arithmetic) and harmonic-arithmetic (or pure harmonic) average techniques, respectively. Being of the analytical type, these methods are very fast and will provide a first quantification of the upscaling uncertainty. Validate the performance of the two composite methods against that of the geo model using a single-phase streamline simulator with the actual well pattern.
If the performance of any of the two methods in Step 1 is within an acceptable range of the geo model, then terminate the procedure and choose the appropriate method. If, on the other hand, the performance of both methods is unsatisfactory because of the complexity of the geo model, then use the upper (open-sides) and lower (sealed-sides) bound diagonal tensor methods to narrow the uncertainty in the flow-model performance. Validate the performance of the two tensor methods against that of the reference solution.
If the performance of any of the methods in Step 2 is within an acceptable range of the reference solution, then terminate the procedure and choose the appropriate method. If, on the other hand, the performance of both methods is unsatisfactory, then the following alternatives may be worth considering:
a. If time allows, refine or coarsen the flow grid (whatever is best) to achieve a better separation of the length scales. Then repeat the upscaling of the outer bounds (in Step 2) and redo the validation to check if the performance of either method has improved.

b. Apply a tensor method with periodic or semi-open-boundary conditions [the semi-open boundary conditions alternative is available in some applications using a multiplier between 0 (sealed) and 1 (open) to the side faces of the computational region]. As previously mentioned, this should result in an intermediate estimate of the effective permeability tensor and therefore provide a flow model performance that lies somewhere in between the two methods in Step 2.

c. Select the best of the two tensor methods in Step 2 and convert the local method to a regional method using a buffer region of modest size. Validate its performance. If necessary, repeat this step using an increasingly larger buffer region until a satisfactory performance of the flow model is obtained.

d. If for some reason none of the previous alternatives are an option, then one needs to apply the method that best satisfies the wanted flow behavior of a given cell or cells in a given region. In other words, a combination of the outer bound techniques within one and the same model may be a fourth alternative.

As previously mentioned, the method using open-sides boundary conditions is a good estimator of sand continuity and quality, whereas the sealed-sides boundary conditions method is better at detecting the presence and effect of barriers.

Consider a long horizontal oil producer in the Troll West Gas Province that is protected against coning from the overlying gas cap by a calcite barrier just above the well. If the vertical grid resolution in the Troll full-field model was such that one could apply the open-sides boundary conditions technique on the cells containing the well, and the sealed-sides boundary conditions technique on the cells containing the barrier, then this would be the optimum local upscaling approach.

However, if both a segment of the well trajectory and a segment of the calcite are present inside the same flow cell, then the open-sides boundary conditions technique will give a good estimate of the well’s PI but result in a much too early gas breakthrough, whereas the sealed-sides boundary conditions method will better capture the effect of the calcite but give a too low estimate of the well’s PI. If this is the case, one needs to consider applying one of the alternatives A or C (alternative B will, in this particular example, give the same result as the sealed-sides method).

The multistep procedure of the absolute upscaling approach is graphically illustrated in Fig. 17.29.

Fig. 17.29 – Flow chart for the absolute upscaling approach.

The Relative Upscaling Approach. This approach acknowledges the fact that an exact validation of the upscaling results cannot be achieved in practice. Therefore, instead of trying to validate the absolute performance of an upscaling method, the approach diagnoses the relative performance of outer bound methods using the actual multiphase finite-volume simulator. This implies that a full black-oil simulation is run on the flow model for each of the upscaling methods that are to be evaluated. The deviation in the simulated performance between outer bound methods will then reflect the part of the model uncertainty that originates from the upscaling process itself. To ensure a high degree of relevance in the diagnostics, it is important that the test simulations contain a representative description of the actual flow model.

If the project is pressed for time, the simulations may be skipped altogether in favor of a faster, although less robust, way of performing the diagnostics. Instead of analyzing simulation profiles, a normalized difference parameter may be computed on a cell by cell basis using the formula

....................(17.18)

The relative upscaling approach that is presented here utilizes a multistep procedure that applies outer bound methods of increasing accuracy until the best possible upscaling method can be identified. The recommended procedure involves the following steps:

Compute the upper and lower bounds of the effective permeability using the arithmetic-harmonic (or pure arithmetic) and harmonic-arithmetic (or pure harmonic) average techniques, respectively. Being of the analytical type, these methods are very fast and will provide a first quantification of the upscaling uncertainty. Run the finite-volume simulator for each of the two composite methods (or compute a grid-based difference parameter) and perform the diagnostics.
If the performance gap between the two methods in Step 1 is acceptable (small), then terminate the procedure and choose either of the two. If the deviation in performance is unsatisfactory because of the complexity of the geo model, then use the upper (open-sides) and lower (sealed-sides) bound diagonal tensor methods to narrow the uncertainty in the flow model performance. Run the finite-volume simulator for each of the two tensor methods (or compute a grid-based difference parameter) and perform the diagnostics.
If this reduces the upscaling uncertainty to within acceptable limits, then either of the two diagonal tensor methods may be used to provide the final permeability field for the flow model. If the deviation in the performance is unsatisfactory (large), then the following alternatives may be worth considering:
1. If time allows, refine or coarsen the flow grid (whatever is the best) to achieve a better separation of the length scales. Then repeat the upscaling of the outer bounds (in Step 2) and rerun the simulations (or recompute the difference parameter) to check if the performance gap (upscaling uncertainty) has narrowed.
2. Apply a tensor method with semi-open or periodic boundary conditions. As previously mentioned, this should result in an intermediate estimate of the effective permeability tensor and hence provide a flow model performance that lies somewhere in between the two methods in Step 2.
3. Convert the two local tensor methods in Step 2 into regional methods using a buffer region of modest size. Rerun the finite-volume simulator (or recompute the difference parameter) and check the performance gap. If necessary, repeat this step using an increasingly larger buffer region until the upscaling uncertainty reaches acceptable limits, at least as far as practically possible.
4. If for some reason none of the previous alternatives are an option, then one needs to apply the method that best satisfies the wanted flow behavior of a given cell or cells in a given region. In other words, a combination of the outer bound techniques within one and the same model may be a fourth alternative.

The multistep procedure of the relative upscaling approach is graphically illustrated in Fig. 17.30.

Fig. 17.30 – Flow chart for the relative upscaling approach.

Streamline Simulation—Rod P. Batycky and Marco R. Thiele

Introduction

Streamline-based flow simulation differentiates itself from cell-based simulation techniques such as finite-differences and finite-elements in that phase saturations and components are transported along a flow-based grid defined by streamlines (or streamtubes) rather than moved from cell-to-cell. This difference allows streamlines to be extremely efficient in solving large, heterogeneous models if key assumptions in the formulation are met by the physical system being simulated (see below). Specifically, large relates to the number of active grid cells.

Streamlines represent a snapshot of the instantaneous flow field and thereby produce data such as drainage/irrigation regions associated with producing/injecting wells and flow rate allocation between injector/producer pairs that are not easily determined by other simulation techniques. The computational speed and novel solution data available have made streamlines an important, complementary approach to traditional simulation approaches to perform sensitivity runs, quantify the impact of upscaling algorithms used to move models from the geomodeling scale to the simulation scale, visualize the flow field, perform more reliable full-field simulations where sector models would normally be used, enable the ranking of predicted field behavior of given multiple production scenarios and input parameters, evaluate the efficiency of injectors and producers, reduce turnaround time in history matching, and perform other established reservoir engineering tasks. A comprehensive overview on streamline-based flow simulation has recently been presented by Thiele.^[82]

Applicability of Streamline Simulation

The power of streamline simulation lies in its simplicity. The main objective is to capture how injected reservoir volumes (usually water and/or gas) displace resident reservoir volumes given well locations, well rates, reservoir geometry, and geological description. One of the key underlying assumptions in streamline simulation is that the system be close to incompressibility. This decouples saturations from the underlying pressure field and allows each streamline to be treated as being independent from the streamlines next to it.

Many fields under waterflooding or other pressure maintenance schemes are excellent candidates for streamline modeling and have been successfully modeled in this way.^[83]^[84] Forecast simulations under the assumption of voidage replacement are another good example where streamlines can be very effective. Even miscible gas injection schemes have been successfully modeled.^[85]^[86] At high pressures, the displacement of resident oil by gas is primarily an issue of simulating local sweep efficiency and channeling, something streamlines are designed to model without incurring numerical difficulties associated with other formulations.^[87] An early application to streamdrive projects was presented by Emanuel.^[88] Crane and Blunt^[89] used streamlines to model solute transport. More recently, streamlines have been shown to be very effective in modeling fractured reservoirs using a dual porosity formulation.^[90]

Historical

Streamlines have been in the petroleum literature as early as Muskat and Wyckoff’s 1934 paper.^[91] In 1937, Muskat presented the governing analytical solutions for the stream function and the potential function for 2D domains using the assumption of incompressible flow.^[92] Since then, streamlines and streamtubes have received repeated attention as a way to numerically predict the movement of fluids, even after the advent of finite-difference methods in the early 1960s. Important early contributions were made by Fay and Pratts,^[93] Higgins and Leighton,^[94]^[95] Bommer and Schechter,^[96] Martin and Wegner,^[97] Lake et al.,^[98] and Emanuel et al.^[99]

In the early 1990s, streamlines were revived because advances in geological modeling techniques were producing models that were too large for finite differences to simulate in an acceptable time frame. For streamlines to be applicable to real field cases, important advances were made that extended streamlines to 3D using a time-of-flight variable^[87]^[100]^[101]^[102]^[103] allowed for streamlines to be periodically updated^[87]^[104]^[105] and included gravity.^[106]

Mathematics of the Streamline Method

The streamline method and the underlying mathematics for incompressible multiphase flow are briefly outlined here. For a detailed discussion as well as additional references describing streamline methods, see Batycky et al.,^[87] Batycky,^[103] and Blunt et al.^[107]

Governing IMPES Equations. The streamline method is an IMPES-type formulation with the pressure field solved for implicitly and the oil/gas/water saturations solved for explicitly along streamlines. The governing equation for pressure, P, for multiphase incompressible flow without capillary or diffusion effects is given by

....................(17.19)

where D is the depth below datum, g is gravitational acceleration constant, k is the permeability tensor, k_rj is the relative permeability, μ_j is viscosity, and ρ_j is the phase density of phase j. The total velocity,

, is derived from the 3D solution to the pressure equation and application of Darcy’s law. The explicit material-balance equation for each incompressible phase j is then given by

....................(17.20)

Each phase fractional flow, f_j, is given by

....................(17.21)

and the phase velocity resulting from gravity effects because of phase density differences is given by

....................(17.22)

The difference between finite-difference simulation and streamline simulation is the way the explicit material balance equations (Eqs. 17.20 through 17.22) are solved. In finite difference, the material balance equations are solved between gridblocks, whereas in streamline simulation the material balance equations are solved along streamlines. How this is done is explained next.

Solution to the Transport Equation. In a standard finite-difference method, Eq. 17.20 is discretized and solved on the underlying grid on which the pressure field is computed. The solution to Eq. 17.20 is governed by the grid CFL condition, which can lead to prohibitively small timestep sizes, particularly for models with high permeability contrasts and/or high local flow velocities. With streamlines, this grid CFL limit is avoided completely by solving Eq. 17.20 along each streamline using a time-of-flight (TOF) coordinate transform.^[106]

Streamlines are traced from sources to sinks based on the underlying total velocity field. As each streamline is traced, compute the TOF along the streamline, which is defined as

....................(17.23)

and leads to the definition

....................(17.24)

Using Eq. 17.19, rewrite Eq. 17.20 as

....................(17.25)

Because the gravity term is not aligned along a streamline direction, Eq. 17.25 is split into two parts (operator splitting), giving two 1D equations. The convective portion of the material-balance equation along streamlines is given by

....................(17.26)

while the portion resulting from phase-density differences solved along gravity lines is given by

....................(17.27)

Both Eqs. 17.26 and 17.27 represent 1D equations that are solved using standard finite-difference numerical techniques. There are still CFL limits that restrict timestep sizes in these equations, but these are local to each streamline or gravity line, rather then at the 3D grid level.

Timestepping. In field-scale displacements, the streamline paths change with time because of the changing fluid distributions and the changing well conditions. As a result, the total velocity field is periodically updated, and new streamlines are recomputed to reflect the nonlinear nature of the displacement.

To move the 3D saturation distribution forward in time between successive streamline distributions from time T_i to T_i+1 = T_i+ dT_i, the algorithm pictured in Fig. 17.31 is used.

Fig. 17.31 – Flow chart showing the basic algorithm for streamline-based flow simulation.

The basic algorithm for streamline-based flow simulation is as follows: (1) Given initial conditions (i.e., pressures and saturations for each active cell in the system) and well conditions, the pressure is solved implicitly for each cell, as is done in conventional finite-difference methods (Eq. 17.19). (2) With the pressures known, the total velocity for each cell interface can be determined using Darcy’s Law. The total velocity is then used to trace streamlines using Pollock’s algorithm.^[100] (3) 1D mass conservation equations are then solved along each streamline, independently of each other (Eq. 17.26). The initial conditions for the streamlines are obtained by a mapping from the underlying 3D grid onto each streamline. The mass-transport problem is marched forward in time along each streamline for a pre-specified global timestep dT_i, and then the solution is mapped back onto the 3D grid. Gravity is included by considering a vertical segregation step along gravity lines after movement along all streamlines (Eq. 17.27). While simple in its approach, important details must be considered. In particular:

The algorithm is similar to an IMPES approach, in that the pressure is solved implicitly for a new time level n+1 assuming saturations at level n. The saturations at time n are given by mapping back solutions from each streamline onto the 3D grid at the previous timestep. Because of the implicit nature of the pressure solution, there is no limitation on the timestep to reach n+1. However, for compressible systems numerical convergence problems might limit the actual size of the timestep. This is no different than in FD simulation.
The tracing of the streamlines using Pollock’s algorithm assumes Cartesian cells. Nonorthogonal corner-point cells require an isoparametric transformation for tracing streamlines.^[108]
For incompressible systems, streamline will start at injection wells and end at production wells. For compressible systems, streamline can start/end anywhere in the system, because any gridblock in the system might act as a source (volume expansion) or a sink (volume contraction). Multiphase gravity effects can give rise to circulation cells for both incompressible and compressible systems.
Initial launching of streamlines from wells can be proportional to the total flux at the wells, though this will in general leave many cells in the system without a streamline passing through them. For missed cells, tracing begins at the center of the missed cell and then traced backward until a source is encountered. If a cell does not have a streamline pass through it, then it is not possible to assign an updated saturation back to that cell.
In practice, it is not possible to have all streamlines carry the same flux and ensure at least one streamline per cell. Thus, streamlines do not carry the same flux. Furthermore, for incompressible problems the flux along each streamline is a constant, while for compressible systems it is not.
The tracing of streamlines using the TOF variable produces a highly irregular 1D grid along each streamline. To numerically solve the 1D problem efficiently, the 1D grid must be regularized, solved using an implicit approach, or regridded in some way to allow for a more efficient solution.
The tracing of the streamlines relies on an accurate solution of the velocity field. Excessive distortions of the grid (nonorthogonal) or a pressure solution that has not been solved to a small enough tolerance can cause problems in tracing streamline paths.

The Computational Efficiency of Streamlines

One advantage of streamline simulation over more traditional approaches is its inherent efficiency, both in terms of memory and computational speed. Specifically, streamline-based simulation can exhibit a near-linear scaling in run times as a function of active cells in the model. Memory efficiency is a result of two key aspects of the formula: streamline-based simulation is an IMPES-type method and therefore involves only the implicit solution of pressure, and tracing of streamlines and solution of the relevant transport problem along each streamline is done sequentially. Only one streamline needs to be kept in memory at any given time.

Computational speed, on the other hand, is achieved because the transport problem is decoupled from the 3D grid and instead solved along each streamline. Because transport along streamlines is 1D, they can be solved efficiently. Because the number of streamlines increases linearly with the number of active cells, and streamlines only need to be updated infrequently, the computational time exhibits a near-linear scaling with increasing number of gridblocks

The number of global timesteps is related to how often the flow field (streamlines) requires updating. Specifically, changing flow paths are a function of heterogeneity, mobility changes, gravity, and changing well conditions. For many practical problems, it is the changing well rates that introduced the greatest impact on a changing flow field and is therefore the limiting factor in deciding on global timestep sizes. Grouping well events into semiyearly or yearly intervals and assuming that the streamlines remain unchanged over each period is reasonable. This is why field simulations with 30- to 40-year histories are successfully and routinely simulated with 1-year timesteps.^[84]

A good example to demonstrate the efficiency of SL simulation is Model 2 of the 10th SPE Comparative solution project.^[12] The total run time, T, of any streamline simulation is approximately proportional to

....................(17.28)

A near-linear scaling arises because:

The number of timesteps (streamline updates) is independent of the model size, heterogeneity, and any other geometrical description of the 3D model. It is mainly a function of the number of well events and the actual displacement physics. For the SPE10 problem in Fig. 17.32, all cases were run with the exact same number of streamline updates—24.
An efficient pressure solver is expected to have a near-linear behavior as well.^[109]
The number of streamlines tends to increase linearly with the number of gridblocks, all else being equal. Fig. 17.32 illustrates this behavior.
The time to solve the transport problem along each streamline can be made efficient by regularizing the underlying TOF grid and choosing the number of nodes to use along each streamline regardless of the size of the underlying 3D grid.

The linear behavior with model size is the main reason why streamline simulation is so useful in modeling large systems. In FDs, finer models not only cause smaller timesteps because of smaller gridblocks but usually face problems because of increased heterogeneity as finer models tend to have wider permeability and porosity distributions. The usual workaround is to use an implicit or adaptive-implicit formulation, but for large problems these solutions can become prohibitively expensive, both in terms of CPU time and memory.

Fig. 17.32 – Example of linear scaling of run time and number of streamlines as a function of active cells for SPE Comparative Solution Project 10 using 3DSL, a commercial streamline simulator, on PIII 866MHz PC.
Fig. 17.33 – Streamlines automatically allow the determination of the allocation of flow between wells by summing the flux of all streamlines associated with a particular well, well pair, or group of wells. Using this information and the visual display of streamlines allows patterns to be balanced correctly and efficiently. From left to right: rates are progressively changed to yield a balanced pattern.

Novel Data Produced by Streamlines

Streamlines produce new data not available with conventional simulators. Because streamlines start at a source and end in a sink, it is possible to determine which injectors are (or which part of an aquifer is) supporting a particular producer, and exactly by how much. A high water cut in a producing well can therefore be traced back to specific injection wells or boundaries with water influx. Conversely, it is possible to determine just how much volume from a particular injection well is contributing to the producers it is supporting—particularly valuable information when trying to balance patterns ( Fig. 17.33) or optimize water injection over a field.

Streamlines can also identify the reservoir volume associated with any well in the system, because a block traversed by a streamline attached to a particular well will belong to that well’s drainage volume. It is therefore possible to divide the reservoir into dynamically defined drainage zones attached to wells (Fig. 17.33). Properties normally associated with reservoir volumes can now be expressed on a per-well basis, such as oil in place, water in place, and average pressure, just to mention a few.

The most successful uses of new data produced by streamlines are in the area of waterflood management and reservoir surveillance,^[110]^[111] and in the area of history matching.^[112]^[113]^[114]^[115]

Applications of Streamlines

Streamlines are a powerful complementary tool to more traditional simulation techniques, and they are expected to play an important part in optimizing field production and management in the future. Specifically, streamlines can be used to:

Validate upscaling techniques by allowing to generate reference solutions of fine-scale models.^[81]
Efficiently perform parametric studies.
Visualize flow.
Balance patterns.
Determine efficiency of injectors and producers using data provided by streamlines.
Aid in history matching.
Enable ranking of production scenarios/geological models.
Optimize and manage field injection/production.
Conduct reservoir surveillance.

It is important to underline that the theory on which streamline simulation rests is firmly rooted in the incompressible formulation of exact voidage replacement. Thus, streamline simulation is particularly powerful for modeling systems that are not a strong function of absolute pressure, but are instead governed by a pressure gradient. In addition, the strong assumptions of independence between streamlines favors modeling displacements that are not a strong function of diffusive phenomena, such as capillary pressure, transverse diffusion, or compressibility. For example, streamline simulation offers little or no advantage over conventional simulation for modeling primary production. This is because the main feature of modeling primary production is to accurately capture the pressure decline over time, not the movement of a saturation front.

The Future of Streamline Simulation

The next few years are expected to bring a further maturing and extended application of streamline-based flow-simulation technology. It is reasonable to expect that most companies using conventional simulation technology today will in one form or another use SL simulation in their future work. What remains uncertain is whether new user groups, such as geologists and geophysicists, will adopt the technology to bring a dynamic flow component to their analysis. Developments in the following areas are currently under way: the use of streamlines in conditioning static reservoir models to production data, extension of streamline simulation to compositional models, tracing of streamlines in structurally complex reservoirs, modeling of dual-porosity/dual-permeability models, and parallelization of streamline numerics for the solution of large models.

Simulation of Geomechanics—Tony Settari

Historically, much of the simulation has accounted for rock mechanics by simple use of a time-invariant rock compressibility c_R , spatially constant or variable. In reality, rock mechanics is intimately coupled with fluid flow in two aspects. First, the porosity changes are a direct result of the deformation of the skeleton, which is a complex function of both pressure and stress (effective stress state). Second, it is well known that the permeability of the media is also a function of effective stresses in the reservoir. Therefore, rigorous reservoir simulation should include simultaneous solution of multiphase flow and stresses as well as the appropriate dependencies between these processes. While these couplings physically exist to some extent in all reservoirs, they can be often ignored or approximated when the reservoir behaves elastically. However, the changes in porosity and permeability are more pronounced when rock failure occurs, such as in compacting reservoirs or in high-pressure injection operations, and these processes require use of more complex, coupled geomechanical modeling.

Types of Coupling

There are essentially two main types of coupling between reservoir flow and stress: pore-volume coupling and flow-properties coupling. In the past, the first led to the development of various simplified "compaction" modeling techniques, while the latter is reflected in the "pressure-dependent permeability" options available in many simulators. In a coupled geomechanical model, both can be treated more rigorously.

Pore-Volume Coupling. The porosity in the reservoir model is traditionally treated as a function of pressure via rock compressibility:

....................(17.29)

and changes in block pore volume V_p are computed as V_p = V_b⁰Ф, where V_b⁰ is the block (bulk) volume. In reality, pore volume changes are a result of complex interaction of fluid (pore) pressure, stresses acting on the element of the rock and temperature. The deformation of the rock solid (also called the skeleton) caused by combination of stress and pressure changes results in changes in the bulk volume of an element V_b, which is computed at any conditions as

....................(17.30)

where ε_v is the volumetric strain at these conditions, and V_b⁰ is the bulk volume at reference conditions at which the volumetric strain = 0. Then the true porosity is given by

....................(17.31)

where V_p is the pore volume of the element. In Eq. 17.31, both pore volume and bulk volume are variable; therefore, true porosity and pore volumes are both functions of pressure, temperature and stress:

....................(17.32)

which shows the coupling between fluid flow and geomechanics (stress modeling).

In stress modeling, the changes in volumetric strain and porosity are calculated from complex constitutive relations of the material, which define both the stress-strain and volumetric behavior.^[116]^[117] To compute pressure changes correctly in the reservoir simulator, it is necessary to force the changes in pore volume to be the same as computed by the stress model, which is the essence of the "volume coupling." This can be achieved in two ways. The rigorous solution is for the reservoir model to recalculate the block sizes based on the stress solution, and use the true porosity. However, reservoir simulators do not allow for modifying the bulk volume. In this case, one can define a pseudoporosity Ф* = V_p / V_b⁰ , which will give the correct pore volume.^[118] In either case, the usual treatment of porosity by rock compressibility must be replaced by the coupling with stress-strain solution. The porosity changes become more complex when failure of the skeleton is reached either in compression (rock compaction) or in shear.

Flow-Properties Coupling. The primary mechanism is the dependence of permeability on stress, usually of the form

....................(17.33)

This relationship is well established for many types of reservoirs^[119] and generally becomes stronger as permeability decreases. It is dominant (and more complex) in fractured reservoirs where stress-dependent fracture aperture and reopening/creating fractures under injection can cause large, anisotropic changes. The tensor character of permeability may be important in these applications. In soft formations and unconsolidated sands, deformations can lead to porosity dilation, which will also lead to permeability enhancement. On reloading, there is a hysteretic effect. Two to three orders of magnitude enhancements in permeability can occur because of injection at fracture pressure (e.g., in microfractured rock^[120] or coal seams).

The permeability changes are a function of some measure of effective stress, but are often simplified (and laboratory data reported) as a function of pressure. The creation of new fractures can cause a transition from initially single-porosity media to dual-porosity, and in such a case will also have an effect on relative permeabilities.

Modeling of Reservoir Compaction and/or Dilation

Modeling reservoir deformation is of considerable importance in soft and/or thick reservoirs where the results of compaction may provide an important production mechanism, cause well failures, and/or cause ground subsidence or heave with environmental consequences. Review of the compaction mechanics and its consequences for field development is found in Settari^[121]. Initial approach to modeling compaction was based on modifications of reservoir models.^[122]^[123]^[124]^[125] The common feature of such reservoir compaction models is that the compaction is treated as a 1D problem (uniaxial strain) by assuming that (a) only vertical deformations take place, and (b) each vertical column of blocks deforms independently. In such models, the porosity changes are calculated by modifying the conventional compressibility c_R as a function of pressure only, in the form of "compaction tables." The tables are based on results of uniaxial strain laboratory experiments, and the stress problem is not solved. The compaction of the reservoir is then obtained analytically assuming uniaxial deformation. The relation between reservoir compaction and surface subsidence can be then obtained by an independent solution of a stress problem using the computed compaction as a boundary condition. Dilation (increase of porosity) is an important geomechanical mechanism occurring during steam injection into unconsolidated sands. This process has been also modeled by the "compaction-dilation" tables.^[126] In chalk reservoirs, increased water saturation in waterflooding causes weakening of the rock and therefore S_w is an additional variable in Eq. 17.32.^[127]

The major drawback of the use of compaction tables is that the dependence on stress indicated in Eq. 17.32 must be either ignored, or the change of stress must be estimated in terms of change in pressure. A more accurate modeling approach is to couple in some fashion the reservoir simulator with stress-strain (geomechanics) solution. Such models typically combine the solution of the multiphase flow in the reservoir and elastoplastic solution of the deformations in a much larger domain including the reservoir, sideburden, underburden, and overburden. The majority of coupled models use the iterative coupling^[128]; the different variations and their shortcomings are discussed next.

Coupled models have much larger computing requirements compared with a reservoir model of the same reservoir, primarily for two reasons: first, the larger number of unknowns per gridblock, and second, the stress-solution grid must be usually much larger laterally than the reservoir grid to eliminate the effect of boundary conditions as well as extend up to surface to provide subsidence solution. Combination of these factors leads to computing times typically of one order of magnitude larger compared with conventional simulation, and even more if elastoplastic solution is required for the stresses.

Modeling of Stress-Dependent Flow Properties

The primary flow-dependent property is permeability, and the problems to model its dependency in Eq. 17.33 are similar to modeling the pore-volume coupling (Eq. 17.32). However, the problem is somewhat easier because stress-dependent permeability (or transmissibility) does not affect mass-balance formulation.

Again, the traditional approach is to use tables of k vs. pressure in an uncoupled model. However, the problem remains one of replacing the dependency on effective stress by one on pressure only. Even in a single-phase, single-porosity gas-flow case, different assumptions about the stress change during depletion can lead to large errors in well decline.^[129] Different strategies for converting the stress-dependent data to pressure tables are based on local constrainment assumptions.^[130]

In coupled models, the permeability dependency can be usually computed explicitly on a timestep basis, and "loose" coupling can be used. In fact, a "coupled" model that deals only with flow properties coupling and ignores the pore-volume coupling can be run successfully even if the stress solution is done in larger intervals of time compared to the reservoir. Such models have been used extensively to study permeability changes in waterfloods, particularly in fractured or jointed media.^[131]^[132] Here, the advantage of coupled modeling is in its capability to predict the permeability changes from the geomechanics of reopening of fractures or failure (dilation) of joints.^[133] The development of anisotropy is dictated by the orientation of fractures or faults, and requires a "full tensor" treatment of transmissibilities in the flow model. In the stress strain model, different methods of pseudoizing the fracture/joint networks into a continuum are used, which include predicting permeability as a function of effective stress and/or strain. While the need for the tensor transmissibilities in such models has been recognized,^[134] in injection processes dual-porosity media can be created. Therefore, reservoir description may be changing in time because of geomechanics; this aspect has been ignored in coupled models to date.

The same principles can be also applied to model hydraulically induced fractures being represented by dynamically changing transmissibility multipliers in the potential fracture plane.^[135] The effective stress dependency (as opposed to pressure in an uncoupled model) allows capturing the changes of the fracture propagation pressure with time, which can be large, in particular in steam injection. Another application is the prediction of production/injection-induced slippage on faults, which can induce communication between reservoir fault blocks and/or seismicity.

Types of Coupled Models

Coupled models can be either fully coupled (i.e., all unknowns solved simultaneously) or modular (reservoir simulator and stress code). In the latter case, different coupling strategies can be used, with consequences for running speed and accuracy. The majority of coupled models use a conventional finite-difference reservoir simulator coupled with a finite-element (FEM) stress simulator. However, attempts have been made to develop fully coupled FEM codes,^[136] and a fully coupled geomechanics was implemented in a commercial model using a finite-difference stress solution.^[137] Considering the proliferation and sophistication of the geomechanics codes available outside the petroleum engineering, the modular approach offers the best solution.^[119]^[127] Generally, the reservoir simulator is the "host" or "master." Commercial stress simulators are, in principle, easy to couple to it (in particular if only permeability coupling is considered).

Because of the extreme complexity and large computing requirements of coupled models, different simplifications have been developed. The main types of models (in the order of increasing rigorousness, but also computing time) are as follows:

One-Way Coupling. Pressure and temperature changes are passed from the reservoir code to the geomechanics module, but no information is passed back on timestep basis. The geomechanics does not improve the flow solution, but the model can be useful for predictions of wellbore stability for infill drilling, fracturing pressures, and so on. Manual adjustment of "compaction tables" is possible manually through restarts. Such manually coupled solutions^[138] of the stress problem (at intervals of time) or one-way coupling^[139] were often used in early coupled modeling. The method can be satisfactory when the reservoir fluid system is highly compressible (i.e., in gas reservoirs), but can lead to errors when the porosity is strongly coupled to flow.

Loose Coupling. Reservoir and geomechanics modules are run sequentially on a timestep basis, passing converged solutions of flow and stress variables to each other. Pore volumes and permeabilities in the flow model are computed as a function of p, T, and σ_avg with stress variables lagged a timestep. However, the relationships are "distilled" into tables similar to the "compaction tables," but now a function of effective stress. The advantages are functional similarity to the uncoupled reservoirs with "compaction tables," no need for iteration during timestep, and the possibility of updating the stress solution less frequently than the reservoir solution. However, complex constitutive models of the solid (e.g., plasticity) may be difficult to represent.

Iterative Coupling. This method is shown schematically in Fig. 17.34. Iteration is carried out between the reservoir and stress solution at every timestep until the pore volumes and permeabilities calculated from the stress model and those used by the reservoir model agree. In each iteration, the previous guess of the V_pⁿ⁺¹ for the end of the timestep is used to converge the flow solution, and the changes of p and T over timestep are then used to solve for the new deformations and stresses, which in turn provide updated estimate of V_pⁿ⁺¹. The changes of permeability are also iterated on. Therefore, each "geomechanical" iteration costs the equivalent of a timestep solution of the previous methods. The original formulation of the coupling iteration^[128] is always convergent, and its efficiency has been recently improved.^[140]

Fig. 17.34 – Schematic of the iterative coupling algorithm.

The iterative coupling, when converged, is equivalent to a fully coupled code while it retains flexibility. In many problems, it is not necessary to fully converge the timestep, and if used with only one iteration per timestep, the method is similar to the loose coupling, except that the porosity is determined directly by the constitutive model of the solid rather than by a simplified relationship.

Full Coupling. This requires simultaneous formulation of the flow and stress variables and therefore results in larger matrices. The advantage is that consistent approach to discretization can be used, and the model is integrated from the point of view of code development. However, it is very costly to redevelop all the features of the physics and numerics now readily available in stress codes. Moreover, in published fully coupled models, the approach for solving the resulting matrix problem has been to partition the matrix in the same fashion as in the iterative coupling^[141] and to apply the geomechanical iteration at the matrix-solution level. Therefore, the fully coupled formulation, which results in larger, strongly nonlinear matrix equations, does not reduce the difficulty of the problem, and it may need to use geomechanical iteration in the solution process as the best strategy. These aspects need further study.

Future Trends and Needs

As a result of much larger computing requirements, coupled models lag behind conventional models in the size of the problems that can be currently handled. Therefore, they are a prime candidate for the use of massively parallel hardware and will require large future development effort in parallelization. Because of the increased understanding of the complexities of the geomechanics, the current trend is also toward more strongly coupled models with fewer simplifications. This further increases the computing requirements.

Given that not all problems require use of geomechanics, and the cost of the study may increase dramatically, it is important to be able to identify when coupled simulation is needed, and what approximations can be made without compromising the answers. There are no simple rules, but there is a growing need to conduct a "screening" process at an early stage of a reservoir study to determine if geomechanics is an issue. This process requires an integration of reservoir, production, and completion engineering data as well as field experience.

Finally, coupled geomechanical modeling is the future tool for truly integrated reservoir management. Conventional reservoir simulation studies ignore numerous constraints placed on the development scenarios from the point of view of drilling, completion, and operations. These constraints can be incorporated into coupled models, and additional modules can be integrated (e.g., long-term wellbore stability and sand production predictions, subsidence management, 4D seismic interpretation, and so on^[142]).

Pressure/Volume/Temperature (PVT) Treatment—Curtis H. Whitson

The PVT treatment of fluids in reservoir simulation describes the phase behavior of gas, oil, and water at reservoir and surface conditions. Phase behavior of a mixture with known composition consists of defining the number of phases, phase amounts, phase compositions, phase properties (molecular weight, density, and viscosity), and the interfacial tension between phases. In addition to defining the phase behavior of mixtures at a specific reservoir pressure, knowing the derivatives of all phase properties with respect to pressure and composition is important in reservoir simulation.

The calculation of phase behavior in a reservoir model can be made in one of two ways—using a "black-oil" approach^[22]^[143]^[144]^[145]^[146] based on simple interpolation of PVT properties as a function of pressure, or using a "compositional" approach based on a thermodynamically-consistent model such as a cubic EOS.^[14]^[15]^[22] With either approach, the PVT quantities required by a reservoir simulator are essentially the same. Modern reservoir simulators^[14]^[22] are usually written with a general compositional formulation, whereas black-oil PVT properties are converted internally to a two-component "compositional" model; the two components are surface gas and surface oil.

A reservoir simulator keeps track of overall composition in each computational grid cell as a function of time. A grid cell’s overall composition changes because of cell-to-cell fluxes and because of production or injection in wells. The phase fluxes and component movement within the reservoir are greatly affected by phase behavior (e.g., the mobility of each phase and which components are carried in each phase). The surface products from wells are dependent on the phase behavior of produced wellstreams, but at conditions far removed from reservoir pressure and temperature.

In most reservoir simulators, the water phase and water component are treated in a simplified manner—namely, that the water component does not partition into the hydrocarbon phases, and the hydrocarbon components do not partition into the water phase; the term "hydrocarbon" also includes N₂, CO₂, and H₂S. Because of relatively high CO₂ solubility in water, and the potential importance of CO₂ "accounting" in CO₂ floods, some compositional models allow CO₂ partitioning in the water phase.^[147]

Conceptually and computationally, it is feasible to allow complete partitioning of all components in all phases in reservoir simulation. The partitioning related to the water phase and the water component could be treated with simple pressure-dependent tables, or even with an EOS model. The water-related K-values have simple composition dependence that would make EOS-based fugacity updates almost trivial. The main problem would be treating the impact of changing salinity on water-related K-values. From a practical point of view, however, modeling water-related component partitioning will have a marginal impact on reservoir performance. Hereafter, only phase behavior of nonaqueous phases will be discussed.

Number of Phases and Phase Type

The starting point for PVT calculations in a model grid cell is to determine if the overall composition is single-phase or two-phase at the current pressure estimate. If a cell is single-phase, the phase "type" (gas or oil) may also be needed to select the proper relative permeability and capillary pressure curves. For a black-oil PVT model (i.e., any model with simple pressure-dependent K -values), determining the number of phases and phase identification is trivial.

For EOS models, two methods can be used for establishing how many phases exist: the Michelsen^[148] phase stability test or a saturation pressure calculation.^[148] The stability test is relatively "slow" because good K-value estimates are not used to initiate the calculation, whereas a saturation pressure can be quite fast because good K-value estimates are often available. A stability test is more reliable in the sense of giving a conclusive result, whereas a near-critical saturation pressure test may not be found^[149] or converge correctly to a nontrivial solution. Another advantage of using saturation pressure is that it gives a consistent method to define phase type for single-phase cells, something that is not provided by the stability test. The choice of which method to use depends on the tradeoff between speed and reliability. Both methods, if detecting two phases, give excellent starting K-values to initiate the two-phase flash calculation.

Two-Phase Flash

Having established that two phases exist in a cell, one must perform a flash calculation. The flash calculation uses the overall moles n and molar composition z_i to determine the molar amounts of each phase (n_g and n_o), and the phase compositions (y_i and x_i). For simple pressure-dependent K-value models, the Rachford-Rice^[150]^[151] (RR) procedure is used to find molar amounts and compositions,

....................(17.34)

solving for α_g = n_g / (n_g + n_o), with x_i = z_i / [1 + α_g (K_i – 1)] and y_i = K_ix_i.^[152]

For an EOS model, using the RR solution can be used, but an additional constraint must be satisfied; the fugacity of each component must be the same in each phase,^[153] f_gi = f_oi. Using an initial set of K_i estimates, the RR equation is solved. Resulting compositions y and x are then used to calculate fugacities f_gi and f_oi. If the constraint f_gi = f_oi is not satisfied for all components, a new set of K-values is estimated using some form of a successive substitution (SS) update K_i new = K_iold (f_oi / f_gi).

For reservoir simulation, a Newton solution of the flash equation should be much faster than SS methods. Including the equal-fugacity constraint within the set of nonlinear equations used to solve for model pressures "automatically" provides the first Newton update of the flash equation.^[153] Usually only one additional Newton iteration of the flash equation is needed to converge the fugacity constraints to required tolerance.

In summary, K-values alone determine the phase amounts and phase compositions (from the Rachford-Rice equation). The EOS model guarantees rigorous thermodynamic consistency through the equal-fugacity constraint, ensuring that K-values properly account for pressure and composition dependence of the phase split and component partitioning.

Having completed the flash calculation in a cell, the following information is known: phase moles, phase masses, phase densities, phase volumes (saturations), and phase compositions. The viscosity and gas/oil IFT can be calculated using compositional correlations,^[152] or interpolated from pressure-dependent tables for a black-oil model. The impact of PVT properties on reservoir performance will consider individual PVT properties and their impact on individual terms in the flow equations.

Density

Phase molar amounts, n_g and n_o, are converted to phase volumes using phase molecular weights and densities: V_g = n_gM_g / ρ_g and V_o = n_oM_o / ρ_o, where S_g = V_g / V_block and S_o = V_o / V_block. Saturations determine relative permeabilities, which can have a dramatic impact on phase fluxes. Typically, k_r ~ Sⁿ , where n ~ 2 to 4. For n = 3, a 5% error in density results in a 15% error in k_r.

For reservoirs with sufficiently high vertical permeability, gravity will often play an important role in recovery. Gravity segregation is caused by vertical fluxes driven by the potential terms ∇ρ_g = d(p_g + ρ_ggz + P_cgo) / dz, ∇ρ_o = d(p_o + ρ_ogz) / dz, and ∇ρ_w = d(p_w + ρ_wgz + P_cwo) / dz . Errors in densities have a direct impact on gravity segregation. It is actually the density difference between flowing phases that determines the magnitude of gravity segregation (i.e., Δρ_wo = ρ_w – ρ_o , Δρ_go = ρ_o – ρ_g , and Δρ_wg = ρ_w – ρ_g). Likewise, the condition for equilibrium initialization ∇ρ_g = ∇ρ_o = ∇ρ_w = 0 results in saturation-depth distributions, which must honor input Δρ_gog = P_cgo(S_g) and Δρ_wog = P_cwo(S_w) relationships.

A +5% error in oil density (e.g., 630 instead of 600 kg/m³) and a –5% in gas density (e.g., 380 instead of 400 kg/m³) results in gravity segregation potential (Δρ_go = ρ_o − ρ_g) in error by 25% (250 instead of 200 kg/m³). Similarly, in the initialization of a water/oil system, the inadvertent use of ρ_w = 62.5 lbm/ft³ instead of the correct value of 70.0 lbm/ft³ for a high-salinity brine could easily cause a 25% error in initial oil in place for a low-permeability (large capillary transition zone), low-API oil.

In black-oil models, it is particularly easy to model densities incorrectly because of the inadvertent use of "arbitrary" surface densities. The reservoir phase densities are not input directly in a black-oil model. Instead, they are calculated using the formation volume factors (FVF) B_o, B_g, and B_w, solution gas/oil and oil/gas ratios R_s and r_s, and surface densities ρ_gs, ρ_os, and ρ_ws, based on the relations ρ_g = (ρ_gs + p_os r_s) / B_g, ρ_o = (ρ_os + R_s ρ_gs) / B_o, and ρ_w = ρ_ws / B_w. Because reservoir densities are not usually tabulated as output by reservoir simulators, the user may not know how "wrong" the densities might be.

A well-tuned EOS using volume shift factors^[154] should always predict reservoir and surface densities within 2 to 4%, and often within 1 to 2%. The author’s experience has found it generally unnecessary to use different volume shift factors for reservoir and surface calculations unless all densities are needed with accuracies of 1 to 2%. The use of either Peng-Robinson^[155] or Soave-Redlich-Kwong^[156] EOS with volume-shift factors provides densities as good or better than the more-complicated ZJRK^[157] EOS model; in this author’s opinion, the ZJRK no longer has a role in reservoir simulation, as supported by decreasing industrial use (almost to nonexistent) during the past 5 years.

Component Partitioning

The partitioning of components in gas and oil phases, as dictated by K-values, is important for properly describing vaporization and condensation. Vaporization is the process of "stock-tank oil" (STO) components (C₆₊) moving from the oil phase to the gas phase, while condensation is the process of intermediate (C_3–C₁₀₊) components moving from the gas phase to the oil phase. Retrograde condensation is particularly important in the depletion of gas condensate reservoirs, where in-situ retrograde condensation results in decreasing surface liquid yields.^[22]^[152]^[158]

For most gas-injection projects, vaporization will play an important role. Proper modeling of vaporization usually requires a compositional model that uses a detailed description of the heavier (C₇₊) components, typically three to five fractions.^[22]^[159]^[160] The lightest fractions (C₇ to C₁₂) will typically vaporize 80 to 100%, while intermediate heavies (C₁₃ to C₂₀) may vaporize in varying amounts from 30 to 90%, and the heaviest components C₂₀₊ may vaporize from 0 to 50%. The degree of vaporization depends primarily on the local pressure and composition of displacing gas. Cook et al.^[144] suggested a modification of the black-oil model that allows saturated properties to change as a function of how much injection gas has contacted a given cell, thereby allowing reasonable description of a vaporization-dominated process. Tang and Zick^[161] propose a limited three-component compositional model that allows accurate description of developed miscibility for use in models in which grid-related dispersion may result in underpredicted conditions of miscibility using an EOS-based simulator.

In some gas-injection processes, the injection gas may be enriched (1) during injection by adding NGLs (C₃ to C₅ components) and/or (2) within the reservoir by multiple contacts with reservoir oil, which gradually enriches the injection gas with intermediate components C₃ through C₁₀. For either type of enrichment, the vaporization process can become extremely efficient and result in a near-miscible or miscible displacement with near-100% recovery of the oil. When such a displacement results, it is often associated with the development of a complex process whereby the near-miscible front consists of upstream vaporization and downstream condensation—the "condensing/vaporizing" gas drive mechanism first described by Zick.^[162]

For an EOS model to properly describe complex phase behavior related to vaporization, condensation, and near-miscible mechanisms, special PVT experimental data should be measured and used to tune the EOS. For immiscible processes with significant vaporization, a single-cell multicontact PVT test is useful, quantifying the degree of vaporization in terms of oil volumetric stripping and gas compositional changes. For near-miscible or miscible processes, a swelling-type test which has the following features is strongly recommended: (1) five to seven mixtures of reservoir oil and injection gas are used, with two to three mixtures being bubblepoints and two to three mixtures being dewpoints; (2) a constant composition expansion is conducted on each mixture, where saturation pressure and the oil relative volumes are measured and reported; and (3) a single equilibrium "flash" somewhat near the critical point is used for obtaining a set of equilibrium gas and oil compositions (K-values).

Proper EOS tuning of the complex phase behavior measured in such an experiment is difficult, and requires the ability to match near-critical volumetric and compositional changes. Such predictions are almost never available a priori with an EOS model that has been tuned only to simple depletion data. Small but important modifications of the heavy-end properties and binary interaction parameters are usually required to obtain a satisfactory match of near-critical PVT data provided by this special "swelling" test. A tuned EOS model that is able to match all PVT data, including near-critical phase behavior, has a good chance of properly predicting true multicontact near-miscible/miscible behavior.^[152]^[162]

Viscosity

Gas viscosities are typically estimated by correlation^[152] within 5 to 10%, and they are almost never measured experimentally. Such accuracy is adequate for most applications, and gas viscosities seldom vary beyond the range of 0.02 to 0.03 cp. For gas-injection processes at high pressure or near-miscible displacements, gas viscosities can range from 0.03 to 0.1 cp. Gas-viscosity correlations are not usually accurate for this higher range of viscosity, and errors up to 20 to 30% may be expected; compositional viscosity correlations also have the same level of accuracy.

Oil viscosities are notoriously inaccurate based on correlations,^[152] at best being within 10 to 30%, but often in error by 50% or more. Oil viscosities should always be measured and used to tune a viscosity model. A minimum requirement would be measurement of stock-tank oil viscosities, and normally live-oil viscosities are available from a differential liberation test.

For gas condensates, oil viscosities are almost never measured. This may be a serious problem if condensate blockage has a significant impact on well deliverability.^[163] Condensate viscosities are difficult to measure because retrograde condensate volumes can be very small. The use of a separator condensate sample is recommended. Bring it to reservoir conditions (T and P) for density and viscosity measurements, then tune the viscosity model to these data. This approach is feasible for any gas condensate; it is not expensive, and it is better than relying on untuned compositional viscosity correlations.

Our experience has been that the compositional viscosity correlation by Lorenz, Bray, and Clark^[164] (LBC) is not predictive, and it is highly dependent on accurate density predictions. Tuning the heavy-component critical volumes and, sometimes, careful modification of the higher-order constants in the LBC equation, provide the required accuracy for gas and oil viscosities in most reservoir processes. Unfortunately, modification of the LBC correlation to match a limited number of viscosities (in a narrow range of pressure and composition typical of depletion experiments) can lead to unphysical viscosity predictions at conditions developed during a gas injection process.

The Pedersen et al.^[165] compositional viscosity correlation, though more complicated and CPU-intensive than the LBC correlation, has quite good predictive capabilities. The correlation is based on a corresponding states formulation in which methane viscosity is the reference compound.

Gas/Oil Interfacial Tension (IFT)

Gas/oil IFTs may be used in a reservoir simulator^[13] to modify input capillary pressures P_cgo , relative permeabilities k_rog and k_rg , and residual oil saturation S_org . As IFT decreases, P_cgo decreases proportionately. Relative permeabilities and residual oil saturation change only marginally until IFT reaches a fairly low value (e.g., 1 mN/m), whereas at lower IFTs, the relative permeability-related properties change approximately as log(IFT).

Practically, IFT impact on capillary pressure is limited to fractured reservoirs^[166] in which displacement efficiencies may be strongly linked to the balance of gravity and capillary forces on a relatively small distance scale (e.g., block heights of 1 to 10 m).

IFT impact on relative permeabilities and residual saturations may have an impact on some gas-injection processes, though near-miscible and miscible processes have minimal IFT dependence because they are usually dominated by the strong compositional effects (vaporization and condensation) that result in near-100% recoveries. In fact, displacements that are miscible should, by definition, be independent of the relative permeabilities and residual oil saturation.

IFT impact on near-well relative permeabilities in gas condensate reservoirs can, together with high velocities, result in large capillary numbers, which have the tendency to "straighten" relative permeabilities and improve flow performance. It has been shown^[167]^[168] in a number of recent publications that this effect can have an important impact on gas-condensate well deliverabilities. Neglecting IFT and velocity dependencies of relative permeability can lead to overly conservative prediction of well deliverability (i.e., overprediction of condensate blockage).

Black-Oil PVT Models

Black-oil PVT properties are generated in one of two ways. For low- to medium-GOR oils (< 150 Sm³ /Sm³ ), a traditional differential liberation experiment (DLE) is used, with corrections for separator flash to calculate oil FVF B_o and solution GOR R_s, as well as the gas FVF B_g.^[152] This approach assumes the reservoir gas contains unsubstantial amounts of condensate in solution, with solution oil/gas ratio r_s ~ 0.

The more common and general approach to generating black-oil PVT properties uses an EOS model to simulate a depletion-type PVT experiment (differential liberation, constant volume depletion, or constant composition expansion), with the equilibrium gas and equilibrium oil at each stage in the depletion being individually processed to surface conditions to provide the four black-oil properties B_o,R_s, B_g, and r_s.^[22]^[156] For highly volatile oils, the EOS method gives substantially different and improved black-oil properties compared with the traditional laboratory DLE/separator-corrected approach.

The conversion of black-oil PVT data (R_s, B_o, r_s, and B_g) to a compositional model^[144] uses K-values of surface gas and oil pseudo "components" K_gs = (R_s + C_os) / (1 + r_sC_os) / R_s and K_os = r_s(R_s + C_os) / (1 + r_sC_os), with C_os = (RT_sc / p_sc )(ρ_os / M_os). The reservoir-phase densities are calculated from ρ_g = (ρ_gs + ρ_osr_s) / B_g and ρ_o = (ρ_os + R_s ρ_gs) / B_o, while phase molecular weights are given by M_g= (M_gs +r_sM_os C_os)/(1 + r_s C_os) and M_o = (R_sM_gs + M_os C_os) / (R_s + C_os). Viscosities and gas/oil IFTs are interpolated directly from input tables.

Coats et al., ^[22] Coats et al.,^[143] and Fevang et al.^[145] have shown that black-oil models can be used for practically any type of reservoir produced by depletion or waterflooding, including reservoirs with large compositional gradients. Some issues require special treatment for complex fluids systems, including fluid initialization and the method for generating black-oil tables. In a nutshell, the recommended procedures^[145] are to generate the black-oil tables with an EOS model using the fluids with the highest saturation pressure (e.g., at the gas/oil contact) and to initialize with solution GOR (R_sand r_s) vs. depth—not saturation pressure vs. depth.

A common problem with black-oil models is the calculation of "negative compressibility," ^[23] meaning that a small pressure drop results in a reduction in total (gas + oil) volume. Another problem is physical extrapolation of saturated PVT properties to saturation pressures higher than given in the input table (e.g., caused by gas injection, gravity segregation in undersaturated reservoirs, or near-well behavior during rate reductions).

When can the black-oil PVT treatment not be used? Basically, for any gas-injection process with significant component partitioning that changes during the displacement.^[145]^[159] This would include processes with high-pressure vaporization using lean gas, condensation from enriched injection gas, and developed-miscibility processes such as the condensing/vaporizing mechanism. Surprisingly, a black-oil treatment is sometimes adequate even for complex gas injection problems, though it is not usually known a priori. To check the validity of a black-oil model in a gas injection project, the reservoir process should first be simulated with a compositional model, and preferably a relevant 3D model that captures all important flow characteristics.

Equation-of-State Models

The most common EOS used in reservoir simulation are the PR and the SRK models. Both models have two constants, a and b. Each constant must be calculated for each component based on component critical properties (T_c and p_c) and acentric factor (ω).

The PR EOS has two versions—the original 1976 version^[155] and the modified 1979 version^[169]; the latter uses a third-degree polynomial expression for the correction term to constant a. For some systems, the difference in equilibrium calculations for the two PR EOS models is significant.

The Peneloux^[154] volume shift factors should always be used with two-constant EOS models to ensure accurate oil and gas densities. The volume shift factors have no impact on calculated K-values or compositions, only densities. As mentioned earlier, the ZJRK^[157] EOS is outdated and was used before the volume-shift method was introduced in 1980, with complex correction functions to constants a and b to improve liquid density predictions.

Binary interaction parameters (BIPs) k_ij are important for adjusting predictions of equilibrium properties (K-values and compositions).^[170] These parameters represent a correction of the form (1 – k_ij) to the a_ia_j term in the quadratic mixing rule for EOS constant a. BIPs can be positive or negative; they are symmetric (k_ij = k_ji); they are usually ~0 for most hydrocarbon-hydrocarbon pairs, except C₁ to C₇₊ pairs which may reach values as high as 0.25; and they are generally close to 0.1 for nonhydrocarbon (N₂, CO₂, H₂S)-hydrocarbon pairs.

Three-Phase PVT Behavior

Three-phase (L1-L2-V) behavior is an occasional but serious problem for EOS-based compositional models. The third phase (L2) is usually a light liquid and typically appears at low temperatures (< 140°F) in gas-injection processes using CO₂ or NGL-enriched gas.^[171] Physically, three phases may actually exist, and the EOS model is correctly predicting the behavior. Sometimes a three-phase system is predicted without one physically existing; this may result for lower temperatures when the heaviest component properties are improperly modified to fit two-phase gas/oil PVT data.

For reservoir simulators, the three-phase problem is caused by the EOS formulation "allowing" only two phases. If three phases actually exist, the two-phase flash may find any of the three possible two-phase combinations: L1-V, L2-V, or L1-L2. These false two-phase solutions may indicate a single-phase condition, or they may actually appear as a valid solution (meeting the equal fugacity constraints). Unfortunately, the reservoir model, in a given cell during a given timestep, may flip-flop between two of the possible solutions, resulting in unstable behavior because the pressure solution is not continuous from one two-phase solution to the other. Models may have to simply give up because of repeated timestep reductions, which result from the inadequacy of the EOS two-phase model handling a three-phase condition.^[172]

Surface Phase Behavior

In compositional simulation, the surface calculations are usually made using multistage separation with an EOS model, with fixed K-values for each separator, or using so-called "gas plant" factors, which define the fraction of a wellstream component that is processed into the stock-tank oil.

For black-oil models, the surface separation is "built in" to the PVT tables. Consequently, if the surface process changes during a model run, all black-oil PVT tables must be reentered at the appropriate time. This also requires that vertical flow performance (VFP) tables be reentered because surface rate and volume ratio nodes change with the process change.

It is difficult to use traditional black-oil models for fields with various well groups that have significantly different processing facilities.^[143]

Thermal Model PVT Requirements

Additional PVT requirements for thermal processes such as steamflooding include quantifying the temperature dependence of K-values, densities, and viscosities. Water-phase behavior of liquid and steam must also be defined in terms of pressure and temperature. Water-hydrocarbon phase behavior is still assumed to be simple, without water/hydrocarbon component partitioning.

An EOS model can be tuned to distillation data for improving the predictive capabilities of K-value dependence; otherwise, a simple correlation of the form K_i = a_i exp(– b_iT) / p may be used for distillable components.^[173] Using distillation data is an indirect approach to defining K-value behavior, and it is used in lieu of high-temperature gas/oil/water phase-behavior experiments, which are not usually available. Oil viscosities in thermal processes may be difficult to correlate with a compositional correlation, so empirical correlations may be used instead.

Fluid Initialization

As with rock and other petrophysical properties such as permeability and porosity, a reservoir simulator model must also initialize the spatial distribution of initial fluids. For an EOS-based model, the initial molar compositions are defined, z_i(x,y,z). For a black-oil model, the initial solution gas-oil ratio R s and solution oil-gas ratio r s ratio are defined, R_s(x,y,z) and r_s(x,y,z); sometimes saturation pressures are used instead, p_b(x,y,z) and p_d(x,y,z), but this is not recommended.^[26] Specifying a saturated gas-oil contact (GOC) is also a means to initialize fluids vertically in a reservoir simulator, where solution GOR R_s (and bubblepoint) are assumed constant below the GOC, while solution OGR r_s decreases upwards from the GOC to honor the model-imposed assumption that reservoir pressure equals dewpoint pressure, p(z) = p_d(z).

Another fluid initialization data might include temperature T(x,y,z). Some black-oil models allow spatial variation of stock-tank oil density, γ_API(x,y,z), where black-oil properties are correlated in multidimensional tables as a function of pressure and γ_API . Across barriers such as faults and sealing shales, discrete PVT model data may be defined, such as EOS parameters or black-oil tables; such "discontinuous" fluid descriptions may cause physical and model incompatibilities if fluids mix in the reservoir or at the surface.

A typical problem with initialization is that the specified fluid distribution, initial pressure specifications, and fluid contacts lead to fluid movement when the model is turned on (without production or injection). Initial fluid movement may be unimportant without significantly changing the user-specified fluid system; serious inconsistencies may lead to time-zero flow that has an impact on model performance.

High-Performance Computing and Reservoir Simulation—John E. Killough

The motivation for high-performance computing in reservoir simulation has always existed. From the earliest simulation models, computing resources have been severely taxed simply because the level of complexity desired by the engineer almost always exceeded the speed and memory of the hardware. The high-speed vector processors such as the Cray of the late 1970s and early 1980s led to orders-of-magnitude improvement in speed of computation and led to production models of several hundred thousand cells. The relief brought by these models, unfortunately, was short-lived. The desire for increased physics of compositional modeling and the introduction of geostatistically/structurally based geological models led to increases in computational complexity even beyond the large-scale models of the vector processors. Tens of millions of cells with complete reservoir parameters now became available for use by the engineer. Although upscaling provided a tool to dramatically reduce model sizes, the inherent assumptions of the upscaling techniques left a strong desire by the engineer to incorporate all of the available data in studies. The only available solution to this problem became the subdivision of the model into small segments and the use of parallel computers for reservoir simulation. The recent introduction of fast, low-cost commodity hardware based on the INTEL architecture has led to a revolution in higher-performance computing based on clusters.

The Exposition of Data and the Chasm of Scale

From the earliest reservoir simulation models of the 1960s with only a few tens of finite-difference cells, reservoir simulation had progressed by several orders of magnitude in the early 1990s to hundreds of thousands of cells. Unfortunately, these fine-scale models still were far from the scale of the data. This chasm of scale was further exacerbated by the introduction of geocellular models with geostatistically derived attributes. With these geostatistically based models, it was possible to generate geological descriptions for the reservoir simulation with literally tens of millions of finite-difference cells. Although upscaling offered some relief to this problem, the resultant assumptions required often left the engineer wondering how far from reality the resultant solutions had been driven. If the additional degree of freedom of uncertainty of the reservoir data is added to this problem, then the number and size of simulations becomes unlimited. This implies that the need for further improvements in high-performance computing will always exist for reservoir simulation. parallel computing offers one approach to overcome the problem of the explosion of information through the use of fine-scale models with limited or no upscaling.

Characteristics of High-Performance Computing

The basis of high-performance computing is the use of specialized hardware to achieve computational speeds that are much faster than conventional computers. This idea of "speedup" is encapsulated in what is known as Amdahl’s Law:

Speedup = 1 / (scalar + special / n),

where "scalar" is the fraction of the program which is scalar and "special" is the fraction of computations performed on specialized hardware (one minus scalar). Simply stated, Amdahl’s law indicates that if there is a piece of specialized hardware that is n times faster than the scalar processor(s) in a computer or group of computers, it must be used a large fraction of the time to have a significant impact on overall performance. The earlier supercomputers primarily used specialized hardware known as vector processors for the speedup. More recently, parallel computers based on connecting tens to hundreds of processors have formed the basis of high-performance computers. The importance of the serial or single-processor computer cannot be overemphasized, however. The introduction of the reduced instruction set computer (RISC) with "superscalar" performance in 1990 spelled the end of the dominance of the Cray supercomputer in the high-performance reservoir-simulation market because of a rapid decrease in the price/performance of the new hardware. Similarly, the recent introduction of low-cost, commodity hardware based on the INTEL chipset has led to significant improvements in price/performance. The application of this hardware to reservoir simulation has been a natural evolution in the application of the lowest-cost, highest-performance computers to reservoir simulation.

The Parallel Reservoir Simulator

The earlier high-performance computing hardware such as the Cray Computer required that modifications be made to the reservoir simulator programs to efficiently use the hardware.^[174]^[175] The modifications known as vectorization led to different approaches in the data organization. One of the best examples of this is given in Killough and Levesque^[176] and Young^[177], in which a compositional model’s data structure is reorganized so that like phases in the finite-difference gridblocks of the reservoir are grouped. The notion behind parallel reservoir simulation is similar; however, in this case, the reservoir grid is subdivided into 3D blocks or domains. The computational work for each of the domains is then performed in parallel by a separate CPU of one of the computers in the parallel hardware. Profiles of computing workload for a typical reservoir simulation model often show tens, if not hundreds, of subroutines that are involved in a substantial portion of the calculations. Because of this, major reprogramming of reservoir simulation models is required to achieve high parallel efficiencies. As Killough^[178] points out, there are numerous obstacles to parallelization for reservoir simulation: recursive nature of the linear equation solver, load imbalance caused by reservoir heterogeneity and/or physics, and diverse data structure of well and facility management routines.

Several papers^[179]^[180]^[181]^[182]^[183]^[184]^[185]^[186]^[187]^[188]^[189] discuss techniques that have been used to bring about efficient parallel reservoir simulations. Each of these addresses the solutions to parallelization and, in particular, the obstacles mentioned previously in various ways. One of these simulators uses local grid refinement for parallelization.^[179] The main concept of parallelization with local grid refinement (LGR) is that LGR is not necessarily used to add additional cells to the model but simply to facilitate the subdivision of the grid. With LGR, the same simulation program can be used to perform simulations either serially on a single processor or in parallel on multiple processors simply through data manipulation. In addition, a great deal of flexibility for domain decomposition exists, which can lead to enhanced load balance for parallel simulations.

Parallelization using LGR involves assigning each grid to a processor. The same processor can be assigned to many grids or in the limit; each grid can be assigned to a separate processor. Variation of processor assignment and the grid refinement can dramatically affect parallel performance. Different preconditioners for the parallel linear equation solvers can also have a dramatic effect on parallel efficiency. Finally, the flexibility of LGR and processor assignment can be used to achieve improved parallel efficiency. For parallel domain decomposition, the base grid may be totally refined into subgrids so that all gridblocks in the base grid become inactive. Each grid is assigned to a processor so that only a portion of the model resides on a given processor. In this manner, the simulator memory requirements become scaleable. That is, as the number of processors is increased, the grid assigned to a given processor becomes smaller. Alternatively, only a portion of the base or root grid may be refined; the remaining unrefined cells in the root grid are then assigned to a processor. Because of this use of local grid refinement, only an extremely small number of global arrays are required on each processor—for example, the pore-volume array.

The existence of horizontal wells in most modern simulations requires that no restriction be placed on the placement of wells relative to grid and/or processor boundaries in a parallel simulator. To alleviate any difficulty with this, one approach is for well information to be passed by a broadcast among all processors.^[179] The overhead for broadcasting completion level information for all wells is included in the following results and appears to have little overall effect on parallel efficiency. Tubing hydraulics calculations for a predictive well-management routine have also been parallelized. These calculations are computationally intensive enough to support gather/scatter-type operations. For example, in a case with 1,500 producing wells, more than 90% of the computational workload for well management was in the tubing hydraulics calculations. parallelization of only these calculations has shown a substantial improvement in parallel efficiency for several test cases for the overall well-management routines.

The Parallel Linear-Equation Solver

A typical parallel linear solver involves composite grid decomposition and is based on the work of Wallis and Nolen.^[190] At each level of the composite grid, a red-black or block Jacobi ordering of the grid is used. Red-black refers to an ordering of the grids or domains. For example, if the grids form an equal areal subdivision of an areally square model, the red and black grids would appear as a checkerboard when viewed from the top of the grid. For a 1D or strip decomposition of the model, the red-black coloring would appear as alternately colored stripes when viewed perpendicular to the decomposition direction. An approximate block LU factorization approach is used for preconditioning and involves, for example, an L-factor for the red-black ordering as follows:

....................(17.35)

The A_i blocks are ordered first and represent the red grids or domains, while the B_i blocks are the black grids and are ordered second. The diagonals of the B_i submatrices are modified from the original coefficients to maintain rowsum constraints (see Wallis and Nolen^[190]). This approximate factorization leads to independent solutions on each processor for the first the red A_i and then followed by the black B_i grid matrices. For the block Jacobi approach, only a single color exists, and all grids at a given level are solved independently. As described in Killough, et al.^[179], parallelization of this preconditioning can be achieved in a manner almost identical to the technique used for the coefficient routines. Burrows, Ponting, and Wood^[47] provides detail of a similar parallel solver developed independently, but almost simultaneously, by Burrows, Ponting, and Wood. As described previously, each grid is assigned to a processor. These may be the same or different processors. If all active grids lie on the same level, the composite grid solver performs all red or black preconditioning solutions independently [i.e., either block Jacobi (all red) or red-black orderings are used]. If grids are assigned to different levels, each grid of the same level and same color (Jacobi being all red) is solved simultaneously, with communication being performed after the entire solution of all grids on a level is completed. Factorization is performed so that rowsum constraints are also imposed. The flexibility of this solver allows the easy implementation of block red-black algorithms, or, as a special case, a red-black linear decomposition ("strip") decomposition could be used.

Parallel Computers and Message-Passing Interface

With the rapid changes in technology for parallel computing, it is difficult to provide the state-of-the-art because this appears to change daily; however, current classification practice divides parallel computers into two distinct categories: shared and distributed memory. Shared-memory systems allow all of the memory associated with the parallel computer to be accessed by all of the processors. For UNIX systems, the memory sizes currently can range up to several tens of billion of bytes with hundreds of processors. For systems based on PC products, the shared-memory systems usually are limited to a few billion bytes of storage and fewer than 10 processors. (This may change with the recently introduced 64-bit processor Itanium.) The distributed memory systems allow access to local memory only by a small number of processors that reside on a "node" of the parallel system. Each of the "nodes" in the distributed memory parallel system often communicates with other nodes by a high-speed switch, although conventional ethernet communications are possible for some applications. Either of these parallel computers is capable of delivering several hundred million floating point operations per second (MFLOPS) for each node, or several billions of floating point operations per second (GFLOPS) for applications that can use the parallel hardware. parallel PC hardware can be further subdivided into operating systems: Windows- and LINUX-based systems. Current trends indicated that it is likely that the Windows-based systems will dominate the engineer’s desktop because of the wide availability of engineering/business applications available for this operating system. LINUX-based systems will likely be limited to clusters of servers to perform computer-intensive work such as parallel simulations.

To perform calculations in parallel, data must be passed in the form of messages among processors of the parallel machine. Message-passing interface (MPI) has been primarily used for parallelization on all classifications of parallel machines because of the wide availability and portability of MPI. For shared-memory computers (Sun and SGI for example), other forms of message passing exist such as the recently popular OpenMP. There are generally three forms of message passing for a distributed memory system: synchronous, asynchronous, and global operations. Synchronous message passing requires that all messages be completed before computations can proceed. Asynchronous, on the other hand, allows some overlap of processing and communications. Global message-passing operations distribute data from single processors to one processor, or vice-versa. Often some form of arithmetic function, such as summation, is performed as part of the global concatenation operation. These functions are particularly useful for evaluation of reduction functions such as dot products, or for distribution of data from a single node to all nodes in the parallel partition. Global communications are normally optimized for a particular architecture and are much more efficient than a simple set of send-receive operations for achieving the same results.

Application of Parallel Computing

The rationalization of parallel reservoir simulation is based on two concepts: the ability to run larger models, or simply the ability to improve turnaround time on existing models. A typical larger-size model would be considered to have more than 1 million gridblocks. Results for a 1-million-cell, black-oil, implicit model show scaling is almost perfect from 1 to 18 processors with only a small degradation up to 36 processors (speedup = 32) on an SGI shared-memory parallel computer. Although this is a somewhat idealized example, it does point out the excellent scaling possible because of low message-passing overhead for a parallel simulator using domain decomposition. Results for a real-world compositional example with 3.4 million cells have been obtained on an IBM SP parallel computer with 108 processors. For this case, with 12 to 27 processors, the problem scales extremely well with a speedup of a factor of 2 for slightly more than doubling the number of processors. The scaling falls off slightly from 27 to 54 processors with only a factor of 1.7, but from 54 to 108 processors, this scaling is maintained at a factor of 1.7. Overall, the speedup from 12 to 108 processors was a very good factor of 5.6. If perfect scaling on 12 processors is assumed, this indicates a factor of 68 overall for 107 processors for a parallel efficiency of approximately 64%. The total elapsed time for the 108 processor case was approximately 12 hours, indicating that overnight turnaround for this 3.4-million-cell model could be achieved. Another example application of parallel simulation with more than 1 million gridblocks was recently presented in a black-oil simulation study of complex water encroachment in a large carbonate reservoir in Saudi Arabia.^[189]

The Future of High-Performance Simulation

It is clear that the uptake of parallel reservoir simulation has been limited by the hurdle of cost and the robustness of the technology. The recent advances in commodity hardware based on the INTEL architecture have provided a significant boost, however. Currently, the use of parallel INTEL-based clusters using inexpensive high-speed switches has lowered the entry point for significant parallel simulation by an order of magnitude. Just as the UNIX workstation caused a revolution in reservoir simulation, this new lost-cost hardware will likely bring parallel computing to the engineer’s desktop in the near future. The remaining hurdles for parallel simulation then are limited to enhancing parallel simulation technology. In particular, emphasis must be placed on load balancing and efficient linear equation solvers. Load balancing refers to the fact that all processors must perform the same amount of work to achieve high parallel efficiency. If only a few processors perform most of the work for a parallel simulation involving large numbers of processors, then poor parallel efficiency results—a case known as load imbalance. The allocation of optimization techniques to solve the load balancing problem offers some promise.^[191] An area of particular importance for load balancing is the coupled surface network/reservoir simulator. For this case, the network often dominates the simulation by almost an order of magnitude; therefore, to achieve reasonable parallel efficiency, the imbalance caused by the serial surface network calculations must be overcome. One approach to this is mentioned previously, the parallelization of the tubing hydraulics calculations; however, considerable additional parallelization effort will be required, especially if more than a few processors are to be used efficiently. Load imbalance is also brought about by complex geologies and complex locally refined grids. It is likely that the introduction of unstructured grids and associated load balancing using techniques such as "Hilbert space-filling curves"^[192] may well lead to the solution of this problem. Finally, the linear equation solver as the heart of the parallel simulator must be enhanced to provide efficient and robust solutions. Solver degradation often results from the use of large numbers of domains (processors) or a poor choice of the decomposition geometry.^[179] Although the previous examples show good parallel performance, future models with large, implicit systems must be routinely solved on large numbers of processors. Multigrid-like solvers may show the best promise for this improvement.^[183]^[193]

Reservoir Simulation Applications—L. Kent Thomas

Reservoir simulation is a widely used tool for making decisions on the development of new fields, the location of infill wells, and the implementation of enhanced recovery projects. It is the focal point of an integrated effort of geosciences, petrophysics, reservoir, production and facilities engineering, computer science, and economics. Geoscientists using seismic, well-log, outcrop analog data and mathematical models are able to develop geological models containing millions of cells. These models characterize complex geological features including faults, pinchouts, shales, and channels. Simulation of the reservoir at the fine geologic scale, however, is usually not undertaken except in limited cases. Generally, the fine-scaled geological model is partially integrated or "upscaled" to a coarse-grid model, which is computationally more tractable. The grid of the upscaled model is designed to capture the principal geological features of the geologic model to simulate the fluid flow. The grid may also be designed to capture the effects of complex wells. In the upscaling process, the laboratory relative-permeability and capillary-pressure functions may be upscaled to "pseudofunctions." These pseudofunctions attempt to capture fluid-flow behavior that is lost because of the integration of fine-scale geologic features in the upscaling process. Phase-behavior treatment can range from simple black-oil PVT to compositional and thermal processes.

The reservoir simulation model may either be used directly to forecast the performance of a new reservoir or adjusted so that it reasonably models the historical behavior of an existing reservoir and wells. This adjustment process is called history matching. Programs called "preprocessors" and "post-processors" enable the engineer to prepare data, manipulate the model, and view results. Once a history-matched model is obtained, then forecasts are made under a variety of operating conditions. These results are combined together with economic models to enable the engineer to make decisions concerning the operation of the reservoir.

Development of the Geological Model

A sound understanding of the structural elements of the reservoir and the depositional environment under which the sediments were deposited is critical to the development of an accurate geologic model. Today, the geologic model is frequently constructed as a numerical representation of the reservoir and adjacent aquifer and is referred to as a static, or geocellular, model. This model provides the vehicle to capture and combine the seismic structural interpretation and well petrophysical data in a numerically consistent way with known depositional characteristics.^[194]^[195] Petrophysical properties such as porosity, permeability, and water saturation can be distributed throughout the interwell 3D volume using various techniques, many of which rely on geostatistics.^[196] Efforts are also underway to condition these numerical, static models with production^[197] and well-test^[198] data to further reduce geologic uncertainty. The construction of a geocellular model represents a significant collaborative effort between geoscientists, petrophysicists, and reservoir engineers.

Geocellular models today may consist of over 25 to 50 million cells on large and/or geologically complex reservoirs. The ability to build static geologic models of this magnitude has outstripped the reservoir engineer’s ability to simulate an equal number of cells in a full physics reservoir simulator (and will continue to do so). Classical development geologic efforts have focused on defining and describing the reservoir geology using 2D maps, which depict the most likely interpretation of the depositional environment and the variability of the reservoir parameters between wells. These interpretations have historically been referred to as "deterministic" reservoir descriptions. With the advent of geocellular models and the application of such technologies as geostatistics, it is now possible for geoscientists to generate multiple reservoir descriptions for the reservoir engineer to simulate. In some cases, one of these descriptions may be selected to represent the "deterministic" model. Regardless if one or several static models are handed over for reservoir simulation, it is generally necessary to reduce the cell count to run the problem with existing reservoir simulators. Significant effort is being spent improving techniques to reduce the number of reservoir cells in the areal and vertical dimension while maintaining the essential geologic character that impacts the recovery process under consideration. This approach is referred to as upscaling, and it will be discussed in greater detail in the following section. To date, the largest reservoir simulators consist of reservoir descriptions of 2 million grid cells and are run using massively parallel processing power.

Upscaling Geological Model to Reservoir Flow Model

Geological models, which contain the complex structural features of large oil and gas reservoirs, commonly have tens of millions of cells. These models, which contain pinchouts, faults, and other significant information including lithology and facies distributions, are upscaled in both the vertical and areal directions to tens or hundreds of thousands of cells for reservoir simulation.^[199]

Several upscaling methods have been developed over the last several years including analytical techniques, local flow-based methods, and global flow-based methods. Analytical methods use arithmetic, harmonic, power law, and geometric averaging to calculate effective properties for each reservoir model gridblock. Local flow-based methods calculate effective gridblock properties by performing single-phase flow simulations in each direction across the upscaled block.^[200] The diagonal permeability tensor is calculated by sealing the boundaries perpendicular to the applied pressure gradient. The full-permeability tensor can be calculated in a similar manner by leaving the boundaries normal to the imposed pressure gradient open and applying periodic boundary conditions. Global flow-based methods use pressure gradients across the entire field subject to a specific set of wells to calculate the permeability tensor. Local and global flow-based techniques can be used to compute upscaled transmissibilities directly.

Inclusion of Faults in Reservoir Flow Model

Faults and pinchouts of geological layers are incorporated in geological models to capture the complex geometry of many reservoirs. This information is then upscaled into the reservoir model, and it results in both neighbor and non-neighbor connections across the faults and non-neighbor connections across the pinchouts. In Cartesian coordinates, the trace of a fault may need to be represented by a "stair-stepped" line, while a somewhat better representation of faults can be made with corner-point grids. PEBI grids, which will be discussed subsequently, are best suited to accurately model fault geometry.

Models for calculating the transmissibility across the fault and parallel to the fault have been developed based on fault type, displacement, geochemical deposition, and whether open joints occur along the fault.^[201] In general, transmissibilities across a fault can be at least an order of magnitude lower than those parallel to the fault. Inclusion of this information in a reservoir model is frequently a key parameter in reservoir description.

A recent paper describes the analysis that was performed to calculate fluid flow through conductive faults in the Khafji oil field in the Arabian Gulf.^[202] Two sandstone reservoirs separated by a thick continuous shale are both connected to the same large aquifer and had the same initial WOC. The top reservoir has edgewater drive, while the deeper reservoir is bottomwater drive. Early water breakthrough in the upper sand was determined to be a function of supplemental water influx from the aquifer of the lower sand through conductive faults.

Development of Pseudofunctions for Multiphase Flow

Pseudorelative permeability curves are developed for upscaled reservoir models to match multiphase fluid flow at the fine-grid level. Several methods for performing these calculations have been presented in the literature.^[203] In the "10th SPE Comparative Solution Project: A Comparison of Upscaling Techniques," the fine-scale geological model was chosen to be sufficiently detailed such that it would be difficult to run the fine grid and use classical pseudoization methods.^[12] Several participants, however, used some level of fine-grid simulation to develop pseudorelative permeability curves, with two of the participants adjusting the exponents of the Corey equations to effect a reasonable match. This approach can be done manually or with an automated history-matching algorithm.

Gridding Techniques

The majority of reservoir simulation studies conducted today use Cartesian or corner-point structured grids with some application of local grid refinement to evaluate infill well locations or to more accurately calculate water and/or gas coning in a well. In a structured grid, cell locations are specified using a 3D, i, j, k, indexing system. This allows for ready access either numerically or visually, using pre- and post-processing software, to multilayer well information or data and calculated results at any physical location in the reservoir model.

A more flexible approach for modeling reservoirs with complex geometries that still relies on structured gridding was presented by Jenny et al.^[63] Here, a hexahedral multiblock grid is used which is unstructured globally, but maintains an i, j, k structure on each subgrid.

PEBI grids ^[58] are now being used on a limited basis to simulate reservoirs with complex geological features that have been developed with nonconventional wells to maximize recovery.^[204] These grids are unstructured and are described internally in a simulator with a 1D index, i, that ranges from one to the number of nonzero pore volume cells in a model. Evaluation of simulator input and results relies heavily on pre- and post-processing software that allows the user to visually look at the model and make changes during the history-matching phase of a study.

Simulation of Nonconventional and Intelligent Wells

Nonconventional wells are routinely used to maximize production rates and ultimate recovery in oil and gas reservoirs. Wells in this category include deviated, horizontal, horizontal and vertical multilaterals, and multilateral fishbone designs. This latter well type is especially effective in low-permeability or heavy-oil reservoirs.

Simulation of nonconventional wells can be approached in several ways. First, the productivity of each perforation in a conventional model can be approximated by applying the appropriate skin and Peaceman’s equation.^[205] Second, simulation grids can be constructed that closely follow the well path and allow a more accurate calculation of well rates.^[206] Another approach, which is quite appealing, is based on a semi-analytical method.^[207] It results in a good approximation for PIs in nonconventional wells and incorporates the near-wellbore skin because of heterogeneity in this region as well as mechanical skin. This method, which is very efficient, can also include wellbore hydraulic effects.

Nonconventional wells coupled with intelligent completions can be used to improve sweep efficiency and optimize recovery.^[208] One example of this technology is the use of surface-controlled, downhole-adjustable chokes, which can be used to apply different pressure drawdowns to separate zones along the well. This allows a more uniform inflow in the well and control of early water or gas breakthrough. Real-time measurements of wellbore pressures and temperatures are being made for use in conjunction with PLT tests for inflow performance analysis.

Integrated Reservoir and Surface Facilities Models

Integration of reservoir and surface facilities simulation can result in improved production forecasts and allows optimization of the surface facilities structure and operating conditions. An integrated reservoir, well flow string, and surface network model of the Prudhoe Bay oil field was built and successfully applied to a facility optimization study.^[209] Production costs as a result of this effort were reduced by defining the optimum number of separator stages and their connections, and defining the optimum separator operating conditions and by using excess capacity in the Prudhoe Bay facilities to process production from satellite fields. Procedures for the simultaneous solution of the reservoir and surface pipeline network flow equations are described in Fang and Lo^[210] and Litvak^[211]

In the Ekofisk field in the Norwegian sector of the North Sea, integrated reservoir and facilities simulations have been made to optimize throughput in existing surface facilities and to forecast production from planned expansion of current facilities.^[212] This optimization project has resulted in sustained high production of approximately 300,000 STB/D over the last several years. Another important aspect in the management of this field is the inclusion of compaction logic in the model based on both stress and water saturation changes during depletion and waterflooding.^[213] Treatment of geomechanical effects in stress-sensitive reservoirs has received increased attention throughout the industry in recent years.

Simulation of Multiple Reservoirs

Simulation of multiple fields producing into a common production facility is routinely practiced today to capture the interplay between well deliverability, water and gas injection, operating constraints, and contract rates. In the J-Block area fields, in the U.K. sector of the North Sea, an integrated reservoir study was conducted that included a gas condensate reservoir with gas cycling that was simultaneously modeled with volatile oil reservoirs.^[214] The fields were developed with a single platform and one subsea manifold completion and a combination of vertical and horizontal wells. Four separate PVT regions were used to describe the fluid behavior. The integrated model used in this study results in an efficient reservoir management tool for making development and operating decisions.

Another example of reservoir management of multiple fields with shared facilities is the Gannet cluster, located in the U.K. sector of the North Sea, which connects four fields.^[215] Wells from one of the fields are directly linked to the production platform, and the other three fields are subsea tiebacks to the platform. Three of the four fields are oil fields and the fourth is a gas field. An integrated model was built to simulate the interaction of the subsurface and surface processes. The well-management objective on this project was to maximize hydrocarbon recovery while simultaneously meeting a long-term gas contract.

Use of Larger Models

The maximum practical model size has increased from tens of thousands to hundreds of thousands of cells at essentially a linear rate vs. time during the last decade. This trend has developed as a result of the dramatic increase in computer hardware speed accompanied with larger memory and cache. Both high-speed UNIX workstations and high-end PCs are used for reservoir simulation, with a close race developing between the two platforms in regard to run times. Additional advances in computing speed for megamodels have been achieved using parallel hardware along with the necessary developments in model software, discussed in a previous section. An example application of this technology was recently presented in a simulation study of complex water encroachment in a large carbonate reservoir in Saudi Arabia.^[189]

History Matching and Production Forecasting

Once a reservoir simulation model has been constructed, the validity of the model is examined by using it to simulate the performance of a field under past operating conditions. This is usually done by specifying historical controlling rates, such as oil rate in an oil reservoir vs. time, and then making a comparison of the nonspecified performance such as GOR, WOR, and reservoir pressure with measured data. If there are significant differences between the calculated performance and the known performance of the well/reservoir system, then adjustments to the reservoir simulation model are made to reduce this difference. This process is called history matching. These adjustments should be made in a geologically consistent manner.^[198] Modification of those parameters that have the highest degree of uncertainty will give the maximum reduction in the error. The history-matching process should be approached in a consistent manner to minimize the effort.^[216] In addition to well rates and bottomhole pressures, and reservoir pressures measured at the time the well is drilled, production logs, long-term pressure gauges, and time-lapse seismic data enable the engineer to better constrain the model during the history-matching process. Time-lapse (4D) seismic^[217]^[218] is becoming an integral part of the field performance monitoring and history matching. Streamline models together with reservoir simulators^[219]^[197] can be used to improve the history-matching process, especially in waterflood operations. Tools to assist in the history-matching process consist of the use of parallel computers, sensitivity analysis, and gradient techniques.^[220]^[221]

Once a history match is obtained, then forecasts of future well/reservoir performance under various operating scenarios are made. Models of multiphase flow in the wellbore and production lines are used to constrain the production rate. These models may include subsea completions^[222] with very long gathering lines or complex surface facilities with reinjection of produced fluids.^[223] Because of the uncertainty in the geological and reservoir simulation models for new fields, often multiple forecasts^[223] with different reservoir parameters are made to determine the uncertainty in the forecasts. Multiple history-matched models based on multiple geological models,^[224] and experimental design^[225] may also be used to characterize the uncertainty in production forecasts.

Nomenclature

A	=	a matrix
B	=	formation volume factor, reservoir volume/surface volume
b	=	right-side vector
C_os	=	equivalent surface gas volume for a unit volume of surface oil
d	=	dimension in direction of flow
D	=	depth below datum, ft
D	=	future time, time units
f	=	fractional flow
f_i	=	fugacity of component i
F	=	flow rate, volume/time unit
F , F_i, f_iI	=	see Eq. 17.2
f_g	=	gas-phase fractional flow, λ_g / (λ_g + λ_o)
f_g′	=	∂f_g / ∂S_g
g	=	acceleration constant (length/length/time) or gas phase, depending on use
G	=	phase velocity
k	=	permeability, md
k_ij	=	binary interaction parameters
k_r	=	relative permeability, fraction
K	=	permeability, md
K_i	=	equilibrium constant for component i
K_r	=	relative permeability, fraction
k_rocw	=	relative permeability of oil at S_wc, S_g = 0, fraction
k_rwro	=	relative permeability of water at residual oil, fraction
L	=	length, feet
L1	=	liquid type 1
L2	=	liquid type 2
m	=	density, mass/unit volume
M	=	mass, mobility, molecular weight, or preconditioner matrix depending on use
n	=	number of components or number of moles
n_st	=	number of streamlines
n_ts	=	number of timesteps
N	=	number of active gridblocks
o	=	oil phase
p	=	polynomial
p_J	=	phase J pressure, psi
P	=	pressure, psi
P_c	=	capillary pressure, psi
P_cgo′	=	dP_cgo / dS_g
P_i	=	vector of primary unknowns for block i
P_ii	=	ith scalar element of P_i
q_ii	=	well production rate of component i from block i, volume/time
q_Iji	=	interblock flow rate of component I from block j to block i, volume/time
q_x	=	total volumetric flow rate in the x direction, volume/time
Q	=	production or injection volume
r_s	=	solution gas-condensate ratio, volume/volume
R_s	=	solution gas, scf/STB oil or nm³ /nm³oil
S	=	saturation
S_on	=	(1.0 – S_w – S_orw)/(1.0 – S_wc – S_orw), fraction
S_orw	=	residual oil saturation after water displacement
S_wc	=	connate water saturation, fraction
S_wn	=	(S_w – S_wc)/(1.0 – S_wc – S_orw), fraction
t^solver	=	time to solve for a global pressure field
t^st	=	time to solve transport equation for each streamline
T	=	time, time units or transmissibility, md•ft/cp
T_ij	=	transmissibility connecting blocks i and j, md•ft/cp
u	=	pore velocity, v/Ф;
u_t	=	total velocity
v	=	Darcy or superficial velocity
V	=	vapor phase or volume
V_p	=	pore volume, bbl
w	=	water phase
x	=	an unknown vector or concentration of a component in the liquid phase
x,y,z	=	Cartesian coordinates
x_i	=	concentration of component i in oil phase
x_iJ	=	mol fraction of component i in phase J
y	=	an unknown vector
y_i	=	concentration of component i in the gas phase
z	=	gas deviation factor
z_i	=	overall mol fraction of component i in a gridblock
Z	=	depth to gridblock center, measured positive downward, ft or m
α	=	mole fraction vapor phase
β	=	see Linear Solver section, Step 6
γ	=	phase density expressed as psi/ft or bar/m
δ	=	δX = X^l+1 – X^l , where l is the Newton iterate index and X is any quantity
Δp	=	pressure difference, psi
ΔQ	=	potential difference
Δt	=	timestep
Δt*	=	maximum stable IMPES timestep
Δx, Δy, Δz	=	gridblock dimensions
ΔZ	=	Z_j – Z_i
θ	=	dip angle of line joining cell centers, degrees
λ	=	phase mobility, k_r / μ
μ	=	viscosity, cp
ρ	=	density, mass/volume or molar density, mols/volume unless noted otherwise
Σ	=	sum
τ	=	time-of-flight along a streamline
Ф	=	porosity, fraction
ψ	=	λ_oλ_g / (λ_o + λ_g)
ω′	=	acentric factor
∇	=	potential term
∂	=	normalized difference parameter

Subscripts

a	=	cell designation
a-h	=	arithmetic-harmonic average
ari	=	arithmetic average
b	=	adjacent cell designation
c	=	component designation
eff	=	effective
g	=	gas phase
h	=	horizontal
h-a	=	harmonic-arithmetic average
har	=	harmonic average
I	=	component number, index counter, or initial condition
j	=	index counter
J	=	phase number
k	=	index counter
lbc	=	linear (open sides) boundary conditions
low	=	lower bound
nbc	=	no-flow (sealed sides) boundary conditions
o	=	oil phase
p	=	phase designation
pb	=	periodic boundary conditions
r	=	residual
R	=	reservoir
s	=	surface
sc	=	standard conditions
true	=	true
upp	=	upper bound
v	=	vertical
w	=	water phase or well
x	=	cell designation
x,y,z	=	Cartesian x,y,z directions
z_eff	=	effective in z direction

Superscripts

0	=	reference (in compressibility equation)
l	=	Newton iterate index
n	=	timestep number

References

↑ The Merriam-Webster Dictionary, New Revised Edition. 2004. New York City: Merriam-Webster.
↑ Coats, K.H. 1987. Reservoir Simulation. Petroleum Engineering Handbook, H.B. Bradley ed., Chap. 48. Richardson, Texas: SPE.
↑ Odeh, A.S. 1981. Comparison of Solutions to a Three-Dimensional Black-Oil Reservoir Simulation Problem. J Pet Technol 33 (1): 13–25. SPE-9723-PA. http://dx.doi.org/10.2118/9723-PA.
↑ Weinstein, H.G., Chappelear, J.E., and Nolen, J.S. 1986. Second Comparative Solution Project: A Three-Phase Coning Study. J Pet Technol 38 (3): 345-353. SPE-10489-PA. http://dx.doi.org/10.2118/10489-PA.
↑ Kenyon, D. 1987. Third SPE Comparative Solution Project: Gas Cycling of Retrograde Condensate Reservoirs. J Pet Technol 39 (8): 981-997. SPE-12278-PA. http://dx.doi.org/10.2118/12278-PA.
↑ ^6.0 ^6.1 Aziz, K., Ramesh, A.B., and Woo, P.T. 1987. Fourth SPE Comparative Solution Project: Comparison of Steam Injection Simulators. J Pet Technol 39 (12): 1576–1584. SPE-13510-PA. http://dx.doi.org/10.2118/13510-PA.
↑ Killough, J.E. and Kossack, C.A. 1987. Fifth Comparative Solution Project: Evaluation of Miscible Flood Simulators. Presented at the SPE Symposium on Reservoir Simulation, San Antonio, Texas, 1–4 February. SPE-16000-MS. http://dx.doi.org/10.2118/16000-MS.
↑ Firoozabadi, A. and Thomas, L.K. 1990. Sixth SPE Comparative Solution Project: Dual-Porosity Simulators. J Pet Technol 42 (6): 710-715, 762-763. SPE-18741-PA. http://dx.doi.org/10.2118/18741-PA.
↑ Nghiem, L., Collins, D.A., and Sharma, R. 1991. Seventh SPE Comparative Solution Project: Modelling of Horizontal Wells in Reservoir Simulation. Presented at the SPE Symposium on Reservoir Simulation, Anaheim, California, 17-20 February 1991. SPE-21221-MS. http://dx.doi.org/10.2118/21221-MS.
↑ Quandalle, P. 1993. Eighth SPE Comparative Solution Project: Gridding Techniques in Reservoir Simulation. Presented at the SPE Symposium on Reservoir Simulation, New Orleans, Louisiana, 28 February-3 March 1993. SPE-25263-MS. http://dx.doi.org/10.2118/25263-MS.
↑ ^11.0 ^11.1 ^11.2 Killough, J.E. 1995. Ninth SPE Comparative Solution Project: A Reexamination of Black-Oil Simulation. Presented at the SPE Reservoir Simulation Symposium, San Antonio, Texas, USA, 12–15 February. SPE 29110. http://dx.doi.org/10.2118/29110-MS.
↑ ^12.0 ^12.1 ^12.2 ^12.3 Christie, M.A. and Blunt, M.J. 2001. Tenth SPE Comparative Solution Project: A Comparison of Upscaling Techniques. SPE Res Eval & Eng 4 (4): 308–317. SPE-72469-PA. http://dx.doi.org/10.2118/72469-PA.
↑ ^13.0 ^13.1 ^13.2 ^13.3 ^13.4 ^13.5 ^13.6 ^13.7 ^13.8 Coats, K.H. 1980. An Equation of State Compositional Model. J Pet Technol 20 (5): 363–376. SPE-8284-PA. http://dx.doi.org/10.2118/8284-PA.
↑ ^14.0 ^14.1 ^14.2 ^14.3 ^14.4 ^14.5 ^14.6 Young, L.C. and Stephenson, R.E. 1983. A Generalized Compositional Approach for Reservoir Simulation. SPE J. 23 (4): 727–742. SPE-10516-PA. http://dx.doi.org/10.2118/10516-PA.
↑ ^15.0 ^15.1 ^15.2 ^15.3 ^15.4 ^15.5 ^15.6 Ács, G., Doleschall, S., and Farkas, E. 1985. General Purpose Compositional Model. SPE J. 25 (4): 543–553. SPE-10515-PA. http://dx.doi.org/10.2118/10515-PA.
↑ Watts, J.W. 1986. A Compositional Formulation of the Pressure and Saturation Equations. SPE Res Eng 1 (3): 243–252. SPE-12244-PA. http://dx.doi.org/10.2118/12244-PA.
↑ Blair, P.M. and Weinaug, C.F. 1969. Solution of Two-Phase Flow Problems Using Implicit Difference Equations. SPE J. 9 (4): 417-424. SPE-2185-PA. http://dx.doi.org/10.2118/2185-PA.
↑ Stone, H. and Garder, A.O.J. 1961. Analysis of Gas-Cap or Dissolved-Gas Drive Reservoirs. SPE J. 1 (2): 92–104. SPE-1518-G. http://dx.doi.org/10.2118/1518-G.
↑ Sheldon, J.W., Harris, C.D., and Bavly, D. 1960. A Method for General Reservoir Behavior Simulation on Digital Computers. Presented at the Fall Meeting of the Society of Petroleum Engineers of AIME, Denver, Colorado, 2-5 October. http://dx.doi.org/10.2118/1521-G.
↑ Fagin, R.G. and Stewart Jr., C.H. 1966. A New Approach to the Two-Dimensional Multiphase Reservoir Simulator. SPE J. 6 (2): 175-182. SPE-1188-PA. http://dx.doi.org/10.2118/1188-PA.
↑ ^21.0 ^21.1 Coats, K.H.: "Reservoir Simulation: A General Model Formulation and Associated Physical/Numerical Sources of Instability," Boundary and Interior Layers—Computational and Asymptotic Methods, J.J.H. Miller (ed.), Boole Press, Dublin, Ireland (1980).
↑ ^22.0 ^22.1 ^22.2 ^22.3 ^22.4 ^22.5 ^22.6 ^22.7 ^22.8 ^22.9 Coats, K.H., Thomas, L.K., and Pierson, R.G. 1998. Compositional and Black Oil Reservoir Simulation. SPE Res Eval & Eng 1 (4): 372-379. SPE-50990-PA. http://dx.doi.org/10.2118/50990-PA.
↑ ^23.0 ^23.1 Coats, K.H. 2000. A Note on IMPES and Some IMPES-Based Simulation Models. SPE J. 5 (3): 245-251. SPE-65092-PA. http://dx.doi.org/10.2118/65092-PA.
↑ MacDonald, R.C. 1970. Methods for Numerical Simulation of Water and Gas Coning. SPE J. 10 (4): 425-436. SPE-2796-PA. http://dx.doi.org/10.2118/2796-PA.
↑ Spillette, A.G., Hillestad, J.G., and Stone, H.L. 1973. A High-Stability Sequential Solution Approach to Reservoir Simulation. Presented at the Fall Meeting of the Society of Petroleum Engineers of AIME, Las Vegas, Nevada, 30 September-3 October 1973. SPE-4542-MS. http://dx.doi.org/10.2118/4542-MS.
↑ ^26.0 ^26.1 Thomas, G.W. and Thurnau, D.H. 1983. Reservoir Simulation Using an Adaptive Implicit Method. SPE J. 23 (5): 759–768. SPE-10120-PA. http://dx.doi.org/10.2118/10120-PA.
↑ Lantz, R.B. 1971. Quantitative Evaluation of Numerical Diffusion (Truncation Error). SPE J. 11 (3): 315–320. SPE-2811-PA. http://dx.doi.org/10.2118/2811-PA.
↑ Todd, M.R., O'Dell, P.M., and Hirasaki, G.J. 1972. Methods for Increased Accuracy in Numerical Reservoir Simulators. SPE J. 12 (6): 515–530. SPE-3516-PA. http://dx.doi.org/10.2118/3516-PA.
↑ Fung, L.S.-K., Collins, D.A., and Nghiem, L.X. 1989. An Adaptive-Implicit Switching Criterion Based on Numerical Stability Analysis (includes associated paper 18727 ). SPE Res Eng 4 (1): 45-51. SPE-16003-PA. http://dx.doi.org/10.2118/16003-PA.
↑ Russell, T.F. 1989. Stability Analysis and Switching Criteria for Adaptive Implicit Methods Based on the CFL Condition. Presented at the SPE Symposium on Reservoir Simulation, Houston, Texas, 6–8 February. SPE-18416-MS. http://dx.doi.org/10.2118/18416-MS.
↑ Grabenstetter, J., Li, Y.-K., Collins, D.A. et al. 1991. Stability-Based Switching Criterion for Adaptive-Implicit Compositional Reservoir Simulation. Presented at the SPE Reservoir Simulation Symposium, Anaheim, California, 17–20 February. SPE-21225-MS. http://dx.doi.org/10.2118/21225-MS.
↑ Watts, J.W. and Rame, M. 1999. An Algebraic Approach to the Adaptive Implicit Method. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 14-17 February 1999. SPE-51900-MS. http://dx.doi.org/10.2118/51900-MS.
↑ Coats, K.H. 2001. IMPES Stability: The CFL Limit. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 11-14 February 2001. SPE-66345-MS. http://dx.doi.org/10.2118/66345-MS.
↑ Coats, K.H. 2001. IMPES Stability: The Stable Step. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 11-14 February 2001. SPE-69225-MS. http://dx.doi.org/10.2118/69225-MS.
↑ Young, L.C. and Russell, T.F. 1993. Implementation of an Adaptive Implicit Method. Presented at the SPE Symposium on Reservoir Simulation, New Orleans, Louisiana, 28 February-3 March 1993. SPE-25245-MS. http://dx.doi.org/10.2118/25245-MS.
↑ ^36.0 ^36.1 Vinesome, P.K.W. 1976. Orthomin, an Iterative Method for Solving Sparse Banded Sets of Simultaneous Linear Equations. Presented at the SPE Symposium on Numerical Simulation of Reservoir Performance, Los Angeles, 19–20 February. SPE-5729-MS. http://dx.doi.org/10.2118/5729-MS.
↑ ^37.0 ^37.1 ^37.2 Appleyard, J.R. 1983. Nested Factorization. Presented at the SPE Reservoir Simulation Symposium, San Francisco, California, USA, 15-18 November. SPE-12264-MS. http://dx.doi.org/10.2118/12264-MS.
↑ Watts III, J.W. 1981. A Conjugate Gradient-Truncated Direct Method for the Iterative Solution of the Reservoir Simulation Pressure Equation. Society of Petroleum Engineers Journal 21 (3): 345-353. SPE-8252-PA. http://dx.doi.org/10.2118/8252-PA.
↑ Forsyth, P.A. and Sammon, P.H. 1986. Local Mesh Refinement and Modeling of Faults and Pinchouts. SPE Form Eval 1 (3): 275-285. SPE-13524-PA. http://dx.doi.org/10.2118/13524-PA.
↑ Verma, S.K. and Aziz, K. 1997. A Control Volume Scheme for Flexible Grids In Reservoir Simulation. Presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, 8–11 June. SPE-37999-MS. http://dx.doi.org/10.2118/37999-MS.
↑ Saad, Y. and Schultz, M.H. 1986. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing 7 (3): 856-869. http://dx.doi.org/10.1137/0907058.
↑ Tan, T.B.S. and Letkeman, J.P. 1982. Application of D4 Ordering And Minimization In An Effective Partial Matrix Inverse Iterative Method. Presented at the SPE Reservoir Simulation Symposium, New Orleans, Louisiana, 31 January-3 Feburary 1982. SPE-10493-MS. http://dx.doi.org/10.2118/10493-MS.
↑ Watts, J.W. 1971. An Iterative Matrix Solution Method Suitable for Anisotropic Problems. Society of Petroleum Engineers Journal 11 (1): 47-51. SPE-2802-PA. http://dx.doi.org/10.2118/2802-PA.
↑ ^44.0 ^44.1 Wallis, J.R. 1983. Incomplete Gaussian Elimination as a Preconditioning for Generalized Conjugate Gradient Acceleration. Presented at the SPE Reservoir Simulation Symposium, San Francisco, California, USA, 15-18 November. SPE-12265-MS. http://dx.doi.org/10.2118/12265-MS.
↑ Behie, A. and Vinsome, P.K.W. 1982. Block Iterative Methods for Fully Implicit Reservoir Simulation. Society of Petroleum Engineers Journal 22 (5): 658-668. SPE-9303-PA. http://dx.doi.org/10.2118/9303-PA.
↑ Holmes, J.A. 1983. Enhancements to the Strongly Coupled, Fully Implicit Well Model: Wellbore Crossflow Modeling and Collective Well Control. Presented at the SPE Reservoir Simulation Symposium, San Francisco, California, 15-18 November 1983. SPE-12259-MS. http://dx.doi.org/10.2118/12259-MS.
↑ ^47.0 ^47.1 Burrows, R., Ponting, D., and Wood, L. 1996. Parallel Reservoir Simulation with Nested Factorization. Presented at the 1996 European Conference on the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.
↑ Killough, J.E. et al. 1996. Parallelization of a General-Purpose Reservoir Simulator. Presented at the 1996 European Conference on the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.
↑ Cuthill, E. and McKee, J. 1969. Reducing the Bandwidth of Sparse Symmetric Matrices. Proc., 24th National Conference of the Assoc. for Computing Machinery, Brandon Press, Princeton, New Jersey, 157–172.
↑ D'Azevedo, E., Forsyth, P.A., and Tang, W.-P. 1991. An Automatic Ordering Method for Incomplete Factorization Iterative Solvers. Presented at the SPE Symposium on Reservoir Simulation, Anaheim, California, 17-20 February 1991. SPE-21226-MS. http://dx.doi.org/10.2118/21226-MS.
↑ D’Azevedo, E.F., Forsyth, P.A., and Tang, W.-P. 1992. Ordering Methods for Preconditioned Conjugate Gradient Methods Applied to Unstructured Grid Problems. SIAM Journal of Matrix Analysis and Applications 13 (3): 944-961. http://dx.doi.org/10.1137/0613057.
↑ ^52.0 ^52.1 Aziz, K. and Settari, A. 1979. Petroleum Reservoir Simulation, ISBN 0–9730614-0-5, reprinted 2002. Calgary: Blitzprint Ltd.
↑ ^53.0 ^53.1 Rozon, B.J. 1989. A Generalized Finite Volume Discretization Method for Reservoir Simulation. Presented at the SPE Symposium on Reservoir Simulation, Houston, Texas, 6-8 February 1989. SPE-18414-MS. http://dx.doi.org/10.2118/18414-MS.
↑ ^54.0 ^54.1 Aavatsmark, I. et al. 2000. MPFA for Faults and Local Refinements With Application to Field Simulations. Proc., 2000 European Conference on the Mathematics of Oil Recovery, Baveno, Italy, 5–8 September.
↑ Forsyth, P.A. 1990. A Control-Volume, Finite-Element Method for Local Mesh Refinement in Thermal Reservoir Simulation. SPE Res Eng 5 (4): 561-566. SPE-18415-PA. http://dx.doi.org/10.2118/18415-PA.
↑ Hegre, T.M., Dalen, V., and Henriquez, A. 1986. Generalized Transmissibilities for Distorted Grids in Reservoir Simulation. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 5-8 October 1986. SPE-15622-MS. http://dx.doi.org/10.2118/15622-MS.
↑ ^57.0 ^57.1 Young, L.C. 1999. Rigorous Treatment of Distorted Grids in 3D. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 14-17 February 1999. SPE-51899-MS. http://dx.doi.org/10.2118/51899-MS.
↑ ^58.0 ^58.1 ^58.2 Aziz, K. 1993. Reservoir Simulation Grids: Opportunities and Problems. J Pet Technol 45 (7): 658-663. SPE-25233-PA. http://dx.doi.org/10.2118/25233-PA.
↑ Aavatsmark, I., Barkve, T., and Mannseth, T. 1998. Control-Volume Discretization Methods for 3D Quadrilateral Grids in Inhomogeneous, Anisotropic Reservoirs. SPE J. 3 (2): 146-154. SPE-38000-PA. http://dx.doi.org/10.2118/38000-PA.
↑ Edwards, M.G. and Rogers, C.F. 1994. A Flux Continuous Scheme for the Full Tensor Pressure Equation. Proc., 1994 European Conference on the Mathematics of Oil Recovery, Roros, Norway, 7–10 June.
↑ Lee, S.H., Tchelepi, H., and DeChant, L.F. 1999. Implementation of a Flux-Continuous Finite Difference Method for Stratigraphic, Hexahedron Grids. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 14-17 February 1999. SPE-51901-MS. http://dx.doi.org/10.2118/51901-MS.
↑ Ponting, D.K. 1989. Corner Point Geometry in Reservoir Simulation. Proc., First European Conference on the Mathematics of Oil Recovery.
↑ ^63.0 ^63.1 Jenny, P., Wolfsteiner, C., Lee, S.H. et al. 2002. Modeling Flow in Geometrically Complex Reservoirs Using Hexahedral Multiblock Grids. SPE J. 7 (2): 149-157. SPE-78673-PA. http://dx.doi.org/10.2118/78673-PA.
↑ Heinemann, Z.E. and Brand, C.W. 1988. Gridding Techniques in Reservoir Simulation. Proc., First Intl. Forum on Reservoir Simulation, Alpbach, Austria, 339.
↑ Heinemann, Z.E., Brand, C.W., Munka, M. et al. 1991. Modeling Reservoir Geometry With Irregular Grids. SPE Res Eng 6 (2): 225–232. SPE-18412-PA. http://dx.doi.org/10.2118/18412-PA.
↑ Palagi, C.L. and Aziz, K. 1994. Use of Voronoi Grid in Reservoir Simulation. SPE Advanced Technology Series 2 (2): 69–77. SPE-22889-PA. http://dx.doi.org/10.2118/22889-PA.
↑ Gunasekera, D., Cox, J., and Lindsey, P. 1997. The Generation and Application of K-Orthogonal Grid Systems. Presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, 8–11 June. SPE-37998-MS. http://dx.doi.org/10.2118/37998-MS.
↑ Verma, S. and Aziz, K. 1996. Two- and Three-Dimensional Flexible Grids for Reservoir Simulation. Proc., 1996 European Conference on the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.
↑ Lasseter, T.J. 2002. A New Approach for the Efficient Construction of 3D Geological Models for Reservoir Applications. Proc., Eighth European Conference on the Mathematics of Oil Recovery, Freiberg, Germany.
↑ Warren, J.E. and Price, H.S. 1961. Flow in Heterogeneous Porous Media. SPE J. 1 (3): 153–169. SPE-1579-G. http://dx.doi.org/10.2118/1579-G.
↑ King, M.J., MacDonald, D.G., Todd, S.P. et al. 1998. Application of Novel Upscaling Approaches to the Magnus and Andrew Reservoirs. Presented at the European Petroleum Conference, The Hague, 20–22 October. SPE-50643-MS. http://dx.doi.org/10.2118/50643-MS.
↑ King, M.J. and Mansfield, M. 1999. Flow Simulation of Geologic Models. SPE Res Eval & Eng 2 (4): 351-367. SPE-57469-PA. http://dx.doi.org/10.2118/57469-PA.
↑ Guerillot, D. et al. 1989. An Integrated Model for Computer Aided Reservoir Description: From Outcrop Study to Fluid Flow Simulations. Paper presented at the 1989 IOR Symposium, Budapest, Hungary.
↑ Samier, P. 1990. A Finite-Element Method for Calculating Transmissibilities. Proc., 1990 European Conference on the Mathematics of Oil Recovery, Arles, France, 11–14 September.
↑ ^75.0 ^75.1 Durlofsky, L.J. and Chung, E.Y. 1990. Effective Permeability of Heterogeneous Reservoir Regions. Proc., 1990 European Conference on the Mathematics of Oil Recovery, Arles, France, 11–14 September.
↑ ^76.0 ^76.1 Durlofsky, L.J. 1991. Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media. Water Resour. Res. 27 (5): 699-708. http://dx.doi.org/10.1029/91wr00107.
↑ Kirkpatrick, S. 1971. Classical Transport in Disordered Media: Scaling and Effective-Medium Theories. Phys. Rev. Lett. 27 (25): 1722-1725. http://dx.doi.org/10.1103/PhysRevLett.27.1722.
↑ Tjølsen, C., Bøe, Ø., and Damsleth, E. 1992. Homogenisation of Absolute Permeability—Does it Affect the Final Results? Proc., Second Lerkendal Petroleum Engineering Workshop.
↑ Bøe, Ø. 1994. Analysis of an upscaling method based on conservation of dissipation. Transport Porous Media 17 (1): 77-86. http://dx.doi.org/10.1007/bf00624051.
↑ ^80.0 ^80.1 ^80.2 Holden, L. and Nielsen, B.F. 2000. Global Upscaling of Permeability in Heterogeneous Reservoirs; The Output Least Squares (OLS) Method. Transport Porous Media 40 (2): 115-143. http://dx.doi.org/10.1023/a:1006657515753.
↑ ^81.0 ^81.1 Samier, P., Quettier, L., and Thiele, M. 2002. Applications of Streamline Simulations to Reservoir Studies. SPE Res Eval & Eng 5 (4): 324–332. SPE-78883-PA. http://dx.doi.org/10.2118/78883-PA.
↑ Thiele, M.R. 2001. Streamline Simulation. Keynote address at the 2001 Intl. Forum on Reservoir Simulation, Schloss Fuschl, Austria, 3–7 September.
↑ Lolomari, T., Bratvedt, K., Crane, M. et al. 2000. The Use of Streamline Simulation in Reservoir Management: Methodology and Case Studies. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 1–4 October. SPE-63157-MS. http://dx.doi.org/10.2118/63157-MS.
↑ ^84.0 ^84.1 Baker, R.O., Kuppe, F., Chug, S. et al. 2001. Full-Field Modeling Using Streamline-Based Simulation: 4 Case Studies. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 11–14 February. SPE-66405-MS. http://dx.doi.org/10.2118/66405-MS.
↑ Thiele, M.R., Batycky, R.P., and Blunt, M.J. 1997. A Streamline-Based 3D Field-Scale Compositional Reservoir Simulator. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 5–8 October. SPE-38889-MS. http://dx.doi.org/10.2118/38889-MS.
↑ Thiele, M.R., Batycky, R.P., and Thomas, L.K. 2002. Miscible WAG Simulations Using Streamlines. Paper presented at the 2002 European Conference on the Mathematics of Oil Recovery, Freiberg, Germany, 3–6 September.
↑ ^87.0 ^87.1 ^87.2 ^87.3 Batycky, R.P., Blunt, M.J., and Thiele, M.R. 1997. A 3D Field-Scale Streamline-Based Reservoir Simulator. SPE Res Eng 12 (4): 246-254. SPE-36726-PA. http://dx.doi.org/10.2118/36726-PA.
↑ Emanuel, A.S. 1993. Development of an Analytical Streamtube Model for Estimating Steam Drive Performance. SPE Advanced Technology Series 1 (1): 81-89. SPE-21756-PA. http://dx.doi.org/10.2118/21756-PA.
↑ Crane, M.J. and Blunt, M.J. 1999. Streamline-based simulation of solute transport. Water Resour. Res. 35 (10): 3061-3078. http://dx.doi.org/10.1029/1999wr900145.
↑ Di Donato, G., Huang, W., and Blunt, M. 2003. Streamline-Based Dual Porosity Simulation of Fractured Reservoirs. Presented at the SPE Annual Technical Conference and Exhibition, 5–8 October. SPE-84036-MS. http://dx.doi.org/10.2118/84036-MS.
↑ Muskat, M. and Wyckoff, R.D. 1934. A Theoretical Analysis of Waterflooding Networks. Trans., AIME, 107: 62–77.
↑ Muskat, M. 1982. Flow of Homogeneous Fluids. Boston, Massachusetts: Intl. Human Resources Development Corp.
↑ Fay, C.H. and Prats, M. 1951. The Application of Numerical Methods to Cycling and Flooding Problems. Proc., Third World Petroleum Congress, The Hague, 28 May–6 June.
↑ Higgins, R.V. and Leighton, A.J. 1962. A Computer Method to Calculate Two-Phase Flow in Any Irregularly Bounded Porous Medium. J Pet Technol 14 (6): 679-683. SPE-243-PA. http://dx.doi.org/10.2118/243-PA.
↑ Higgins, R.V. and Leighton, A.J. 1962. Computer Prediction of Water Drive of Oil and Gas Mixtures Through Irregularly Bounded Porous Media Three-Phase Flow. J Pet Technol 14 (9): 1048-1054. SPE-283-PA. http://dx.doi.org/10.2118/283-PA.
↑ Bommer, P.M. and Schechter, R.S. 1979. Mathematical Modeling of In-Situ Uranium Leaching. Society of Petroleum Engineers Journal 19 (6): 393-400. SPE-7533-PA. http://dx.doi.org/10.2118/7533-PA.
↑ Martin, J.C. and Wegner, R.E. 1979. Numerical Solution of Multiphase, Two-Dimensional Incompressible Flow Using Stream-Tube Relationships. Society of Petroleum Engineers Journal 19 (5): 313-323. SPE-7140-PA. http://dx.doi.org/10.2118/7140-PA.
↑ Lake, L.W., Johnston, J.R., and Stegemeier, G.L. 1981. Simulation and Performance Prediction of a Large-Scale Surfactant/Polymer Project. Society of Petroleum Engineers Journal 21 (6): 731-739. SPE-7471-PA. http://dx.doi.org/10.2118/7471-PA.
↑ Emanuel, A.S., Alameda, G.K., Behrens, R.A. et al. 1989. Reservoir Performance Prediction Methods Based on Fractal Geostatistics(includes associated papers 20011 and 20158 ). SPE Res Eng 4 (3): 311-318. SPE-16971-PA. http://dx.doi.org/10.2118/16971-PA.
↑ ^100.0 ^100.1 Pollock, D.W. 1988. Semianalytical Computation of Path Lines for Finite-Difference Models. Ground Water 26 (6): 743.
↑ King, M.J., Blunt, M.J., Mansfield, M.M. et al. 1993. Rapid Evaluation of the Impact of Heterogeneity on Miscible Gas Injection. Presented at the SPE Western Regional Meeting, Anchorage, Alaska, 26-28 May 1993. SPE-26079-MS. http://dx.doi.org/10.2118/26079-MS.
↑ Datta-Gupta, A. and King, M.J. 1995. A semianalytic approach to tracer flow modeling in heterogeneous permeable media. Adv. Water Resour. 18 (1): 9-24. http://dx.doi.org/http://dx.doi.org/10.1016/0309-1708(94)00021-V.
↑ ^103.0 ^103.1 Batycky, R.P. 1997. A Three-Dimensional Two-Phase Field Scale Streamline Simulator. PhD dissertation, Stanford U., Stanford, California.
↑ Renard, G. 1990. A 2D reservoir streamtube EOR model with periodical automatic regeneration of streamlines. In Situ 14 (2): 175-200.
↑ Thiele, M.R., Batycky, R.P., Blunt, M.J. et al. 1996. Simulating Flow in Heterogeneous Systems Using Streamtubes and Streamlines. SPE Res Eng 11 (1): 5–12. SPE-27834-PA. http://dx.doi.org/10.2118/27834-PA.
↑ ^106.0 ^106.1 Bratvedt, F., Gimse, T., and Tegnander, C. 1996. Streamline computations for porous media flow including gravity. Transport Porous Media 25 (1): 63-78. http://dx.doi.org/10.1007/bf00141262.
↑ Blunt, M.J., Liu, K., and Thiele, M.R. 1996. A Generalized Streamline Method to Predict Reservoir Flow. Petroleum Geoscience 2 (3): 259-269. http://dx.doi.org/10.1144/petgeo.2.3.259.
↑ Prévost, M., Edwards, M.G., and Blunt, M.J. 2002. Streamline Tracing on Curvilinear Structured and Unstructured Grids. SPE J. 7 (2): 139–148. SPE-78663-PA. http://dx.doi.org/10.2118/786637-PA.
↑ Thiele, M.R. and Batycky, R.P. 2003. Water Injection Optimization Using a Streamline-Based Workflow. Presented at the SPE Annual Technical Conference and Exhibition, Denver, 5-8 October. SPE 84080. http://dx.doi.org/10.2118/84080-MS.
↑ Batycky, R.P., Thiele, M.R., Baker, R.O. et al. 2005. Revisiting Reservoir Flood-Surveillance Methods Using Streamlines. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 9-12 October 2005. SPE-95402-MS. http://dx.doi.org/10.2118/95402-MS.
↑ Wang, Y. and Kovscek, A.R. 2000. Streamline Approach for History Matching Production Data. SPE J. 5 (4): 353-362. SPE-58350-PA. http://dx.doi.org/10.2118/58350-PA.
↑ Agarwal, B. and Blunt, M.J. 2003. Streamline-Based Method With Full-Physics Forward Simulation for History-Matching Performance Data of a North Sea Field. SPE J. 8 (2): 171-180. SPE-84952-PA. http://dx.doi.org/10.2118/84952-PA.
↑ Caers, J., Krishnan, S., Wang, Y. et al. 2002. A Geostatistical Approach to Streamline-Based History Matching. SPE J. 7 (3): 250-266. SPE-73144-PA. http://dx.doi.org/10.2118/73144-PA.
↑ Fenwick, D., Thiele, M., Agil, M. et al. 2005. Reconciling Prior Geologic Information With Production Data Using Streamlines: Application to a Giant Middle-Eastern Oil Field. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 9–12 October. SPE-95940-MS. http://dx.doi.org/10.2118/95940-MS.
↑ Stüben, K. 2000. Algebraic Multigrid (AMG): An Introduction with Applications. Multigrid, U. Trottenberg, C.W. Oosterlee, and A. Schüller eds., appendix. New York City: Academic Press.
↑ Darve, F. 1990. Geomaterials: Constitutive Equations and Modeling. New York City: Elsevier Applied Science.
↑ Coussy, O. 1995. Mechanics of Porous Continua. Chichester, England: John Wiley and Sons.
↑ Settari, A. and Walters, D.A. 2001. Advances in Coupled Geomechanical and Reservoir Modeling With Applications to Reservoir Compaction. SPE J. 6 (3): 334–342. SPE-51927-PA. http://dx.doi.org/10.2118/51927-MS.
↑ ^119.0 ^119.1 Lorenz, J.C. 1999. Stress Sensitive Reservoirs. J Pet Technol (January): 61.
↑ Settari, A., Warren, G.M., Jacquemont, J. et al. 1999. Brine Disposal Into a Tight Stress-Sensitive Formation at Fracturing Conditions: Design and Field Experience. SPE Res Eval & Eng 2 (2): 186-195. SPE-56001-PA. http://dx.doi.org/10.2118/56001-PA.
↑ Settari, A. 2002. Reservoir Compaction. Distinguished Author Series, J Pet Technol 54 (8): 62–69. SPE-76805-PA. http://dx.doi.org/10.2118/76805-MS.
↑ Finol, A. and Farouq Ali, S.M. 1975. Numerical Simulation of Oil Production With Simultaneous Ground Subsidence. Society of Petroleum Engineers Journal 15 (5): 411-424. SPE-4847-PA. http://dx.doi.org/10.2118/4847-PA.
↑ Merle, H.A., Kentie, C.J.P., van Opstal, G.H.C. et al. 1976. The Bachaquero Study - A Composite Analysis of the Behavior of a Compaction Drive/Solution Gas Drive Reservoir. J Pet Technol 28 (9): 1107-1115. SPE-5529-PA. http://dx.doi.org/10.2118/5529-PA.
↑ Rattia, A.J. and Farouq Ali, S.M. 1981. Effect of Formation Compaction on Steam Injection Response. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 4-7 October 1981. SPE-10323-MS. http://dx.doi.org/10.2118/10323-MS.
↑ Chase, C.A.J. and Dietrich, J.K. 1989. Compaction Within the South Belridge Diatomite. SPE Res Eng 4 (4): 422–428. SPE-17415-PA. http://dx.doi.org/10.2118/17415-PA.
↑ Beattie, C.I., Boberg, T.C., and McNab, G.S. 1991. Reservoir Simulation of Cyclic Steam Stimulation in the Cold Lake Oil Sands. SPE Res Eng 6 (2): 200–206. SPE-18752-PA. http://dx.doi.org/10.2118/18752-PA.
↑ ^127.0 ^127.1 Chin, L.Y. and Thomas, L.K. 1999. Fully Coupled Analysis of Improved Oil Recovery by Reservoir Compaction. Presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, 3–6. October. SPE-56753-MS. http://dx.doi.org/10.2118/56753-MS.
↑ ^128.0 ^128.1 Settari, A. and Mourits, F.M. 1998. A Coupled Reservoir and Geomechanical Simulation System. SPE J. 3 (3): 219–226. SPE-50939-PA. http://dx.doi.org/10.2118/50939-PA
↑ Settari, A. and Price, H.S. 1984. Simulation of Hydraulic Fracturing in Low-Permeability Reservoirs. Society of Petroleum Engineers Journal 24 (2): 141-152. SPE-8939-PA. http://dx.doi.org/10.2118/8939-PA.
↑ Settari, A.T., Bachman, R.C., and D.A.Walters. 2005. How To Approximate Effects of Geomechanics in Conventional Reservoir Simulation. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, 9-12 October. SPE 97155. http://dx.doi.org/10.2118/97155-MS.
↑ Heffer, K.J. et al. 1992. The Influence of Natural Fractures, Faults and Earth Stresses on Reservoir Performance—Geomechanical Analysis by Numerical Modelling. Proc., Third Intl. Conference on North Sea Oil and Gas Reservoirs, Trondheim, Norway, 201–211.
↑ Jin, M., Somerville, J., and Smart, B.G.D. 2000. Coupled Reservoir Simulation Applied to the Management of Production Induced Stress-Sensitivity. Presented at the International Oil and Gas Conference and Exhibition in China, Beijing, China, 7–10 November. SPE-64790-MS. http://dx.doi.org/10.2118/64790-MS.
↑ Coupled Thermo-Hydro-Mechanical Processes of Fractured Media. 1996. Developments in Geotechnical Engineering, O. Stephansson, L. Jing, and S.-F. Tsang eds., 575. New York City: Elsevier.
↑ Koutsabeloulis, N.C. and Hope, S.A. 1998. “Coupled” Stress/Fluid/Thermal Multi-Phase Reservoir Simulation Studies Incorporating Rock Mechanics. Presented at the SPE/ISRM Rock Mechanics in Petroleum Engineering, Trondheim, Norway, 8–10 July. SPE-47393-MS. http://dx.doi.org/10.2118/47393-MS.
↑ Settari, A. 2000. Coupled Fracture and Reservoir Modeling. Proc., Advanced Hydraulic Fracturing Workshop, 2000 North American Rock Mechanics Conference, Seattle, Washington.
↑ Lewis, R.W. and Schrefler, B.A. 1987. The Finite Element Method in the Deformation and Consolidation of Porous Media, 344. New York City: John Wiley and Sons.
↑ Stone, T., Bowen, G., Papanastasiou, P. et al. 2000. Fully Coupled Geomechanics in a Commercial Reservoir Simulator. Presented at the SPE European Petroleum Conference, Paris, France, 24–25 October. SPE-65107-MS. http://dx.doi.org/10.2118/65107-MS.
↑ Sulak, R.M., Thomas, L.K., and Boade, R.R. 1991. 3D Reservoir Simulation of Ekofisk Compaction Drive (includes associated papers 24317 and 24400 ). J Pet Technol 43 (10): 1272-1278. SPE-19802-PA. http://dx.doi.org/10.2118/19802-PA.
↑ Fredrich, J.T. et al. 2001. One-Way Coupled Reservoir-Geomechanical Modeling of the Lost Hills Oil Field, California. Proc., 38th U.S. Rock Mechanics Symposium, Washington, DC.
↑ Tran, D., Settari, A., and Nghiem, L. 2002. New Iterative Coupling Between a Reservoir Simulator and a Geomechanics Module. Presented at the SPE/ISRM Rock Mechanics Conference, Irving, Texas, 20-23 October 2002. SPE-78192-MS. http://dx.doi.org/10.2118/78192-MS.
↑ Dean, R.H. 2000. A Poroelastic Multicomponent Reservoir Simulator. Paper SPE 62806 available from SPE, Richardson, Texas.
↑ Settari, A., Walters, D.A., and Behie, G.A. 2000. Use of Coupled Reservoir and Geomechanical Modelling for Integrated Reservoir Analysis and Management. Presented at the Canadian International Petroleum Conference, Calgary, Alberta, Jun 4 - 8, 2000 2000. PETSOC-2000-078. http://dx.doi.org/10.2118/2000-078.
↑ ^143.0 ^143.1 ^143.2 Coats, B.K., Nolen, J.S., and Sherman, A.H. 1988. Applicability and Accuracy of Two Component Characterizations of Reservoir Oils for Numerical Modeling Purposes. Presented at the International Meeting on Petroleum Engineering, Tianjin, China, 1-4 November 1988. SPE-17614-MS. http://dx.doi.org/10.2118/17614-MS.
↑ ^144.0 ^144.1 ^144.2 Cook, R.E., Jacoby, R.H., and Ramesh, A.B. 1974. A Beta-Type Reservoir Simulator for Approximating Compositional Effects During Gas Injection. Society of Petroleum Engineers Journal 14 (5): 471-481. SPE-4272-PA. http://dx.doi.org/10.2118/4272-PA.
↑ ^145.0 ^145.1 ^145.2 ^145.3 Fevang, Ø. and Whitson, C.H. 1996. Modeling Gas-Condensate Well Deliverability. SPE Res Eng 11 (4): 221–230. SPE-30714-PA. http://dx.doi.org/10.2118/30714-PA.
↑ Spivak, A. and Dixon, T.N. 1973. Simulation of Gas-Condensate Reservoirs. Presented at the SPE Symposium on Numerical Simulation of Reservoir Performance, Houston, Texas, 11-12 January 1973. SPE-4271-MS. http://dx.doi.org/10.2118/4271-MS.
↑ Chang, Y.-B., Coats, B.K., and Nolen, J.S. 1998. A Compositional Model for CO2 Floods Including CO2 Solubility in Water. SPE Res Eval & Eng 1 (2): 155-160. SPE 35164. http://dx.doi.org/10.2118/35164-PA.
↑ ^148.0 ^148.1 Michelsen, M.L. 1985. Saturation point calculations. Fluid Phase Equilib. 23 (2–3): 181-192. http://dx.doi.org/http://dx.doi.org/10.1016/0378-3812(85)90005-6.
↑ Michelsen, M.L. 1982. The isothermal flash problem. Part I. Stability. Fluid Phase Equilib. 9 (1): 1-19. http://dx.doi.org/http://dx.doi.org/10.1016/0378-3812(82)85001-2.
↑ Muskat, M. and McDowell, J.M. 1949. An Electrical Computer for Solving Phase Equilibrium Problems. J Pet Technol 1 (11): 291-298. http://dx.doi.org/10.2118/949291-G.
↑ Rachford Jr., H.H. and Rice, J.D. 1952. Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium. J Pet Technol 4 (10): 19, 3. http://dx.doi.org/10.2118/952327-G.
↑ ^152.0 ^152.1 ^152.2 ^152.3 ^152.4 ^152.5 ^152.6 Whitson, C.H. and Brule, M.R. 2000. Phase Behavior, Vol. 20. Richardson, Texas: Monograph Series, SPE.
↑ ^153.0 ^153.1 Michelsen, M.L. 1982. The isothermal flash problem. Part II. Phase-split calculation. Fluid Phase Equilib. 9 (1): 21-40. http://dx.doi.org/http://dx.doi.org/10.1016/0378-3812(82)85002-4.
↑ ^154.0 ^154.1 Péneloux, A., Rauzy, E., and Fréze, R. 1982. A consistent correction for Redlich-Kwong-Soave volumes. Fluid Phase Equilib. 8 (1): 7-23. http://dx.doi.org/http://dx.doi.org/10.1016/0378-3812(82)80002-2.
↑ ^155.0 ^155.1 Peng, D.Y. and Robinson, D.B. 1976. A New-Constant Equation of State. Ind. Eng. Chem. Fund. 15 (1): 59-64. http://dx.doi.org/10.1021/i160057a011.
↑ ^156.0 ^156.1 Soave, G. 1972. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci. 27 (6): 1197-1203. http://dx.doi.org/http://dx.doi.org/10.1016/0009-2509(72)80096-4.
↑ ^157.0 ^157.1 Zudkevitch, D. and Joffe, J. 1970. Correlation and prediction of vapor-liquid equilibria with the redlich-kwong equation of state. AIChE J. 16 (1): 112-119. http://dx.doi.org/10.1002/aic.690160122.
↑ Whitson, C.H. and Torp, S.B. 1983. Evaluating Constant-Volume Depletion Data. J Pet Technol 35 (3): 610-620. SPE-10067-PA. http://dx.doi.org/10.2118/10067-PA.
↑ ^159.0 ^159.1 Cook, A.B., Spencer, G.B., and Bayazeed, A.F. 1967. The Role of Vaporization in High Percentage Oil Recovery by Pressure Maintenance. J Pet Technol 19 (2): 245-250. SPE-1646-PA. http://dx.doi.org/10.2118/1646-PA.
↑ Whitson, C.H. 1983. Characterizing Hydrocarbon Plus Fractions. SPE J. 23 (4): 683–693. SPE-12233-PA. http://dx.doi.org/10.2118/12233-PA.
↑ Tang, D.E. and Zick, A.A. 1993. A New Limited Compositional Reservoir Simulator. Presented at the SPE Symposium on Reservoir Simulation, New Orleans, Louisiana, 28 february–3 March. SPE-25255-MS. http://dx.doi.org/10.2118/25255-MS.
↑ ^162.0 ^162.1 Zick, A.A. 1986. A Combined Condensing/Vaporizing Mechanism in the Displacement of Oil by Enriched Gases. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, 5–8 October. SPE 15493. http://dx.doi.org/10.2118/15493-MS.
↑ Fevang, Ø., Singh, K., and Whitson, C.H. 2000. Guidelines for Choosing Compositional and Black-Oil Models for Volatile Oil and Gas-Condensate Reservoirs. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 1-4 October 2000. SPE-63087-MS. http://dx.doi.org/10.2118/63087-MS.
↑ Lohrenz, J., Bray, B.G., and Clark, C.R. 1964. Calculating Viscosities of Reservoir Fluids From Their Compositions. J Pet Technol 16 (10): 1171–1176. SPE-915-PA. http://dx.doi.org/10.2118/915-PA.
↑ Pedersen, K.S. and Fredenslund, A. 1987. An Improved Corresponding States Model for the Prediction of Oil and Gas Viscosities and Thermal Conductivities. Chem. Eng. Sci. 42: 182.
↑ Saidi, A.M. 1996. Twenty Years of Gas Injection History into Well-Fractured Haft Kel Field (Iran). Presented at the International Petroleum Conference and Exhibition of Mexico, Villahermosa, Mexico, 5–7 March. SPE 35309. http://dx.doi.org/10.2118/35309-MS.
↑ Whitson, C.H., Fevang, Ø., and Sævareid, A. 1999. Gas Condensate Relative Permeability for Well Calculations. Presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, 3–6 October. SPE-56476-MS. http://dx.doi.org/10.2118/56476-MS.
↑ Mott, R.E., Cable, A.S., and Spearing, M.C. 2000. Measurements of Relative Permeabilities for Calculating Gas-Condensate Well Deliverability. SPE Res Eval & Eng 3 (6): 473–479. SPE-68050-PA. http://dx.doi.org/10.2118/68050-PA.
↑ Robinson, D.B. and Peng, D.Y. 1978. The Characterization of the Heptanes and Heavier Fractions. Research Report 28. Tulsa, Oklahoma: Gas Producers Assn.
↑ Katz, D.L. and Firoozabadi, A. 1978. Predicting Phase Behavior of Condensate/Crude-Oil Systems Using Methane Interaction Coefficients. J Pet Technol 30 (11): 1649–1655. SPE-6721-PA. http://dx.doi.org/10.2118/6721-PA.
↑ Shelton, J.L. and Yarborough, L. 1977. Multiple Phase Behavior in Porous Media During CO2 or Rich-Gas Flooding. J Pet Technol 29 (9): 1171-1178. SPE-5827-PA. http://dx.doi.org/10.2118/5827-PA.
↑ Nghiem, L.X. and Li, Y.K. 1986. Effect of Phase Behavior on CO2 Displacement Efficiency at Low Temperatures: Model Studies With an Equation of State. SPE Res Eng 1 (4): 414–422. SPE-13116-PA. http://dx.doi.org/10.2118/13116-PA.
↑ Coats, K.H. 1976. Simulation of Steamflooding With Distillation and Solution Gas. SPE J. 16 (5): 235–247. SPE-5015-PA. http://dx.doi.org/10.2118/5015-PA.
↑ Killough, J.E. 1979. The Use of Vector Processors in Reservoir Simulation. Presented at the SPE Reservoir Simulation Symposium, Denver, Colorado, 31 January-2 Feburary 1979. SPE-7673-MS. http://dx.doi.org/10.2118/7673-MS.
↑ Nolen, J.S., Kuba, D.W., and Jr., M.J.K. 1981. Application of Vector Processors To Solve Finite Difference Equations. Society of Petroleum Engineers Journal 21 (4): 447-453. SPE-7675-PA. http://dx.doi.org/10.2118/7675-PA.
↑ Killough, J.E. and Levesque, J.M. 1982. Reservoir Simulation and the In-house Vector Processor: Experience for the First Year. Presented at the SPE Reservoir Simulation Symposium, New Orleans, Louisiana, 31 January-3 Feburary 1982. SPE-10521-MS. http://dx.doi.org/10.2118/10521-MS.
↑ Young, L.C. 1987. Equation of State Compositional Modeling on Vector Processors. Presented at the SPE Symposium on Reservoir Simulation, San Antonio, Texas, 1-4 February 1987. SPE-16023-MS. http://dx.doi.org/10.2118/16023-MS.
↑ Killough, J.E. 1993. Is Parallel Computing Ready for Reservoir Simulation? A Critical Analysis of the State of the Art. Presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, 3-6 October 1993. SPE-26634-MS. http://dx.doi.org/10.2118/26634-MS.
↑ ^179.0 ^179.1 ^179.2 ^179.3 ^179.4 Killough, J.E. et al. 1996. A General-Purpose parallel Reservoir Simulator. Presented at the 1996 European Conference on the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.
↑ van Daalen, D.T. et al. 1990. The parallelization of BOSIM, Shell’s Black/Volatile Oil Reservoir Simulator. Proc. of the First IMA/SPE European Conference on the Mathematics of Oil Recovery, Oxford U.
↑ Killough, J.E. and Bhogeswara, R. 1991. Simulation of Compositional Reservoir Phenomena on a Distributed-Memory Parallel Computer. J Pet Technol 43 (11): 1368-1374. SPE-21208-PA. http://dx.doi.org/10.2118/21208-PA.
↑ Wheeler, J.A. and Smith, R.A. 1990. Reservoir Simulation on a Hypercube. SPE Res Eng 5 (4): 544-548. SPE-19804-PA. http://dx.doi.org/10.2118/19804-PA.
↑ ^183.0 ^183.1 Rutledge, J.M., Jones, D.R., Chen, W.H. et al. 1992. Use of a Massively Parallel SIMD Computer for Reservoir Simulation. SPE Comp App 4 (4): 16-21. SPE-21213-PA. http://dx.doi.org/10.2118/21213-PA.
↑ Gautam S. et al. 1996. Parallel Computing Alters Approaches, Raises Integration Challenges in Reservoir Modeling. Oil & Gas J. (20 May): 48.
↑ Tchelepi, H.A., Durlofsky, L.J., Chen, W.H. et al. 1997. Practical Use of Scale Up and Parallel Reservoir Simulation Technologies in Field Studies. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 5-8 October 1997. SPE-38886-MS. http://dx.doi.org/10.2118/38886-MS.
↑ Chien, M.C.H. et al.: "A Scalable parallel Multi-Purpose Reservoir Simulator," paper SPE 37976 presented at the 1997 SPE Reservoir Simulation Symposium, Dallas, 8–11 June.
↑ Shiralkar, G.S., Stephenson, R.E., Joubert, W. et al. 1997. Falcon: A Production Quality Distributed Memory Reservoir Simulator. Presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, 8-11 June 1997. SPE-37975-MS. http://dx.doi.org/10.2118/37975-MS.
↑ Wang, P., Yotov, I., Wheeler, M. et al. 1997. A New Generation EOD Compositional Reservoir Simulator: Part 1—Formulation and Discretization. Presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, 8–11 June. SPE-37979-MS. http://dx.doi.org/10.2118/37979-MS.
↑ ^189.0 ^189.1 ^189.2 Pavlas Jr., E.J. 2001. MPP Simulation of Complex Water Encroachment in a Large Carbonate Reservoir in Saudi Arabia. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71628-MS. http://dx.doi.org/10.2118/71628-MS.
↑ ^190.0 ^190.1 Wallis, J.R. and Nolen, J.S. 1993. Efficient Iterative Linear Solution of Locally Refined Grids Using Algebraic Multilevel Approximate Factorizations. Presented at the SPE Symposium on Reservoir Simulation, New Orleans, Louisiana, 28 February-3 March 1993. SPE-25239-MS. http://dx.doi.org/10.2118/25239-MS.
↑ Song, T. 1996. A Load-Balancing Technique for Reservoir Simulation Based on Dantzig’s Transportation Model. MS thesis, U. of Houston, Houston.
↑ Reference Manual for Zoltan. 1998. Los Alamos Natl. Laboratory, Los Alamos, New Mexico.
↑ Bhogeswara, R. and Killough, J.E. 1994. Parallel Linear Solvers for Reservoir Simulation: A Generic Approach for Existing and Emerging Computer Architectures. SPE Comp App 6 (1): 5-11. SPE-25240-PA. http://dx.doi.org/10.2118/25240-PA.
↑ Guerberoff, D., Zucchi, H., Victoria, M. et al. 2001. Lateral Delineation of Sandstone Bodies Guided by Seismic and Petrophysical Data Using Geocellular Model: Geocellular Model: Canadon Secon Formation, San Jorge Basin, Argentina. Presented at the SPE Latin American and Caribbean Petroleum Engineering Conference, Buenos Aires, Argentina, 25-28 March 2001. SPE-69487-MS. http://dx.doi.org/10.2118/69487-MS.
↑ Agarwal, B., Hermansen, H., Sylte, J.E. et al. 2000. Reservoir Characterization of Ekofisk Field: A Giant, Fractured Chalk Reservoir in the Norwegian North SeaHistory Match. SPE Res Eval & Eng 3 (6): 534-543. SPE-68096-PA. http://dx.doi.org/10.2118/68096-PA.
↑ Caers, J., Avseth, P., and Mukerji, T. 2001. Geostatistical integration of rock physics, seismic amplitudes and geological models in North-Sea turbidite systems. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71321-MS. http://dx.doi.org/10.2118/71321-MS.
↑ ^197.0 ^197.1 Qassab, H.M.A., Rahmeh, B.A., Khalifa, M.A.A. et al. 2001. Conditioning Integrated Geological Models to Dynamic Flow Data of Giant Saudi Arabian Reservoir. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71319-MS. http://dx.doi.org/10.2118/71319-MS.
↑ ^198.0 ^198.1 Raghavan, R., Dixon, T.N., Robinson, S.W. et al. 2000. Integration of Geology, Geophysics, and Numerical Simulation in the Interpretation of A Well Test in a Fluvial Reservoir. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 1-4 October 2000. SPE-62983-MS. http://dx.doi.org/10.2118/62983-MS.
↑ Durlofsky, L.J. 2005. Upscaling and gridding of fine-scale geological models for flow simulation. Proc., 2005 Intl. Forum on Reservoir Simulation, Stresa, Italy, 20–25 June.
↑ Peaceman, D.W. 1997. Effective Transmissibilities of a Gridblock by Upscaling - Comparison of Direct Methods with Renormalization. SPE J. 2 (3): 338-349. SPE-36722-PA. http://dx.doi.org/10.2118/36722-PA.
↑ Flodin, E.A., Aydin, A., Durlofsky, L.J. et al. 2001. Representation of fault zone permeability in reservoir flow models. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September–3 October. SPE-71617-MS. http://dx.doi.org/10.2118/71617-MS.
↑ Nishikiori, N. and Hayashida, Y. 2000. Investigation of Fluid Conductive Faults and Modeling of Complex Water Influx in the Khafji Oil Field, Arabian Gulf. SPE Res Eval & Eng 3 (5): 401-407. SPE-66223-PA. http://dx.doi.org/10.2118/66223-PA.
↑ Barker, J.W. and Thibeau, S. 1997. A Critical Review of the Use of Pseudorelative Permeabilities for Upscaling. SPE Res Eng 12 (2): 138–143. SPE-35491-PA. http://dx.doi.org/10.2118/35491-PA.
↑ Beckner, B.L., Hutfilz, J.M., Ray, M.B. et al. 2001. EMpower: New Reservoir Simulation System. Presented at the SPE Middle East Oil Show, Bahrain, 17–20 March. SPE-68116-MS. http://dx.doi.org/10.2118/68116-MS.
↑ Peaceman, D.W. 1983. Interpretation of Well-Block Pressures in Numerical Reservoir Simulation With Nonsquare Grid Blocks and Anisotropic Permeability. SPE J. 23 (3): 531–543. SPE-10528-PA. http://dx.doi.org/10.2118/10528-PA.
↑ Mlacnik, M.J. and Heinemann, Z.E. 2003. Using Well Windows in Full Field Reservoir Simulation. SPE Res Eval & Eng 6 (4): 275–285. SPE-85709-PA. http://dx.doi.org/10.2118/85709-PA.
↑ Wolfsteiner, C., Durlofsky, L.J., and Aziz, K. 2000. Approximate Model for Productivity of Nonconventional Wells in Heterogeneous Reservoirs. SPE J. 5 (2): 218–226. SPE-56754-PA. http://dx.doi.org/10.2118/62812-PA.
↑ Valvatne, P.H., Durlofsky, L.J., and Aziz, K. 2001. Semi-Analytical Modeling of the Performance of Intelligent Well Completions. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 11-14 February 2001. SPE-66368-MS. http://dx.doi.org/10.2118/66368-MS.
↑ Litvak, M.L., Clark, A.J., Fairchild, J.W. et al. 1997. Integration of Prudhoe Bay Surface Pipeline Network and Full Field Reservoir Models. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, 5-8 October. SPE 38885. http://dx.doi.org/10.2118/38885-MS.
↑ Fang, W.Y. and Lo, K.K. 1996. A Generalized Well-Management Scheme for Reservoir Simulation. SPE Res Eng 11 (2): 116–120. SPE-29124-PA. http://dx.doi.org/10.2118/29124-PA.
↑ Litvak, M.L. and Wang, C.H. 2000. Simplified Phase-Equilibrium Calculations in Integrated Reservoir and Surface-Pipeline-Network Models. SPE J. 5 (2): 236-241. SPE-64498-PA. http://dx.doi.org/10.2118/64498-PA.
↑ Hermansen, H., Thomas, L.K., Sylte, J.E. et al. 1997. Twenty Five years of Ekofisk Reservoir Management. Presented at the Annual Technical Conference and Exhibition, San Antonio, Texas, 5–8 October. SPE-38927-MS. http://dx.doi.org/10.2118/38927-MS.
↑ Sylte, J.E., Thomas, L.K., Rhett, D.W. et al. 1999. Water Induced Compaction in the Ekofisk Field. Presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, 3–6 October. SPE-56426-MS. http://dx.doi.org/10.2118/56426-MS.
↑ Ashton, K., Sylte, J.E., Thomas, L.K. et al. 1998. Judy/Joanne Field Development. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 27-30 September 1998. SPE-49128-MS. http://dx.doi.org/10.2118/49128-MS.
↑ Pieters, J. and Por, G.J.A. 1995. Total system modelling - a tool for effective reservoir management of multiple fields with shared facilities. Presented at the Offshore Europe, Aberdeen, United Kingdom, 5-8 September 1995. SPE-30442-MS. http://dx.doi.org/10.2118/30442-MS.
↑ Williams, M.A., Keating, J.F., and Barghouty, M.F. 1998. The Stratigraphic Method: A Structured Approach to History-Matching Complex Simulation Models. SPE Res Eval & Eng 1 (2): 169-176. SPE-38014-PA. http://dx.doi.org/10.2118/38014-PA.
↑ Ditzhuijzen, R.v., Oldenziel, T., and Kruijsdijk, C.P.J.W.v. 2001. Geological Parameterization of a Reservoir Model for History Matching Incorporating Time-Lapse Seismic Based on a Case Study of the Statfjord Field. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71318-MS. http://dx.doi.org/10.2118/71318-MS.
↑ Bogan, C., Johnson, D., Litvak, M. et al. 2003. Building Reservoir Models Based on 4D Seismic & Well Data in Gulf of Mexico Oil Fields. Presented at the SPE Annual Technical Conference and Exhibition, Denver, Colorado, 5-8 October 2003. SPE-84370-MS. http://dx.doi.org/10.2118/84370-MS.
↑ Milliken, W.J., Emanuel, A.S., and Chakravarty, A. 2000. Applications of 3D Streamline Simulation to Assist History Matching. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 1-4 October 2000. SPE-63155-MS. http://dx.doi.org/10.2118/63155-MS.
↑ Schulze-Riegert, R.W., Axmann, J.K., Haase, O. et al. 2001. Optimization Methods for History Matching of Complex Reservoirs. Presented at the SPE Reservoir Simulation Symposium, Houston, 11-14 February. SPE 66393. http://dx.doi.org/10.2118/66393-MS.
↑ Grussaute, T. and Gouel, P. 1998. Computer Aided History Matching of a Real Field Case. Presented at the European Petroleum Conference, The Hague, Netherlands, 20-22 October 1998. SPE-50642-MS. http://dx.doi.org/10.2118/50642-MS.
↑ Zapata, V.J., Brummett, W.M., Osborne, M.E. et al. 2001. Advances in Tightly Coupled Reservoir/Wellbore/Surface-Network Simulation. SPE Res Eval & Eng 4 (2): 114–120. SPE-71120-PA. http://dx.doi.org/10.2118/71120-PA.
↑ ^223.0 ^223.1 Gorell, S. and Bassett, R. 2001. Trends in Reservoir Simulation: Big Models, Scalable Models? Will you Please Make up Your Mind? Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71596-MS. http://dx.doi.org/10.2118/71596-MS.
↑ Zabalza-Mezghani, I., Mezghani, M., and Blanc, G. 2001. Constraining Reservoir Facies Models to Dynamic Data - Impact of Spatial Distribution Uncertainty on Production Forecasts. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71335-MS. http://dx.doi.org/10.2118/71335-MS.
↑ Manceau, E., Mezghani, M., Zabalza-Mezghani, I. et al. 2001. Combination of Experimental Design and Joint Modeling Methods for Quantifying the Risk Associated With Deterministic and Stochastic Uncertainties—An Integrated Test Study. Presented at the SPE Annual Conference and Technical Exhibition, New Orleans, Louisiana, 30 September–3 October. SPE-71620-MS. http://dx.doi.org/10.2118/71620-MS.

SI Metric Conversion Factors

bbl	×	1.589 873	E – 01	=	m³
cp	×	1.0*	E – 03	=	Pa•s
dyne	×	1.0*	E – 02	=	mN
ft	×	3.048*	E – 01	=	m
ft³	×	2.831 685	E − 02	=	m³
°F		(°F − 32)/1.8		=	°C
lbm	×	4.535 924	E – 01	=	kg

*

Conversion factor is exact.

[r1-1] The Merriam-Webster Dictionary, New Revised Edition. 2004. New York City: Merriam-Webster.

[r2-2] Coats, K.H. 1987. Reservoir Simulation. Petroleum Engineering Handbook, H.B. Bradley ed., Chap. 48. Richardson, Texas: SPE.

[r3-3] Odeh, A.S. 1981. Comparison of Solutions to a Three-Dimensional Black-Oil Reservoir Simulation Problem. J Pet Technol 33 (1): 13–25. SPE-9723-PA. http://dx.doi.org/10.2118/9723-PA.

[r4-4] Weinstein, H.G., Chappelear, J.E., and Nolen, J.S. 1986. Second Comparative Solution Project: A Three-Phase Coning Study. J Pet Technol 38 (3): 345-353. SPE-10489-PA. http://dx.doi.org/10.2118/10489-PA.

[r5-5] Kenyon, D. 1987. Third SPE Comparative Solution Project: Gas Cycling of Retrograde Condensate Reservoirs. J Pet Technol 39 (8): 981-997. SPE-12278-PA. http://dx.doi.org/10.2118/12278-PA.

[r6-6] 6.0 ^6.1 Aziz, K., Ramesh, A.B., and Woo, P.T. 1987. Fourth SPE Comparative Solution Project: Comparison of Steam Injection Simulators. J Pet Technol 39 (12): 1576–1584. SPE-13510-PA. http://dx.doi.org/10.2118/13510-PA.

[r7-7] Killough, J.E. and Kossack, C.A. 1987. Fifth Comparative Solution Project: Evaluation of Miscible Flood Simulators. Presented at the SPE Symposium on Reservoir Simulation, San Antonio, Texas, 1–4 February. SPE-16000-MS. http://dx.doi.org/10.2118/16000-MS.

[r8-8] Firoozabadi, A. and Thomas, L.K. 1990. Sixth SPE Comparative Solution Project: Dual-Porosity Simulators. J Pet Technol 42 (6): 710-715, 762-763. SPE-18741-PA. http://dx.doi.org/10.2118/18741-PA.

[r9-9] Nghiem, L., Collins, D.A., and Sharma, R. 1991. Seventh SPE Comparative Solution Project: Modelling of Horizontal Wells in Reservoir Simulation. Presented at the SPE Symposium on Reservoir Simulation, Anaheim, California, 17-20 February 1991. SPE-21221-MS. http://dx.doi.org/10.2118/21221-MS.

[r10-10] Quandalle, P. 1993. Eighth SPE Comparative Solution Project: Gridding Techniques in Reservoir Simulation. Presented at the SPE Symposium on Reservoir Simulation, New Orleans, Louisiana, 28 February-3 March 1993. SPE-25263-MS. http://dx.doi.org/10.2118/25263-MS.

[r11-11] 11.0 ^11.1 ^11.2 Killough, J.E. 1995. Ninth SPE Comparative Solution Project: A Reexamination of Black-Oil Simulation. Presented at the SPE Reservoir Simulation Symposium, San Antonio, Texas, USA, 12–15 February. SPE 29110. http://dx.doi.org/10.2118/29110-MS.

[r12-12] 12.0 ^12.1 ^12.2 ^12.3 Christie, M.A. and Blunt, M.J. 2001. Tenth SPE Comparative Solution Project: A Comparison of Upscaling Techniques. SPE Res Eval & Eng 4 (4): 308–317. SPE-72469-PA. http://dx.doi.org/10.2118/72469-PA.

[r13-13] 13.0 ^13.1 ^13.2 ^13.3 ^13.4 ^13.5 ^13.6 ^13.7 ^13.8 Coats, K.H. 1980. An Equation of State Compositional Model. J Pet Technol 20 (5): 363–376. SPE-8284-PA. http://dx.doi.org/10.2118/8284-PA.

[r14-14] 14.0 ^14.1 ^14.2 ^14.3 ^14.4 ^14.5 ^14.6 Young, L.C. and Stephenson, R.E. 1983. A Generalized Compositional Approach for Reservoir Simulation. SPE J. 23 (4): 727–742. SPE-10516-PA. http://dx.doi.org/10.2118/10516-PA.

[r15-15] 15.0 ^15.1 ^15.2 ^15.3 ^15.4 ^15.5 ^15.6 Ács, G., Doleschall, S., and Farkas, E. 1985. General Purpose Compositional Model. SPE J. 25 (4): 543–553. SPE-10515-PA. http://dx.doi.org/10.2118/10515-PA.

[r16-16] Watts, J.W. 1986. A Compositional Formulation of the Pressure and Saturation Equations. SPE Res Eng 1 (3): 243–252. SPE-12244-PA. http://dx.doi.org/10.2118/12244-PA.

[r17-17] Blair, P.M. and Weinaug, C.F. 1969. Solution of Two-Phase Flow Problems Using Implicit Difference Equations. SPE J. 9 (4): 417-424. SPE-2185-PA. http://dx.doi.org/10.2118/2185-PA.

[r18-18] Stone, H. and Garder, A.O.J. 1961. Analysis of Gas-Cap or Dissolved-Gas Drive Reservoirs. SPE J. 1 (2): 92–104. SPE-1518-G. http://dx.doi.org/10.2118/1518-G.

[r19-19] Sheldon, J.W., Harris, C.D., and Bavly, D. 1960. A Method for General Reservoir Behavior Simulation on Digital Computers. Presented at the Fall Meeting of the Society of Petroleum Engineers of AIME, Denver, Colorado, 2-5 October. http://dx.doi.org/10.2118/1521-G.

[r20-20] Fagin, R.G. and Stewart Jr., C.H. 1966. A New Approach to the Two-Dimensional Multiphase Reservoir Simulator. SPE J. 6 (2): 175-182. SPE-1188-PA. http://dx.doi.org/10.2118/1188-PA.

[r21-21] 21.0 ^21.1 Coats, K.H.: "Reservoir Simulation: A General Model Formulation and Associated Physical/Numerical Sources of Instability," Boundary and Interior Layers—Computational and Asymptotic Methods, J.J.H. Miller (ed.), Boole Press, Dublin, Ireland (1980).

[r22-22] 22.0 ^22.1 ^22.2 ^22.3 ^22.4 ^22.5 ^22.6 ^22.7 ^22.8 ^22.9 Coats, K.H., Thomas, L.K., and Pierson, R.G. 1998. Compositional and Black Oil Reservoir Simulation. SPE Res Eval & Eng 1 (4): 372-379. SPE-50990-PA. http://dx.doi.org/10.2118/50990-PA.

[r23-23] 23.0 ^23.1 Coats, K.H. 2000. A Note on IMPES and Some IMPES-Based Simulation Models. SPE J. 5 (3): 245-251. SPE-65092-PA. http://dx.doi.org/10.2118/65092-PA.

[r24-24] MacDonald, R.C. 1970. Methods for Numerical Simulation of Water and Gas Coning. SPE J. 10 (4): 425-436. SPE-2796-PA. http://dx.doi.org/10.2118/2796-PA.

[r25-25] Spillette, A.G., Hillestad, J.G., and Stone, H.L. 1973. A High-Stability Sequential Solution Approach to Reservoir Simulation. Presented at the Fall Meeting of the Society of Petroleum Engineers of AIME, Las Vegas, Nevada, 30 September-3 October 1973. SPE-4542-MS. http://dx.doi.org/10.2118/4542-MS.

[r26-26] 26.0 ^26.1 Thomas, G.W. and Thurnau, D.H. 1983. Reservoir Simulation Using an Adaptive Implicit Method. SPE J. 23 (5): 759–768. SPE-10120-PA. http://dx.doi.org/10.2118/10120-PA.

[r27-27] Lantz, R.B. 1971. Quantitative Evaluation of Numerical Diffusion (Truncation Error). SPE J. 11 (3): 315–320. SPE-2811-PA. http://dx.doi.org/10.2118/2811-PA.

[r28-28] Todd, M.R., O'Dell, P.M., and Hirasaki, G.J. 1972. Methods for Increased Accuracy in Numerical Reservoir Simulators. SPE J. 12 (6): 515–530. SPE-3516-PA. http://dx.doi.org/10.2118/3516-PA.

[r29-29] Fung, L.S.-K., Collins, D.A., and Nghiem, L.X. 1989. An Adaptive-Implicit Switching Criterion Based on Numerical Stability Analysis (includes associated paper 18727 ). SPE Res Eng 4 (1): 45-51. SPE-16003-PA. http://dx.doi.org/10.2118/16003-PA.

[r30-30] Russell, T.F. 1989. Stability Analysis and Switching Criteria for Adaptive Implicit Methods Based on the CFL Condition. Presented at the SPE Symposium on Reservoir Simulation, Houston, Texas, 6–8 February. SPE-18416-MS. http://dx.doi.org/10.2118/18416-MS.

[r31-31] Grabenstetter, J., Li, Y.-K., Collins, D.A. et al. 1991. Stability-Based Switching Criterion for Adaptive-Implicit Compositional Reservoir Simulation. Presented at the SPE Reservoir Simulation Symposium, Anaheim, California, 17–20 February. SPE-21225-MS. http://dx.doi.org/10.2118/21225-MS.

[r32-32] Watts, J.W. and Rame, M. 1999. An Algebraic Approach to the Adaptive Implicit Method. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 14-17 February 1999. SPE-51900-MS. http://dx.doi.org/10.2118/51900-MS.

[r33-33] Coats, K.H. 2001. IMPES Stability: The CFL Limit. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 11-14 February 2001. SPE-66345-MS. http://dx.doi.org/10.2118/66345-MS.

[r34-34] Coats, K.H. 2001. IMPES Stability: The Stable Step. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 11-14 February 2001. SPE-69225-MS. http://dx.doi.org/10.2118/69225-MS.

[r35-35] Young, L.C. and Russell, T.F. 1993. Implementation of an Adaptive Implicit Method. Presented at the SPE Symposium on Reservoir Simulation, New Orleans, Louisiana, 28 February-3 March 1993. SPE-25245-MS. http://dx.doi.org/10.2118/25245-MS.

[r36-36] 36.0 ^36.1 Vinesome, P.K.W. 1976. Orthomin, an Iterative Method for Solving Sparse Banded Sets of Simultaneous Linear Equations. Presented at the SPE Symposium on Numerical Simulation of Reservoir Performance, Los Angeles, 19–20 February. SPE-5729-MS. http://dx.doi.org/10.2118/5729-MS.

[r37-37] 37.0 ^37.1 ^37.2 Appleyard, J.R. 1983. Nested Factorization. Presented at the SPE Reservoir Simulation Symposium, San Francisco, California, USA, 15-18 November. SPE-12264-MS. http://dx.doi.org/10.2118/12264-MS.

[r38-38] Watts III, J.W. 1981. A Conjugate Gradient-Truncated Direct Method for the Iterative Solution of the Reservoir Simulation Pressure Equation. Society of Petroleum Engineers Journal 21 (3): 345-353. SPE-8252-PA. http://dx.doi.org/10.2118/8252-PA.

[r39-39] Forsyth, P.A. and Sammon, P.H. 1986. Local Mesh Refinement and Modeling of Faults and Pinchouts. SPE Form Eval 1 (3): 275-285. SPE-13524-PA. http://dx.doi.org/10.2118/13524-PA.

[r40-40] Verma, S.K. and Aziz, K. 1997. A Control Volume Scheme for Flexible Grids In Reservoir Simulation. Presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, 8–11 June. SPE-37999-MS. http://dx.doi.org/10.2118/37999-MS.

[r41-41] Saad, Y. and Schultz, M.H. 1986. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing 7 (3): 856-869. http://dx.doi.org/10.1137/0907058.

[r42-42] Tan, T.B.S. and Letkeman, J.P. 1982. Application of D4 Ordering And Minimization In An Effective Partial Matrix Inverse Iterative Method. Presented at the SPE Reservoir Simulation Symposium, New Orleans, Louisiana, 31 January-3 Feburary 1982. SPE-10493-MS. http://dx.doi.org/10.2118/10493-MS.

[r43-43] Watts, J.W. 1971. An Iterative Matrix Solution Method Suitable for Anisotropic Problems. Society of Petroleum Engineers Journal 11 (1): 47-51. SPE-2802-PA. http://dx.doi.org/10.2118/2802-PA.

[r44-44] 44.0 ^44.1 Wallis, J.R. 1983. Incomplete Gaussian Elimination as a Preconditioning for Generalized Conjugate Gradient Acceleration. Presented at the SPE Reservoir Simulation Symposium, San Francisco, California, USA, 15-18 November. SPE-12265-MS. http://dx.doi.org/10.2118/12265-MS.

[r45-45] Behie, A. and Vinsome, P.K.W. 1982. Block Iterative Methods for Fully Implicit Reservoir Simulation. Society of Petroleum Engineers Journal 22 (5): 658-668. SPE-9303-PA. http://dx.doi.org/10.2118/9303-PA.

[r46-46] Holmes, J.A. 1983. Enhancements to the Strongly Coupled, Fully Implicit Well Model: Wellbore Crossflow Modeling and Collective Well Control. Presented at the SPE Reservoir Simulation Symposium, San Francisco, California, 15-18 November 1983. SPE-12259-MS. http://dx.doi.org/10.2118/12259-MS.

[r47-47] 47.0 ^47.1 Burrows, R., Ponting, D., and Wood, L. 1996. Parallel Reservoir Simulation with Nested Factorization. Presented at the 1996 European Conference on the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.

[r48-48] Killough, J.E. et al. 1996. Parallelization of a General-Purpose Reservoir Simulator. Presented at the 1996 European Conference on the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.

[r49-49] Cuthill, E. and McKee, J. 1969. Reducing the Bandwidth of Sparse Symmetric Matrices. Proc., 24th National Conference of the Assoc. for Computing Machinery, Brandon Press, Princeton, New Jersey, 157–172.

[r50-50] D'Azevedo, E., Forsyth, P.A., and Tang, W.-P. 1991. An Automatic Ordering Method for Incomplete Factorization Iterative Solvers. Presented at the SPE Symposium on Reservoir Simulation, Anaheim, California, 17-20 February 1991. SPE-21226-MS. http://dx.doi.org/10.2118/21226-MS.

[r51-51] D’Azevedo, E.F., Forsyth, P.A., and Tang, W.-P. 1992. Ordering Methods for Preconditioned Conjugate Gradient Methods Applied to Unstructured Grid Problems. SIAM Journal of Matrix Analysis and Applications 13 (3): 944-961. http://dx.doi.org/10.1137/0613057.

[r52-52] 52.0 ^52.1 Aziz, K. and Settari, A. 1979. Petroleum Reservoir Simulation, ISBN 0–9730614-0-5, reprinted 2002. Calgary: Blitzprint Ltd.

[r53-53] 53.0 ^53.1 Rozon, B.J. 1989. A Generalized Finite Volume Discretization Method for Reservoir Simulation. Presented at the SPE Symposium on Reservoir Simulation, Houston, Texas, 6-8 February 1989. SPE-18414-MS. http://dx.doi.org/10.2118/18414-MS.

[r54-54] 54.0 ^54.1 Aavatsmark, I. et al. 2000. MPFA for Faults and Local Refinements With Application to Field Simulations. Proc., 2000 European Conference on the Mathematics of Oil Recovery, Baveno, Italy, 5–8 September.

[r55-55] Forsyth, P.A. 1990. A Control-Volume, Finite-Element Method for Local Mesh Refinement in Thermal Reservoir Simulation. SPE Res Eng 5 (4): 561-566. SPE-18415-PA. http://dx.doi.org/10.2118/18415-PA.

[r56-56] Hegre, T.M., Dalen, V., and Henriquez, A. 1986. Generalized Transmissibilities for Distorted Grids in Reservoir Simulation. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 5-8 October 1986. SPE-15622-MS. http://dx.doi.org/10.2118/15622-MS.

[r57-57] 57.0 ^57.1 Young, L.C. 1999. Rigorous Treatment of Distorted Grids in 3D. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 14-17 February 1999. SPE-51899-MS. http://dx.doi.org/10.2118/51899-MS.

[r58-58] 58.0 ^58.1 ^58.2 Aziz, K. 1993. Reservoir Simulation Grids: Opportunities and Problems. J Pet Technol 45 (7): 658-663. SPE-25233-PA. http://dx.doi.org/10.2118/25233-PA.

[r59-59] Aavatsmark, I., Barkve, T., and Mannseth, T. 1998. Control-Volume Discretization Methods for 3D Quadrilateral Grids in Inhomogeneous, Anisotropic Reservoirs. SPE J. 3 (2): 146-154. SPE-38000-PA. http://dx.doi.org/10.2118/38000-PA.

[r60-60] Edwards, M.G. and Rogers, C.F. 1994. A Flux Continuous Scheme for the Full Tensor Pressure Equation. Proc., 1994 European Conference on the Mathematics of Oil Recovery, Roros, Norway, 7–10 June.

[r61-61] Lee, S.H., Tchelepi, H., and DeChant, L.F. 1999. Implementation of a Flux-Continuous Finite Difference Method for Stratigraphic, Hexahedron Grids. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 14-17 February 1999. SPE-51901-MS. http://dx.doi.org/10.2118/51901-MS.

[r62-62] Ponting, D.K. 1989. Corner Point Geometry in Reservoir Simulation. Proc., First European Conference on the Mathematics of Oil Recovery.

[r63-63] 63.0 ^63.1 Jenny, P., Wolfsteiner, C., Lee, S.H. et al. 2002. Modeling Flow in Geometrically Complex Reservoirs Using Hexahedral Multiblock Grids. SPE J. 7 (2): 149-157. SPE-78673-PA. http://dx.doi.org/10.2118/78673-PA.

[r64-64] Heinemann, Z.E. and Brand, C.W. 1988. Gridding Techniques in Reservoir Simulation. Proc., First Intl. Forum on Reservoir Simulation, Alpbach, Austria, 339.

[r65-65] Heinemann, Z.E., Brand, C.W., Munka, M. et al. 1991. Modeling Reservoir Geometry With Irregular Grids. SPE Res Eng 6 (2): 225–232. SPE-18412-PA. http://dx.doi.org/10.2118/18412-PA.

[r66-66] Palagi, C.L. and Aziz, K. 1994. Use of Voronoi Grid in Reservoir Simulation. SPE Advanced Technology Series 2 (2): 69–77. SPE-22889-PA. http://dx.doi.org/10.2118/22889-PA.

[r67-67] Gunasekera, D., Cox, J., and Lindsey, P. 1997. The Generation and Application of K-Orthogonal Grid Systems. Presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, 8–11 June. SPE-37998-MS. http://dx.doi.org/10.2118/37998-MS.

[r68-68] Verma, S. and Aziz, K. 1996. Two- and Three-Dimensional Flexible Grids for Reservoir Simulation. Proc., 1996 European Conference on the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.

[r69-69] Lasseter, T.J. 2002. A New Approach for the Efficient Construction of 3D Geological Models for Reservoir Applications. Proc., Eighth European Conference on the Mathematics of Oil Recovery, Freiberg, Germany.

[r70-70] Warren, J.E. and Price, H.S. 1961. Flow in Heterogeneous Porous Media. SPE J. 1 (3): 153–169. SPE-1579-G. http://dx.doi.org/10.2118/1579-G.

[r71-71] King, M.J., MacDonald, D.G., Todd, S.P. et al. 1998. Application of Novel Upscaling Approaches to the Magnus and Andrew Reservoirs. Presented at the European Petroleum Conference, The Hague, 20–22 October. SPE-50643-MS. http://dx.doi.org/10.2118/50643-MS.

[r72-72] King, M.J. and Mansfield, M. 1999. Flow Simulation of Geologic Models. SPE Res Eval & Eng 2 (4): 351-367. SPE-57469-PA. http://dx.doi.org/10.2118/57469-PA.

[r73-73] Guerillot, D. et al. 1989. An Integrated Model for Computer Aided Reservoir Description: From Outcrop Study to Fluid Flow Simulations. Paper presented at the 1989 IOR Symposium, Budapest, Hungary.

[r74-74] Samier, P. 1990. A Finite-Element Method for Calculating Transmissibilities. Proc., 1990 European Conference on the Mathematics of Oil Recovery, Arles, France, 11–14 September.

[r75-75] 75.0 ^75.1 Durlofsky, L.J. and Chung, E.Y. 1990. Effective Permeability of Heterogeneous Reservoir Regions. Proc., 1990 European Conference on the Mathematics of Oil Recovery, Arles, France, 11–14 September.

[r76-76] 76.0 ^76.1 Durlofsky, L.J. 1991. Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media. Water Resour. Res. 27 (5): 699-708. http://dx.doi.org/10.1029/91wr00107.

[r77-77] Kirkpatrick, S. 1971. Classical Transport in Disordered Media: Scaling and Effective-Medium Theories. Phys. Rev. Lett. 27 (25): 1722-1725. http://dx.doi.org/10.1103/PhysRevLett.27.1722.

[r78-78] Tjølsen, C., Bøe, Ø., and Damsleth, E. 1992. Homogenisation of Absolute Permeability—Does it Affect the Final Results? Proc., Second Lerkendal Petroleum Engineering Workshop.

[r79-79] Bøe, Ø. 1994. Analysis of an upscaling method based on conservation of dissipation. Transport Porous Media 17 (1): 77-86. http://dx.doi.org/10.1007/bf00624051.

[r80-80] 80.0 ^80.1 ^80.2 Holden, L. and Nielsen, B.F. 2000. Global Upscaling of Permeability in Heterogeneous Reservoirs; The Output Least Squares (OLS) Method. Transport Porous Media 40 (2): 115-143. http://dx.doi.org/10.1023/a:1006657515753.

[r81-81] 81.0 ^81.1 Samier, P., Quettier, L., and Thiele, M. 2002. Applications of Streamline Simulations to Reservoir Studies. SPE Res Eval & Eng 5 (4): 324–332. SPE-78883-PA. http://dx.doi.org/10.2118/78883-PA.

[r82-82] Thiele, M.R. 2001. Streamline Simulation. Keynote address at the 2001 Intl. Forum on Reservoir Simulation, Schloss Fuschl, Austria, 3–7 September.

[r83-83] Lolomari, T., Bratvedt, K., Crane, M. et al. 2000. The Use of Streamline Simulation in Reservoir Management: Methodology and Case Studies. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 1–4 October. SPE-63157-MS. http://dx.doi.org/10.2118/63157-MS.

[r84-84] 84.0 ^84.1 Baker, R.O., Kuppe, F., Chug, S. et al. 2001. Full-Field Modeling Using Streamline-Based Simulation: 4 Case Studies. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 11–14 February. SPE-66405-MS. http://dx.doi.org/10.2118/66405-MS.

[r85-85] Thiele, M.R., Batycky, R.P., and Blunt, M.J. 1997. A Streamline-Based 3D Field-Scale Compositional Reservoir Simulator. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 5–8 October. SPE-38889-MS. http://dx.doi.org/10.2118/38889-MS.

[r86-86] Thiele, M.R., Batycky, R.P., and Thomas, L.K. 2002. Miscible WAG Simulations Using Streamlines. Paper presented at the 2002 European Conference on the Mathematics of Oil Recovery, Freiberg, Germany, 3–6 September.

[r87-87] 87.0 ^87.1 ^87.2 ^87.3 Batycky, R.P., Blunt, M.J., and Thiele, M.R. 1997. A 3D Field-Scale Streamline-Based Reservoir Simulator. SPE Res Eng 12 (4): 246-254. SPE-36726-PA. http://dx.doi.org/10.2118/36726-PA.

[r88-88] Emanuel, A.S. 1993. Development of an Analytical Streamtube Model for Estimating Steam Drive Performance. SPE Advanced Technology Series 1 (1): 81-89. SPE-21756-PA. http://dx.doi.org/10.2118/21756-PA.

[r89-89] Crane, M.J. and Blunt, M.J. 1999. Streamline-based simulation of solute transport. Water Resour. Res. 35 (10): 3061-3078. http://dx.doi.org/10.1029/1999wr900145.

[r90-90] Di Donato, G., Huang, W., and Blunt, M. 2003. Streamline-Based Dual Porosity Simulation of Fractured Reservoirs. Presented at the SPE Annual Technical Conference and Exhibition, 5–8 October. SPE-84036-MS. http://dx.doi.org/10.2118/84036-MS.

[r91-91] Muskat, M. and Wyckoff, R.D. 1934. A Theoretical Analysis of Waterflooding Networks. Trans., AIME, 107: 62–77.

[r92-92] Muskat, M. 1982. Flow of Homogeneous Fluids. Boston, Massachusetts: Intl. Human Resources Development Corp.

[r93-93] Fay, C.H. and Prats, M. 1951. The Application of Numerical Methods to Cycling and Flooding Problems. Proc., Third World Petroleum Congress, The Hague, 28 May–6 June.

[r94-94] Higgins, R.V. and Leighton, A.J. 1962. A Computer Method to Calculate Two-Phase Flow in Any Irregularly Bounded Porous Medium. J Pet Technol 14 (6): 679-683. SPE-243-PA. http://dx.doi.org/10.2118/243-PA.

[r95-95] Higgins, R.V. and Leighton, A.J. 1962. Computer Prediction of Water Drive of Oil and Gas Mixtures Through Irregularly Bounded Porous Media Three-Phase Flow. J Pet Technol 14 (9): 1048-1054. SPE-283-PA. http://dx.doi.org/10.2118/283-PA.

[r96-96] Bommer, P.M. and Schechter, R.S. 1979. Mathematical Modeling of In-Situ Uranium Leaching. Society of Petroleum Engineers Journal 19 (6): 393-400. SPE-7533-PA. http://dx.doi.org/10.2118/7533-PA.

[r97-97] Martin, J.C. and Wegner, R.E. 1979. Numerical Solution of Multiphase, Two-Dimensional Incompressible Flow Using Stream-Tube Relationships. Society of Petroleum Engineers Journal 19 (5): 313-323. SPE-7140-PA. http://dx.doi.org/10.2118/7140-PA.

[r98-98] Lake, L.W., Johnston, J.R., and Stegemeier, G.L. 1981. Simulation and Performance Prediction of a Large-Scale Surfactant/Polymer Project. Society of Petroleum Engineers Journal 21 (6): 731-739. SPE-7471-PA. http://dx.doi.org/10.2118/7471-PA.

[r99-99] Emanuel, A.S., Alameda, G.K., Behrens, R.A. et al. 1989. Reservoir Performance Prediction Methods Based on Fractal Geostatistics(includes associated papers 20011 and 20158 ). SPE Res Eng 4 (3): 311-318. SPE-16971-PA. http://dx.doi.org/10.2118/16971-PA.

[r100-100] 100.0 ^100.1 Pollock, D.W. 1988. Semianalytical Computation of Path Lines for Finite-Difference Models. Ground Water 26 (6): 743.

[r101-101] King, M.J., Blunt, M.J., Mansfield, M.M. et al. 1993. Rapid Evaluation of the Impact of Heterogeneity on Miscible Gas Injection. Presented at the SPE Western Regional Meeting, Anchorage, Alaska, 26-28 May 1993. SPE-26079-MS. http://dx.doi.org/10.2118/26079-MS.

[r102-102] Datta-Gupta, A. and King, M.J. 1995. A semianalytic approach to tracer flow modeling in heterogeneous permeable media. Adv. Water Resour. 18 (1): 9-24. http://dx.doi.org/http://dx.doi.org/10.1016/0309-1708(94)00021-V.

[r103-103] 103.0 ^103.1 Batycky, R.P. 1997. A Three-Dimensional Two-Phase Field Scale Streamline Simulator. PhD dissertation, Stanford U., Stanford, California.

[r104-104] Renard, G. 1990. A 2D reservoir streamtube EOR model with periodical automatic regeneration of streamlines. In Situ 14 (2): 175-200.

[r105-105] Thiele, M.R., Batycky, R.P., Blunt, M.J. et al. 1996. Simulating Flow in Heterogeneous Systems Using Streamtubes and Streamlines. SPE Res Eng 11 (1): 5–12. SPE-27834-PA. http://dx.doi.org/10.2118/27834-PA.

[r106-106] 106.0 ^106.1 Bratvedt, F., Gimse, T., and Tegnander, C. 1996. Streamline computations for porous media flow including gravity. Transport Porous Media 25 (1): 63-78. http://dx.doi.org/10.1007/bf00141262.

[r107-107] Blunt, M.J., Liu, K., and Thiele, M.R. 1996. A Generalized Streamline Method to Predict Reservoir Flow. Petroleum Geoscience 2 (3): 259-269. http://dx.doi.org/10.1144/petgeo.2.3.259.

[r108-108] Prévost, M., Edwards, M.G., and Blunt, M.J. 2002. Streamline Tracing on Curvilinear Structured and Unstructured Grids. SPE J. 7 (2): 139–148. SPE-78663-PA. http://dx.doi.org/10.2118/786637-PA.

[r109-109] Thiele, M.R. and Batycky, R.P. 2003. Water Injection Optimization Using a Streamline-Based Workflow. Presented at the SPE Annual Technical Conference and Exhibition, Denver, 5-8 October. SPE 84080. http://dx.doi.org/10.2118/84080-MS.

[r110-110] Batycky, R.P., Thiele, M.R., Baker, R.O. et al. 2005. Revisiting Reservoir Flood-Surveillance Methods Using Streamlines. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 9-12 October 2005. SPE-95402-MS. http://dx.doi.org/10.2118/95402-MS.

[r111-111] Wang, Y. and Kovscek, A.R. 2000. Streamline Approach for History Matching Production Data. SPE J. 5 (4): 353-362. SPE-58350-PA. http://dx.doi.org/10.2118/58350-PA.

[r112-112] Agarwal, B. and Blunt, M.J. 2003. Streamline-Based Method With Full-Physics Forward Simulation for History-Matching Performance Data of a North Sea Field. SPE J. 8 (2): 171-180. SPE-84952-PA. http://dx.doi.org/10.2118/84952-PA.

[r113-113] Caers, J., Krishnan, S., Wang, Y. et al. 2002. A Geostatistical Approach to Streamline-Based History Matching. SPE J. 7 (3): 250-266. SPE-73144-PA. http://dx.doi.org/10.2118/73144-PA.

[r114-114] Fenwick, D., Thiele, M., Agil, M. et al. 2005. Reconciling Prior Geologic Information With Production Data Using Streamlines: Application to a Giant Middle-Eastern Oil Field. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 9–12 October. SPE-95940-MS. http://dx.doi.org/10.2118/95940-MS.

[r115-115] Stüben, K. 2000. Algebraic Multigrid (AMG): An Introduction with Applications. Multigrid, U. Trottenberg, C.W. Oosterlee, and A. Schüller eds., appendix. New York City: Academic Press.

[r116-116] Darve, F. 1990. Geomaterials: Constitutive Equations and Modeling. New York City: Elsevier Applied Science.

[r117-117] Coussy, O. 1995. Mechanics of Porous Continua. Chichester, England: John Wiley and Sons.

[r118-118] Settari, A. and Walters, D.A. 2001. Advances in Coupled Geomechanical and Reservoir Modeling With Applications to Reservoir Compaction. SPE J. 6 (3): 334–342. SPE-51927-PA. http://dx.doi.org/10.2118/51927-MS.

[r119-119] 119.0 ^119.1 Lorenz, J.C. 1999. Stress Sensitive Reservoirs. J Pet Technol (January): 61.

[r120-120] Settari, A., Warren, G.M., Jacquemont, J. et al. 1999. Brine Disposal Into a Tight Stress-Sensitive Formation at Fracturing Conditions: Design and Field Experience. SPE Res Eval & Eng 2 (2): 186-195. SPE-56001-PA. http://dx.doi.org/10.2118/56001-PA.

[r121-121] Settari, A. 2002. Reservoir Compaction. Distinguished Author Series, J Pet Technol 54 (8): 62–69. SPE-76805-PA. http://dx.doi.org/10.2118/76805-MS.

[r122-122] Finol, A. and Farouq Ali, S.M. 1975. Numerical Simulation of Oil Production With Simultaneous Ground Subsidence. Society of Petroleum Engineers Journal 15 (5): 411-424. SPE-4847-PA. http://dx.doi.org/10.2118/4847-PA.

[r123-123] Merle, H.A., Kentie, C.J.P., van Opstal, G.H.C. et al. 1976. The Bachaquero Study - A Composite Analysis of the Behavior of a Compaction Drive/Solution Gas Drive Reservoir. J Pet Technol 28 (9): 1107-1115. SPE-5529-PA. http://dx.doi.org/10.2118/5529-PA.

[r124-124] Rattia, A.J. and Farouq Ali, S.M. 1981. Effect of Formation Compaction on Steam Injection Response. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 4-7 October 1981. SPE-10323-MS. http://dx.doi.org/10.2118/10323-MS.

[r125-125] Chase, C.A.J. and Dietrich, J.K. 1989. Compaction Within the South Belridge Diatomite. SPE Res Eng 4 (4): 422–428. SPE-17415-PA. http://dx.doi.org/10.2118/17415-PA.

[r126-126] Beattie, C.I., Boberg, T.C., and McNab, G.S. 1991. Reservoir Simulation of Cyclic Steam Stimulation in the Cold Lake Oil Sands. SPE Res Eng 6 (2): 200–206. SPE-18752-PA. http://dx.doi.org/10.2118/18752-PA.

[r127-127] 127.0 ^127.1 Chin, L.Y. and Thomas, L.K. 1999. Fully Coupled Analysis of Improved Oil Recovery by Reservoir Compaction. Presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, 3–6. October. SPE-56753-MS. http://dx.doi.org/10.2118/56753-MS.

[r128-128] 128.0 ^128.1 Settari, A. and Mourits, F.M. 1998. A Coupled Reservoir and Geomechanical Simulation System. SPE J. 3 (3): 219–226. SPE-50939-PA. http://dx.doi.org/10.2118/50939-PA

[r129-129] Settari, A. and Price, H.S. 1984. Simulation of Hydraulic Fracturing in Low-Permeability Reservoirs. Society of Petroleum Engineers Journal 24 (2): 141-152. SPE-8939-PA. http://dx.doi.org/10.2118/8939-PA.

[r130-130] Settari, A.T., Bachman, R.C., and D.A.Walters. 2005. How To Approximate Effects of Geomechanics in Conventional Reservoir Simulation. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, 9-12 October. SPE 97155. http://dx.doi.org/10.2118/97155-MS.

[r131-131] Heffer, K.J. et al. 1992. The Influence of Natural Fractures, Faults and Earth Stresses on Reservoir Performance—Geomechanical Analysis by Numerical Modelling. Proc., Third Intl. Conference on North Sea Oil and Gas Reservoirs, Trondheim, Norway, 201–211.

[r132-132] Jin, M., Somerville, J., and Smart, B.G.D. 2000. Coupled Reservoir Simulation Applied to the Management of Production Induced Stress-Sensitivity. Presented at the International Oil and Gas Conference and Exhibition in China, Beijing, China, 7–10 November. SPE-64790-MS. http://dx.doi.org/10.2118/64790-MS.

[r133-133] Coupled Thermo-Hydro-Mechanical Processes of Fractured Media. 1996. Developments in Geotechnical Engineering, O. Stephansson, L. Jing, and S.-F. Tsang eds., 575. New York City: Elsevier.

[r134-134] Koutsabeloulis, N.C. and Hope, S.A. 1998. “Coupled” Stress/Fluid/Thermal Multi-Phase Reservoir Simulation Studies Incorporating Rock Mechanics. Presented at the SPE/ISRM Rock Mechanics in Petroleum Engineering, Trondheim, Norway, 8–10 July. SPE-47393-MS. http://dx.doi.org/10.2118/47393-MS.

[r135-135] Settari, A. 2000. Coupled Fracture and Reservoir Modeling. Proc., Advanced Hydraulic Fracturing Workshop, 2000 North American Rock Mechanics Conference, Seattle, Washington.

[r136-136] Lewis, R.W. and Schrefler, B.A. 1987. The Finite Element Method in the Deformation and Consolidation of Porous Media, 344. New York City: John Wiley and Sons.

[r137-137] Stone, T., Bowen, G., Papanastasiou, P. et al. 2000. Fully Coupled Geomechanics in a Commercial Reservoir Simulator. Presented at the SPE European Petroleum Conference, Paris, France, 24–25 October. SPE-65107-MS. http://dx.doi.org/10.2118/65107-MS.

[r138-138] Sulak, R.M., Thomas, L.K., and Boade, R.R. 1991. 3D Reservoir Simulation of Ekofisk Compaction Drive (includes associated papers 24317 and 24400 ). J Pet Technol 43 (10): 1272-1278. SPE-19802-PA. http://dx.doi.org/10.2118/19802-PA.

[r139-139] Fredrich, J.T. et al. 2001. One-Way Coupled Reservoir-Geomechanical Modeling of the Lost Hills Oil Field, California. Proc., 38th U.S. Rock Mechanics Symposium, Washington, DC.

[r140-140] Tran, D., Settari, A., and Nghiem, L. 2002. New Iterative Coupling Between a Reservoir Simulator and a Geomechanics Module. Presented at the SPE/ISRM Rock Mechanics Conference, Irving, Texas, 20-23 October 2002. SPE-78192-MS. http://dx.doi.org/10.2118/78192-MS.

[r141-141] Dean, R.H. 2000. A Poroelastic Multicomponent Reservoir Simulator. Paper SPE 62806 available from SPE, Richardson, Texas.

[r142-142] Settari, A., Walters, D.A., and Behie, G.A. 2000. Use of Coupled Reservoir and Geomechanical Modelling for Integrated Reservoir Analysis and Management. Presented at the Canadian International Petroleum Conference, Calgary, Alberta, Jun 4 - 8, 2000 2000. PETSOC-2000-078. http://dx.doi.org/10.2118/2000-078.

[r143-143] 143.0 ^143.1 ^143.2 Coats, B.K., Nolen, J.S., and Sherman, A.H. 1988. Applicability and Accuracy of Two Component Characterizations of Reservoir Oils for Numerical Modeling Purposes. Presented at the International Meeting on Petroleum Engineering, Tianjin, China, 1-4 November 1988. SPE-17614-MS. http://dx.doi.org/10.2118/17614-MS.

[r144-144] 144.0 ^144.1 ^144.2 Cook, R.E., Jacoby, R.H., and Ramesh, A.B. 1974. A Beta-Type Reservoir Simulator for Approximating Compositional Effects During Gas Injection. Society of Petroleum Engineers Journal 14 (5): 471-481. SPE-4272-PA. http://dx.doi.org/10.2118/4272-PA.

[r145-145] 145.0 ^145.1 ^145.2 ^145.3 Fevang, Ø. and Whitson, C.H. 1996. Modeling Gas-Condensate Well Deliverability. SPE Res Eng 11 (4): 221–230. SPE-30714-PA. http://dx.doi.org/10.2118/30714-PA.

[r146-146] Spivak, A. and Dixon, T.N. 1973. Simulation of Gas-Condensate Reservoirs. Presented at the SPE Symposium on Numerical Simulation of Reservoir Performance, Houston, Texas, 11-12 January 1973. SPE-4271-MS. http://dx.doi.org/10.2118/4271-MS.

[r147-147] Chang, Y.-B., Coats, B.K., and Nolen, J.S. 1998. A Compositional Model for CO2 Floods Including CO2 Solubility in Water. SPE Res Eval & Eng 1 (2): 155-160. SPE 35164. http://dx.doi.org/10.2118/35164-PA.

[r148-148] 148.0 ^148.1 Michelsen, M.L. 1985. Saturation point calculations. Fluid Phase Equilib. 23 (2–3): 181-192. http://dx.doi.org/http://dx.doi.org/10.1016/0378-3812(85)90005-6.

[r149-149] Michelsen, M.L. 1982. The isothermal flash problem. Part I. Stability. Fluid Phase Equilib. 9 (1): 1-19. http://dx.doi.org/http://dx.doi.org/10.1016/0378-3812(82)85001-2.

[r150-150] Muskat, M. and McDowell, J.M. 1949. An Electrical Computer for Solving Phase Equilibrium Problems. J Pet Technol 1 (11): 291-298. http://dx.doi.org/10.2118/949291-G.

[r151-151] Rachford Jr., H.H. and Rice, J.D. 1952. Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium. J Pet Technol 4 (10): 19, 3. http://dx.doi.org/10.2118/952327-G.

[r152-152] 152.0 ^152.1 ^152.2 ^152.3 ^152.4 ^152.5 ^152.6 Whitson, C.H. and Brule, M.R. 2000. Phase Behavior, Vol. 20. Richardson, Texas: Monograph Series, SPE.

[r153-153] 153.0 ^153.1 Michelsen, M.L. 1982. The isothermal flash problem. Part II. Phase-split calculation. Fluid Phase Equilib. 9 (1): 21-40. http://dx.doi.org/http://dx.doi.org/10.1016/0378-3812(82)85002-4.

[r154-154] 154.0 ^154.1 Péneloux, A., Rauzy, E., and Fréze, R. 1982. A consistent correction for Redlich-Kwong-Soave volumes. Fluid Phase Equilib. 8 (1): 7-23. http://dx.doi.org/http://dx.doi.org/10.1016/0378-3812(82)80002-2.

[r155-155] 155.0 ^155.1 Peng, D.Y. and Robinson, D.B. 1976. A New-Constant Equation of State. Ind. Eng. Chem. Fund. 15 (1): 59-64. http://dx.doi.org/10.1021/i160057a011.

[r156-156] 156.0 ^156.1 Soave, G. 1972. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci. 27 (6): 1197-1203. http://dx.doi.org/http://dx.doi.org/10.1016/0009-2509(72)80096-4.

[r157-157] 157.0 ^157.1 Zudkevitch, D. and Joffe, J. 1970. Correlation and prediction of vapor-liquid equilibria with the redlich-kwong equation of state. AIChE J. 16 (1): 112-119. http://dx.doi.org/10.1002/aic.690160122.

[r158-158] Whitson, C.H. and Torp, S.B. 1983. Evaluating Constant-Volume Depletion Data. J Pet Technol 35 (3): 610-620. SPE-10067-PA. http://dx.doi.org/10.2118/10067-PA.

[r159-159] 159.0 ^159.1 Cook, A.B., Spencer, G.B., and Bayazeed, A.F. 1967. The Role of Vaporization in High Percentage Oil Recovery by Pressure Maintenance. J Pet Technol 19 (2): 245-250. SPE-1646-PA. http://dx.doi.org/10.2118/1646-PA.

[r160-160] Whitson, C.H. 1983. Characterizing Hydrocarbon Plus Fractions. SPE J. 23 (4): 683–693. SPE-12233-PA. http://dx.doi.org/10.2118/12233-PA.

[r161-161] Tang, D.E. and Zick, A.A. 1993. A New Limited Compositional Reservoir Simulator. Presented at the SPE Symposium on Reservoir Simulation, New Orleans, Louisiana, 28 february–3 March. SPE-25255-MS. http://dx.doi.org/10.2118/25255-MS.

[r162-162] 162.0 ^162.1 Zick, A.A. 1986. A Combined Condensing/Vaporizing Mechanism in the Displacement of Oil by Enriched Gases. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, 5–8 October. SPE 15493. http://dx.doi.org/10.2118/15493-MS.

[r163-163] Fevang, Ø., Singh, K., and Whitson, C.H. 2000. Guidelines for Choosing Compositional and Black-Oil Models for Volatile Oil and Gas-Condensate Reservoirs. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 1-4 October 2000. SPE-63087-MS. http://dx.doi.org/10.2118/63087-MS.

[r164-164] Lohrenz, J., Bray, B.G., and Clark, C.R. 1964. Calculating Viscosities of Reservoir Fluids From Their Compositions. J Pet Technol 16 (10): 1171–1176. SPE-915-PA. http://dx.doi.org/10.2118/915-PA.

[r165-165] Pedersen, K.S. and Fredenslund, A. 1987. An Improved Corresponding States Model for the Prediction of Oil and Gas Viscosities and Thermal Conductivities. Chem. Eng. Sci. 42: 182.

[r166-166] Saidi, A.M. 1996. Twenty Years of Gas Injection History into Well-Fractured Haft Kel Field (Iran). Presented at the International Petroleum Conference and Exhibition of Mexico, Villahermosa, Mexico, 5–7 March. SPE 35309. http://dx.doi.org/10.2118/35309-MS.

[r167-167] Whitson, C.H., Fevang, Ø., and Sævareid, A. 1999. Gas Condensate Relative Permeability for Well Calculations. Presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, 3–6 October. SPE-56476-MS. http://dx.doi.org/10.2118/56476-MS.

[r168-168] Mott, R.E., Cable, A.S., and Spearing, M.C. 2000. Measurements of Relative Permeabilities for Calculating Gas-Condensate Well Deliverability. SPE Res Eval & Eng 3 (6): 473–479. SPE-68050-PA. http://dx.doi.org/10.2118/68050-PA.

[r169-169] Robinson, D.B. and Peng, D.Y. 1978. The Characterization of the Heptanes and Heavier Fractions. Research Report 28. Tulsa, Oklahoma: Gas Producers Assn.

[r170-170] Katz, D.L. and Firoozabadi, A. 1978. Predicting Phase Behavior of Condensate/Crude-Oil Systems Using Methane Interaction Coefficients. J Pet Technol 30 (11): 1649–1655. SPE-6721-PA. http://dx.doi.org/10.2118/6721-PA.

[r171-171] Shelton, J.L. and Yarborough, L. 1977. Multiple Phase Behavior in Porous Media During CO2 or Rich-Gas Flooding. J Pet Technol 29 (9): 1171-1178. SPE-5827-PA. http://dx.doi.org/10.2118/5827-PA.

[r172-172] Nghiem, L.X. and Li, Y.K. 1986. Effect of Phase Behavior on CO2 Displacement Efficiency at Low Temperatures: Model Studies With an Equation of State. SPE Res Eng 1 (4): 414–422. SPE-13116-PA. http://dx.doi.org/10.2118/13116-PA.

[r173-173] Coats, K.H. 1976. Simulation of Steamflooding With Distillation and Solution Gas. SPE J. 16 (5): 235–247. SPE-5015-PA. http://dx.doi.org/10.2118/5015-PA.

[r174-174] Killough, J.E. 1979. The Use of Vector Processors in Reservoir Simulation. Presented at the SPE Reservoir Simulation Symposium, Denver, Colorado, 31 January-2 Feburary 1979. SPE-7673-MS. http://dx.doi.org/10.2118/7673-MS.

[r175-175] Nolen, J.S., Kuba, D.W., and Jr., M.J.K. 1981. Application of Vector Processors To Solve Finite Difference Equations. Society of Petroleum Engineers Journal 21 (4): 447-453. SPE-7675-PA. http://dx.doi.org/10.2118/7675-PA.

[r176-176] Killough, J.E. and Levesque, J.M. 1982. Reservoir Simulation and the In-house Vector Processor: Experience for the First Year. Presented at the SPE Reservoir Simulation Symposium, New Orleans, Louisiana, 31 January-3 Feburary 1982. SPE-10521-MS. http://dx.doi.org/10.2118/10521-MS.

[r177-177] Young, L.C. 1987. Equation of State Compositional Modeling on Vector Processors. Presented at the SPE Symposium on Reservoir Simulation, San Antonio, Texas, 1-4 February 1987. SPE-16023-MS. http://dx.doi.org/10.2118/16023-MS.

[r178-178] Killough, J.E. 1993. Is Parallel Computing Ready for Reservoir Simulation? A Critical Analysis of the State of the Art. Presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, 3-6 October 1993. SPE-26634-MS. http://dx.doi.org/10.2118/26634-MS.

[r179-179] 179.0 ^179.1 ^179.2 ^179.3 ^179.4 Killough, J.E. et al. 1996. A General-Purpose parallel Reservoir Simulator. Presented at the 1996 European Conference on the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.

[r180-180] van Daalen, D.T. et al. 1990. The parallelization of BOSIM, Shell’s Black/Volatile Oil Reservoir Simulator. Proc. of the First IMA/SPE European Conference on the Mathematics of Oil Recovery, Oxford U.

[r181-181] Killough, J.E. and Bhogeswara, R. 1991. Simulation of Compositional Reservoir Phenomena on a Distributed-Memory Parallel Computer. J Pet Technol 43 (11): 1368-1374. SPE-21208-PA. http://dx.doi.org/10.2118/21208-PA.

[r182-182] Wheeler, J.A. and Smith, R.A. 1990. Reservoir Simulation on a Hypercube. SPE Res Eng 5 (4): 544-548. SPE-19804-PA. http://dx.doi.org/10.2118/19804-PA.

[r183-183] 183.0 ^183.1 Rutledge, J.M., Jones, D.R., Chen, W.H. et al. 1992. Use of a Massively Parallel SIMD Computer for Reservoir Simulation. SPE Comp App 4 (4): 16-21. SPE-21213-PA. http://dx.doi.org/10.2118/21213-PA.

[r184-184] Gautam S. et al. 1996. Parallel Computing Alters Approaches, Raises Integration Challenges in Reservoir Modeling. Oil & Gas J. (20 May): 48.

[r185-185] Tchelepi, H.A., Durlofsky, L.J., Chen, W.H. et al. 1997. Practical Use of Scale Up and Parallel Reservoir Simulation Technologies in Field Studies. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 5-8 October 1997. SPE-38886-MS. http://dx.doi.org/10.2118/38886-MS.

[r186-186] Chien, M.C.H. et al.: "A Scalable parallel Multi-Purpose Reservoir Simulator," paper SPE 37976 presented at the 1997 SPE Reservoir Simulation Symposium, Dallas, 8–11 June.

[r187-187] Shiralkar, G.S., Stephenson, R.E., Joubert, W. et al. 1997. Falcon: A Production Quality Distributed Memory Reservoir Simulator. Presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, 8-11 June 1997. SPE-37975-MS. http://dx.doi.org/10.2118/37975-MS.

[r188-188] Wang, P., Yotov, I., Wheeler, M. et al. 1997. A New Generation EOD Compositional Reservoir Simulator: Part 1—Formulation and Discretization. Presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, 8–11 June. SPE-37979-MS. http://dx.doi.org/10.2118/37979-MS.

[r189-189] 189.0 ^189.1 ^189.2 Pavlas Jr., E.J. 2001. MPP Simulation of Complex Water Encroachment in a Large Carbonate Reservoir in Saudi Arabia. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71628-MS. http://dx.doi.org/10.2118/71628-MS.

[r190-190] 190.0 ^190.1 Wallis, J.R. and Nolen, J.S. 1993. Efficient Iterative Linear Solution of Locally Refined Grids Using Algebraic Multilevel Approximate Factorizations. Presented at the SPE Symposium on Reservoir Simulation, New Orleans, Louisiana, 28 February-3 March 1993. SPE-25239-MS. http://dx.doi.org/10.2118/25239-MS.

[r191-191] Song, T. 1996. A Load-Balancing Technique for Reservoir Simulation Based on Dantzig’s Transportation Model. MS thesis, U. of Houston, Houston.

[r192-192] Reference Manual for Zoltan. 1998. Los Alamos Natl. Laboratory, Los Alamos, New Mexico.

[r193-193] Bhogeswara, R. and Killough, J.E. 1994. Parallel Linear Solvers for Reservoir Simulation: A Generic Approach for Existing and Emerging Computer Architectures. SPE Comp App 6 (1): 5-11. SPE-25240-PA. http://dx.doi.org/10.2118/25240-PA.

[r194-194] Guerberoff, D., Zucchi, H., Victoria, M. et al. 2001. Lateral Delineation of Sandstone Bodies Guided by Seismic and Petrophysical Data Using Geocellular Model: Geocellular Model: Canadon Secon Formation, San Jorge Basin, Argentina. Presented at the SPE Latin American and Caribbean Petroleum Engineering Conference, Buenos Aires, Argentina, 25-28 March 2001. SPE-69487-MS. http://dx.doi.org/10.2118/69487-MS.

[r195-195] Agarwal, B., Hermansen, H., Sylte, J.E. et al. 2000. Reservoir Characterization of Ekofisk Field: A Giant, Fractured Chalk Reservoir in the Norwegian North SeaHistory Match. SPE Res Eval & Eng 3 (6): 534-543. SPE-68096-PA. http://dx.doi.org/10.2118/68096-PA.

[r196-196] Caers, J., Avseth, P., and Mukerji, T. 2001. Geostatistical integration of rock physics, seismic amplitudes and geological models in North-Sea turbidite systems. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71321-MS. http://dx.doi.org/10.2118/71321-MS.

[r197-197] 197.0 ^197.1 Qassab, H.M.A., Rahmeh, B.A., Khalifa, M.A.A. et al. 2001. Conditioning Integrated Geological Models to Dynamic Flow Data of Giant Saudi Arabian Reservoir. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71319-MS. http://dx.doi.org/10.2118/71319-MS.

[r198-198] 198.0 ^198.1 Raghavan, R., Dixon, T.N., Robinson, S.W. et al. 2000. Integration of Geology, Geophysics, and Numerical Simulation in the Interpretation of A Well Test in a Fluvial Reservoir. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 1-4 October 2000. SPE-62983-MS. http://dx.doi.org/10.2118/62983-MS.

[r199-199] Durlofsky, L.J. 2005. Upscaling and gridding of fine-scale geological models for flow simulation. Proc., 2005 Intl. Forum on Reservoir Simulation, Stresa, Italy, 20–25 June.

[r200-200] Peaceman, D.W. 1997. Effective Transmissibilities of a Gridblock by Upscaling - Comparison of Direct Methods with Renormalization. SPE J. 2 (3): 338-349. SPE-36722-PA. http://dx.doi.org/10.2118/36722-PA.

[r201-201] Flodin, E.A., Aydin, A., Durlofsky, L.J. et al. 2001. Representation of fault zone permeability in reservoir flow models. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September–3 October. SPE-71617-MS. http://dx.doi.org/10.2118/71617-MS.

[r202-202] Nishikiori, N. and Hayashida, Y. 2000. Investigation of Fluid Conductive Faults and Modeling of Complex Water Influx in the Khafji Oil Field, Arabian Gulf. SPE Res Eval & Eng 3 (5): 401-407. SPE-66223-PA. http://dx.doi.org/10.2118/66223-PA.

[r203-203] Barker, J.W. and Thibeau, S. 1997. A Critical Review of the Use of Pseudorelative Permeabilities for Upscaling. SPE Res Eng 12 (2): 138–143. SPE-35491-PA. http://dx.doi.org/10.2118/35491-PA.

[r204-204] Beckner, B.L., Hutfilz, J.M., Ray, M.B. et al. 2001. EMpower: New Reservoir Simulation System. Presented at the SPE Middle East Oil Show, Bahrain, 17–20 March. SPE-68116-MS. http://dx.doi.org/10.2118/68116-MS.

[r205-205] Peaceman, D.W. 1983. Interpretation of Well-Block Pressures in Numerical Reservoir Simulation With Nonsquare Grid Blocks and Anisotropic Permeability. SPE J. 23 (3): 531–543. SPE-10528-PA. http://dx.doi.org/10.2118/10528-PA.

[r206-206] Mlacnik, M.J. and Heinemann, Z.E. 2003. Using Well Windows in Full Field Reservoir Simulation. SPE Res Eval & Eng 6 (4): 275–285. SPE-85709-PA. http://dx.doi.org/10.2118/85709-PA.

[r207-207] Wolfsteiner, C., Durlofsky, L.J., and Aziz, K. 2000. Approximate Model for Productivity of Nonconventional Wells in Heterogeneous Reservoirs. SPE J. 5 (2): 218–226. SPE-56754-PA. http://dx.doi.org/10.2118/62812-PA.

[r208-208] Valvatne, P.H., Durlofsky, L.J., and Aziz, K. 2001. Semi-Analytical Modeling of the Performance of Intelligent Well Completions. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 11-14 February 2001. SPE-66368-MS. http://dx.doi.org/10.2118/66368-MS.

[r209-209] Litvak, M.L., Clark, A.J., Fairchild, J.W. et al. 1997. Integration of Prudhoe Bay Surface Pipeline Network and Full Field Reservoir Models. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, 5-8 October. SPE 38885. http://dx.doi.org/10.2118/38885-MS.

[r210-210] Fang, W.Y. and Lo, K.K. 1996. A Generalized Well-Management Scheme for Reservoir Simulation. SPE Res Eng 11 (2): 116–120. SPE-29124-PA. http://dx.doi.org/10.2118/29124-PA.

[r211-211] Litvak, M.L. and Wang, C.H. 2000. Simplified Phase-Equilibrium Calculations in Integrated Reservoir and Surface-Pipeline-Network Models. SPE J. 5 (2): 236-241. SPE-64498-PA. http://dx.doi.org/10.2118/64498-PA.

[r212-212] Hermansen, H., Thomas, L.K., Sylte, J.E. et al. 1997. Twenty Five years of Ekofisk Reservoir Management. Presented at the Annual Technical Conference and Exhibition, San Antonio, Texas, 5–8 October. SPE-38927-MS. http://dx.doi.org/10.2118/38927-MS.

[r213-213] Sylte, J.E., Thomas, L.K., Rhett, D.W. et al. 1999. Water Induced Compaction in the Ekofisk Field. Presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, 3–6 October. SPE-56426-MS. http://dx.doi.org/10.2118/56426-MS.

[r214-214] Ashton, K., Sylte, J.E., Thomas, L.K. et al. 1998. Judy/Joanne Field Development. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 27-30 September 1998. SPE-49128-MS. http://dx.doi.org/10.2118/49128-MS.

[r215-215] Pieters, J. and Por, G.J.A. 1995. Total system modelling - a tool for effective reservoir management of multiple fields with shared facilities. Presented at the Offshore Europe, Aberdeen, United Kingdom, 5-8 September 1995. SPE-30442-MS. http://dx.doi.org/10.2118/30442-MS.

[r216-216] Williams, M.A., Keating, J.F., and Barghouty, M.F. 1998. The Stratigraphic Method: A Structured Approach to History-Matching Complex Simulation Models. SPE Res Eval & Eng 1 (2): 169-176. SPE-38014-PA. http://dx.doi.org/10.2118/38014-PA.

[r217-217] Ditzhuijzen, R.v., Oldenziel, T., and Kruijsdijk, C.P.J.W.v. 2001. Geological Parameterization of a Reservoir Model for History Matching Incorporating Time-Lapse Seismic Based on a Case Study of the Statfjord Field. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71318-MS. http://dx.doi.org/10.2118/71318-MS.

[r218-218] Bogan, C., Johnson, D., Litvak, M. et al. 2003. Building Reservoir Models Based on 4D Seismic & Well Data in Gulf of Mexico Oil Fields. Presented at the SPE Annual Technical Conference and Exhibition, Denver, Colorado, 5-8 October 2003. SPE-84370-MS. http://dx.doi.org/10.2118/84370-MS.

[r219-219] Milliken, W.J., Emanuel, A.S., and Chakravarty, A. 2000. Applications of 3D Streamline Simulation to Assist History Matching. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 1-4 October 2000. SPE-63155-MS. http://dx.doi.org/10.2118/63155-MS.

[r220-220] Schulze-Riegert, R.W., Axmann, J.K., Haase, O. et al. 2001. Optimization Methods for History Matching of Complex Reservoirs. Presented at the SPE Reservoir Simulation Symposium, Houston, 11-14 February. SPE 66393. http://dx.doi.org/10.2118/66393-MS.

[r221-221] Grussaute, T. and Gouel, P. 1998. Computer Aided History Matching of a Real Field Case. Presented at the European Petroleum Conference, The Hague, Netherlands, 20-22 October 1998. SPE-50642-MS. http://dx.doi.org/10.2118/50642-MS.

[r222-222] Zapata, V.J., Brummett, W.M., Osborne, M.E. et al. 2001. Advances in Tightly Coupled Reservoir/Wellbore/Surface-Network Simulation. SPE Res Eval & Eng 4 (2): 114–120. SPE-71120-PA. http://dx.doi.org/10.2118/71120-PA.

[r223-223] 223.0 ^223.1 Gorell, S. and Bassett, R. 2001. Trends in Reservoir Simulation: Big Models, Scalable Models? Will you Please Make up Your Mind? Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71596-MS. http://dx.doi.org/10.2118/71596-MS.

[r224-224] Zabalza-Mezghani, I., Mezghani, M., and Blanc, G. 2001. Constraining Reservoir Facies Models to Dynamic Data - Impact of Spatial Distribution Uncertainty on Production Forecasts. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September-3 October 2001. SPE-71335-MS. http://dx.doi.org/10.2118/71335-MS.

[r225-225] Manceau, E., Mezghani, M., Zabalza-Mezghani, I. et al. 2001. Combination of Experimental Design and Joint Modeling Methods for Quantifying the Risk Associated With Deterministic and Stochastic Uncertainties—An Integrated Test Study. Presented at the SPE Annual Conference and Technical Exhibition, New Orleans, Louisiana, 30 September–3 October. SPE-71620-MS. http://dx.doi.org/10.2118/71620-MS.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

[39]

[40]

[41]

[42]

[43]

[44]

[45]

[46]

[47]

[48]

[49]

[50]

[51]

[52]

[53]

[54]

[55]

[56]

[57]

[58]

[59]

[60]

[61]

[62]

[63]

[64]

[65]

[66]

[67]

[68]

[69]

[70]

[71]

[72]

[73]

[74]

[75]

[76]

[77]

[78]

[79]

[80]

[81]

[82]

[83]

[84]

[85]

[86]

[87]

[88]

[89]

[90]

[91]

[92]

[93]

[94]

[95]

[96]

[97]

[98]

[99]

[100]

PEH:Reservoir Simulation

Overview

The Generalized Model

Model Formulations

Advances in Model Forumlations

Stable Timestep and Switching Criteria

The Linear Solver

Cartesian Grids and Reservoir Definition

Linear Solver—John Wallis and J.W. Watts

Gridding in Reservoir Simulation

Introduction

Regular Cartesian Grids

Hexahedral Grids

Multiple-Domain Hexahedral Grids

Grid Refinement

Unstructured Grids

Truncated Regular Grids[69]

Other Gridding Systems

Future Directions

Upscaling of Grid Properties—Alan Grindheim

Definition

Can Upscaling Be Avoided?

Upscaling Techniques for Absolute Permeability

Upscaling Schemes for Absolute Permeability

Local Upscaling

Regional Upscaling

Global Upscaling

Best Practice Guidelines

Streamline Simulation—Rod P. Batycky and Marco R. Thiele

Introduction

Applicability of Streamline Simulation

Historical

Mathematics of the Streamline Method

The Computational Efficiency of Streamlines

Novel Data Produced by Streamlines

Applications of Streamlines

The Future of Streamline Simulation

Simulation of Geomechanics—Tony Settari

Types of Coupling

Modeling of Reservoir Compaction and/or Dilation

Modeling of Stress-Dependent Flow Properties

Types of Coupled Models

Future Trends and Needs

Pressure/Volume/Temperature (PVT) Treatment—Curtis H. Whitson

Number of Phases and Phase Type

Two-Phase Flash

Density

Component Partitioning

Viscosity

Gas/Oil Interfacial Tension (IFT)

Black-Oil PVT Models

Equation-of-State Models

Three-Phase PVT Behavior

Surface Phase Behavior

Thermal Model PVT Requirements

Fluid Initialization

High-Performance Computing and Reservoir Simulation—John E. Killough

The Exposition of Data and the Chasm of Scale

Characteristics of High-Performance Computing

The Parallel Reservoir Simulator

The Parallel Linear-Equation Solver

Parallel Computers and Message-Passing Interface

Application of Parallel Computing

The Future of High-Performance Simulation

Reservoir Simulation Applications—L. Kent Thomas

Development of the Geological Model

Upscaling Geological Model to Reservoir Flow Model

Inclusion of Faults in Reservoir Flow Model

Development of Pseudofunctions for Multiphase Flow

Gridding Techniques

Simulation of Nonconventional and Intelligent Wells

Integrated Reservoir and Surface Facilities Models

Simulation of Multiple Reservoirs

Use of Larger Models

History Matching and Production Forecasting

Nomenclature

Subscripts

Superscripts

References

SI Metric Conversion Factors

Truncated Regular Grids^[69]