A numerical simulation of channel reservoirs containing vertical and horizontal wells*

527

Journal of Hydrodynamics Ser.B, 2006,18(5): 527-536

sdlj.chinajournal.net.cn

A NUMERICAL SIMULATION OF CHANNEL RESERVOIRS CONTAINING VERTICAL AND HORIZONTAL WELLS* LIU Fu-ping Beijing Institute of Graphic Communication, Beijing 102600, China，E-mail:[email protected] WANG Xu-song，WANG Jun Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China WANG An-ling Beijing Institute of Graphic Communication, Beijing 102600, China

(Received April. 14, 2005) ABSTRACT: In this article, the bounding surfaces of channels were modeled by Bayesian stochastic simulation, which is a boundary-valued problem with observed valley erosion thickness at the locations of wells (hard data). In this study, it was assumed that the cross-section of the channel shows a parabolic shape, and the case that the vertical well and the horizontal well are located in the channel was considered. Peaceman’s equations were modified to simultaneously solve both the vertical well problem and the horizontal well problem. In porous media, a 3D fluid equation was solved with iteration in the spatial domain, which had channels, vertical wells, and horizontal wells. As an example, the spatial distributions of pressure were calculated for channel reservoirs containing vertical and horizontal wells. KEY WORDS: fluid equations in porous media, porosity, channel reservoirs, numerical simulation, heterogeneous media, hrizontal well, vertical well

1. INTRODUCTION Stochastic models have been developed to describe the different types of sequence stratigraphic bounding surfaces. In particular, the modeling surfaces with incised valleys have drawn considerable attention. Such surfaces are typically complex and their geometries between the wells are generally intrinsically related to significant uncertainties. In fluvial reservoir, the surface geometries are important as they may control the lateral extent of shale barriers or the spatial distribution of reservoir sandstone. In

recent years, several authors [1-3] have considered the estimation of a channel’s location and their size and shape from field-specific observations and general geological information. One of the most general article on this problem is by Hektoen et al. [1]. The channel geometric parameters are modeled as random variables or fields. The authors presented a method to generate realization of all the model parameters. Bi et al. [2, 4] also proposed a method of generating the channel surfaces by stochastic simulation, with simple geometric shapes, i.e., the width of the channel does not vary in the vertical direction. As a horizontal well significantly increases the contact area between the well-bore and stratum, the oil output can be enhanced. This problem has been of great concern to many researchers ever since the first horizontal well has been drilled. Several researches have been carried out, sometimes only in theory and sometimes for a simplified model [1], and at times for a specific model, however, all these methods resolved the problems of only a simple or a special model [5-8]. The above-mentioned methods could not be used to treat the problem and meet the needs of oil production because the horizontal well is usually very long and the distribution of flow in the horizontal well is heterogeneous and various complex strata may be encountered. When numerical simulation was used to resolve the problem with regard to the heterogeneous medium, it has its own advantages, as a result of which no other method can be used as a substitute

* Project supported by the Scientific Research Common Program of Beijing Municipal Commission of Education (Grant No: KM200510015003). Biography: LIU Fu-ping (1960-), Ph.D., Professor

528

for it. In this study, by assuming that the cross-section of the channel is of parabolic shape, Bayesian stochastic simulation is used to simulate the bounding surfaces of the channel, which is constrained by observing valley erosion thickness in wells [9,10] and it is assumed that there are two homogeneous facies in the reservoir, i.e., the floodplain facies (outside the channel) and the channel fill facies (channel interior). If the vertical and horizontal wells are located in the channel, for solving the complex geological problem, the iteration algorithm of a 3-D fluid equation is employed to obtain satisfactory results. It is crucial to resolve the complex geological problem radically using numerical simulation of a 3-D fluid equation to simulate the model with channel reservoirs containing vertical and horizontal wells.

2. FLUID EQUATIONS AND INITIAL BOUNDARY VALUE PROBLEM IN HETEROGENEOUS POROUS MEDIA 2.1 Fluid equations in heterogeneous porous media For the three-dimensional single-phase fluid problem considered here, the gravity effects are neglected, and thus the governing equation can be written as[11,12]

p ( x, y, z , t ) = pb , ( x, y, z ) ∈ Σ

