Upscaling and Simulation of Waterflooding in Heterogeneous ...

Transp Porous Med (2008) 72:311–338 DOI 10.1007/s11242-007-9152-1

Upscaling and Simulation of Waterflooding in Heterogeneous Reservoirs Using Wavelet Transformations: Application to the SPE-10 Model Mohammad Reza Rasaei · Muhammad Sahimi

Received: 26 September 2006 / Accepted: 13 June 2007 / Published online: 20 July 2007 © Springer Science+Business Media B.V. 2007

Abstract We have developed an accurate and highly efficient method for upscaling and simulation of immiscible displacements in three-dimensional (3D) heterogeneous reservoirs, which is an extension of the technique that we developed previously for 2D systems. The method utilizes wavelet transformations (WTs) to upscale the geological model of a reservoir, based on the spatial distribution of the single-phase permeabilities and the locations of the wells in the reservoir. It generates a non-uniform grid in which the resolved structure of the fine grid around the wells, as well as in the high-permeability sectors, are preserved, but the rest of the grid is upscaled. A robust uplayering procedure is used to reduce the number of the layers, and the WTs are used to upscale each layer areally. To demonstrate the method’s accuracy and efficiency, we have applied it to the geological model of a highly heterogeneous reservoir put forward in the tenth Society of Petroleum Engineers comparative solution project (the SPE-10 model), and carried out simulation of waterflooding in the upscaled model. Various upscaling scenarios were examined, and although some of them resulted in efficient simulations and accurate predictions, the results when non-uniform upscaling is used based on the WT technique are in excellent agreement with the solution of the same problem in the fine grid of the SPE-10 model. Most importantly, the speed-up factors that we obtain are several orders of magnitude. Hence, the method renders it unnecessary to use massively parallel computations for such problems. Keywords Two-phase flows · Geological model · Upscaling · Wavelet transformations · SPE-10 model

M. Reza Rasaei Institute for Petroleum Engineering, University of Tehran, Tehran 11365-4563, Iran M. Reza Rasaei · M. Sahimi (B) Mork Family Department of Chemical Engineering & Materials Science, University of Southern California, Los Angeles, CA 90089-1211, USA e-mail: [email protected]

123

312

M. R. Rasaei, M. Sahimi

1 Introduction Progress in developing various techniques for measuring or estimating the properties of oil reservoirs, and advances in the development of methods for their characterization, when combined with modern geostatistical techniques, enable one to generate a highly resolved computational grid, with up to a few million grid blocks, that represents the geological model (GM) of the reservoirs. Computer simulation of multiphase flows in such models entails carrying out highly intensive computations (Peacemen 1977; Aziz and Settari 1979; Thomas 1982; Sahimi 2007), so much so that it is usually not feasible to carry out the simulations using such grids. It is, therefore, necessary to upscale the properties of the grid blocks of the GM in order to develop a coarsened grid that can be used in reservoir simulation with an affordable amount of computation time, while preserving the accuracy of the results. For simulation of single-phase flows, the most important parameter is the distribution of the absolute permeabilities of the grid blocks, for which upscaling methods are well-established. Simulation of multiphase flows, on the other hand, involves, in addition to computing the equivalent absolute permeabilities of the upscaled blocks, adjustments to the flow of the fluid phases through the connected blocks of the upscaled grid. To accomplish this task many methods have been proposed, a brief review of most of which was given in a recent article (Sahimi et al. 2005) and, therefore, are not described here. Significant progress was made when Durlofsky et al. (1997) introduced a method whereby finer resolution is used in the regions of high-fluid velocities, and upscaled, homogenized description is utilized for the rest of the domain. In their approach no upscaling scheme is used for the relative permeabilities, as the original rock curves are used for the upscaled grid blocks, hence making the technique process-independent. Non-uniform grid coarsening offers significant computational speed-up and general applicability, but at the cost of increasing numerical dispersion and decreasing accuracy for the nonseparable scales. Such upscaling techniques are suitable for the sectors that are far from the wells. Since the pressure field in the near-well regions usually changes rapidly in the radial direction, the above upscaling approaches are not suitable. For such situations, well pseudo-functions to account for the pressure changes were suggested by Chappelear and Hirasaki (1976) and King et al. (1991). Durlofsky et al. (2000) developed a method for calculating the transmissibility and well index for single-phase flow based on the solution to the local well-driven flow. Grid selection methods that rely on dynamic responses, such as streamline simulation, in order to identify the locations of the grid blocks in the GM through which most of the fluids pass, were also proposed by Verma and Aziz (1996) and Castellini et al. (2000). Static methods, suggested by Garcia et al. (1992), Li and Beckner (2000), and Younis and Caers (2001, 2002), that rely on the spatial distribution of the permeability for developing upscaled grids, are also robust. Recently, we proposed a new approach for grid coarsening and upscaling of single-phase permeabilities, based on the use of wavelet transformations (WTs) (Mehrabi and Sahimi 1997; Ebrahimi and Sahimi 2002). The method is a multiscale or multiresolution approach that acts as an automatic grid generator at all the relevant length scales that are incorporated in the GM. In three recent articles we demonstrated the accuracy and efficiency of the technique by applying it to simulations of unstable miscible displacements (Ebrahimi and Sahimi 2004) and two-phase flows in two-dimensional (2D) models of heterogeneous reservoirs (Sahimi et al. 2005; Rasaei and Sahimi 2007). The goal of the present article is to further demonstrate the method’s virtues by applying it to a 3D model of a highly heterogeneous reservoir used by the Society of Petroleum Engineers (Christie and Blunt 2001; see below) in a competition for upscaling techniques. The model is considered a benchmark for testing the efficiency and accuracy any upscaling or coarsening technique. Since the model is 3D, application of our

123

Upscaling and Simulation of Waterflooding in Heterogeneous Reservoirs

313

method to it also entails extending it significantly. As before, we show that the method results in several orders of magnitude speed-up in the computations without any loss of accuracy. Applications of the WTs to problems in oil reservoirs and other types of porous media have recently been increasing. In addition to our own work mentioned above, and also analyzing seismic data which was one of the prime motivations for developing the WTs, such applications also include numerical simulation of flow and transport equations (Moridis et al. 1996), pressure transient analysis (Kikani and He 1998; Athichanagorn et al. 1999; Soliman et al. 2001; Ebrahimi 2002; Ebrahimi and Sahimi 2006), reservoir characterization (Brewer and Wheatcraft 1994; Panda et al. 1996; Mehrabi et al. 1997; Lu and Horne 2000; Neupauer and Powell 2005), determination of the fracture density in oil and gas reservoirs (Sahimi and Hashemi 2001), as well as problems involving turbulent transport and chemical reactions in the atmosphere (Heidarinasab et al. 2004). We point out that Chu et al. (1996) also attempted to use the WTs for upscaling using a method different from what we have developed in the past and extend in the present article. The rest of this article is organized as follows. In the next section we describe the details of the geological model that we upscale and use in the simulation of the waterflooding problem. Section 3 describes the wavelet technique for upscaling of a given grid, while Sect. 4 describes briefly the coarsening of the computational grid in the vertical direction and calculation of the effective permeabilities of the coarsened blocks. Numerical simulation of the waterflooding problem that we carry out in the upscaled grids is described in Sect. 5, while representation of the wells in the grid is briefly discussed in Sect. 6. The results are described and discussed in Sect. 7, while the efficiency of the method is discussed in Sect. 8. The article is summarized in Sect. 9, where we also discuss further refinement of the method.

