advanced horizontal well model

SPE 102801

ADVANCED HORIZONTAL WELL MODEL Imaging Seismic Deformation Induced by Hydraulic Fracture Complexity

Alexey N. Cheremisin, Alexander A. Chesnokov, Sergey V. Golovin, Dmitry S. Kuznetsov S.C. Maxwell, C.K. Waltman, N.R. Warpinski, M.J. Mayerhofer, and N. Boroumand, Pinnacle Technologies. single-phase numerical simulator. Those two values are used injection and awellbore pre-existing fractureGeneral network2-6solutions . Even between to calculate the effective radius. are neighbouring wells, the geometry of the stimulated fracture developed for horizontal wells with multiple 2D and 3D fracIntroduction network shows a high degree of variability due to localized This paper was selected for presentation by an SPE Program Committee following review of The reservoir is homogenous with constant tures. porosity differences in the fracture network2. The Barnett has veryand low information contained in an abstract submitted by the author(s). Contents of the paper, as presented, have not been reviewed by the Society of Petroleum Engineers and are subject to intrinsic matrix and the permeability constant permeabilities k k has impermeThe standard analytical method of solving the diffusion equax , ky ,permeability, z . The reservoir correction by the author(s). The material, as presented, does not necessarily reflect any associated with thereservoir fracture stimulation results in of the Society of Petroleum Engineers, officers, or members. Papers able presentedtop at and enhancement bottom; infinite acting or box-shaped tion is Laplace transformationposition with respect to time withitsfurther SPE meetings are subject to publication review by Editorial Committees of the Society of permeable fracture networks sufficient for economic gas Petroleum Engineers. Electronic reproduction, distribution, or storage of any part reservoir of this paper is considered. Multiple fractures with infinite or fiFourier transformation with respect to spatial variables, then, for commercial purposes without the written consent of the Society of Petroleum Engineers is recovery in the field. Previous studies have shown a Permission to reproduce in print is restricted to an abstract of notnite more than conductivity can be arbitrary oriented to the wellbore. according to particular shapeprohibited. of reservoir and boundary concorrelation between the volume of the reservoir that is 300 words; illustrations may not be copied. The abstract must contain conspicuous of where and by whom the paper was presented. Write Librarian, SPE, P.O. The inner boundary wellbore) total flow ditions, the image of solutionacknowledgment is obtained. Inverse Fourier and stimulated,condition as measured(at by the the volume of theis reservoir that Box 833836, Richardson, TX 75083-3836, U.S.A., fax 01-972-952-9435. emits microseisms during the stimulation, and the production rate q(t). Laplace transformations give the formula for solution. ultimately realized from well2,4,6. element The correlation is Abstract Fully analytical model based onthe boundary method The described procedure has many restrictions in use. The attributed to larger fracture networks being stimulated in wells Microseismic mapping is extensively utilized in the Barnett combined with Laplace are active used volume to perform a fracdomain (reservoir) must be of special shape (ellipse, rectanguwhere a largetransform microseismically of the reservoir Shale, to map hydraulic fracture complexity associated with ture flow model reservoirresulting flow model [3]. permeable Based onfracture solular), the external and internalinteractions boundaryofconditions are usually has and beena realized, in more the stimulation with pre-existing fractures. pathways connected to the well and hence a higher potential tion built for vertical wells in a closed reservoir, a solution for constant values of pressure orPrevious flow rate. The properties of the studies have indicated a fair correlation between the for intersected gas flow to the Recently, many operators in the wells bywell. multiple hydraulic fractures is well performance of theline seismically active horizontal volume. wellbore are simplified, its deviation from and theextent straight is Barnett have attempted horizontal completions, which have However, in addition to this measure of the extentdeveloped. of the The fluid is supposed to be slightly compressible, often omitted. allowed large volumes of the reservoir to be stimulated with stimulated fracture network, the characteristics of thisthe fracture reservoir anisotropy intoMany account. Thecompletions wellbore has Taking into account all restrictions above, numerical methlarge fractureis taken networks. of these use network is also expected to impact the well performance. In infinite conductivity, and a flow rate q(t) of the well is known. perforated, cemented liners and the microseismic images allow ods for solving reservoir problems are in the vanguard now. particular, the fracture spacing is believed to be important for indentification of improved perforation Itstaging to Earlyapproximations are derived. can be Nevertheless, numerical methods are powerless against factor controlling the potential gas high flow. In this paper, we and long-time maximize the stimulated reservoir volume (SRV)4. utilize the density of the total seismic moment release (a used in welltest analysis for evaluation of reservoir properties. gradient solutions, which occur in singular points of the domeasure of or the in microseism strength) as an indicationAn of interesting combination of the numerical fracture model main such as source (or sink)robust for example, the neighborMany of the Barnett stimulations are water fracs, where high the seismic deformation that may correlate to the fracture with model wasat developed [4].possible The hood of faults. So, one problem caused by high gradients is volumes reservoir of water are injected high rate7. One density. The study uses a set of microseismic maps ofan analytical mechanism for the success of waterfracs that increased fluid reservoir is anisotropic (kx 6= ky 6= kz ) andis infinite acting. accuracy and stability of numerical Anotherincluding prob- cases where hydraulic methods. fracture stimulations, the pressure in natural fractures induced shear failure, resulting in stimulated reservoir volume measured by the extent of the The solution is obtained for one well with one fracture, but the lem is time resources and machine memory needed for accurate fracture dilation associated with mismatched surfaces on seismically active region poorly correlated with the well of fractures number can be extended by the method of superpocalculation. opposite sides of the fracture. Within this conceptual performance. Incorporating the seismic moment density to sition. The fracture hasthe varying properties conductivity), Common features of the existing analytical models are dicframework, microseismic events (e.g. correspond to the actual assess the fracture density with the network extent, an fracture movement. earlier investigations of the SRV and its direction might not beThe straight. The boundary conditated by methods of solution. These are: constant reservoir improved correlation with the well performance was observed. measured the totalflow volume the It microseismically tion at the wellbore is known rateofq(t). is shown thatactive the properties with respect to spatial variables (x, y, z), type of region. However, this measure of the stimulated volume does approximationnotoftake complex fracture geometry with pseudo-skin external boundary conditionsIntroduction (constant pressure or no-flow). into account the properties of the fracture network, of hydraulic fracture stimulations givehas non-satisfactory the early stageswell of production. Also, the pressure drop in theMicroseismic wellbore ismapping neglected. which hasmistake also beenatindicated to impact performance6. become a common technique to map the fracture growth and 1-5solution, but it also simThe article [5] compares different methods of estimating Furthermore, the permeability enhancement of the fracture Of course, all of this simplify the geometry . Microseismic images provide details of the related to deformation with fracturing. thefrom horizontalmay wellbe productivity. These associated are well-known Joshi plifies the description of the fluid flow in the reservoir. fracture azimuth, height, length, and complexity resulting Beyond the basic hypocentral of the microseisms with pre-existing The resulting images can its augmented method, version, a new locations analytical solution preIn the review below onlyinteraction models with multiple fratures. transverse used to calculate the SRV, additional seismic signal be orientation used to calibrate numerical simulations of the sented fracture in paper, and a numerical simulation. The constantfractures are considered. The of the fractures to characteristics allow investigation of the source of the growth, allowing more confident modeling of other pressure outermechanical boundarydeformation ellipse is resulting used in the model. The nuthe wellbore is arbitrary. in the microseisms. In stimulations in the field, and a better identification of the 8 particular, the seismic momentwith , a robust measure ofone. the simulation showed a deviation an analytical The inner boundary condition (atregion the that wellbore) is known stimulated may ultimately be drained by themerical well. of an earthquake or microearthquake, be used to The cause of strength this deviation was identified, and ancan augmentaflow rate q(t). This is due to fact that the majority of results quantify the seismic deformation9. Arguably, the Barnett Shale is the field that has had the most tion to Joshi’s solution was done. Also a new more appropriate are supposed to be used in welltest analysis. fracs mapped over the last several years. Microseismic analytical solution to the simulaThe well is usually positioned at the geometry center of the In this provided paper, we results examine much severalcloser published microseismic mapping in the Barnett has repeatedly demonstrated extreme projects in the Barnett for correlation between the production tion.the drainage region (reservoir), and the well is treated as afrom straight fracture complexity resulting interaction between Based on the information above, and business needs, a new line, no deviation. semi-analytical method for 3-D simulation of reservoir flow toA short overview of most relevant existing models is done wards the horizontal well was developed. Our model is a combelow. bination of discrete approach in one direction, and analytical An analytical solution for diffusivity equation was perway of solution in the rest 2 directions plus time. formed [1]. The reservoir is homogenous, external boundary Following requirements for the model were made: conditions are either of Dirichlet or Neumann type. Boundary condition at the wellbore is flow rate q(t). The fluid is • Permeabilities are step functions along the direction of slightly compressible. The conductivity of the wellbore is inthe horizontal wellbore; finite. There can exist multiple transverse fractures along the • Wellbore has finite conductivity; wellbore. The conductivity of fractures is finite, non-Darcy flow is considered inside the fracture. A solution for multiple • Trajectory of the horizontal wellbore is not a straight wells with fractures in an isolated reservoir is obtained. line; Another semi-analytical model of horizontal well with mul• There can exist multiple transverse hydraulic fractures. tiple fractures in a brick-shaped drainage volume was perEach fracture has its own size and conductivity; all fracformed in [2]. 2D Fourier transform is used to obtain anatures can be positioned at any point of the wellbore. lytical solution. At the same time pressure is calculated using Copyright 2009, Society of Petroleum Engineers Inc.

Copyright 2006, Society of Petroleum Engineers

This paper was prepared for presentation at the 2006 SPE Annual Technical Conference and Exhibition held in San Antonio, Texas, U.S.A., 24–27 September 2006.

2

A.Cheremisin, A.Chesnokov, S.Golovin, D.Kuznetsov

The Background The differential equation for pressure distribution in reservoir is a consequence of the Darcy law k w ~ = − ∇p µ

(0.1)

and continuity equation ∂(ρϕ) + div(ρw) ~ = Fˆ . ∂t

(0.2)

System (0.1), (0.2) of differential equations for unknown velocity w ~ and pressure p is closed if functions ϕ = ϕ(p, ~ x),

ρ = ρ(p, ~x),

µ = µ(p, ~ x),

k = k(p, ~ x)

defining porosity ϕ, liquid density ρ, dynamic liquid viscosity µ, and permeability coefficient k as functions of pressure p(t, ~ x) and coordinates ~ x = (x, y, z), and the right-hand side Fˆ (t, x, y, p) of Eq. (0.2) are specified. The following assumptions are valid:

SPE number

For our current purpose the reservoir permeability is a function of special type: kx = kx (x), ky = ky (x), kz = kz (x). The equation (0.5) reduces to kx (x) ∂P ky (x) ∂ 2 P kz (x) ∂ 2 P ∂p ∂ β∗ = + + + ∂t ∂x µ ∂x µ ∂y 2 µ ∂z 2 + q(t, x)δ y − yw (x), z − yw (x) + F. (0.7) The main physical parameters have the following dimensions m kg kg kg [w] ~ = , [p] = , [ρ] = 3 , [µ] = , [k] = m2 , s m · s2 m m·s 1 m · s2 m2 m3 [δ(yw , zw )] = 2 , [β ∗ ] = , [q] = 3 , [Q] = , m kg s s porosity ϕ is non-dimensional. Piezoconductivity equation (0.7) equipped with suitable initial and boundary conditions determines pressure distribution in the reservoir and fluid inflow to the wellbore.

1

Heterogeneous Reservoir Without Fracture

• Dynamic viscosity µ is constant;

In this part the principal algorithm for pressure distribution calculation in the reservoir is described.

• Permeability tensor k is diagonal and depends on (x, y, z);

1.1

• Liquid is almost incompressible:

Ω = {(x, y, z) ∈ R3 | x ∈ (0, X), y ∈ (0, Y ), z ∈ (0, H)}.

ρ(p) = ρ0 (1 + βl (p − p0 ));

(0.3)

• Porous media is almost incompressible: ϕ(p) = ϕ0 + βp (p − p0 ).

(0.4)

Here volume liquid compressibility βl and volume rock compressibility βp are defined as follows βl = −

1 dVl , Vl dp

βp =

dVi V dp

(Vl is the liquid volume, Vi is the pore volume, V is the stratum volume). Relations (0.3) and (0.4) reduce the system (0.1), (0.2) to one differential equation for unknown pressure. To do this, one must represent the term ρϕ as follows: ρϕ = ρ0 ϕ0 + (ρ0 ϕ0 βl + ρ0 βp )(p − p0 )+ + ρ0 βp βl (p − p0 )2 ≈ ρ0 ϕ0 + ρ0 β ∗ (p − p0 ) (quantity ρ0 βp βl (p − p0 )2 is negligible in comparison with the remaining terms). Here β ∗ = ϕ0 βl + βp is the storage coefficient. Next, the following approximate formulae are derived: ∂(ρϕ) ∂p = ρ0 β ∗ , ∂t ∂t ρ0 ρ0 div (1 + βl (p − p0 ))k∇p ≈ − div(k∇p). µ µ Thus, a parabolic equation for pressure is obtained: k ∂p β∗ = div ∇p + F. (0.5) ∂t µ div(ρw) ~ =−

ˆ Here F = ρ−1 0 F . The horizontal wellbore is introduced to the modelby means of the additional term q(t, x)δ y − yw (x), z − yw (x) . This term models a horizontal wellbore with the trajectory W = {(x, yw (x), zw (x)) |0 ≤ x ≤ X}, (0.6) and time-dependent specific flow rate q(t, x). The function R Q(t) = q(t, x) dx is a total flow rate. W

Formulation of the Problem The equation (0.7)

is considered in domain

(see Fig. 1.1) Viscosity µ and storage coefficient β ∗ = ϕct are constants. The term of equation (0.7) with Delta-function δ is responsible for the fluid inflow q(t, x) to the wellbore. The bottom and the top of the reservoir are impermeable: ∂P = 0. (1.1) ∂z z=0,H At the remaining part of the outer boundary Γ either pressure P Γ = pe (t)

(1.2)

∂P =0 ∂n Γ

(1.3)

or zero flow condition

is specified. Here n is the outer normal vector to the boundary Γ. At the initial time t = 0 the pressure distribution is uniform: P |t=0 = pe (0). (1.4) The sink, which models a wellbore, is distributed along the curve W (0.6). Pressure along the wellbore is distributed according to the Hagen-Poiseuille law: ∂ ∂P γπR4 Ψ(x) = q(t, x), where Ψ(x) = (1.5) ∂x ∂x 8µ Functions R(x), γ(x) are radius and empirical conductivity of the wellbore. Sections with horizontal well

Sections with fracture

y z

x

Sections without horizontal well

Fig. 1.1. Schematic view of the reservoir

SPE number

Advanced horizontal well model

The inner boundary condition is specified either as given pressure pw (t) at the start point of the horizontal wellbore P = pw (t),

(x = xl , y = yw (xl ), z = zw (xl ))

(1.6)

or as given total flow rate Q(t) of the well Z q(t, x) dx = Q(t).

(1.7)

W

At the end point of the wellbore the flow rate equals zero: Ψ(x)

∂P = 0, ∂x

(x = xr , y = yw (xr ), z = zw (xr )).

(1.8)

The mathematical statement of the problem is the following: it is required to satisfy equation (0.7) in the domain Ω with given • initial condition (1.4); • outer boundary condition (1.1), and either (1.2) or (1.3); • inner boundary condition either (1.6), or (1.7);

3

All functions are assumed to be piecewise constant with respect to x, equal to their mean values over intervals [xi−1/2 , xi+1/2 ] (see Fig. 1.2). Average on x over the interval [xi−1/2 , xi+1/2 ] value of function is denoted by the same symbol as the function with upper index i. Discretization of the second x-derivative is performed as ∂ kx (x) ∂p = ∂x µ ∂x x=xi " # i+1/2 i+1 i−1/2 i p − pi kx p − pi−1 1 kx = − , hi µ ∆i µ ∆i−1 where kxi−1 + kxi . 2 In the case of the first-type boundary condition (p Γ = 0), it is formally assumed that kxi+1/2 =

p0 = 0,

• closure boundary condition (1.5), (1.8). The reservoir pressure and either flow rate or bottomhole flowing pressure should be determined. Replacing pressure P with a function p(t, x, y, z) = P (t, x, y, z) − pe (t)

Nr X

q i (t)hi = Q(t)

Laplace Image of Solution The idea of solving for-

mulated problem is following:

stands for (1.7). 1.2.2 Pressure Drop at the Wellbore. Let us transform the formula (1.5) describing the pressure drop along the wellbore. Discretization on x yields pj+1 − pj j+1/2 pj − pj−1 j−1/2 Ψ − Ψ = hj q j . ∆j ∆j−1

• Laplace transform with respect to time; • Discretization along Ox;

Making sum of these equations with indexes j = n + 1, . . . , j = Nr , taking into account (1.8) lead us to expression

• Fourier transform in y, z • Linear algebraic system for Fourier coefficients; • Inverse Fourier and Laplace transformations for solution images.

pn+1 − pn = −

Sections with

Sections with fracture horizontal well Wellbore The Ox-axis is 1.2.1 Discretization Along the splitted into N subintervals by the set of points x1 , . . . , xN , y where xi < xi+1 . Thez domain Ω is accordingly divided into N layers Ωi having the shape of right parallelepipeds: horizontal well Ωi = {~x ∈ Ω | x ∈ [xi−Sections 1 , xwithout i+ 1 ]}, i = 1, . . . , N. x

2

2

+ xi±1 ) denote boundaries of each layer