for ( x, y , z ) at Σ and t > 0 ,where pb is the pressure at the boundary. Assume that the initial pressure, p i , is uniform . The initial condition is then specified by

p ( x, y, z , t = 0) = pi

（4）

for ( x, y , z ) in Ω . Within Ω , the inner boundary condition can be written as Nm

Nw

m =1

m =1 k = m1

m2

q = ∑ qm = ∑ ∑ qm ( x, y, z )δ ( x − xi )i

δ ( y − y j )δ ( z − zk )

（5）

where N w is the number of wells, q is the total flow rate, q m is the total flow rate of the well m , and

∂ ⎛ k x ∂p ⎞ ∂ ⎛ k y ∂p ⎞ ∂ ⎟+ i ⎜ ⎟+ ⎜ ∂x ⎝ μ ∂x ⎠ ∂y ⎝ μ ∂y ⎠ ∂z

（3）

δ ( x − xi ) , δ ( y − y j ) and δ ( z − z k ) denote

the Dirac delta functions.

where Ω is the solution domain, p is the pressure,

3. STOCHASTIC MODEL OF THE CHANNEL The model is constructed by simulating a number of surfaces from the base of the reservoir. The surfaces are boundary-valued at observed depths in the wells. Each surface is modeled as a 2-D Gaussian field

t is the time, k x is the permeability in the x-direction, k x = k x ( x, y, z ) and similarly for

S i ( χ ) = ui ( χ ) + ε i ( χ )

⎛ k z ∂p ⎞ ∂p , ( x, y , z ) ∈ Ω ⎜ ⎟ − q = φ ct ∂t ⎝ μ ∂z ⎠

(1)

k y and k z , μ is the fluid viscosity, φ is the porosity of the rock， ct is the total compressibility, and q is the source or sink term at time t . 2.2 Initial boundary value problem For no-flow boundary conditions,

∂p ( x, y , z , t ) = 0, ( x, y , z ) ∈ Σ ∂n

（2）

for ( x, y , z ) at Σ and t > 0 , where Σ is the boundary of Ω . In case of constant pressure boundary, the boundary condition is

（6）

where u i ( χ ) is the expectation of surface i , and ε i ( χ ) is a Gaussian field with zero expectationand a specified covariance function. 3.1 Model for expected surface geometry The expectation surface ui ( χ ) is assumed to consistute a certain number of valleys at a given depth with a 2-D surface and is completely described using i the direction lines l i , the shape parameters ρ , and the corresponding vectors V i ( s ) , which describe the valley geometry along the lines， i = 1,2, …… n ， where n is the number of valleys. The direction line l i can be parameterized by a reference point and two

529

direction vectors. V i ( s ) , defining the position and size of the valley, is a correlated Gaussian function with a one-dimensional reference along l i :

V i ( s ) = [u y ( s), u z ( s ), uw ( s ), uT ( s ), u DP ( s )] (7) where u y ( s ) and u z (s ) are deviations from the projection line of the direction line in the x-y and x-z planes, respectively, u w (s ) is the valley width,

uT (s) is the valley erosion thickness, and u DP (s) is the deviation of the deepest point from the central valley at the position s along l i .The univariate properties are defined by an expectation and a spatial covariance function. Let u y , u z , uW , uT , u DP and σ u y , σ u z , σ uW , σ uT , 2

2

2

2

σ u2 denote the means and Dp

variances of the corresponding normal variables, respectively. The cross-sectional form of a valley is defined as c 0 r β + c1 , where r is measured normal to the line, with r = 0 at the position of the deepest point. Two sets of constants (c 0− , c1− ) and (c 0+ , c1+ ) ，for r < 0 and r > 0, respectively, are computed from the vector V i ( s ) to define an asymmetric valley cross-section. 3 .2 Direction line of the channel The direction line of the channel is a spatial straight-line. It controls the average alignment or the main tendency of the channel. The direction line of the channel is denoted as l p , which is expressed as

l p = l ( x p , y p , z p ) = ( x0 , y0 , z0 ) + t (1, s xy , s xz ) (8) where x 0 is the x -coordinate of the starting point of