2 Details of the SPE-10 Model We apply our WT technique to Model 2 of the 10th Society of Petroleum Engineers Comparative Solution Project (Christie and Blunt 2001), hereafter referred to as the SPE-10 model. The model represents a very heterogeneous reservoir, and is represented by a regular Cartesian grid. Its dimensions are 1,200 × 2,200×170 ft. It consists of part of a Brent sequence, with the top 70 ft (35 layers) representing the Tarbert formation, which is a representation of a prograding near shore environment, and the bottom 100 ft (constituting 50 layers) representing the Upper Ness, which is fluvial. The fine-scale (geological) model contains 60×220×85 blocks or cells (for a total of 1,122,000 cells), with the blocks’ size at the fine scale being 20 × 10 × 2 ft. The porosity and permeability maps were downloaded from the Project’s web site at: http://www.spe.org/csp In the model the porosity varies between 0 and 0.5 (see Fig. 1). The horizontal permeability K h varies over eight orders of magnitude, between 6.65 × 10−4 to 2 × 104 md. The reservoir contains channeling paths, with K v /K h = 0.3 in the channels, and K v /K h = 10−4 in the background, where K v is the vertical permeability. The model has no top structure or fault, to allow maximum flexibility in selecting the upscaled grids. The water properties used in the simulations are, Bw = 1.01, Cw = 3 × 10−6 psi−1 , and µw = 0.3 cp. The dead oil pressure/volume/temperature (PVT) data are shown in Table 1. The relative permeability functions are given by krw =

S − Swc 1 − Swc − Sor

2 ,

(1)

123

314


Channeling structure

Fig. 1 The porosity distribution of the SPE-10 model Table 1 PVT data for the dead oil

P (psi)

Bo

µ0 (cp)

300 800 8000

1.05 1.02 1.01

2.85 2.99 3.0

kro =

1 − S − Sor 1 − Swc − Sor

2 ,

(2)

with the residual and irreducible saturations of the water and oil phases being, Swc = Swi = Sor = 0.2. All the wells are vertical and completed throughout the formation. The central injector has an injection rate of 5,000 bbl/day (under reservoir conditions) and a maximum injection bottom hole pressure of 10,000 psi. There are four producers in the four corners of the model, each producing at 4,000 psi bottom-hole pressure. The initial reservoir pressure is set at 6,000 psi, and the simulation period is 2,000 days. All the data were utilized in the waterflooding simulations, when necessary. This very heterogeneous GM is now upscaled via the WT technique that we have developed, and extend in this article to 3D.

3 Wavelet Coarsening of a Computational Grid For mathematical foundations of the WTs see, Mallat (1999), while Nievergelt (1999) provides an excellent introduction. Wavelet transformation of the spatially varying single-phase permeability k(x) (also called its wavelet detail coefficient), is defined by ∞ ∞ 1 D(a, b) = k(x)ψab (x)dx = √ k(x)ψ[(x − b)/a]dx, (3) a −∞ −∞ where ψ(x)—the mother wavelet—is selected from a family of such functions, hence providing one with great flexibility in analyzing the spatial distribution of k(x) at many length scales. Here, a > 0 is a rescaling parameter and b represents translation of the wavelet. Equation 3 indicates that, using WT of k(x), one can analyze the permeability distribution

123


315

at increasingly coarser (a > 1) or finer (a < 1) length scales. D(a, b) contains information on the differences between two approximations of k(x) in two successive length scales. The most accurate approximation of k(x) at a fixed scale is obtained by computing the wavelet approximate or scale coefficient, defined by +∞ S(a, b) = φab (x)k(x)dx, (4) −∞

where the wavelet scale-function φ(x) is orthogonal to ψ(x), and the definition of φab (x) is similar to that of ψab (x). An important property of the WTs is that they are recursive, so that they can be applied in succession to any set of properties produced by using the wavelets, to generate another level of averages and another level of details. These properties are utilized in our method (see below). In the computations described below, upscaled porosities were calculated as simple arithmetic averages of those of the blocks in the fine-grid model of the GM. Upscaled saturations were computed as porosity-weighted averages. We utilize the wavelet technique to coarsened the grid areally. That is, each horizontal layer is coarsened (see Sects. 4 and 7 for its practical implementation) in order to reduce the number of the blocks and make the computations efficient. Fully 3D wavelet technique can also be used (Pazhoohesh et al. 2006), but we do not utilize it in the present article (see also below). We consider the 2D square grid in each layer or horizontal plane. Each of the grid’s square blocks is assigned the permeabilities k x and k y . We then apply a one-level discrete WT to the permeabilities, i.e., we attempt to upscale the grid by a factor of 2, but the coarsening is not done uniformly. Associated with the WT of the permeabilities of a block with its center at x = (i 1 , i 2 ) is a set of four wavelet coefficients that consists of one wavelet scale coefficient and three wavelet detail coefficients. Following Eqs. 3 and 4, the four wavelet coefficients are given, in more precise forms, by (d) (d) D j (i 1 , i 2 ) = k(x, y)ψ j,i1 ,i2 (x, y)d xd y, (5)

S j (i 1 , i 2 ) =

k(x, y)φ j,i1 ,i2 (x, y)d xd y,

(6)

where j is the level of coarsening ( j = 1 for the GM), and represents the grid. As mentioned above, S j (i 1 , i 2 ) contains information about k(x) associated with a block centered (d) at x in the coarser grid, while D j (i 1 , i 2 ) measure the contrast between k(x) in the coarser scale and those of the block’s neighbors in the previous finer scale grid, with d = 1, 2, and 3 measuring the contrasts between the blocks in the x, y, and the diagonal directions, respectively. Two thresholds, s and d , are then introduced. s is a measure of the blocks’ permeabilities, with which is associated a corresponding wavelet scale coefficient, and is set as a fraction of the largest scale coefficient in the entire 2D grid. The second threshold, d , measures the contrast between the permeabilities of the neighboring blocks and, similar to s , is set as a fraction of the largest detail coefficient in the grid. To carry out the upscaling process, we compare the threshold s with the scale coefficient associated with each block represented by its center at (i 1 , i 2 ). If S j (i 1 , i 2 ) > s , implying that the block’s permeability (k x or k y ) is large enough, we move on to the next block. If, however, S j (i 1 , i 2 ) < s , we check the (d) block’s associated detail coefficients D j (i 1 , i 2 ) and set to zero all such coefficients that are smaller than d , implying that the neighbor of the block centered at (i 1 , i 2 ) corresponding to the direction (d), which in the finer-scale grid is just one block (or one diagonal block) away

123

316


from (i 1 , i 2 ), is merged with the block centered at (i 1 , i 2 ) to form a larger block. Therefore, depending on the structure of the spatial distribution of k(x) and the numerical values of the two thresholds, a number of blocks in the original fine-scale grid are coarsened. Since, as pointed out earlier, the WTs have recursive properties, the resulting grid is coarsened further by applying again the discrete WT to the scale coefficients obtained at the previous level which contain information about k(x) in the present upscaled grid, and calculating a new set of four wavelet coefficients for each block of the grid in its current state. The new detail coefficients are again set to zero if they are smaller than the threshold d , and the corresponding blocks in the grid are merged. In effect, at upscaling level j, each block is compared with those that are at a distance 2 j−1 from it, where the distance is measured in units of the blocks’ length in the original grid. The process is repeated again until no significant number of the blocks is upscaled. Typically, and for fixed d and s , after 3 or 4 levels of coarsening the grid is no longer effectively upscaled, thus yielding very efficiently the final upscaled grid. Clearly, the higher the two thresholds d and s , the larger the number of the blocks that can be upscaled. The thresholds are fixed by the amount of computation time that one can afford. The effective permeabilities of the coarsened blocks are then computed by the methods described below. To carry out these computations, one must use 2D wavelets. One way of constructing 2D wavelets, which is used in the present article, is by tensor products of the 1D wavelets, j j which gives rise to one 2D wavelet scaling function, φ j,i1 ,i2 (x, y) = φi1 (x)φi2 (y), and three (1)

j

j

(2)

j

j

(3)

wavelets, ψ j,i1 ,i2 (x, y) = φi1 (x)ψi2 (y), ψ j,i1 ,i2 (x, y) = ψi1 (x)φi2 (y), and ψ j,i1 ,i2 (x, y) = j

j