x1/2 = 0, xNl −1/2 = xl , x

(1.11)

i=Nl

in the case of condition (1.3), we arrive to the similar problem, but with zero outer boundary conditions.

1 (xi 2

pN +1 = 0.

Integrals with respect to x along the wellbore are replaced by sums. Hence, the sum

p(t, x, y, z) = P (t, x, y, z) − pe (0)

Points xi±1/2 = such that

kxi−1/2 =

(1.9) The boundary condition of the second type (∂p/∂n Γ=0 = 0) leads to p0 = p1 , pN +1 = pN . (1.10)

in the case of condition (1.2), or

1.2

kxi + kxi+1 , 2

=x , x

Schematic view of reservoir r Nr +1/2 with horizontal well N +1/2

= X.

Nr X

∆n Ψn+1/2

hj q j .

j=n+1

Summation of this equations for n = Nl , . . . , i − 1 gives pi − pNl = −

i−1 X n=Nl

∆n Ψn+1/2

Nr X

hj q j .

(1.12)

j=n+1

Granting (1.6), one can derive that flux through the start point of the wellbore and pressure is connected via equation pNl = pw − F −

Nr ∆Nl X hj q j . 2ΨNl j=N l

The following notations are introduced ∆i = xi+1 − xi ,

hi = xi+1/2 − xi−1/2 =

∆i + ∆i−1 . 2

The wellbore area is splitted into Nr − Nl layers, the interval between the left xl (right xr ) end of the wellbore and the left x = 0 (right x = X) boundary of the reservoir is splitted into Nl (N − Nr ) layers. hi x1/2 =0 xNl -1/2 =xl

x i-1

xi

i

x i+1

xN+1/2 =xr r

x N+1/2 =X

This equation together with (1.12) gives an equation for flow rates q i : Nr X pi = pw − F − bij hj q j , (1.13) j=Nl

where   ∆Nl /(2ΨNl ),    k−1  P  ∆Nl /(2ΨNl ) + ∆n /Ψn+1/2 , i bj = n=Nl   i−1  P  N  ∆n /Ψn+1/2 , ∆Nl /(2Ψ l ) +

i = Nl ; Nl + 1 6 j 6 i; i + 1 6 j 6 Nr .

n=Nl

(1.14) Fig. 1.2. Discretization along Ox

These expressions will be used below.

4


1.2.3 Laplace Transformation. The Laplace images [6] of ¯ F¯ consequently. main functions are denoted by q¯i , p¯w , Q, An application of Laplace transformation procedure to the system of equations from Section 1.1 gives us the following system of linear equations, which will be solved by Fourier analysis: kyi hi ∂ p¯i kzi hi ∂ p¯i εi,1 p¯i+1 −(εi,0 +hi β ∗ s)¯ pi +εi,−1 p¯i−1 + + + µ ∂y 2 µ ∂z 2 i i + q¯i hi δ y − yw , z − zw − hi F¯0 = 0, for i = 1, . . . , N ; = p¯w − F¯ −

p¯i

W

Nr X

p¯Nl −1/2

bij hj q¯j ;

W

= p¯w − F¯ ;

j=Nl Nr X

∂ p¯i = 0; ∂z z=0;H

¯ q¯j hj = Q;

j=Nl

p¯i y=0;Y

= 0;

0

N +1

p¯ = p¯

i

∂ p¯ = 0; ∂y y=0;Y

=0

p¯0 = p¯1 ;

i+1/2

i+1/2

p¯N +1 = p¯N .

i−1/2

(1.15)

i−1/2

Coefficients are given by (1.14). Functions q¯i equal zero in the layers without a wellbore:

1.3

i = 0, . . . , Nl − 1

and

1 νˆ00 = 0,

1 νˆn0 =

2(1 − (−1)n ) , for n > 1; πn

1.3.2 Tridiagonal Matrix for Fourier Coefficients. Plugging the obtained decompositions into the first equation of the system (1.15) and collection of coefficients at the same products of basic functions give an algebraic system for unknown coefficients c¯inm :

i = Nr + 1, . . . , N. (1.16)

problem for the Laplace images of function is solved in terms of trigonometric Fourier series. The problem (1.15) is reduced to a linear system of algebraic equations for the coefficients of Fourier series for the unknown pressure function. 1.3.1 Fourier Decomposition. Fourier series are performed with respect to y and z. The Laplace image of the pressure is expressed as i

1 1 2 νnm = −ˆ νnm hi β ∗ (pe (0) − s¯ pe ), νnm =0 i i ky nπ 2 kz mπ 2 Θ2 = + + β ∗ s. µ Y µ H

p¯ (s, y, z) =

c¯0nm = 0,

where constants κnm are κ00 =

1 , 4

κ10 = κ01 =

1 , 2

κnm = 1

(m, n > 1).

(1.18)

In the case of the first-type boundary condition (p|Γ = 0) basic functions ϕ1n (y), ψm (z) are sines and cosines: ϕ1n (y) = sin

πny , Y

ψm (z) = cos

πmz . H

In the case of the second-type boundary (∂p/∂n|Γ=0 = 0) basic functions are only cosines: ϕ2n (y) = cos

πny , Y

ψm (z) = cos

πmz . H

(1.19) conditions

(1.20)

Throughout the paper the upper index α = 1 corresponds to sine-decomposition over the interval (0, Y ), whereas index α = 2 signals the cosine-decomposition. i1 i2 By symbols γnm , γnm we denote coefficients of the Fourier i i decomposition for the Dirac delta function at point (yw , zw ), the coordinates of the wellbore in the i-th layer: iα γnm

4 α i i = ϕn (yw )ψm (zw ), YH

α = 1, 2.

+1 c¯N nm = 0

(1.22)

in case of the first-type boundary conditions, and c¯0nm = c¯1nm ,

+1 c¯N ¯N nm = c nm

(1.23)

for the boundary conditions of the second type. Application of the standard sweep procedure to the system (1.21) gives expressions for the coefficients c¯inm :

(1.17)

n,m=0

(α = 1, 2),

For fixed indices n and m, equations (1.21) form a tridiagonal system of algebraic equations for unknown c¯inm , with the property of the diagonal predominance. Index i takes integer values from 1 to N . According to boundary conditions, for the limiting values of the index one has

c¯inm =

κnm c¯inm (s)ϕα n (y)ψm (z),

(1.21)

where

Construction of the Solution The solution of the

∞ X

1 νˆnm = 0 for m > 0.

iα α = −¯ q i hi γnm − νnm ,

kx kx kx kx ; εi,0 = + ; εi,−1 = . µ∆i µ∆i µ∆i−1 µ∆i−1

for

The decomposition of unity with respect to sines, which is needed for the first-type of boundary condition has the form

or

bij

qi = 0

Given function F¯0 from (1.15) is determined as follows: ( β ∗ [s¯ pe (s) − pe (0)], α = 1, 0 ¯ F = 0, α = 2.

2 εi,1 c¯i+1 cinm + εi,−1 c¯i−1 nm − (Θ hi + εi,0 )¯ nm =

Here the following notations are used for i = 1, . . . , N εi,1 =

SPE number

Nr X

ij jα α ωnm q¯j hj γnm + νnm .

(1.24)

j=Nl

In (1.24) formulae (1.16) were taken into account. In case of α the second-type boundary conditions terms νnm in (1.21) vanish, hence this system of equations becomes homogenous with respect to c¯inm and q i . 1.3.3 Closure Relations and Determination of Flow Rates. To determine the flow rates q¯i we use the equality condition between reservoir pressure at the wellbore and wellbore pressure. Functions p¯i are plugged into the second equation of the system (1.15). A pressure coming from the reservoir is calculated at points √ √ i i i i + Ri / 2), (ˆ yw , zˆw ) = (yw + Ri / 2, zw where Ri is the average wellbore radius in i-th layer. So, one arrives to the following system Nr n X ∞ X j=Nl

o i ij jα α i ϕn yˆw ψm zˆw κnm ωnm γnm + bij hj q¯j =

n,m=0

= p¯w − F¯ −

Nr X ∞ X

ij α i κn0 ωn0 νn0 ϕα ˆw . n y

(1.25)

j=Nl n=0

Here α = 1, 2 depending on the type of the boundary problem, and i = Nl , . . . , Nr .

SPE number


Relations (1.25) form the system of Nr − Nl + 1 algebraic equations for Nr −Nl +1 unknowns q¯j . System for the determination of local inflows q¯j is closed if pressure pw (t) at the start point of the wellbore is specified. Otherwise, this pressure pw is treated as an extra unknown, and the linear algebraic system for q¯i must be extended by the equation (1.11). Column vectors ~ q , f~ and matrix A = {Aij } are: f~ = (f Nl +¯ pw , ..., f Nr +¯ pw )T ,

~ q = (¯ q Nl , ..., q¯Nr )T ,

A = {Aij },

5

2.1

The Mathematical Model Hydraulic fractures are taken into account by adding several terms to the right hand side of the main equation (0.7). These terms are responsible for the inflow to the wellbore from fractures: kx (x) ∂p ky (x) ∂ 2 p kz (x) ∂ 2 p ∂p ∂ β∗ = + + + 2 ∂t ∂x µ ∂x µ ∂y µ ∂z 2 + q(t, x)δ y − yw (x), z − zw (x) + Nr X

+ Nr X

f = −F¯ − i

∞ X

ij α κn0 ωn0 νn0 ϕα n

Aij = bij hj +

ij jα i i κnm ωnm γnm hj ϕα ˆw ψm zˆw . n y

n,m=0

Here i and j are indexes of the rows and the columns of the matrix A; upper index T is matrix transposition. These notations allow one to represent equations (1.25) in the standard form A~ q = f~. In case of the given total flux Q(t) the unknown vector ~ q, and function p¯w determines by the system (1.25), extended by the relation (1.11). This leads to the similar system A0 q~0 = f~0 for “extended” vectors q~0 , f~0 and matrix A0 : ¯ T, ~ q = (¯ q Nl , ..., q¯Nr , p¯w )T , f~ = (f Nl , ..., f Nr , Q)   ANl Nl . . . ANl Nr −1  .. .. ..   . . .  A0 =    AN N . . . AN N −1  l r l l h Nl ... h Nr 0 Solution of this system provides one with functions q¯j , which give expression of pressure p¯ inside the reservoir in the form of the Fourier series with coefficients (1.24). The inverse Laplace transformation restores p as the function of time t.

2

(2.1)

k=Nl

i yˆw ,

j=Nl n=0 ∞ X

Φk (t, y)δ(x − xk ) − F 0 (t).

Heterogenous Reservoir With Multiple Hydraulic Fractures

The algorithm described in the previous section can be adopted to the case with multiple transverse hydraulic fractures along the horizontal wellbore. It is assumed that each layer along Ox axis (except for the first and the last one) may contain a hydraulic fracture (HF), located at the center point x = xi of the corresponsing layer (see Fig. 2.1). We assume that the height of each HF equals the reservoir thickness Z. In horizontal direction along Oy-axis each fracture has its individual size. Hydraulic fracture (HF)

Fluid inflow to the fracture in i-th layer is determined by the discharge intensity Φi (t, y). For simplicity and due to the thickness of the reservoir in z-direction we assume that functions Φi do not depend on z. In the absence of the fracture in i-th layer, the corresponding discharge intensity is equal to zero: Φi ≡ 0. The finite conductivity of the fracture is taken into account in the following way: pressure drop in the HF satisfies the equation ∂p ∂ CF = Φ(t, y). (2.2) ∂y ∂y Here the empirical function CF determines the conductivity of the HF. Differential operators in (2.2) are taken only with respect to variables (y, z). The following boundary conditions for equation (2.2) are taken: ∂p = 0, p yi = pi W (2.3) CF 0 ∂y yi ,yi r

l

yli

yri

Here and are the left and right limiting values of y along the HF in i-th layer; y0i is the position of the wellbore; pi W is the wellbore pressure in i-th layer.

2.2

Discretization and Laplace Transformation

Using the same approach as in Section 1, we rewrite the equation (2.1) in Laplace images as follows i+1/2

kx

µ +

i−1/2

p¯i+1 − p¯i kx − ∆i µ

kyi hi ∂ 2 p¯i p¯i − p¯i−1 + + ∆i−1 µ ∂y 2

kzi hi ∂ 2 p¯i i i − sβ ∗ hi p¯i + hi q¯N1 (i) (s)δ(y − yw , z − zw )+ µ ∂z 2 ¯ i (s, y) − hi F¯ 0 (s) = 0. (2.4) +Φ

Definition of integers N1 (i) will be given below. Note that • each layer along Ox-axis contains not more than one HF; • if presented, HF is located at the some mesh point xi . 2.3 Fourier Decomposition Next, we represent functions p¯i by the Fourier series (1.17). Basic functions are defined in (1.19), (1.20). Coefficients c¯inm (s) are unknown values. Substituting (1.17) into (2.4) we obtain tridiagonal algebraic system of equations: εi,1 c¯i+1 ¯inm + εi,−1 c¯i−1 nm − εi,0 c nm = iα iα α = −γnm hi q¯N1 (i) − β¯nm − νnm .

y

z

Here i+1/2

εi,1 = wellbore

x i-1/2

xi

x i+1/2

x

Fig. 2.1. Hydraulic fracture in i-th layer

kx , µ∆i

i−1/2

εi,−1 =

kx , µ∆i−1

εi,0 = εi,1 + εi,−1 + hi Θ2 ,

kyi n2 π 2 kzi m2 π 2 4 α i iα i + + β ∗ s, γnm = ϕn (yw )ψm (zw ). µY 2 µH 2 YH n 2(1 − (−1) ) 1 1 2 νn0 =− hi β ∗ (pe (0)−s¯ pe ), νnm = 0, νnm = 0. πn iα i ¯ Here βnm are Fourier coefficients of function Φ . Index α specifies the type of the boundary conditions as previously (see section 1.3.2). Θ2 =

xi-1

(2.5)

6


It is worth to note, that the second and the third terms in the right-hand side of (2.5) are nonzero only for m = 0. Again, using the sweep method, we obtain formulaes for unknowns c¯inm , which contain unknown flow rates q¯i : Nr X

c¯inm =

jα ij ωnm γnm hj q¯N1 (j) ,

m > 0;

j=Nl Nr X

(2.6)

ij jα iα α ωn0 γn0 hj q¯N1 (j) + β¯n0 + νn0

c¯in0 =

j=Nl

2.4

Closure Relations Closure relations are needed to ¯ i . The foldetermine unknown flow rates q¯i and intensities Φ lowing conditions will be used: • The reservoir pressure at the wellbore coincide with the wellbore pressure; • The reservoir pressure at the HF coincide with the HF Hydraulic fracture (HF) pressure. At least three approaches can be used to determine the reservoir pressure at the HF:

z

• Pressure at the HF coincide with the average layer pressure: pf = pi . (2.7)

y

• Pressure at the HF is calculated by linear approximation (see Appendix): pf =

x

∆i−1 (4pi − pi−1 ) + ∆i (4pi − pi+1 ) . 3(∆i−1 + ∆i )

(2.8)

x

z

x

in the case of the first-type boundary conditions, and i2 β¯00 =

4 YH

4 i2 β¯n0 = πn

N2 (i)

X

q¯k ,

k=N1 (i)+1 N2 (i)

X k=N1 (i)+1

πn∆ik q¯k 2 i ϕ ek−1/2 sin , i n 2Y Dk

x

k=Nl

The flow rate from i-th layer into each vertical strip of the hydraulic fracture is numbered by q¯k , where k = N1 (i) + 1, . . . , N2 (i). So, there are N2 (i) − N1 (i) strips and the same number of flow rates. If N2 (i) = N1 (i) the fracture in i-th layer is absent. The total flow rate from i-th layer into the fracture NP 2 (i) q¯N1 (k) . Hence, the overall flow rate to the is equal to k=N1 (i)+1

wellbore from the reservoir and from all fractures is given by N Pr i q¯total , where i=Nl N2 (i) i q¯total (s) = hi q¯N1 (i) +

N1(i)+2

DNi (i)+1 DNi (i)+2 1 1 e

e

...

q

N2 (i)

DNi (i) 2

...

N1(i)+1

N1(i)

e

N2(i)

y

Fig. 2.2. Representation of the hydraulic fracture at section x = xi (in (y, z)-plane)

The region occupied by the HF in i-th layer is: HFi = {x = xi ,

eN1 (i) 6 y 6 eN2 (i) ,

0 6 z 6 H}.

¯i

For simplicity we consider function Φ (s, y) piecewise constant with respect to y (see Fig. 2.2). Intervals of constancy, so-called strips, are y ∈ [eik−1 , eik ] with function values ¯ i (s, y) = q¯k /Dki . Here Dki = (eik − eik−1 )H is the area of the Φ corresponding strip; q¯k (k = N1 (i)+1, . . . , N2 (i)) are flow rates to each strip. Hence, ¯ i (s, y) = Φ

∞ X

iα α β¯n0 ϕn (y)

n=0

where for n > 0 4 i1 β¯n0 = πn

N2 (i)

X k=N1 (i)+1

N2 X

k bik q¯total .

(2.10)

k=N1

no HF

no HF

q

q¯k .

The wellbore pressure in i-th layer is determined from equation (1.13), where instead of q¯k we use the total inflow to the k wellbore q¯total . This gives a linear system of Nr − Nl + 1 algebraic equations for unknown flow rates q¯N1 (k) (k = Nl , . . . , Nr ): √ i √ i + Ri / 2, zw p¯i W = p¯i (s, yw + Ri / 2) = = p¯w − F¯ −

well N1(i)+1

X k=N1 (i)+1

Hydraulic fracture (HF)

q

n>0

in the case of the second-type boundary conditions. Hereafter the notations eii−1/2 = (eik + eik−1 )/2, ∆ik = eik − eik−1 are used. The following numbering scheme of flow rates q¯j is adopted in the sequel. Each i-th layer is supplied with two numbers: N1 (i) and N2 (i). The flow rate q¯ directly to the wellbore from i-th layer is numbered by N1 (i). So, the total flow rate to the wellbore directly from the reservoir (not taking into account N Pr hk q¯N1 (k) . inflows from fractures) is given by

wellbore

• Pressure at the HF is taken from the solution for 1D nonstationary fluid motion in i-th layer (see Appendix): i-1 i-1/2 i+1/2 i λ∆ 1 3λ∆ i p ¯ λ∆ ch p¯f = + ch − 4 4 2sh 3λ∆ 4 λ∆ −(¯ pi−1 + p¯i+1 )sh (2.9) 4

x

SPE number

πn∆ik q¯k 1 i ϕ ek−1/2 sin , i n 2Y Dk

Coefficients bik are determined by (1.14). The second part of the system, which determines flow rates q¯k , follows from either of conditions (2.7)–(2.9). In all cases we equate the reservoir pressure in the middle of each strip y ∈ [eik−1 , eik ] at the hydraulic fracture with the fluid pressure in the fracture. The closing relations are p¯f (s, eik−1/2 , H/2) = p¯ff liud (s, eik−1/2 ),

(2.11)

k = N1 (i) + 1, . . . , N2 (i). At that, the left-hand side is taken from either of formulae (2.7), (2.8), or (2.9). The right-hand side of (2.11) is determined by equation (2.2) with boundary conditions (2.3), where yli = eN1 (i) , p¯i W

yri = eN2 (i) , √ i √ i = p¯i (s, yw + Ri / 2, zw + Ri / 2).

Discretization of the equation (2.2) along Oy-axis and integration in the same way as it was done earlier for the wellbore pressure equation (1.13) allows one to express pressure at each stripe of the HF in terms of the wellbore pressure in the i-th layer. Finally, the conditions (2.10), (2.11) provide with the required linear algebraic equations for unknowns q¯i . Solution of this system gives functions q¯j , which finally produces expression of pressure p¯ inside the whole reservoir in the form of the Fourier series (1.17) with coefficients (2.6). The inverse Laplace transformation restores functions p as the function of time t.

SPE number

3


Verification

A reliable finite-difference simulator Eclipse–100 was used to validate this new model, which can have application not only for slightly compressible fluid flow, but also for gas flows. The usual approach to switch from liquid to gas is a use of pseudofunctions, e.g. pseudo-pressure and pseudo-time [?]. Standard SPE tests for horizontal well simulators are not applicable here, because they are developed for multi-phase flow instead of 1-phase flow considered in our model. Tests have started from simple reservoir models, which complexity was increased with each step. All types of tests are listed in the tree below. • Homogenous Reservoir

3.1.3 Parameters to Compare. The AHW model permits to track the following reservoir and production parameters: • Pressure as a function of time at the starting point of the wellbore; • Flow rates from each layer as a function of time; • Total flow rate of the wellbore; • Pressure as a function of y at any fixed time t = const, and any section x = const. 3.1.4 Stability of the Models. Before running the verification tests, both models in Eclipse–100 and AHW model were put on trials with different grid block dimensions. Eclipse–100 gives stable result with an inessential distinction when:

– No fractures

• the number of grid blocks in Ox, Oy directions extends 3, at least (with central block containing well);

– Fractured

• the number of grid blocks in Oz direction is greater than 1.

• Heterogenous (along Ox) reservoir • Deviated horizontal well

Models of the Reservoirs

3.1.1 Numerical Hydrodynamic Model. Productivity layer is supposed to have right parallelepiped form with dimensions 1100 × 100 × 10 m. Cell dimensions of Cartesian grid in Eclipse are selected according to the case tested. Several simulations were done on different grids to be convinced that the result is not sensitive to the grid selection. Simulation of the vertical hydraulic fracture in Oyz–plane was done on refined grid with several buffer blocks in front of and behind the fracture (see Fig. 3.9 for explanation). Permeability of grid cells is chosen depending upon case considered (type of reservoir, fractured or not). Porosity ϕ and total compressibility ct are constant values in all tests. Only 1-phase fluid flow is considered. A numerical aquifer is used in Eclipse–100 to simulate cases with constant outer boundary pressure condition. Fractures in Eclipse–100 are modeled via local permeability increase of the cells in refined grid: the grid in Ox–direction was successively decreased from 100 m to 1 cm. The permeability of the cells containing fracture is chosen so that the equivalent fracture conductivity kf wf , where kf is real fracture permeability, and wf is fracture width, has the same value both in model and reservoir.

Note that the permeability of the reservoir can be a function of only x variable (see Section 1). The AHW model showed good convergence on different grids1 . Typical plots are present on Fig. 3.1. This model requires to have at least 5 sections in Ox direction with mandatory 3 mid-sections containing horizontal well. Such requirement is imposed by the solution method. 30 25 Flow rate, m3/day

3.1

– Number of layers of heterogeneity along the horizontal well;

7(3) Sections

15 10

0

100 Time, days

150

14,8 3 14,7 2,5 14,6 2 14,5 1,5 14,4 1 14,3 0,5 14,2

7(3) Sections 15(13) Sections 7(5) Sections 7(5) Sections 7(3) Sections

0 14,1 140

– Hydraulic fractures, its properties;

25,0

– Wellbore radius and wellbore conductivity;

50 25,2

100

150

Time, days 25,6 25,4

200 25,8

Time, days

– Boundary conditions both at the wellbore and at the outer boundary;

Fig. 3.2. Zoomed plots from Fig. 3.1

3.2

Homogenous Reservoirs 4,5 4

term “grid” for semi-analytical AHW model means how detailed is the discretization along Ox axis. 3,5 error, %

1 The

200

415 3,5 14,9

– Permeability of each layer in x, y and z directions;

– Time period for simulation.

50

Fig. 3.1. Flow rate curves for AHW. Varying discretization along Ox.

4,5

error, % FlowRelative rate, m3/day

– Start and end point of the horizontal well, its deviation survey;

7(5) Sections

0

– Dimensions of the reservoir; – Total compressibility ct

15(13) Sections 20

5

3.1.2 AHW Model. The model is defined by the following parameters:

– Fluid viscosity µ;

7

3 2,5 7(3) Sections

26,0


Flow rate, m3/day

200

150

Eclipse AHW

100

50

0 0

50

100 Time, days

150

200

Fig. 3.3. Flow rates. Case #1. 10000 10,0

Eclipse

9,0

8000

Eclipse AHW

100

0 0

100

200

300 Time, days

400

500

600

Fig. 3.5. Flow rates. Case #2. 8000

Case #1 (see Appendix B for details) shows good correspondence between Eclipse Eclipse and AHW (Fig. 3.3, Fig. 3.4). The overall disagreementAHW for cumulative production at late times is 6000 even less than 1%. Of course, the mismatch in flow rates during first several days is not taken into account. The behavior of production curves make this case nearly a reference one. 4000 Case #2 simulates a depletion drive type of production. Besides this, the time of forecast increased up to 500 days (Fig. 3.5). 2000 A minor disagreement between Eclipse–100 and AHW cumulative production curves has appeared from day 200 of production, after some oscillations in AHW flow rates occur due to transverse Laplace transform procedure. This disappoint0 ing singularity of the semi-analytical approach can 600 be killed 0 100 200 300 400 500 by using another procedure, e.g. Stephest algorithm instead of Time, days transform with Legendre polynomials. Latter the mismatch in flow rates (oscillations around zero) has generated a divergence of 6.5% in cumulative production on late time stages. 3.2.2 Multiple Transverse Fractures Case #3 (see Appendix B) is aimed to make sure that AHW and Eclipse–100 give similar results on refined grid Fig. 3.9, prepared to simulate vertical transverse fractures (Fig. 3.6). It is seen that on refined grid without fracture both simulators perform equally, which shows a good stability of numerical methods used. 10,0

Flow Rate

7,0

6000 Error, %

150

50

AHW

8,0

9,0

Cumulative

6,0

8,0

5,0

7,0

Flow Rate

4000 4,0

6,0

Cumulative

Error, %

Cumulative, m m3

200

Flow raate, m3/day

3.2.1 No Fractures We start from simple homogenous reservoirs with straight (not deviated) horizontal well positioned at the symmetry axis (along Ox) of the rectangular reservoir. All physical parameters of the rock and the fluid remain the same in all these tests. The conductivity of the well in new model was set manually while Eclipse–100 calculates pressure drop at the wellbore due to friction according to engineering formulas used in pipeline flow calculations. The matching with pressure drop along the wellbore in both simulators was done, though a minor difference existed: semi-analytical model uses a Poiseuille equation (1.5) for laminar flows only, while Eclipse–100 differentiates laminar and turbulent types of flow and automatically switches between them. Nevertheless, the divergence between these two approaches is negligible and does not influence on overall production characteristics. Such parameters as porosity, compressibility, viscosity are just multiplicative constants for one-phase flows; permeability influences on duration of different flow regimes in the reservoir and flow regime in the wellbore. Thus, the storage coefficient ϕct in all tests is 4.05 · 10−5 bar−1 , the permeability in all directions is 0.1 mD, the fluid viscosity µ = 0.88 cP.

SPE number

Cumulative, m3

8

3,0

2000 2,0 1,0

5,0 4,0 3,0

0 0,0

2,0

00

50 50

100 100

Time, days

150 150

Time, days

Fig. 3.4. Relative Error. Case #1.

200 200

1,0 0,0 0

50

100 Time, days

150

200


Hydraulic fractures with significantly different conductivities were simulated in cases #4–7. Case #5 gave minor oscillations of flow rate (Fig. 3.8) due to very close placement of fracture tip to the outer boundary of the reservoir, and it has caused in pressure gradient approximately 100 bar within 5 meters of porous media. To reduce the relative error in this case (7.5% between blue and red lines) an increase of number of basic functions is needed. Of course, the time of simulation grows much faster than the rate of convergency is: the number

SPE number

Advanced horizontal well model 10,0

of basic functions in both directions (y, z) was increased by factor 1.5, but the accuracy has increased only by 2.4% (green–red lines, Fig. 3.8).

9,0 8,0

Error, %

9,0 8,0

Error, %

10,0

7,0

Flow Rate

6,0

Cumulative

9

7,0

Flow Rate

6,0

Cumulative

5,0 4,0 3,0 2,0

5,0

1,0

4,0

0,0

3,0

0

50

2,0

100 Time, days

150

200


1,0 0,0 0

50

100 Time, days

150

200

Fig. 3.7. Relative Error. Case #4. 250

Cases #4, #8 with more or less real fracture conductivity (1000–3000 mD · m) shows good agreement (approximately 2%) with numerical simulator (Fig. 3.7, Fig. 3.12), while other cases with extremely low fracture conductivity (Fig. 3.10, Fig. 3.11) gave a relative error of 3–4% in flow rate, and has resulted in total error of 2–4% in cumulative production. 10,0

Eclipse

9,0

AHW short series

150

8,0

AHW extended series 100

Error, %

Flow w rate, m3/day

200

50

7,0

Flow Rate

6,0

Cumulative

5,0 4,0 3,0

0

2,0

0

50

100 Time, days

150

200

1,0 0,0 0


5*100m

1*25m

2*10m

2*1m

1*2.5m

AHW short series

1*0.25m

2*0.1m

2*0.1m

10*0.01m

2*1m

1*0.25m

2*10m

1*2.5m

1*25m

5*100m

12000

200

10,0 9,0

AHW externed series

10000

8,0

5*100m - 5 ggrid blocks with dx =100m

8000

Fig. 3.9. Progressive grid refinement for fracture simulation

4000 10,0

2000

Error, %

Cumulative, m33

150


fracture in one cell of dx =0.01m

7,0

Flow Rate

6,0

Cumulative

5,0 4,0 3,0

9,0

2,0

0 8,0 7,0 0 Error, %

100 Time, days

Eclipse

14000

6000

50

50

6,0

100

Flow 150 Rate

Time, days

Cumulative

200

1,0 0,0 0

5,0 4,0

50

100 Time, days

150

200


3,0 2,0 1,0 0,0 0

50

100 Time, days

150


200

Advanced Horizontal Well model has an option to simulate non-symmetry vertical fractures, which extreme case is when only one wing of fracture exist. Such experiment was launched using a grid for symmetry fractures, and the result is quite good: disagreement between AHW and Eclipse is small (Fig. 3.13), and the diagram for pressure distribution at the section with fracture show how the solution is sensitive to fracture parameters (Fig. 3.14).

10


SPE number

10,0

200

9,0 8,0 7,0

T=0,6

100

6,0

Error, %

Presssure, bar

150

T=3,5

5,0 4,0

T=9,1 50

3,0

T=21,3

2,0

T=148

Flow Rate

1,0

0 0

10

20

30

40

50

60

70

80

90

Cumulative

0,0

100

0

Y, meters

Fig. 3.14. Pressure distribution along the section with non-symmetry fracture. Case #9.

Before making multiple fractures, the stability of calculations was checked on refined grid. As in Case #3, each fracture in Eclipse was simulated via local increase of permeability so that the equivalent permeability (in mD · m) both in AHW and Eclipse was similar. Thus, on Fig. 3.15 there is a production forecast on fine grid with original reservoir properties; plots on Fig. 3.16 are product of 3 fractures simulation at X0 = 550m, X0 = 650m, and X0 = 750m. Again, some oscillations in production curves occur, which has further influence on cumulative production, but the relative error is satisfactory (Fig. 3.17). As discussed before, the modification of transverse Laplace transform procedure is required to overcome this issue.

Flow w rate, m3/day

200

150

Eclipse AHW

100

50

50

100 Time, days

150

200

Fig. 3.17. Relative error. Case #11.

3.2.3 Straight Well Shifted to the Boundary All previous tests were conducted with a well positioned at the symmetry axis of the reservoir. Here (Appendix B, cases #12–14) the well is shifted towards one or two boundaries of the reservoir, no fractures. Extreme cases are considered, when the well is only 0.5 m off the external boundary. In case #12 the well is drawn up to the impermeable upper boundary Z = const (top of the reservoir). An overall error in flow rate is approximately 3% (Fig. 3.18), which shows that AHW and numerical model of Eclipse–100 are still in satisfactory agreement. A disagreement between AHW and Eclipse– 100 lay between 2% and 4% depending upon well deviation from the center line of the reservoir. A better coincidence between Eclipse–100 and AHW (only 0.7% of relative error) is when the well is moved close to open boundary, i.e. Y = const (case #13, Fig. 3.19). Such difference in two previous cases is due to the fact that the lesser number of basic functions is used in Oz–direction in comparison with Oy–direction (100 against 200). Increase in series length in Oz–direction causes in significant increase of time resources for execution with only several fraction of percent of accuracy increase. Thus, the golden mean is preferable in this case: both satisfactory accuracy and running time of the program.

0

150

0

50

100 Time, days

150

200

Flow w rate, m3/day

Fig. 3.15. Flow rates. Case #10. 12000 350 Eclipse

10000 300

AHW

Eclipse

Cumulative, m33 Flow w rate, m3/day

250 8000

Eclipse 100

AHW

50

AHW

200 6000 0

150 4000 100

0

2000 50

100 Time, days

150

200

Fig. 3.18. Flow rates. Case #12.

00

12000

50 50

100 100 Time,days days Time,

150 150

Fig. 3.16. Flow rates. Case #11. 14000 Eclipse AHW

200 200

Case #14 is the superposition of two previous cases. The number of harmonicsEclipse in Fourier series was optimized, and the 10000 total relative error has AHWdropped to 0.6%. This was done just to show that the adjustment between models is possible, but 8000 prior to forecast execution one must set the appropriate mistake 6000 according to common sense and select the series length. This will remain as an option of AHW, similar to Eclipse–100 time 4000 step choice. Cumulative, m33

00

12000

50

Cumulative, m33

10000 8000 6000

2000 0 0

50

100

150

200

SPE number


10,0

#16 (Fig. 3.22) shows this disagreement between Eclipse–100 and AHW. After the conductivity of the well in AHW was manually decreased, the results of simulations nearly coincided (green curves on referred plots).

9,0 8,0 7,0 Error, %

11

6,0

Flow Rate

5,0

Cumulative

5000 4500

4,0

4000 Flow w rate, m3/day

3,0 2,0 1,0 0,0 0

50

100 Time, days

150

200

Eclipse

3500 3000

AHW

2500

AHW low conduct.

2000 1500 1000


500 0

150

0

50

Flow w rate, m3/day

Eclipse 100

100 Time, days

150

200


AHW

1000000 10,0

900000

Eclipse

9,0

800000

0 0

50

100 Time, days

150

200

AHW

7,0

AHW low conduct. conduct

600000 6,0

Flow Rate

500000 5,0

Cumulative

400000 4,0

300000 3,0

200000 2,0


100000 1,0

0

18000 Cases with Extreme Permeability Values Three 3.2.4

Cumulative, m33

cases with reservoirEclipse permeability 0.01 mD, 10 mD and 10 D 16000 were14000 studied (see table below). Extremely low permeability AHW (case #15, Fig. 3.21, Appendix B) causes in inefficient dis12000 agreement between AHW and Eclipse–100 at early-time flow 10000 rates. The flow rates on steady-state regime are nearly equal: ±0.34% 8000 of relative error. It is known that Eclipse–100 is not a good tool for reservoirs with extremely low permeability, and 6000 probably this is the case. 4000 30 2000

0,0

00

50 50

100 100 Time, days Time, days

150 150

200 200


Cases #17, 18 are similar except wellbore friction option. Once increased reservoir permeability up to 10 D has resulted in sufficiently turbulent flow in the wellbore, and the disagreement between two simulators became visibly high (Fig. 3.24). After the wellbore friction option was disabled, the agreement between simulators has been restored (Fig. 3.25).

0 25

180000

0

50

20

100 Time, days

150 Eclipse

200

160000

AHW

140000 Flow w rate, m3/day

Flow w rate, m3/day

8,0

700000

Error, %

Cumulative, m33

50

15 10 5

120000 100000

Eclipse

80000

AHW

60000 40000

0

20000

0

100

200 Time, days

300

400

0 0

50


150

Fig. 3.24. Flow rates. Case #17. 35000000 30000000

Eclipse AHW

25000000 Cumulative, m33

Cumulative, m33

3500 Increase in reservoir permeability to 10 mD results in higher well productivity. The flow in the wellbore switches to turbuEclipse 3000 lent type, and the pressure gradient along the horizontal segAHW ment2500 of the well increases. Eclipse–100 automatically tracks the type of fluid flow in the well: as Reynolds number exceeds2000 4000, the friction factor is calculated using another formula 1500 (for details see Eclipse Technical Description Manual), while AHW model still uses the Poiseuille equation (??). Case

100 Time, days

20000000 15000000

1000

10000000

500

5000000

0

0

200

12

A.Cheremisin, A.Chesnokov, S.Golovin, D.Kuznetsov 10,0

500

9,0

450

8,0

Error, %

Flow w rate, m3/day

Flow Rate

7,0

Cumulative

6,0 5,0 4,0 3,0

400

Eclipse

350

AHW

300 250 200 150

2,0

100

1,0

50

0,0

0 0

50

100 Time, days

150

200

0


50

100 Time, days

150

200


3.3.1 Constant Permeability The reservoirs simulated here have constant, but different permeability in x, y, and z directions. It means that kx , ky , kz are constant values, though they are not equal. The results of Eclipse–100 runs with different number of cells in each coordinate axis direction were analyzed before comparison with AHW model. During these “adjustment” runs the optimal number of cells was selected: the minimal number of cells in each direction, when its further increment has negligible effect on result.

The behavior of production curves described above does Eclipse of cells in Eclipse or number of laynot depend upon number 60000 ers and series lengthsAHW in AHW. 50000 Cumulative, m33

Heterogeneous Reservoirs

10,0

40000 9,0

30000 8,0

Flow Rate

7,0

Cumulative

20000 6,0

Error, %

3.3

SPE number

10000 5,0 4,0

0

3,0

35

0

50

2,0

Flow w rate, m3/day

30 25

Eclipse

1,0

AHW

0,0

100

150

200

150

200

Time, days

20

0

100 Time, days

15


10

50

200

5

Eclipse 0

50

100 Time, days

150

Flow w rate, m3/day

0 200


Cases #19–22 see Appendix B are typical with respect to 3500 the variety of tested Eclipse in this section. The overall trend is as AHW follows: in Oy direction dominates perme3000when permeability ability in both Ox and Oz directions, then AHW model pro2500 duction curves lay slightly below the Eclipse–100 production curves2000 (Fig. 3.26–Fig. 3.28). Then, if kx values are dominant, the rates calculated by 1500 AHW model exceed those obtained via Eclipse–100 (Fig. 3.29).

150 AHW 100

50

Cumulative, m33

0 0

50

100 Time, days

150

200


1000

12000

0 0

50

100 Time, days

150

200

Cumulative, m33

500

3.3.2 Varying Permeability The heterogeneity of the Eclipse can be “magnified” in comparison reservoir in AHW model 10000 with section 3.3.1, meaning that the components permeability AHW tensor k , k , and k can be step functions of x variable (along x y z 8000 the main direction of the wellbore). New 6000 advanced model was tested on reservoirs with the following properties (see table on page 16). 4000 2000 0 0

50

100 Time, days

150

200

SPE number


10,0

10,0

9,0

9,0

8,0

Flow Rate

7,0

Cumulative

6,0

Error, %

Error, %

8,0

Flow Rate

7,0

13

5,0 4,0

Cumulative

6,0 5,0 4,0

3,0

3,0

2,0

2,0

1,0

1,0 0,0

0,0 0

50

100 Time, days

150

0

200

Cases #23–26 Appendix B cover several extreme scenarios of reservoir structure. It is known, that numerical algorithms may have difficulties (instabilities) when functions in the design domain are highly discrete. Special methods of calculations are required to overcome such problems. Listed cases were executed using several grids with different mesh size (for AHW the term “mesh of the grid” mean the number of layers along Ox bounded by planes x = const). Final time of simulation T = 200 days is enough to reach steadystate production period. Typical relative errors are shown on plots Fig. 3.30–Fig. 3.33. The agreement in flow rates and cumulative production is excellent except that AHW rates have oscillations at some time segments. Nevertheless, AHW is a good and fast tool for modeling non-homogeneous pore rocks. 10,0 9,0 Flow Rate

Error, %

7,0

150

200

Cumulative

6,0 5,0

3.4

Deviated Wellbore Survey One more option is available in AHW model: the wellbore can be positioned at any place of the reservoir, and its deviation survey can be not a straight line, but a step function determined at every layer separately. Thus, the well becomes a broken line with infinite conductivity connections at the boundaries of layers with x = const. The reservoir under consideration is still rectangular shaped with original dimensions: 1100 m along Ox, 100 m along Oy, and 10 m along Oz axis. The permeability is constant: 0.1 mD in all directions, fluid and wellbore properties are as described in section 3.1. The well starts at point x0 = 100 m and finishes at the point x1 = 1000 m. The fluid in the wellbore flows to the left, e.g. towards the starting point x0 . In case #27 (see Appendix B) the well changes its direction in Oxz plane. Vertical coordinate of the well center Zc 6= const and jumps from 0.5 m (distance to the bottom of the reservoir) to Zc = 5.5 , which is approximately the center of the reservoir. Next, in case #28 the wellbore lays in Oyz plane, and switches between 16.7 and 83.3 meters stand-off the side reservoir boundary y = 0. Case #29 is the superposition of two previous cases.

4,0

10,0

3,0

9,0

2,0

8,0

1,0

Flow Rate

0

50

100 Time, days

150

200


Error, %

7,0

0,0

Cumulative

6,0 5,0 4,0 3,0 2,0

10,0

1,0

9,0

0,0

8,0

Flow Rate

7,0 Error, %

100 Time, days



8,0

50

0

Cumulative

6,0

50

100 Time, days

150


5,0 4,0 3,0 2,0 1,0 0,0 0

50

100 Time, days

150


200

200

14

A.Cheremisin, A.Chesnokov, S.Golovin, D.Kuznetsov 10,0

5

SPE number

Conclusions and Notes

9,0 8,0

Flow Rate

Error, %

7,0

Cumulative

6,0 5,0 4,0 3,0 2,0 1,0 0,0 0

50

100 Time, days

150

200

Fig. 3.35. Relative error. Case #28. 10,0 9,0 8,0

Flow Rate

Error, %

7,0

Cumulative

6,0 5,0 4,0 3,0 2,0 1,0 0,0 0

50

100 Time, days

150

200


From the plots above one can see how wellbore placement influences production rates at transient and steady-state regimes. Relative error between AHW and Eclipse–100 is less than 1% in all considered cases.

4

Area of Applicability of Advanced Horizontal Well Model

Based on more than 50 comparative runs of AHW and Eclipse– 100, one can determine the boundaries of new model applicability.

4.1

Good Performance • Homogenous and heterogenous reservoirs with nonconstant permeability along horizontal well; • Reservoirs with multiple hydraulic fractures, each with its own physical and geometry properties; • Reservoirs with deviated horizontal wells, minimum recommended distance between well and outer boundary of the reservoir is 0.5 meters (2 feet); • Laminar flow in the wellbore.

4.2

Poor Performance • Reservoirs with high permeability and long horizontal wells cause in turbulent flow inside the wellbore. Then pressure distribution equation (1.5) needs additional adjustment via special type of Ψ–function, or the model should be equipped with turbulent flow regime in the wellbore (Fig. 3.24); • Long time simulation on depletion drive type reservoirs due to minor oscillations of the flow rate curve near zero values (Fig. 3.5);

The Advanced Horizontal Well Model is developed to capture objectives listed in Introduction. According to specific of the task, we used finite-difference approximation in Ox direction combined with Fourier series in both Oy, Oz directions. Also, as governing equations are linear with respect to time t, the Laplace transform technique was applied. A new original approach was used for analytical description of transverse hydraulic fractures, which number and properties can be different. The whole reservoir was divided into sections by planes, parallel to Oyz-plane. The model allows to describe cross-flows between these sections, and moreover, the cross-flows behind the fractures (as tips of fractures do not reach the boundary of the reservoir). Next complication was adding the Poiseuille equation to take into account the pressure drop along the horizontal segment of the wellbore. Together with non-trivial analytics, an original numerical algorithm was created. During debugging tests it has showed perfect stability, and due to this fact the whole AHW model can be evaluated as successful. The verification was started from very simple cases, and models of reservoirs were complicated with each passed test. Below there are brief comments on all new functionalities of the model. Varying Permeability In spite of the fact that proposed model is semi-analytical, the components of permeability tensor are not constants, though they depend only upon x variable, directed along the horizontal wellbore. Comparison tests with Eclipse–100 have showed coincidence, even when the permeability of adjacent sections differed by factor 10 (see section 3.3.2). Hydraulic Fractures Transverse hydraulic fractures can be placed within any section that contains well. All fractures can have non-symmetry wings, the conductivity of each wing may have its own value. Tests have showed that even in extreme conditions AHW gives results very close to Eclipse–100. There are no limitations in fracture conductivity, but there are some limitation of fracture length: the tip of the fracture is undesirable to be close to outer boundary (see Fig. 3.8 and description). Finite Wellbore Conductivity As a first step only laminar flow was considered, and equation (1.5) was used. When the flow rates become extremely high (it always happens with long wells in high-permeable reservoirs), then the type of flow in the wellbore switches from laminar to turbulent. The pressure gradient grows, and well productivity decreases in comparison with the one without friction inside the wellbore. Eclipse–100 automatic tracks the flow type (in reference to Reynolds number), while AHW model has no such option (see Fig. 3.24). It means that the end user must be careful with this option. Some adjustment can be manually done by a Ψ function selection (as it is realized now), but full control over the production flow will happen only after adding new turbulent flow regime into the AHW. Parallelepiped-type Shaped Reservoir A reservoir dimensions have three values: X, Y , and Z. Playing with different ratios between those parameters it was found that the algorithm and model itself produce stable results anywhere.

SPE number


Deviated Well The wellbore in AHW can have non-straight survey. Coordinates of the well center in every section containing well are specified independently from other sections. Thus, the wellbore becomes a broken line, which non-horizontal segments have infinite conductivity. The most complicated cases are when the well is positioned close to outer boundary, when the pressure gradient is high More horizontali.e. well drilled within several meters (or even centimeters). Cases #12–14 and Detailedeven Description Needed #27–29 show that if the wellbore “jumps” from section to section from top of the reservoir to the bottom of the reservoir Reservoir Complexity Well Properties Stimulation Events with only 0.5 m stand-off from the boundary, the disagreement between AHW Adequate and Eclipse–100 is only 2–3%. Well Production Nevertheless itCharacterization is not recommended to put the well closer than 0.5 m to the boundary, and this requirement also fits the physical sense.

Appendix A. Pressure at the Hydraulic Fracture Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Review of existing approaches and methods

The flow in i-th layer is approximated by one-dimensional flow along Ox-axis driven by the pressure drop from periphery to the center of the layer. In the steady state case the pressure Development of Math Model depends linearly on x, which leads to the formula (2.8). In Development of Numerical Algorithm the more realistic non-stationary flow regime the formula (2.9) Testing, Debugging, Verification is obtained. Note that in asymptotic t → ∞ this expression Software Development: Galveston Platform limits (2.8). Technical Assignment, Problem Statement, Possible Solutions

Reports

Linear approximation The pressure distribution in the Mar Apr May Jun Jul Aug Sep Oct Nov Dec fractured layer isFebsupposed to be linear (this corresponds to Phase Routine Phase Activity Stopped theActive exact solution for the 1D stationary fluid flow). Graph of the dependence p¯(x) is presented on Fig. 0.1. Average pressure p¯i in i-th layer is determined by averaging along x of the x-dependent linear function p¯lin .

pi-1 pi+1 pi pf

Area xi-1

xi

=

Area

xi+1

Fig. 0.1. Linear approximation of pressure in the layer with the hydrtaulic fracture

Analytically it looks like this:  f  p¯ − pi−1   (x − xi ) + p¯f , lin ∆ i−1 p¯ =  p¯ − pf   i+1 (x − xi ) + p¯f , ∆i

Averaging this function along i−th layer one obtain p¯f =

xi 6 x 6 xi+1 .

∆i−1 (4¯ pi − p¯i−1 ) + ∆i (4¯ pi − p¯i+1 ) . 3(∆i−1 + ∆i )

(A–1)

Time-dependent one-dimensional approximation Suppose that in the neighborhood of the hydraulic fracture the main flow is one-dimensional, directed along Ox-axis. Then the Laplace image of pressure satisfies the following equation −β ∗ s¯ p+

kx d2 p¯ = 0. µ dx2

(A–2)

The boundary conditions for the pressure in i-th and i + 1-st layers are p¯ = p¯i , p¯ = p¯f , x=xi−1

x=xi

and p¯ x=x = p¯f , i

p¯ x=x

i+1

= p¯i+1

correspondingly. The equation (A–2) with upper boundary conditions can be solved analytically. Again, averaging the solution with respect to x over the interval x ∈ (xi−1 , xi+1 ) we derive at: λi ∆i λi 3∆i λi h p¯f = ch + ch (∆i + ∆i−1 )¯ pi λi−1 × 4 4 4 3∆i−1 λi−1 ∆i−1 λi−1 ∆i−1 λi−1 i−1 + ch − 2¯ p sh − × ch 4 4 4 ∆i−1 λi−1 ∆i−1 λi−1 ∆i λi − p¯i+1 λi−1 ch ch sh × 4 2 4 ∆i λi ∆i λi 3∆i−1 λi−1 ch sh + × λi ch 4 2 4 −1 ∆i−1 λi−1 ∆i−1 λi−1 3∆i λi +λi−1 ch ch sh (A–3) 4 2 4 q q i−1/2 i+1/2 Here λi−1 = µβ ∗ s/kx , λi = µβ ∗ s/kx . Note, that in case t → ∞ (s → 0) and λi → 0, λi−1 → 0 an expression (A–3) coincide with (A–1). The simpler version of (A–3) can be obtained in assumption λi = λi−1 = λ, and ∆i−1 = ∆i = ∆: p¯f =

xi−1 6 x 6 xi ;

15

1 2sh 3λ∆ 4

λ∆ 3λ∆ p¯i λ∆ ch + ch − 4 4 −(¯ pi−1 + p¯i+1 )sh

λ∆ 4

16


SPE number

Appendix B. Cases # X 1 2 3 4 5 6 7 8 9 10 11

Reservoir, m Y Z

1100 1100 1100 1100 1100 1100 1100 1100 1100 1100 1100

100 100 101 101 101 101 101 101 101 101 101

X

10 10 11 11 11 11 11 11 11 11 11

Eclipse Cells Y Z

11 11 482 48 48 48 48 48 48 1504 150

3 3 3 3 3 3 3 3 3 3 3

1 1 1 1 1 1 1 1 1 1 1

Boundary Cond. Well Outer BHP BHP BHP BHP BHP BHP BHP BHP BHP BHP BHP

Fracture Y (m)/N

CP NF CP CP CP CP CP CP CP CP CP

Conductivity, mD · m

N N N 40m 90m 40m 40m 40m 20m3 N 3 × 40m5

1000 1000 10 1 3000 3000 3000

Well standoff6 Y = const Z = const 12 13 14

1100 1100 1100

15 16 17 18

1100 1100 1100 1100

101 101 101

100 100 100 100

11 11 11

11 11 11

10 10 10 10

101 101 101

11 11 11 11

3 3 3 3

100 11 100

1 1 1 1

BHP BHP BHP

BHP BHP BHP BHP

CP CP CP

50.5 0.5 0.5

0.5 5.5 0.5

Permeability, mD

Wellbore friction (Y/N)

0.01 10 10 000 10 000

Y Y Y N

CP CP CP CP

Permeability, mD kx ky kz 19 20 21 22

1100 1100 1100 1100

100 100 100 100

10 10 10 10

11 11 11 11

3 3 3 3

1 10 10 10

BHP BHP BHP BHP

CP CP CP CP

0.1 1 0.1 1

0.1 1 1 0.1

0.01 0.1 0.1 0.1

Y Y Y Y

BHP: Bottomhole pressure control; CP: Constant pressure; NF: No flow.

# kx

0–100 ky

kz

23 24 25 26 # 23 24 25 26

0.01 1 1 1

0.01 0.01 1 1 1 0.1 0.1 0.5 600–700 0.07 0.07 0.07 1 1 1 1 1 0.1 1 0.1 0.5

#

0–100 Yc Zc

27 28 29

N N N

N N N

kx

kz

0.02 0.02 0.1 0.1 0.1 1 1 5 700–800 0.08 0.08 0.08 0.1 0.1 0.1 0.1 0.1 1 0.1 1 5

0.5 5.5 0.5

200–300 Yc Zc 50 83.3 83.3

5.5 5.5 5.5

Permeability along Ox, m in mD 200–300 300–400 ky kz kx ky

kx

0.02 0.1 0.1 0.1

100–200 Yc Zc 50 16.7 16.7

100–200 ky

0.03 1 1 1

0.03 0.03 1 1 1 0.1 0.1 0.5 800–900 0.09 0.09 0.09 1 1 1 1 1 0.1 1 0.1 0.5

300–400 Yc Zc 50 16.7 16.7

0.5 5.5 0.5

0.04 0.1 0.1 0.1 0.1 0.1 0.1 0.1

0.04 0.1 0.1 1 900–1000 0.1 0.1 0.1 1

kz

kx

0.04 0.1 1 5

0.05 1 1 1

0.1 0.1 1 5

400–500 ky

0.05 1 1 0.1 1000–1100 0.11 0.11 1 1 1 1 1 0.1

Wellbore Coordinates in the Layer (Yc , Zc )7 , m 400–500 500–600 600–700 Yc Zc Yc Zc Yc Zc 50 83.3 83.3

5.5 5.5 5.5

50 16.7 16.7

0.5 5.5 0.5

50 83.3 83.3

5.5 5.5 5.5

kz

kx

500–600 ky

kz

0.05 1 0.1 0.5

0.06 0.1 0.1 0.1

0.06 0.1 0.1 1

0.06 0.1 1 5

0.11 1 0.1 0.5

700–800 Yc Zc 50 16.7 16.7

0.5 5.5 0.5

800–900 Yc Zc 50 83.3 83.3

5.5 5.5 5.5

900–1000 Yc Zc 50 16.7 16.7

0.5 5.5 0.5

1000-1100 Yc Zc N N N

N: No well in current layer.

References 1. Brusswell G., Banerjee R., Thambynayagam R.K.M., Spath J.: “Generalized Analytical Solution for Reservoir Problems With Multiple Wells and Boundary Conditions”, SPE 99288 2. Wan J., Aziz K. “Semi-Analytical Well Model of Horizontal Wells With Multiple Hydraulic Fractures” // SPE 81190 3. Zernar A., Bettam Y. “Interpretation of Multiple Hydraulically Fractured Horizontal Wells in Closed Sys-

2 Progressively

tems” // SPE 84888 4. Al-Kobaisi M., Ozkan E., Kazemi H. “A Hybrid Numerical/Analytical Model of a Finite-Conductivity Vertical Fracture Intercepted by a Horizopntal Well” // SPE 92040 5. Economides M., Deimbacher F., Brand C., Heinemann Z. “Comprehensive Simulation of Horizontal Well Performance” // SPE Formation Evaluation, December 1991 6. Krylov V.I., Skoblia N.S. “Handbook in numerical inversion of the Laplace trnsformation” Izdatelstvo “Nauka i tehnika”, Minsk, 1968.

refined towards the layer with fracture Fig. 3.9 fracture: only one wing l = 20m exist. 4 Progressively refined grid analogous to Fig. 3.9 for 3 separate fractures in different Ox-sections 5 3 fractures, each 40m of length (20m — half length) 6 In meters from the corresponding boundary 3 Non-symmetry

N N N