l p and will be set to equal to zero, y 0 and z 0 are the y and z coordinates of the starting point, i.e., the y and z coordinates at the point where the line intersects the x = x 0 plane, t is the argument along l p , which will be set to equal to x , and s xy and

assumed that y 0 , z 0 , s xy , and s xz are independent normal random variables and that the estimates of the means and the variances of these variables are available from the geological data or their interpretation. Let y 0 , z 0 , s xy , s xz and

σ y2 , σ z2 , σ s2 , σ s2 denote the means and variances 0

xy

xz

channel are denoted as

η ，distributed

according to

f (η ) , and η and σ are the prior means and the 2 i

prior variances of η , respectively. It is assumed that the parameters of the model are independent normal random variables so that the matrix of the prior variances is a diagonal matrix. Based on these assumptions, the prior distribution has a probability density function satisfying the following relation:

1 f μ [ μ ( χ )] = A1 exp[− (η − η ) T Cη−1 (η − η )] 2 （9） where A1 is a normalizing constant, Cη is the prior variance matrix,

Cη−1 is the inverse matrix of Cη ,

and (η − η ) represents the transpose of (η − η ) . 3 .4 Boundary valued by hard data In this study, hard data, whose measurement errors are modeled as independent Gaussian random variables with zero mean and prescribed variances, refer to the measurement of the valley erosion thickness at the location of wells. The generated stochastic model of the channel must satisfy the hard data at the locations of the wells. Let d h ,ob represent T

the observed data (hard data) for the valley erosion thickness, and d h represent the hard data function, which is related to the model by a linear operator Gh such that

d h = Ghη

s xz represent the l p projection of the slopes onto the x-y and the x-z planes, respectively. It is noted that this is just a parametric representation of a line ( x0 , y 0 , z 0 ) with the direction containing

0

of the corresponding normal variables, respectively. 3.3 Prior model To simplify the notation in the following sections, i the parameters l i , ρ , and V i ( s ) describing a

（10）

Assume the dimension of d h and d h ,ob as N h , let

d h ,i and d h ,ob ,i , respectively, represent the

components of d h and d h ,ob , i = 1, 2,

, N h such

numbers (1, s xy , s xz ) . It should also be noted

that Gh is an N h × N P matrix, where N p is the

that y 0 , z 0 , s xy , and s xz are random variables. It is

number of total model parameters to be estimated. Moreover, let d h ,ob ,i represent the measured value of

530

η ri ， 1 ≤ i ≤ N h ， 1 ≤ ri ≤ N p , Eq.(10) can then be written as

⎡ d h ,ob ,1 ⎤ ⎡ g h ,1,1 ⎢d ⎥ ⎢g ⎢ h ,ob ,2 ⎥ = ⎢ h ,2,1 ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢⎣ d h ,ob , Nh ⎥⎦ ⎢⎣ g h, Nh ,1

g h ,1,2 g h ,2,2 g h, Nh ,2

g h ,1, N p ⎤ ⎥ g h ,2, N p ⎥ ⎥i ⎥ g h, Nh , N p ⎥ ⎦

⎡ η1 ⎤ ⎢η ⎥ ⎢ 2⎥ ⎢ ⎥ ⎢ ⎥ η ⎣⎢ N p ⎦⎥

(11)

and hard data to be significant is obtained by maximizing f D (η d h,ob ) .

4. FINITE DIFFERENCE EQUATIONS 4.1 Finite difference equations Following are the difference equations for the initial-boundary-value problem specified by Eqs. (1) to (4). A purely implicit seven-point difference is used. Suppose the spatial meshes are N x × N y × N z . In the region there are N x grid-blocks in the x-direction, N y grid-blocks in the y-direction, and N z grid-blocks in the z-direction with the

grid-block

i = 1, 2, where,

i = 1, 2

j = 1, 2,

Nx ,

i = 0,1,

Let xi +1 2 for

⎧1, j = ri g h ,i , j = ⎨ , ⎩0, j ≠ ri N h , j = 1, 2,

centers

boundaries in

at

( xi , y j , z k ) for

N y , and k = 1, 2,

Nz .

N x denote the grid-block

the x-direction with x1 2 = 0 and

x N X +1 2 = L x , where L x is the dimension in the N p ,1 ≤ ri ≤ N p

Similar to Eq. (9), the hard data satisfy the following probability density function:

C (Ghη − d h,ob )]

definition,

,

Δt n = t n − t n −1 n i , j ,k

,

and

denote the pressure

obtained at the grid-block center ( x i , y j , z k ) by （12）

where C h is the diagonal covariance matrix of the hard-data measurement errors, i.e., the errors in the valley erosion thickness measurement at the locations of the wells, and A2 is a normalizing constant. For this model, according to Bayes’s theorem, the a posteriori probability density function is given as

solving the system of finite difference equations at time t n . With the preceding notations, the implicit finite difference equation at any grid-block can be written as

Txi −1/ 2, j ,k pin−1, j ,k + Tyi , j −1/ 2,k pin, j −1,k + Tzi , j ,k −1/ 2 pin, j ,k −1 − Ti , j ,k pin, j ,k + Txi +1/ 2, j ,k pin+1, j ,k + Tyi , j +1/ 2,k pin, j +1.k +

f D (η d h,ob ) = fη (d h,ob η ) f μ [ μ ( χ )] =

Tzi , j ,k +1/ 2 pin, j ,k +1 − qin, j ,k =

1 A exp{− [(η − η )T Cη−1 (η − η ) + 2

(Ghη − d h ,ob )T C h−1 (Ghη − d h ,ob )]}

t0 = 0

t n = t n −1 + Δt n . Let p