ψi1 (x)ψi2 (y). The 1D discrete wavelets themselves are constructed by setting a = 2 j and b = 2 j i in Eqs. 3 and 4, where i and j are both integer. By rescaling and translating ψ(x) j j and φ(x) using ψi (x) = 2− j/2 ψ(2− j x − i) and φi (x) = 2− j/2 φ(2− j x − i), one can also construct orthonormal wavelets that are nonzero over only small intervals of x. An important set of wavelets, developed by Daubechies (1988), consists of orthonormal wavelets of order M (usually referred to as DBM) which possess the properties that their first M moments are zero, and that φ(x) =

L−1 √ 2 h i φ(2x − i),

(7)

i=0

ψ(x) =

L−1 √ 2 gi φ(2x − i),

(8)

i=0

where L = 2M. Here, h i and gi —called the filter coefficients—are related by, gi = (−1)i h L−i−1 , with i = 0, 1, . . . , L − 1; they are usually nonzero for only √ a few√values of for the DB2 wavelet, one has, (h 0 , h 1 , h 2 , h 3 ) = (1/4 2)(1 + 3, 3 + √ √ √ i. For example, 3, 3 − √ 3, 1 − 3), and (g0 , g1 , g2 , g3 ) = (−h 3 , h 2 , −h 1 , h 0 ) (in both cases the extra factor 1/ 2 is included to ensure the orthonormality of φ and ψ). Press et al. (1992) list values of h k for many other wavelets. Based on our previous simulations and experience for 2D models (Mehrabi and Sahimi 1997; Ebrahimi and Sahimi 2002, 2004; Sahimi et al. 2005; Rasaei and Sahimi 2007), the results reported in this article were obtained by using the DB4 wavelet. Figure 2 presents wavelet√and the√scaling Its √ DB4 √ √ functions. √ √ filter coefficients √the 1D (h 0 , h 1 , h 2 , h 3 ) are, 18 [( 2 + 6), (3 2 + 6), (3 2 − 6), ( 2 − 3)]. As a test, we also used the Haar wavelet—also known as the DB1 wavelet—for which, ψ0,1 (x) = 1 and

123


317

Fig. 2 The DB4 scaling (top) and wavelet (bottom) functions

−1 for 0 ≤ x ≤ 1/2 and 1/2 ≤ x < 1, respectively, and ψ0,1 (x) = 0 for all other values of √ x. Its filter coefficients are, (h 0 , h 1 ) = (1/ 2)(1, 1).

4 Uplayering of the Geological Model To carry out a thoroughly integrated upscaling, we must try to utilize all the information available about the reservoir, such as its structure, the distributions of its petrophysical parameters, the thermodynamic properties of the reservoir’s fluids, and the wells’ locations and operating conditions. The first operation to be carried out is uplayering; that is, reducing the number of layers (and, hence, the number of horizontal planes in the computational grid), which can be performed both uniformly and non-uniformly, depending on the initial number of the geological layers and the final desired number of layers in the upscaled model (determined by the computation times that can be afforded). For example, in uniform uplayering, every few layers may be commingled together to form one thicker layer in the reservoir model, without any consideration for the geological and/or petrophysical data; the uplayering

123

318


is done just based on the level of the details that one wishes to achieve in the computational grid. This is often not the best practice since layers with very different properties might mix together, resulting in over-smoothing of the layers’ heterogeneities which, in turn, may have strong negative influence on the true fluid flow pattern in the reservoir. In this article we carry out non-uniform uplayering based on the directional effective permeabilities of each layer, as well as the way they vary laterally within the geological layers. For the sake of comparison, we also carry out uniform uplayering. Thus, we consider three different non-uniform uplayering scenarios: (1) The permeabilities in each layer in the x, y, and z directions are considered, and the layers with extreme permeability values are considered either separately or only low levels of uplayering are carried out for them. (2) One may keep the layers with the most lateral changes intact, and compute the first-level wavelet coefficients of each layer (see above for details of the calculations). Then, layers with the maximum wavelet coefficients (i.e., maximum lateral changes in the layer’s permeability distribution; see above) are considered separately, or only a low level of uplayering is carried out for them. (3) One may combine and upscale the layers based on their injectivity indices (see below). Carrying out the uplayering by one of such scenarios, we end up with fewer new layers that have been upscaled vertically. They are, therefore, ready to be upscaled areally one by one. Areal upscaling of each layer is then done by the 2D wavelet technique, as described above. This yields a non-uniform upscaled grid which has preserved its resolved structure around the wells (see below) and the high-permeability zones in each layer, but upscaled to different degrees in the rest of the reservoir. The next step after construction of the coarsened layers and identification of blocks’ sizes everywhere is calculating the equivalent properties of each coarsened block. While the procedures for calculating the upscaled porosity, fluids’ saturation, and block pressures are well-established based on volume averaging, the same is not as easy and straightforward for the blocks’ permeabilities. Many approaches have been proposed for computing the upscaled absolute permeabilities, ranging from simple harmonic and geometric averages to complex global pressure-solver methods (see, for example, Wen and Gomez-Hernandez 1996; Renard and de Marsily 1997). We used two different approaches for calculating the effective permeabilities. One was based on a combination of the arithmetic-harmonic averages. For this purpose, the analytical relations given by Henriette et al. (1989) and Li et al. (2001) can be used to calculate layers’ effective permeabilities in each direction. The horizontal effective permeabilities are approximately equal to the harmonic mean of the arithmetic averages of the fine-grid permeabilities: Lx

Kx = L y Lz

Nx

xi

,

(9)

.

(10)

N y Nz i=1 y j z i, j,k k x,i, j,k j=1 k=1

Ly

Ky = Lx Lz

Ny j=1

yi Nz Nx i=1 k=1

123

xi z i, j,k k y,i, j,k


319

The vertical effective permeability for each layer is approximately equal to the arithmetic average of the harmonic means of the fine-grid vertical permeabilities:

Kz =

Ny Nx xi yi z i, j 1 . Nz Lx Ly z i, j,k i=1 j=1

k=1

(11)

k z,i, j,k

Here, L x , L y , and L z are the layers’ dimensions, while x, y, and z represent the corresponding values for each block in the fine grid, with k x , k y , and k z being their permeabilities. N x , N y , and Nz are the number of blocks in the fine grid in the three directions in each layer (Nz = 1 for a single layer). The second method, utilizing the inverse of WT, is based on the reconstruction of the permeability field in each horizontal layer, i.e., computing its inverse WT after some of its detail coefficients are set to zero (i.e., after some of the blocks have been merged to form larger upscaled blocks; see above). The horizontal permeabilities of the upscaled blocks are then assigned based on reconstructed permeability field. The blocks’ vertical permeabilities, located in the combined layers, are calculated using Eq. 11; see below for further explanations. We found that the reconstruction method yields more accurate results. Accordingly, it is this method that we used to calculate the effective permeabilities of the upscaled blocks in the horizontal planes. Let us point out that one can use the wavelet technique for upscaling the entire 3D grid, rather than the method that we use in this article, namely, using 2D WTs to upscale the horizontal planes, and apply a separate uplayering in the vertical direction. To do so, one must use fully 3D wavelet functions, and define separate wavelet thresholds for the vertical direction and the horizontal planes. Elsewhere, we have utilized (Pazhoohesh et al. 2006) such a method for upscaling of a computational grid used for computing the effective hopping conductivity of amorphous semiconductors. As we show below, for the present problem, use of 2D wavelets for coarsening the grid in the horizontal planes, together with an appropriate scheme for vertical uplayering, suffice for accurate and efficient computations while, at the same time, being simpler than a full 3D wavelet method. However, use of a method based on 3D wavelets is completely general, and should be applicable to most models of oil reservoirs. Let us also point out an important point regarding the efficiency of the wavelet method: the broader the permeability distribution, the more efficient is the method. The reason is due to the fact that, if the heterogeneities are relatively mildly distributed, then, significant fluid flow takes place through many sectors of the reservoir. If, on the other hand, the medium is highly heterogeneous, then only a small fraction of the reservoir contributes significantly to fluid flow. This important concept, known in physics of heterogeneous media as the critical-path analysis, has proven to be very powerful for analyzing the properties of a wide variety of highly disordered media (Sahimi 2003). As a result, the number of grid blocks in the upscaled model, when the heterogeneities are mildly distributed, is larger than those in a highly heterogeneous medium. We recently demonstrated (Pazhoohesh et al. 2006) this aspect of the method for simulation of frequency-dependent conduction in 3D amorphous semiconductors, in which the local conductivities are distributed over up to seventeen orders of magnitude. Such a distribution is broader, by orders of magnitude, than the permeability distribution in the SPE-10 model that we work with in the present article. The wavelet method reduced the number of grid blocks from 2,097,152 in the fine-grid model of the material to only 630. The difference between the effective conductivities of the two models, when the frequencies were varied over ten orders of magnitude, was at most 5%.

123

320


5 Numerical Simulations The problem to be solved in the SPE-10 model reservoir is a 3D waterflood. Water is injected into a well at the center of the reservoir, and oil is produced from four wells at the corners. Thus, combining the Darcy’s law and the mass conservation equation yields the following equations for two-phase flow of oil and water in the reservoir: qo ∂ φ So −1 , (12) = Bo ∇ · (K λ0 ∇o ) − ρo ∂t Bo Bw−1 ∇

qw ∂ · (K λw ∇w ) − = ρw ∂t

φ Sw Bw

,

(13)

where, w = Po − Pc − ρw gh and o = Po − ρo gh, with Pc = Po − Pw , Po and Pw being the pressures in the oil and water phases, λo and λw are, respectively, the mobility of the oil and water, (Bo , qo ) and (Bw , qw ) their corresponding formation volume factor and flow rate, respectively, and h is a block’s (center) depth from the sea level. We then combine Eqs. 12 and 13 to obtain a single equation for Po , containing no explicit time derivatives of the phase saturations: λo λw qo qw + Bw ∇ · K ∇ Po − ρ0 g − (∇ Po − ∇ Pc ) − ρw g − Bo ∇ · K Bo ρo Bw ρw ∂ Po = φC , (14) ∂t where Ci = −Bi−1 ∂ Bi /∂ Pi (i = o, w) represents the phase compressibility, Cr = −φ −1 ∂φ/ ∂ Po is the rock’s compressibility, C = Cr + Co So + Cw Sw , and the terms ρo g and ρw g, representing the gravitational terms, operate only in the vertical direction. The governing equations for Po and So are solved by a finite-volume (FV) method, which is most suitable for the type of computational grid that we generate by the upscaling method, which is one with unequal blocks. To utilize the FV method, we first integrate both sides of Eq. 14 over a control volume in 3D to obtain, ∂ Po λo qo Cφ dV = Bo ∇ · K ∇ Po − ρo g − ∂t Bo ρo qw λw (∇ Po − ∇ Pc ) − ρw g − d V. (15) +Bw ∇ · K Bw ρw The left side of Eq. 15 is evaluated by assuming that the oil pressure in any block is represented by an average value. Likewise, the flow rates qo and qw are assumed to represent the average rates in any gridblock. The time derivative is then approximated by a suitable finite-difference (FD) form, such as, CφV (Pin+1 − Pin )/t, if we use a forward FD approximation, where V is the volume of the block, Pin is the pressure after n time steps at grid point i representing a block, and t is the time step. One then invokes the divergence theorem to convert the volume integrals on the right side of Eq. 15 to surface integrals. Since the upscaled or coarsened grid contains blocks of various sizes, all the terms on the right side of Eq. 15, after their conversion to surface integrals via the divergence theorem, must be evaluated carefully. Consider, for example, the configuration of the blocks shown in Fig. 3 (the third dimension is not shown). For the western side of the base block b we write the

123

Upscaling and Simulation of Waterflooding in Heterogeneous Reservoirs Fig. 3 A block b and its “western” neighbors j in the finite-volume formulation

321

lj

j

lb

b

lb

j

lj lb

j

lj

surface integral as, K w

N λo,up m P j − Po

j + b λo Sbj . ∇ Po d S = Bo

j /K j + b /K b Bo wn j + b

(16)

j=1

Here, Sbj is the common surface between two adjacent blocks j and b, j and b are, respectively, the cell-to-face distances for the neighbor cell j and the base cell b, N ≥ 1 is the number of neighboring blocks on the western side of the base cell b (N = 1 if the two blocks have the same size), (P j , K j ) and (Pb , K b ) are the average pressures and permeabilities representing blocks j and b, and λo,up is the mobility of the oil phase in the front block. The superscript m is the time step number at which the quantities are evaluated. The discretization represented by Eq. 16 results in a local leading truncation error on the order of O(1/ h), when applied to h-adaptive non-uniform grids (Quandalle and Besset 1985; Forsyth and Sammon 1985). The analysis of Edwards (1996) indicated that in grids with large aspect ratios the truncation error depends on the grid interface and aspect ratios, the permeabilities, and the pressure gradient acting tangentially on the grid interface (between two unequal blocks). Thus, in order to eliminate such errors (or reduce them to very small values), we followed the approach suggested by Edwards (1996). Consider, for example, the configuration of the blocks shown in Fig. 4 (the third dimension is not shown). To calculate the flux between a and c in the figure, we use the standard expression based on the harmonic averaging of the blocks’ permeabilities: Jac = 2

h x Ka Kc (Pa − Pc ) , h y Ka + Kc

(17)

whereas to compute, for example, the flux between a and b we use, Jab = 4

Ka Kb Kc hx 1 Pb − (Pa + Pc ) , h y 4K a K c + K a K b + K b K c 2

(18)

with a similar expression for Jbc , where h x and h y are the blocks’ linear dimensions at the finer scale. Consequently, Eq. 16 is used only when equal-size blocks are neighbors (N = 1), but for the case in which unequal blocks (with an interface ratio of 2) are neighbors, Eq. 16

123

322


Fig. 4 Three blocks of unequal sizes, used for computing the fluxes between them (see Eqs. 17–19)

is modified to, K w

N λo,up m Sbj Ka Kb Kc λo ∇ Po d S = 4 Bo 4K a K c + K a K b + K b K c Bo wn lb j=1 1 × Pb − (Pa + Pc ) . 2

(19)

As described by Edwards (1996), the extension to grids of larger interface ratios is straightforward. The simulator first solves the governing equations for the pressure Po of the oil phase, and then computes the oil saturation So by solving Eq. 12. To obtain the solution, we use a combination of the IMPES (implicit pressure, explicit saturation) and the fully implicit methods. The limits m = n and m = n + 1 (where n is the current time step number) of Eq. 16 correspond, respectively, to the IMPES and fully implicit procedures. Expressions similar to Eqs. 16 and 19 are also written down for all the neighbors of the base block b, and the resulting terms are all added up to represent the corresponding integral. The IMPES method is conditionally stable and converges to the correct solution if the time step t is selected carefully. We imposed the condition that the saturation and pressure changes in any grid block between two consecutive time steps must remain tightly bounded, by adjusting t such that the saturation and pressure changes between t and t +t in any grid block were, respectively, no more than 0.05 and 100 psi. In addition, an adaptive time-step method was used. The limitations of the IMPES method are that, (a) because a large portion of the reservoir far from the wells may experience very slow changes in Po and So , it is not efficient to use very small t, and (b) the method is not accurate if large variations in the dependent variables occur rapidly. The time-truncation errors are generally larger in implicit simulators. A fully implicit method provides the required stability, but at a considerably higher-computational cost. The largest variations in So and Po occur in the near-well regions, which also control the maximum allowed time step t when the flood front is close to them. Therefore, at the beginning of the injection, t must be small, but when the front is no longer near the injection well, it can be much larger. When the front reaches the vicinity of the production well, t must again be small. Therefore, in a small zone of grid blocks (of size 8 × 8 × 8) around each well we use the fully implicit procedure for discretizing the governing equations, while the IMPES method is used in the remaining part of the grid. For the fully implicit part, we guess the So distribution at time t and solve the discretized equations for Po . The So distribution is then computed based on the newly calculated Po distribution. The procedure is iterated several times until converged solutions are obtained. For the IMPES part, after computing the Po distribution after each time step, we solve Eq.