1 fη (d h ,ob η ) = A2 exp[− (Ghη − d h ,ob )T i 2 −1 h

x-direction. Similarly, the grid-block boundaries in the y- and z-directions are defined. In finite difference techniques, the discrete times, t n , where n = 0,1, 2 , are determined by the

−Vi ,nj ,k pin, −j 1,k （13）

where A is a normalizing constant. The most probable model (the maximum a posteriori estimate) that considers prior information

where

Txi +1/ 2, j ,k =

k x ( xi +1/ 2 , y j , zk )Δy j Δzk

μ ( xi +1 − xi )

（14）

,

531

Txi −1/ 2, j ,k =

k x ( xi −1/ 2 , y j , zk )Δy j Δzk

（15a）

μ ( xi − xi −1 )

Tx1 2, j ,k = TxN x +1 2, j , k = 0

Tyi , j +1/ 2,k =

Tyi , j −1/ 2, k =

N x − 1 and all j and k ,

k y ( xi , y j +1/ 2 , zk )Δxi Δzk

k x ,i −1/ 2, j ,k =

k y ( xi , y j −1/ 2 , zk )Δxi Δzk

k y ,i , j +1/ 2,k =

μ ( y j +1 − y j )

μ ( y j − y j −1 )

(16a)

Tyi ,1 2,k = Txi , N y +1 2, j ,k = 0 for j = 1,2,

Tzi , j , k +1/ 2 =

Tzi , j ,k −1/ 2 =

(16b)

k z ( xi , y j , zk +1/ 2 )Δxi Δy j

μ ( zk +1 − zk )

k z ,i , j ,k +1/ 2 =

,

k z ( xi , y j , zk −1/ 2 )Δxi Δy j

(17a)

μ ( zk − zk −1 )

k z ,i , j ,k −1/ 2 =

2 ( xi +1 − xi ) k x ,i +1, j ,k k x ,i , j ,k k x ,i , j ,k Δxi +1 + k x ,i +1, j ,k Δxi

,

2 ( xi − xi −1 ) k x ,i , j ,k k x ,i −1, j ,k

（20）

k x ,i −1, j ,k Δxi + k x ,i , j ,k Δxi −1

2 ( y j +1 − y j ) k y ,i , j +1,k k y ,i , j ,k

,

k y ,i , j ,k Δy j +1 + k y ,i , j +1,k Δy j 2 ( y j − y j −1 ) k y ,i , j ,k k y ,i , j −1,k

（21）

k y ,i , j −1,k Δy j + k y ,i , j ,k Δy j −1 2 ( zk +1 − zk ) k z ,i , j ,k +1k z ,i , j ,k k z ,i , j ,k Δzk +1 + k z ,i , j ,k +1Δzk

2 ( zk − zk −1 ) k z ,i , j ,k k z ,i , j ,k −1 k z ,i , j ,k −1Δzk + k z ,i , j ,k Δzk −1

,

（22）

n

In Eq. (14), qi , j ,k represents the internal sink or

(17b)

N z − 1 and all i and j . Here, Δxi ,

Δy j , and Δz k are the dimensions of the grid-block, and its center is at ( x i , y j , z k ) . In Eq. (14),

Vi ,nj ,k =

k y ,i , j −1/ 2,k =

N y − 1 and all i and k ,

Tz i , j ,1 2 = Tz i , j , N z +1 2 = 0 for k = 1,2,

grid-block interfaces are computed using the standard harmonic average as shown below:

k x ,i +1/ 2, j ,k =

(15b) for i = 1, 2,

k y ,i , j ,k = k y ( xi , y j , z k ) . The permeabilities at the

φ ( xi , y j , z k )ct Δxi Δy j Δz k

（18）

Δt n

the source term at t n in the grid-block (i, j , k ) and is nonzero only if the grid-block is penetrated by a well. n When a grid-block is penetrated by a well, qi , j ,k represents the total sand-face flow rate into or out of the well over the interval z k −1 2 ≤ z ≤ z k +1 2 . 4.2 Peaceman’s equations 4.2.1 Peaceman’s equations Consider a well located at a position, with the coordinates ( x, y ) given by ( x i , y j ) , and denote the n

Ti , j ,k = Txi +1/ 2, j ,k + Txi −1/ 2, j ,k + Tyi , j +1/ 2,k +

total sand-face flow rate at time t n by qi , j , with the n

individual source or sink terms given by qi , j ,k for

Tyi , j −1/ 2,k + Tzi , j ,k +1/ 2 + Tzi , j ,k −1/ 2 + Vi ,nj ,k

k = m1 , m1 + 1, (19)

for all (i, j , k ) , i = 1, 2

k = 1, 2

, Nx ,

j = 1, 2

, N y ,and

, N z . Throughout the article, the following

notations are used: k x ,i , j , k = k x ( x i , y j , z k ) and

source

or

the

, m2 , i.e., qin, j ,k provides the sink

term

for

the

grid-block

( xi , y j , z k ) . The individual source or sink terms are related to the well-bore pressure by applying Peaceman’s equations at each location, i.e.,

532

n qin, j ,k = WI i , j ,k ( pin, j ,k − pwf ,i , j )

where p

n wf ,i , j

WI i , j ,k =

（23）

0.28073Δy j roi =

is the wellbore pressure, and

2πΔzk k x ,i , j , k k y ,i , j ,k

k x , i , j , k ⎛ Δy j ⎞ ⎜ ⎟ k y ,i , j , k ⎝ Δxi ⎠ 1 + k x ,i , j , k k y ,i , j , k

q

2

n j ,k

m2

∑ WI ( p

k = m1

∑q

k = m1

n i , j ,k

i, j ,k

n i , j ,k

= q ( y j , zk ) = ∑ qin, j ,k = i = i1

for i = i1 , i1 + 1,

n − pwf ,i , j )

n − pwf , j ,k )

（26）

, i2 .

5.

MODIFICATIONS OF THE FLUID EQUATION AND PEACEMAN’S EQUATIONS IN POROUS MEDIA By substituting Eqs.（23）and（25）in（14）, the following equation is obtained:

Tzi , j ,k −1/ 2 pin, j ,k −1 − (Ti , j ,k + WI i , j ,k + （24）

Wki , j ,k ) pin, j ,k + Txi +1/ 2, j ,k pin+1, j ,k +

4.2.2 Extension of Peaceman’s equations It is known that Peaceman’s equations can be used to solve only the vertical well problems not the horizontal well problems. If the gravity effects are neglected, Peaceman’s equations can be extended not only to solve the vertical well problem but also to solve the horizontal well problem. Peaceman’s equations is expanded as