123


323

12 in which the time-derivative term is discretized explicitly. In both cases the discretized equations are solved by a combination of the Newton–Raphson and biconjugate-gradient methods. Improved approximations for estimating the upstream flow functions are required, if one is to obtain accurate solutions for Po and So . This requires accurate evaluations of the phase mobilities. Harten and co-workers (Harten 1979; Harten and Osher 1987; Yee and Harten 1987) introduced the total variation diminishing (TVD) method for evaluating the phase mobilities, which is free of the instabilities. It is increasingly being used in reservoir simulation (Rubin and Blunt 1991; Wattenbarger et al. 1997), especially when numerical dispersion dominates the physical dispersion (Peddibhotla et al. 1997; Sammon et al. 2001; Thiele and Edwards 2001; Shrivastava et al. 2002; Mallison et al. 2003). The method requires a front detection technique based on calculating successive mobility ratios in the neighboring blocks. Higher-order TVD methods are constructed by introducing a variable mobility at each block face given by, 1 λo,face = λo,up + ηi+1/2 (λo,down − λo,up ), 2

(20)

where ηi+1/2 is called the mobility limiter. Here, λo,face = λi+1/2 , λo,up = λi , and λo,down = λi+1 . Moreover, ri+1/2 =

λo,up − λo,up−1 . λo,down − λo,up

(21)

Here, λo,up−1 denotes the oil mobility of up-side block, attached to the upstream block of the base cell. ri+1/2 varies in such a way that the scheme is, at least spatially, second order. The results obtained by using Eq. 20 are usually free of spurious oscillations. Several expressions for ηi+1/2 have been developed (van Leer 1979; Sweby 1985; Roe 1985). One of the most accurate limiters (van Leer 1979), used in the present work, is given by, 1 ηi+1/2 = max min 2ri+1/2 , (ri+1/2 + 1), 2 , 0 . (22) 2 According to Eq. 22, in addition to the upstream and downstream mobilities, second-neighbor mobilities λo,up−1 enter the formulation as well. Since in a non-uniform grid the base block might have several second-neighbor upstream blocks, the maximum mobility among such blocks is used as the representative of the upside mobility; that is, λo,up−1 = λi−1,max . Such a choice guarantees that the approaching front will be detected by at least one of the neighboring blocks of the base block. Let us point out that, due to its multiscale nature of the wavelets, no assumption needs to be made regarding the significance, or lack thereof, of the capillary pressure. Clearly, since one obtains a grid with various block sizes, Pc can be (and is) neglected in the large ones.

6 Well Representation It is well-known that, except in the case of a central well in a problem involving cylindrical coordinates, it is impractical to represent a well with an internal boundary, because the ratio of well’s radii and the desired grid block size can be of the order of 10−3 or smaller. One procedure, used by almost all the numerical simulators, is to represent the well by a source and express the relationship between the source strength, the wellbore flow, and the flow in the surrounding grid blocks by a well model. The source/sink representation of a well is a

123

324


local, approximate, steady, singular solution of the flow equations. Separating a singularity of this type for special treatment is an old idea in applied mathematics (Collatz 1960; Carslaw and Jaeger 1963; Fox 1971). We use the ideas described by Crichlow (1977) and Williamson and Chappelear (1981). There are essentially two methods for representing a well in a simulator: by a rate or by a pressure constraint (Fanchi 2001). Employing the concept of productivity index Ip , the oil production rate qo is written as (Peaceman 1977; Aziz and Settari 1979) qo = I p

λo Bo

(Pr − Pwf ) ,

(23)

where Pr and Pwf are, respectively, the reservoir’s and the well’s flowing pressures, with Pr > Pwf . The pressure P at the well grid block is usually used for Pr , so that the variables are qo and P, while the remaining terms are considered parameters. For a rate-constrained well, qo is specified and the simulator solves the pressure equation for P. If a well is pressureconstrained, Eq. 23 is used in place of qo in the pressure equation. The subsequent well pressure is then substituted into Eq. 23 in order to determine qo . The same procedure is used for the injection wells, except that the productivity index Ip is replaced by the injectivity index in Eq. 23.

7 Results and Discussion All the computations were carried out on a Pentium-4 computer at 3.2 GHz. We examined three upscaling scenarios, but at the end of the discussion we identify what we believe to be the most accurate and efficient scenario. (1) The grid was uniformly upscaled in both the vertical direction and the horizontal planes. We refer to this scheme as the UVUH scenario. (2) The grid was coarsened uniformly in the vertical direction, but non-uniformly in the horizontal planes. We refer to this scheme as the UVNUH scenario. (3) The grid was upscaled non-uniformly in both the vertical direction and the horizontal planes. We call this scheme the NUVNUH scenario. The five participants in the original SPE-10 project reported the results of their simulations with the fine grid of the GM model described above and, as expected, they were all similar. In particular, Landmark carried out the full fine-grid simulations using the VIP simulator on a parallel computer. Thus, we did not simulate the waterflooding problem in the fine grid of the GM of the SPE-10 model, but in what follows we compare the results that we obtained with each upscaling scenario to those obtained with the fine-grid model reported by Landmark. We note that the original article describing the project (Christie and Blunt 2001) also used the fine-grid solutions reported by Landmark as the basis for comparison between different upscaling strategies. In what follows, we first provide some details of how the GM was upscaled, and then present the results obtained according to each upscaling scenario. The results that are presented and compared with Landmark’s solution are for production wells 1 and 3 in the model. They are the wells (and the only ones) for which the results obtained with the fine grid of the SPE-10 were presented in the original article (Christie and Blunt 2001). The results that we have computed for the other two wells are given elsewhere (Rasaei 2005).

123


325

7.1 The UVUH Scenario This scheme represents the simplest upscaling scenario that we considered. We first carried out a uniform uplayering procedure to reduce the number of layers. As for the grid in the horizontal planes, as pointed out earlier, we coarsened it using the WT technique. Uniform coarsening of the grid in such planes corresponds to setting both the wavelet scale and detail thresholds to be, s = d = 1. Four levels of wavelet coarsening in the horizontal planes were implemented for all the upscaling scenarios, including UVUH scenario that we describe in this section. Thus, the largest block in the upscaled grid was 16 times larger than the blocks in the fine-grid model. As described above, only in the zones around the wells the structure of the fine grid for the GM was kept intact. In this way, we generated four upscaled models by uniformly uplayering the layers in groups of 2, 3, 4, and 5 layers of the GM. In the fifth model, no uplayering was carried out (all the original 85 layers were kept intact). Thus, the five upscaled models contained 85 (no uplayering), 43, 29, 22, and 17 layers, with the corresponding total number of blocks being 40,035, 19,782, 13,188, 9,981, and 8,007. Note that, for example, the number of blocks in the most coarsened model—8,007—is smaller than that of the original fine-grid model (1.22 × 106 ) by a factor of 140. The total computation times for generating the upscaled grids and their upscaled properties, as well as carrying out simulation of the waterflood problem, were 15, 2.6, 1.4, 1.1, and 0.82 CPU hours. The speed-up factors of the models will be described in Sect. 8. Since the coarsening was done uniformly, both in the horizontal planes and the vertical direction, we did not expect the results to be highly accurate. Despite this, the scenario did yield predictions that, in some cases, were in good agreement with those obtained with the fine grid for the SPE-10 model. Figure 5 presents the reservoir’s total oil production, computed using the five upscaled grids. The maximum difference between the predictions of all the five upscaled models and those of the original SPE-10 presented by Landmark is 8%. Figure 6 compares the results for the oil production from well 3 and, once again, the agreement is very good, implying that even the most coarsened grid with only 8,007 blocks yields accurate results in such instances. The difference between the results obtained with this scenario and those obtained with the UVNUH and NUVNUH schemes (see below) is manifested in the computed breakthrough times and watercuts at wells 1 and 3. Figures 7 and 8 present, respectively, the computed watercuts at these wells, using the five upscaled models, and compare them with those calculated by Landmark using the original fine grid for the reservoir’s GM. The most accurate breakthrough times and watercuts for well 1 are obtained with the two densest grids with 85 and 43 layers (containing, respectively, 40,035 and 19,782 blocks), which are also in excellent agreement with Landmark’s results for the fine grid. The same is true for well 3, except that the next two coarser grids also yield relatively accurate results. Only the most coarsened grid does not yield results that are in quantitative agreement with those of Landmark. As we show shortly, the UVNUH and NUVNUH scenarios yield equally, or more, accurate results (see below), but with much coarser grids and, hence, at much lower computation cost. 7.2 The UVNUH Scenario This scenario is similar to the UVUH scheme described above except that, since the horizontal planes are not upscaled uniformly with the wavelet technique, the thresholds s and d are not, in general, equal to 1. As in the UVUH scheme, we carried out four wavelet upscaling levels. For the first level we set, s = d = 1, while for the next three levels we used, s = d = 0.6, selected based on our simulations of various problems in 2D models.

123

326


Fig. 5 Comparison of the reservoir’s total oil production, as computed with the upscaled models in the UVUH scenario with various number of layers, with the fine-grid results of Landmark

Fig. 6 Same as in Fig. 5, but for well 3 of the reservoir

The uplayering in the vertical direction was the same as in the UVUH scheme. Thus, we generated five coarsened grids with 85, 43, 29, 22, and 17 layers, with the corresponding total number of grid blocks being 48,513, 24,603, 17,049, 13,041, and 10,719. Note that, due to non-uniform coarsening of the horizontal planes, the number of square blocks in the upscaled layers is not the same for all the layers. The corresponding total computation times were 20.5, 4.5, 2.3, 1.85, and 1.6 CPU hours. The resulting speed-up factors will be discussed in Sect. 8. Figure 9 presents a comparison between the reservoir’s total oil production using the five upscaled model and the results presented by Landmark using the reservoir’s original

123


327

Fig. 7 Comparison of the watercut produced at well 1 of the reservoir, as computed with the upscaled model in the UVUH scenario with various number of layers, with the fine-grid results of Landmark


fine-grid model. The agreement is excellent. Figure 10 compares the reservoir’s average pressure, computed with the five coarsened model, with the results reported by Landmark. The two densest grids provide very accurate results for the pressure. This was not the case with the UVUH scenario described above, hence indicating the effect of non-uniform wavelet upscaling of the horizontal planes. Oil production results for wells 1 and 3 and their comparison with those of Landmark are presented in Figs. 11 and 12. In the case of well 3 (and similar to the UVUH scheme), all the results are in excellent agreement with Landmark’s,

123

328


Fig. 9 Comparison of the reservoir’s total oil production, as computed with the upscaled models in the UVNUH scenario with various number of layers, with the fine-grid results of Landmark

Fig. 10 Comparison of the average reservoir pressure, as computed with the upscaled models in the UVNUH scenario with various number of layers, with the fine-grid results of Landmark

while for well 1 the two densest grids provide accurate predictions. Let us mention that the greatest difference between the results of the participants in the SPE-10 project was also for those associated with well 1 (Christie and Blunt 2001). The computed watercuts for wells 1 and 3 are compared in Figs. 13 and 14 with the results of Landmark. Figure 13 should also be compared with Fig. 7 that was obtained with the UVUH scenario. Clearly, the UVNUH scheme improves the predictions. For well 1, the three desnsest grids provide quantitative predictions for the breakthrough times and the watercuts, while in the case of well 3, only the most coarsened grid with the fewest blocks

123


329


Fig. 12 Same as in Fig. 9, but for well 3

cannot do so. Let us again mention that, even the densest grid that we generate by our method requires only about 21 CPU hours of computation on a single, relatively fast workstation. 7.3 The NUVNUH Scenario The main problem with the UVUH and UVNUH scenarios is that, the vertical layers are upscaled uniformly. In practice, this is not the best strategy, as the layers usually have widely different permeabilities. Therefore, for the NUVNUH scenario we first established some criteria for vertical upscaling based on the geological features of the reservoir, so as to improve the accuracy of the method and its efficiency.

123

330


Fig. 13 Comparison of the watercut produced at well 1 of the reservoir, as computed with the upscaled models in the UVNUH scenario with various number of layers, with the fine-grid results of Landmark

Fig. 14 Same as in Fig. 13, but for well 1

We should emphasize again that in the wavelet method, consistent with the theory of WTs (Daubechies 1988), it is in the wavelet space that the significance of the blocks’ permeabilities to the upscaling procedure is best manifested. Therefore, as a preliminary step, we first apply a one-level WT to each layer to obtain its maximum wavelet scale and detail coefficients. Layers with the largest detail coefficients (representing the contrasts between neighboring permeabilities) have the strongest variation in their permeability distribution, while those with the maximum wavelet scale coefficients have the largest (absolute)

123


331

Fig. 15 Comparison between the reservoir’s total oil production, as computed with the upscaled models in the NUVNUH scenario with various number of layers, and the fine-grid results of Landmark

permeability values. In this way, we identify layers with high-lateral permeability changes and can treat them individually in developing the upscaled models. Horizontal upscaling of each layer is carried out similar to the UVNUH scenario described above with the same wavelet thresholds, d = s = 0.6 Hence, five upscaled models were generated: (1) Layers with the maximum effective permeabilities (K x , K y , K z ) do not commingle with other layers. This resulted in an upscaled model with 22 layers. (2) Layers with the maximum and minimum effective permeabilities (K x , K y , K z ) do not combine with other layers. This yielded a model with 24 layers. (3) Layers with maximum wavelet detail coefficients (maximum contrasts between the neighboring blocks) do not mix with other layers. The resulting model had 30 layers. Note the increase in the number of layers when the choice of layers is done based on the wavelet coefficients. (4) Layers with maximum wavelet scale coefficients (maximum permeabilities in the wavelet space) and wavelet detail coefficients (maximum contrasts between the neighboring blocks) do not commingle with other layers. This generated an upscaled model with 36 layers. (5) Uplayering in the vertical direction was carried out based on the well index of the injection well in each layer. To make the water injection rates comparable in all the layers in the upscaled model, layers with high-injectivity indices were considered separately and those with low-injectivity indices were combined together to improve the combined layer injectivity. This resulted in an upscaled model with 36 layers. The total numbers of blocks in the five upscaled models were, respectively, 13,806, 14,907, 18,522, 21,711, and 21,303. The corresponding total computation times for the five models were, 1.9, 2.5, 2.7, 4.0, and 3.7 CPU hours, which are very small. The corresponding speed-up factors will be discussed Sect. 8. Figure 15 shows a comparison between the reservoir’s total oil production, computed with the five upscaled model, and the results reported by Landmark. The agreement is excellent, with all the five upscaled models producing nearly identical results. Figure 16 compares

123

332


Fig. 16 Comparison of the reservoir’s average pressure, as computed with the upscaled models in the NUVNUH scenario with various number of layers, with the fine-grid simulation of Landmark


the average reservoir pressure that the five upscaled models predict with Landmark’s finegrid simulations. The most accurate predictions are provided by the fourth model described above, which is completely based on the wavelet method described in Sect. 3. Note that the differences between our results, shown in Fig. 15, and those reported by Landmark for the fine grid are smaller than those reported by almost all the participants in the SPE-10 project (see Figs. 29–32 in Christie and Blunt 2001). Figure 17 presents the computed oil production for well 1 and compares them with those reported by Landmark. Once again, the fourth model described above (which is based on the WTs) provides the most accurate results, with the differences between our results and Landmark’s being smaller than what all the participants in the SPE-10 project reported.