(

n i , j ,k

Txi −1/ 2, j ,k pin−1, j ,k + Tyi , j −1/ 2,k pin, j −1,k +

=

n qin, j ,k = Wk i , j ,k pin, j ,k − p wf , j ,k

i , j ,k

i =i1

layer k . As the total rate of a well equates the summation of each individual grid-block rate, the following equation is obtained: m2

i2

n

∑ Wk ( p

s k represents an effective skin for the well grid-block centered at ( x i , y j , z k ) , s k is the skin factor for

qin, j = q n ( xi , y j ) =

2

1 + k y ,i , j , k k z ,i , j , k

i2

rw is the well-bore radius of the well,

and

⎞ ⎟⎟ ⎠

The total rate of a horizontal well is given as

μ [ln(rok rw ) + sk ]

0.28073Δxi 1 + rok =

⎛ Δz k 1 + y ,i , j , k ⎜ k k z ,i , j ,k ⎜⎝ Δy j

)

（25）

Tyi , j +1/ 2,k pin, j +1.k + Tzi , j ,k +1/ 2 pin, j ,k +1 + n WI i , j ,k pwfn ,i , j + Wki , j ,k pwf , j ,k =

−Vi ,nj ,k pin, −j 1,k

（27）

Equations（24）and（26）are rewritten as

⎛ m2 ⎞ n n WI p − ⎜ ∑ WI i , j , k ⎟ pwf ,i , j = ∑ i , j ,k i , j ,k k = m1 ⎝ k = m1 ⎠ m2

where

p

n wf , j , k

is the well-bore pressure of the

horizontal well, and

Wki , j ,k =

2πΔxi k y ,i , j ,k k z ,i , j ,k

μ [ln(roi rw ) + sk ]

q n ( xi , y j ) ,

i2

∑Wk i =i1

i , j ,k

n i , j ,k

p

（28）

⎛ i2 ⎞ − ⎜ ∑Wki , j ,k ⎟ pwfn , j ,k = ⎝ i =i1 ⎠ q n ( y j , zk )

（29）

533

Equations (27) to (29) are called the modifications of fluid equation and Peaceman’s equations. If the well rate is specified, then by solving the n equations, the grid-block pressure pi , j , k and the

σ u2 = 0.01 , σ u2 = 0.01 z

y

The width:

well-bore pressure p wf ,i , j and p wf , j ,k can be obtained.

uW = 175m, σ u2W = 0.01

With known grid-block pressure and well-bore pressure, the layer rates can be computed using Eqs. （23）and（25）. For a problem with the well-bore pressures specified, Eq. (14) can be rewritten as

The valley erosion thickness:

n i −1, j , k

Txi −1/ 2, j ,k p

n

+ Tyi , j −1/ 2,k p

n i , j −1, k

n i , j , k −1

Tzi , j ,k −1/ 2 p

+

− (Ti , j ,k + WI i , j ,k +

Tyi , j +1/ 2,k p

+ Tzi , j ,k +1/ 2 p

u DP = 0 and σ u2 = 0.01 Dp

d h ,ob ,1 = 15m, at (x =200m, y =500m)

n i , j , k +1

=

n n −WI i , j ,k pwf ,i , j − Wki , j , k pwf , j , k −

Vi ,nj ,k pin, −j 1,k

The deviation of the deepest point from the central valley at the position s along l i :

The observed valley erosion thickness in the wells:

Wki , j ,k ) pin, j ,k + Txi +1/ 2, j ,k pin+1, j ,k + n i , j +1.k

uT = 27.5m , σ u2T = 0.01

d h ,ob ,2 = 18m, at (x =800m, y =520m) d h ,ob ,3 = 16m, at (x =1600m, y =480m)

（30）

Equation (30) can be used to solve the grid-block pressures individually. Equatins（23）, (24), (25), and (26) can then be applied to compute the flow rates.

6. NUMERICAL EXAMPLES Figures 1, 2, and 3 show a realization of a 3-D single channel on a simulation grid with the direction line of the channel along the x-direction. The dimensions of the grid are 2000 × 1000 × 300 m3. In generating this realization, the following values of the model parameters were used: The direction line:

y0 = 500m, z0 = 162m, sxy = 0.00075,

Figure 1 shows the cross-sectional shape of the bounding surfaces of the channel in the x=700m plane. Figure 2 illustrates the banks of the stochastic channel. Figure 3 shows the contours of the channel depths, and, the black dots in the figure represent the observed well locations. 160

z/m

n

140

120 350

420

490

560

y/m

Fig.1 Cross-section of bounding surfaces of the channel

sxz = 0.0025, σ

2 y0

= 1.0, σ = 1.0, 2 z0

σ s2 = 0.0001, σ s2 = 0.0001 xy

xz

The vertical and horizontal deviation from the direction line:

u y = 0.0 , u z = 0.0 ,

Thus, the modified fluid equation and the modified Peaceman’s equations in porous media have been solved by the iteration. All reservoir boundaries are no-flow boundaries. The spatial domain is 2000 × 1000 × 300m3. The vertical well is located at （x=250m，y=500m）, and the perforated interval

534

μ = 5 ×10−4 Pai s (5cp)

750

y/m

Compressibility: 500

ct = 2.41× 10 −9 Pa −1

250 0

600

1200 x/m

1800

Fig.2 Banks of the stochastic channel

Fig.3 Contour of the channel depths

ranges from z=140 to z=240. The horizontal well is located at（y=500m，z=140m）, and the perforated interval ranges from x=250 to x=750. The other reservoir and fluid properties are as follows: it is assumed that there are two homogeneous facies in the reservoir, the floodplain facies (outside the channel), and the channel-fill facies (channel interior). System permeability:

Initial pressure:

pi = 4.144 × 107 Pa All examples are based on the above parameters. The numerical results of the flow rates that demonstrates the accuracy of this method are obtained as follows: The total flow rate computed using Eqs.(27), (28), and (29) is 86.39 m3/d, whereas the total flow rate specified in the known condition is 86.40 m3/d. By comparing the computed result with the specified result, it is seen that the flow rate obtained using this method is in excellent agreement with the specified flow rate. Hence, the numerical results obtained are reliable. Figure 4(a) shows the pressure distribution in the horizontal plane (z=138m), which incises the horizontal well and the channel and is perpendicular to the vertical well. Figure 4(b) shows the contours of the pressure distribution in the same plane as in Fig.4(a). It is seen from the figure that the pressure inside the channel is lower than that outside the channel. In the figure, the horizontal perforated interval is clearly shown. Because the geometric shape of the channel is asymmetric, the pressure distribution should also be asymmetric in the plane. Indeed, this point can be clearly seen in Fig. 4.

k x = k y = k z = 0.965 ×10−15 m2（0.001darcy）, outside the channel

k x = k y = k z = 1.93 ×10−13 m2（0.2darcy）, inside the channel Porosity:

φ = 0.01 , Outside the channel φ = 0.2, Channel interior The total flow rate:

q = 0.001m3 s −1 （ 86.4 m3/d） Fig.4 Pressure distribution in the horizontal plane z=138m

Fluid viscosity:

535

vertical well and horizontal well and the channel are obviously identified. The vertical well is located at x=250m.

Fig.7 Pressure distribution in the vertical plane y=300m Fig.5 Pressure distribution in the horizontal plane z=102m

Figure 5(a) shows the pressure distribution in the horizontal plane z=102m, which is parallel to the horizontal well and perpendicular to the vertical well. The distance is 36 m from the plane to the horizontal well. The plane does not incise the channel and the horizontal well. Figure 5(b) shows the contours of the pressure distribution in the same plane. The effect of the channel and the horizontal well is shown in Fig. 5. The variation of the pressure distribution in Fig. 5 is less sharp than that in Fig.4.

Figure 7(a) shows the pressure distribution in the vertical plane y=300m, which is parallel to the vertical well and the horizontal well. The distance is 200 m from the plane to the horizontal well. Figure 7(b) shows the contours of the pressure distribution in the same plane. In Fig.7, the perforated intervals of the vertical well and the horizontal well are not identified, however, the asymmetry of the pressure distribution is shown.

Fig.6 Pressure distribution in the vertical plane y=500m

Fig.8 Pressure distribution in the vertical plane x=460m

Figure 6(a) shows the pressure distribution in the vertical plane y=500m, which incises the vertical and horizontal wells and the channel. Figure 6(b) shows the contours of the pressure distribution in the same plane. In the figure, the perforated intervals of the

Figure 8(a) shows the pressure distribution in the vertical plane x=460m, which incises the vertical well and the channel and is perpendicular to the horizontal well and also shows the funneled distribution of pressure, which is attributed to the effect of the

536

horizontal well. Figure 8(b) shows the contours of the pressure distribution in the plane. In Fig .8, the perforated interval of the vertical well and the cross-section of the channel are identified. The vertical well is located at y=500m. The pressure inside the channel is obviously lower than that outside the channel.

solve the problem of channel reservoirs containing vertical and horizontal wells is a good option.

ACKNOWLEDGEMENTS This work was supported by the signal Processing Basic Research Grants of Beijing Institute of Graphic Communication, Recruiting Personnel Research Grants of Beijing Institute of Graphic Communication.

REFERENCES

Fig.9 Pressure distribution in the vertical plane x=1280m

Figure 9(a) shows the pressure distribution in the vertical plane x=1280m, which is parallel to the vertical well and is perpendicular to the horizontal well and also shows the funneled distribution of pressure. Figure 9(b) shows the contours of the pressure distribution on the plane. In Fig. 9, the cross-section of the channel is clearly identified. The pressure change inside the channel is obviously higher than that outside the channel.

7. CONCLUSIONS It is known that the stochastic model has been developed to describe the sequence stratigraphic bounding surfaces. Gausian fields have appeared to be well suited for the modeling of the complex geometry in these surfaces because of the large degree of flexibility associated with the trend surfaces. By the use of Bayesian theorem and by combining the specific hard data, realizations from a posteriori model have been produced. The model can meet the needs of the complex boundary conditions. The vertical and the horizontal wells are incorporated into the complex model. The problem itself with universality and difficulty, whereas this problem cannot be resolved using traditional methods. For the complex geologic problem, the finite difference iteration has been applied to solve the problem. The numerical examples prove that the results are reliable. Hence, it is seen that using numerical simulation of a 3-D fluid equation to

[1] HEKTOEN Anne Lise, HOLDEN Lars. Bayesian modelling of sequence stratigraphic bounding surfaces [C]. Fifth International Geostatistics Congress. Wollongong, Australia, 1996, 22-27. [2] BI Z. X. Conditioning three dimensional stochastic channels to well-test pressure data[D]. Ph. D. Thesis, Tulsa, USA: University of Tulsa, 1999, 6-21. [3] ZHANG F. J., REYNOLDS A. C. and OLIVER D. S. Model errors inherent in conditioning a stochastic channel to pressure data[C]. SPE 62987. 2000, Dallas, USA, 315-325. [4] BI Z. X., OLIVER D. S. and REYNOLDS A. C. Conditioning 3D stochastic channels to pressure data[J]. SPE Journal, 2000, 5(4): 474-484. [5] LI Tang et al. Analysis of horizontal well production formulas[J]. Petroleum Exploration and Development, 1997, 24 (5):76-79. (in Chinese) [6] DIKKEN B. J. Pressure drops in horizontal wells and its effect on their production performance[C]. SPE 19824. Houston, USA,1989, 1-7. [7] ZHAO Chun-sen, CUI Hai-qing, SONG Wen -ling. The productivity of well pattern with horizontal wells and vertical wells[J]. Journal of Hydrodynamics, Ser. B, 2003,15(5):60-64. [8] YIN Hong-jun, HE Ying-fu, FU Chun-quan. Pressure transient analysis of heterogeneous reservoirs with impermeability barrier using perturbation boundary element method[J]. Journal of Hydrodynamics, Ser. B, 2005,17(1):102-109. [9] LIU Fu-ping, ZHANG Yu, YANG Chang-chun. Adaptive nonuniform grid integration equation upscaling method of channel reservoir[J]. Journal of Hydrodynamics, Ser.A, 2003,18(6):775-782. (in Chinese) [10] LIU Fu-ping, KONG Fan-qun, LIU Li-feng et al. Adaptive nonuniform grid upscaling method of channel reservoir[J].Chinese Journal of Computational Mechanics, 2003,20(4): 456-461. (in Chinese) [11] KONG Xiang-yan. Advanced mechanics of fluids in porous media [M]. Hefei ： Press of University of Science and Technology of China, 1999, 18-57. (in Chinese) [12] HE N. Q. Three dimensional reservoir by inverse problem theory using well-test pressure and geostatistical data[D]. Ph. D. Thesis, Tulsa, USA: University of Tulsa, 1997, 63-93.