123


333


Fig. 19 Comparison of the watercut produced at well 1 of the reservoir, as computed with the upscaled models in the NUVNUH scenario with various number of layers, and the fine-grid results of Landmark

Figure 18 makes the same comparison as in Fig. 17, but for well 3. The agreement between the upscaled results and those reported by Landmark is excellent. The computed watercuts for wells 1 and 3 are shown, respectively, in Figs. 19 and 20. For well 3, all the upscaled models provide predictions that are in excellent agreement with Landmark’s results. As for well 1, the fourth model described above provides, once again, the most accurate results, as its prediction for the breakthrough time is also in excellent agreement with Landmark’s.

123

334



Finally, let us point out that, for all the quantites that we computed using the five upscaled models in the NUVNUH scenario, the accuracy of the predictions, when compared with the fine-grid simulation of Landmark, were at least as good as, and in many cases better than, the results presented by the participants in the SPE-10 project competition.

8 Speed-up Factors of the Computations The SPE-10 model that we upscaled was designed in such a way that one could not carry out simulation of the waterflooding process on the fine-grid model on a single workstation. The computation time would be too long for the simulations to be feasible. Thus, the participants in the project carried out their fine-grid simulations on parallel computers with many processors. Since we used only a Pentium-4 computer for all of our computations, there is no direct comparison between the times that it took to carry out the simulations with the upscaled models and those spent by the participants in the SPE-10 project for their fine-grid simulations. It was also not feasible for us to carry out the fine-scale simulations on the Pentium-4 computer, as we estimated that it would take more than two years to complete the simulations (see below). Thus, to provide the reader with some idea about the speed-up factors in our upscaling and simulations, we offer the followings: (1) As mentioned earlier, the largest computation time for any of the upscaled models was only about 21 CPU hours on a Pentium-4 computer at 3.2 GHz. The computation time for model 4 in the NUVNUH scenario, which we identified as the most accurate and reliable upscaled model, was less than 4 CPU hours. Given the accuracy of the computations, the modest speed of the computer, and the size of the fine grid of the SPE-10 model, such computation times indicate the very high efficiency of the method. (2) We attempted to obtain an estimate for the computation time of the fine-grid model of the SPE-10 model on our Pentium-4 computer, by carrying out simulations of several

123


335

Fig. 21 Correlation between the computation time t of the grids with the number N of blocks that they contain

models that we constructed based on the fine grid of the SPE-10 model, all upscaled the same way but coarsened at different levels, so that they had different number of blocks. We then plotted the CPU times of the upscaled models versus the total number of blocks that they contained, and extrapolated the resulting correlation to the number of blocks in the fine grid of the SPE-10 model which is 1.122 × 106 . To do so, we started with the one-level upscaled model in which the wavelet thresholds d and s were set to 0.6 (see above), in order to upscale the grid in the horizontal planes, and no uplayering was carried out in the vertical direction (i.e., the 85 layers in the fine grid of the SPE-10 model were preserved). We then carried out uniform uplayering first, and applied the second-level wavelet upscaling with the same thresholds to the horizontal planes thereafter, in order to obtain the second upscaled model. Repeating this procedure, we generated five upscaled models with 85, 43, 29, 22, and 17 layers, with their total number of blocks N being, respectively, 68,928, 38,043, 27,150, 21,476, and 16,371. The upscaled models were then used in the simulation of the waterflooding process. Figure 21 presents the resulting correlation between the total computation times of the upscaled models (that is, generating the models, assigning their effective properties, and carrying out simulation of waterflooding) and the total number of grid blocks N that they contain. Based on this correlations, we estimate that, had we utilized our Pentium-4 computer, carrying out simulation of waterflooding in the fine grid of the SPE-10 model would have taken about 863 days (or 20,712 CPU hours). Admittedly, this is only a rough estimate, but is consistent with the main premise of the SPE-10 model that it should not be feasible to carry out the simulations with its fine-grid model on a single processor. Based on this estimate, we can now compute speed-up factors of all the upscaling scenarios that we considered. They vary from three to more than four orders of magnitude, huge savings in computation times, while preserving the accuracy of the results, as we demonstrated above. For example, the speed-up factor for the fourth model in the NUVNUH scenario—the upscaled model that we identified as being the most accurate—is about 5,000, more than three orders magnitude savings in the computation time.

123

336


9 Summary The various wavelet upscaling scenarios that we presented in this article yielded predictions with an accuracy that, when compared with the fine-grid simulations of the SPE-10 model, ranged from good to excellent. The computation times of all the models never exceeded 21 CPU hours on a relatively fast, but modest, computer. At the same time, however, the work described in this article suggests that a grid coarsening technique based on an appropriate non-uniform vertical uplayering and wavelet coarsening of the horizontal planes—model 4 in the NUVNUH scheme described above—generates a computational grid that is a faithful representation of the extremely broad heterogeneities of the SPE-10 model, and results in accurate and highly efficient simulations. The model can certainly be improved. One way is to use a more accurate scheme for uplayering of the layers. The second and most promising way is a scheme based on 3D wavelet transformations, in order not to rely on a separate uplayering technique for the vertical direction. As mentioned earlier, we have already used (Pazhoohesh et al. 2006) such a method for a simpler problem, namely, computations of the effective conductivity of amorphous semiconductors at low temperatures, in which the local, fine-scale conductivities varied over up to seventeen orders of magnitude—a conductivity distribution much broader than the permeability distribution of the SPE-10 model. We are currently developing the 3D wavelet method for reservoir simulation. The results will be reported in a future article. Acknowledgement The work of M.R.R. was supported in part by the NIOC.

References Athichanagorn, S., Horne, R.N., Kikani, J.: Processing and interpretation of long-term data acquired from permanent pressure gages, SPE paper 56419 (1999) Aziz, K., Settari, A.: Petroleum Reservoir Simulation. Applied Sci. Publishing, London (1979) Brewer, K.E., Wheatcraft, S.W.: Including multi-scale information in the characterization of hydraulic conductivity distributions. In: Wavelets in Geophysics, p. 213. Academic Press, San Diego (1994) Carslaw, H.S., Jaeger, J.C.: Operational Methods in Applied Mathematics. Dover, New York (1963) Castellini, A., Edwards, M.G., Durlofsky, L.J.: Flow based modules for grid generation in two and three dimensions. Proceedings of Seventh European Conference on the Mathematics of Oil Recovery (2000) Chappelear, A., Hirasaki, G.J.: A model of oil-water coning for two dimensional areal reservoir simulation, SPE paper 4980 (1976) Christie, M.A., Blunt, M.J.: Tenth SPE comparative solution project: A comparison of upscaling techniques, SPE paper 66599 (2001) Chu, L., Schatzinger, R.A., Tham, M.K.: Application of wavelet analysis to upscaling of rock properties, SPE paper 36517 (1996) Collatz, L.: The Numerical Treatment of Differential Equations. Springer-Verlag, Berlin (1960) Crichlow, H.B.: Modern Reservoir Engineering – A Simulation Approach. Prentice-Hall, Englewood Cliffs (1977) Daubechies, I.: Orthonormal basis of compactly supported wavelets. Commun. Pure Appl. Math. 41, 901 (1988) Durlofsky, L.J., Jones, R.C., Milliken, W.J.: A non-uniform coarsening approach for the scale up of displacement processes in heterogeneous porous media. Adv. Water Resour (1997) Durlosfky, L.J., Milliken, W.J., Bemath, A.: Scaleup in the near-well region, SPE paper 61855 (2000) Ebrahimi, F.: Ph.D. Dissertation, Ferdowsi University of Mashhad, Mashhad, Iran (2002) Ebrahimi, F., Sahimi, M.: Multiresolution wavelet coarsening and analysis of transport in heterogeneous porous media. Physica A 316, 160 (2002) Ebrahimi, F., Sahimi, M.: Multiresolution wavelet scale up of unstable miscible displacements in flow through porous media. Trans. Porous Med. 57, 75 (2004) Ebrahimi, F., Sahimi, M.: Grid coarsening, simulation of transport processes in, and scale-up of heterogeneous media: Application of multiresolution wavelet transformations. Mech. Mater. 38, 772 (2006)

123


337

Edwards, M.G.: Elimination of adaptive grid interface errors in the discrete cell centered pressure equation. J. Comput. Phys. 126, 356 (1996) Fanchi, J.R.: Principles of Applied Reservoir Simulation. 2nd ed. Gulf Publishing, Boston (2001) Forsyth, P.A., Sammon, P.H.: Local mesh refinement and modeling of faults and pinchouts, SPE paper 13534 (1985) Fox, L.: Some experiments with singularities in linear elliptic partial differential equations. Proc. R. Soc. Lond. A 323, 179 (1971) Garcia, M.H., Journel, A.G., Aziz, K.: An automatic grid generation and adjustment method for modeling reservoir heterogeneity, SPE paper 21471 (1992) Harten, A.: On a class of highly accurate convective modeling procedure based on quadratic upstream interpolation. Comput. Methods Appl. Mech. Eng. 19, 59 (1979) Harten, A., Osher, S.: Uniformly high-order accurate non-oscillatory schemes. SIAM. J. Numer. Anal. xx, 279 (1987) Heidarinasab, A., Dabir, B., Sahimi, M.: Multiresolution wavelet-based simulation of transport and photochemical reactions in the atmosphere. Atmos. Envir. 38, 6381 (2004) Henriette, A., Jacquin, C.G., Adler, P.M.: The effective permeability of heterogeneous porous media. PhysicoChem. Hydro. 11, 63 (1989) Kikani, J., He, M.: Multiresolution analysis of long-term pressure transient data using wavelet methods, SPE paper 48996 (1998) King, P.R., Snyder, D.E., Obut, T.S., Perkins, R.L.: A case study of the full-field simulation of a reservoir containing bottomwater, SPE paper 21203 (1991) Li, D., Beckner, B.: Optimal uplayering for scaleup of multimillion-cell geologic models, SPE paper 62927 (2000) Li, D., Beckner, B., Kumar, A.: A new efficient averaging for scaleup of multimillion-cell geologic models, SPE paper 72599 (2001) Lu, P., Horne, R.N.: A multiresolution approach to reservoir parameter estimation using wavelet analysis, SPE paper 62985 (2000) Mallat, S.: A Wavelet Tour of Signal Processing. Academic Press, San Diego (1999) Mallison, B., Gerritsen, M., Jessen, K., Orr, F.M.: High order upwind schemes for two-phase multicomponent flow, SPE paper 79671 (2003) Mehrabi, A.R., Rassamdana, H., Sahimi, M.: Characterization of long-range correlation in complex distributions and profiles. Phys. Rev. E 56, 712 (1997) Mehrabi, A.R., Sahimi, M.: Coarsening of heterogeneous media: application of wavelets. Phys. Rev. Lett. 79, 4385 (1997) Neupauer, R.M., Powell, K.L.: A fully-anisotropic Morlet wavelet to identify dominant orientations in a porous medium. Comput. Geosci. 31, 465 (2005) Neivergelt, Y.: Wavelets Made Easy. Birkhäser, Boston (1999) Panda, M.N., Mosher, C., Chopra, A.K.: Application of wavelet transforms to reservoir data analysis and scaling, SPE paper 3656 (1996) Pazhoohesh, E., Hamzehpour, H., Sahimi, M.: Numerical simulation of ac conduction in three-dimensional heterogeneous materials. Phys. Rev. B 73, 174206 (2006) Peaceman, D.W.: Fundamentals of Numerical Reservoir Simulations. Elsevier, New York (1977) Peddibhotla, S., Gupta, A.D., Xue, G.: Multiphase streamline modeling in three-dimensions: further generalizations and a field application, SPE paper 38003 (1997) Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P.: Numerical Recipes. 2nd ed. Cambridge University Press, London (1992) Quandalle, P., Besset, P.: Reduction of grid effects due to local subgriding in simulation using a composite grid, SPE paper 13527 (1985) Rasaei, M.R.: Ph.D. Dissertation, Institute of Petroleum Engineering. University of Tehran, Iran (2005) Rasaei, M.R., Sahimi, M.: Upscaling of two-dimensional models of heterogeneous reservoirs by multiresolution wavelet transformations and its use in the simulation of multiphase flows. Comput. Geosci, submitted (2007) Renard, P., de Marsily, G.: Calculating equivalent permeability: a review. Adv. Water Resour. 20, 253 (1997) Roe, P.L.: Some contribution to the modeling of discontinuous flows. Lect. Appl. Math. 22, 163 (1985) Rubin, B., Blunt, M.J.: Higher order implicit flux limiting schemes for black oil simulation, SPE paper 21222 (1991) Sahimi, M.: Heterogeneous Materials I, Linear Transport and Optical Properties. Springer, New York (2003) Sahimi, M.: Flow and Transport in Porous Media and Fractured Rock, Second Revised and Expanded Edition, Wiley-VCH, Weinheim (2007)

123

338


Sahimi, M., Hashemi, M.: Wavelet identification of the spatial distribution of fractures. Geophys. Res. Lett. 28, 611 (2001) Sahimi, M., Rasaei, M.R., Ebrahimi, F., Haghighi, M.: Upscaling of unstable displacements and multiphase flows using multiresolution wavelet transformation, SPE paper 93320 (2005) Sammon, P.H., Kurihara, M., Liang, J.: Applying high-resolution numerical schemes in reservoirs described by complex corner-point grids, SPE paper 66344 (2001) Settari, A., Aziz, K.: Treatment of nonlinear terms in the numerical solution of partial differential solutions for multiphase flow in porous media. Int. J. Multiphase Flow 1, 817 (1975) Shrivastava, V.K., Nghiem, L.X., Moore, R.G.: A novel approach for incorporating physical dispersion in miscible displacement, SPE paper 77724 (2002) Soliman, M.Y., Ansah, J., Stephenson, S., Mandal, B.: Application of wavelet transform to analysis of pressure transient data, SPE paper 71571 (2001) Sweby, P.K.: High resolution TVD schemes using flux limiters. Lect. Appl. Math. 22, 289 (1985) Thiele, M., Edwards, M.G.: Physically based higher order Godunov schemes for compositional simulation, SPE paper 66403 (2001) Thomas, C.W.: Principles of Hydrocarbon Reservoir Simulation. International Human Resources Development Corporation, Boston (1982) van Leer, B.: Towards the ultimate conservative difference scheme, V. A second order sequal to Godunov’s method. J. Comp. Phys. 32, 101 (1979) Verma, S., Aziz, K.: Two and three dimensional flexible grids for reservoir simulation, Proceedings of the Fifth European Conference on the Mathematics of Oil Recovery, Leoben, Austria (1996) Wattenbarger, R.C., Aziz, K., Orr, F.M.: High-throughput TVD-based simulation of tracer flow, SPE paper 29097 (1997) Wen, X.-H., Gomez-Hernandez, J.J.: Upscaling hydraulic conductivities in heterogeneous media: an overview. J. Hydro. 183, 1996 (1996) Williamson, A.S., Chappelear, J.E.: Representing wells in numerical reservoir simulation: part 1-Theory, SPE paper 7697 (1981) Yee, H.C., Harten, A.: Implicit TVD schemes for hyperbolic conservative laws in curvilinear coordinates. AIAA J. 25, 266 (1987) Younis, R.M., Caers, J.: A method for reservoir scale-up by static-based non-uniform gridding, Stanford Center for Reservoir Forcasting Annual Report (2001) Younis, R.M., Caers, J.: A method for static-base upgridding, Proceedings of the Eighth European Conference on the Mathematics of Oil Recovery, Germany (2002)

123