Determining multilayer formation properties from transient temperature and pressure measurements in gas wells with commingled zones

Journal of Natural Gas Science and Engineering 9 (2012) 60e72

Contents lists available at SciVerse ScienceDirect

Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse

Determining multilayer formation properties from transient temperature and pressure measurements in gas wells with commingled zones Weibo Sui a, *, Ding Zhu b a b

China University of Petroleum (Beijing), C. Ehlig-Economides, Beijing, China A.D. Hill, Texas A&M University, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 27 June 2011 Received in revised form 21 August 2011 Accepted 4 May 2012 Available online 29 June 2012

With the evolution of downhole permanent monitoring techniques, transient temperature and pressure data can play an important role in reservoir description due to their inherent real-time characteristics. Previous studies presented a completely new analysis technique for quantifying permeability and altered zone permeability and radius for multiple commingled layers. However, the previous model mainly applies for single-phase oil flow. A new wellbore/reservoir coupled flow model has been developed for multilayer commingled gas reservoirs including both damage and non-Darcy skin in each commingled layer. The non-Darcy effects are considered as permeability alteration and are incorporated to the reservoir flow model by using Forchheimer equation. Additionally, this coupled flow model can consider the pressure drop due to friction and kinetic energy changes in wellbore over the producing layers, which yields more accurate transient layer flow rate allocation. This coupled flow model is used to provide the wellbore pressure distribution and the radial reservoir pressure gradient for the coupled wellbore/reservoir temperature model. The temperature model is formulated using wellbore and reservoir energy balance equations considering subtle thermal factors such as JouleeThomson effect and also using fluid properties which are dependent on in-situ pressure and temperature. The inverse method is adopted from previous study directly and is used for determining formation properties by doing nonlinear regression. The mathematical model is solved numerically and used to study the sensitivity of transient temperature behavior to formation properties. The results show that transient temperature behavior in the wellbore at strategic locations is very sensitive to formation property values and has some interesting characteristics. However, due to the non-Darcy effects, each producing layer in multilayer gas reservoirs has non-Darcy skin more or less, which makes the transient temperature changes in gas reservoirs show more complex behavior. In the end, two hypothetical examples are presented to show the performance of the inverse method. The regression results show that the damage skin location and magnitude can be determined correctly using the proposed testing method. Ó 2012 Elsevier B.V. All rights reserved.

Keywords: Transient temperature analysis Temperature sensor Formation damage

1. Introduction Since the first generation fiber-optic sensing system was installed in Shell’s subsea wells successfully in early 1990’s, fiberoptic downhole sensing technology has made a significant progress in the past twenty years. With the development of technology and the reduction of installation and maintenance cost, downhole sensing technology is gaining increasing interest today in petroleum industry, and it brings evidently additional assets to reservoir and production management. Usually, downhole sensors include * Corresponding author. Tel.: þ86 13601332040. E-mail address: [email protected] (W. Sui). 1875-5100/$ e see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jngse.2012.05.010

pressure/temperature gauges, distributed or array temperature sensors (DTS or ATS), and multiphase flowmeters. Downhole pressure measurement is a relatively mature technology, which is proved by the well-established well testing theory and interpretation technology. While real-time downhole temperature measurement is a developing technology and attracts people’s interests due to its inherent characteristics. In the past five years, the industry reported various field applications of downhole temperature sensors from conventional flow profiling (Achnivu and Zhu, 2008), SAGD temperature observations (Wang, 2008) to stimulation treatments (Glasbergen et al., 2009; Huckabee, 2009). Meanwhile, researchers contributed many modeling works to assist various field applications of downhole temperature sensing,

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

and to discover some other potentials of temperature measurements. Kabir et al. (2008) improved Wang’s model (2008) to estimate total flow rate using DTS measurements. Huebsch et al. (2008) presented a developed thermal simulator to interpret inflow profile in velocitystring gas wells. Their models can be used for steady-state for the whole well depth or transient flow over the non reservoir intervals of the well. Li and Zhu, 2009 replenished Yoshioka’s model (2007) to a 3D, fully transient wellbore/reservoir model and investigated the water conning problem for infinite water-drive reservoir. In Li’s work, streamline simulation technique was used in forward modeling to speed up forward model simulation. Some other researchers focus on perfecting the wellbore/reservoir thermal model to simulate the reality and complexity of well completions (Almutairi and Davies, 2008; Muradov and Davies, 2009). All the models mentioned above are helpful for the industry to interpret downhole temperature measurements, and most of them addressed the usefulness of the distributed temperature data in flow profiling. Sui et al. (2008) a new testing approach and interpretation technique, which was aimed at determining multilayer formation properties (individual layer permeability, damage permeability, and damage radius). The working mechanism of the proposed testing approach is recording the transient wellbore temperature and pressure behavior at some strategic locations in wellbore during a transient flow test. Then the temperature and pressure measurement data are used for determining formation properties by using an inverse method. The general data acquisition configuration is shown in Fig. 1. In our previous work, a 2D (rez) single-well model was developed for multilayer single-phase oil reservoir, which is extended to gas reservoir case in this paper. Considering the great impact of non-Darcy effect on high-rate gas wells, we include the non-Darcy effect in this model and investigate the influence of non-Darcy effects on transient temperature behavior. The non-Darcy effect in gas reservoirs have been realized and studied by researchers for many years. Research results show that gas flow in the near well region cannot be adequately represented by Darcy law due to high flow velocities. In conventional well testing studies, the non-Darcy effect is usually treated as a flowrate-dependent skin factor (Smith, 1961; Swift and Kiel, 1962) that causes an instant pressure drop at the wellbore position. Since

61

the proposed testing method is aimed at determining damage radius and damage permeability instead of a single skin value, the simplified treatment of non-Darcy effect cannot meet our needs. Therefore, we applied the Forchheimer Equation (Forchheimer, 1901) in our reservoir flow model to account for the spatial distribution of non-Darcy effect. By using the new developed model, the transient temperature behavior in gas reservoirs is checked in this work, which shows a great potential for applying the transient temperature test in commingled gas reservoirs. 2. Methodology As an indirect testing method of reservoir characterization, we are facing an inverse problem similar to conventional well transient pressure test. The different point is that the formation properties cannot be obtained from analytical solutions explicitly due to the complexity of the heat transfer processes involved in this situation. To interpret multilayer formation properties from transient temperature and pressure measurements in gas reservoirs, a forward simulation model is established first to simulate the transient temperature and pressure behavior during transient flow tests with varying formation properties. With the established forward model, the inverse problem can be solved by using a nonlinear regression technique that could calculate the true formation properties from the observed temperature and pressure data by minimizing the objective function. For the forward simulation, a 2D (rez) single-phase gas wellbore/reservoir coupled model is established for this study. The forward model consists of four parts, i.e. reservoir flow model, reservoir thermal model, wellbore flow model, and wellbore thermal model. The reservoir flow model in previous study was an analytical model, which is replaced by a numerical 2D reservoir flow model with accounting near well non-Darcy effects. Other pieces of the model are adopted from previous study, thus only major equations are listed here. 2.1. Reservoir flow model The reservoir 2D flow model can be derived based on mass balance equation in porous media. The mass balance equation for single-phase gas flow in porous media is given by


 v
rg f ¼ V$ rg vg vt

(1)

From the definition of FVF (Formation Volume Factor) and Darcy’s law, the Darcy gas flow equation can be written as

" #   kg v f ¼ V$ Vp m g Bg vt Bg

(2)

In this study, the “effective permeability” technology (Pereira et al., 2006) is adopted to incorporate the near-well non-Darcy flow effect due to high flow rates. For single-phase gas flow, the original Forchheimer equation is given by

Vp ¼

mg kg

  vg þ bg rg vg vg 

(3)

By rearranging Eq. (3), the non-Darcy velocity can be represented by

v g ¼  1 þ bg rg Fig. 1. Data acquisition configuration (Sui et al., 2008).

kg 

mg

 vg 

!1 "

kg

mg

# Vp

(4)

62

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

The effective permeability, b k g , is defined by the following relationships,

b kg ¼

kg 1þa

(5)

kg 

mg

 vg 

6:15  1010 ; k1:55

(7)

(8)

 vp ¼ 0 vr r¼re

pjr¼re ¼ pe :

Re ¼

rvk0:5 ; m

the relationship between effective permeability and Reynolds number can be derived using Eqs. 6 through 10, and Eq. (5) becomes

b kg ¼

kg sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi: Re 0:5 þ 0:25 þ 1:03283  1013 1:05 kg

(11)

Eq. (11) shows a strong relationship between Reynolds number and effective permeability. k g yields the gas reservoir flow Replacing kg in Eq. (2) with b equation accounting for non-Darcy flow,

# "   b kg v f Vp : ¼ V$ mg Bg vt Bg



! 2p b k g dz vp : r mg Bg vr

(18)

Using control volume finite difference (CVFD) method, Eq. (12) can be discretized into radial logarithmic grid blocks and solved implicitly. The logarithmic correlation coefficient alg is determined by the reservoir outer radius re reservoir inner radius rw, and the number of grid blocks NR in the radial direction.

 alg ¼

re rw



1 NR1

(19)

At the same time, the interblock transmissibilities can be determined by considering pressure drop relationships for steadystate flow between neighbor grid blocks. For cylindrical coordinate system, the transmissibilities are defined as follows (Pedrosa and Aziz, 1986)

Tri1=2;j

2pDzj krg ¼ alg lnalg alg  1 mg Bg ln ln alg  1 lnalg þ n n b b k ri1;j k ri;j

! (20) ri1=2;j

(12)

For a cylindrical coordinate system, initial and boundary conditions of the reservoir flow model are given as follows. The initial pressure distribution is set to be the initial reservoir pressure distribution. Different initial layer pressure distribution can be considered.

pðr; z; tÞjt¼0 ¼ pri ðr; zÞ;

Zmax X z ¼ Zmin

(9)

(10)

(17)

The inner boundary condition is given by the constant surface flow rate,

qg;sc ¼

Since the Reynolds number is usually defined as the function of permeability for flow in porous media,

(16)

or

where b is in 1/ft, k is in md. For SI units, Eq. (8) becomes

1:03283  1013 b¼ : k1:55

(15)

where Zmin and Zmax are the depth of upper and lower boundaries of the model computation region. The outer boundary can be either no-flow boundary or constant-pressure boundary,

where b is the non-Darcy coefficient, which is evaluated by using Jones’ correlation (1987) in this study. Although it is very hopeful to evaluate the non-Darcy coefficient together with the other formation property parameters using the proposed testing approach, this study didn’t go that far away at this stage. Here b is evaluated by using the correlation by Jones (1987),

b¼

and

(6)

and

b ¼ bg rg a

(14)

 vp ¼ 0; vz z¼Zmax

where

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b 1 1 þ 4a a ¼ 2

 vp ¼ 0 vz z¼Zmin

Triþ1=2;j

2pDzj krg ¼ alg  1 alg lnalg mg Bg ln ln alg  1 lnalg þ n n b b k ri;j k riþ1;j

(13)

where pri is the initial reservoir pressure. Since the computation region is much larger than the producing zones, the upper and lower boundaries of the model are set to be sealed to flow,

Tzi;j1=2 ¼


  2 Vbi;j Dzj

Dzj1=2 n b k zi;j

þ

Dzj1=2 n b k zi;j1

krg mg Bg

! (21) riþ1=2;j

! (22) riþ1=2;j

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

63

Therefore, the discretized flow equation is given by

0 ðvÞ nþ1

Tzi;j1=2

¼

ðvþ1Þ B nþ1 @pi;j1



f Dt Bg

Vbi;j

0 i;j

1

0

ðvþ1Þ ðvÞ ðvþ1Þ nþ1 C nþ1 B nþ1  pi;j A þ Tri1=2;j @pi1;j

2 ðvþ1Þ 6 nþ1 4pi;j

1

0

ðvþ1Þ ðvÞ ðvþ1Þ nþ1 C nþ1 B nþ1  pi;j A þ Triþ1=2;j @piþ1;j

1

0

1

ðvþ1Þ ðvþ1Þ ðvÞ nþ1 C nþ1 B nþ1  pi;j A þ Tzi;jþ1=2 @pi;jþ1

ðvþ1Þ nþ1 C  pi;j A

3

(23)

7  pni;j 5;

where Vbi;j is the volume of grid block (i, j), the subscripts n and v represent time step and iteration step respectively. During the simulation, the transmissibility coefficients are linearized by simple iteration method. The accumulation term is evaluated by the chord slope method (Abou-Kassem et al., 2006). Another complexity of solving the reservoir flow model comes from the flow rate allocation among multiple producing zones. According to common reservoir simulation strategy, the layer flow rate allocation in the commingled multilayer reservoirs is determined by layer flow capacities,

krg h j qg;scj ¼ PN

qg;sc p krg h j j¼1

error iteration method to achieve the correct flow rate allocation in this study. The relative error is controlled within 1% to satisfy our accuracy requirement. Since the wellbore pressure opposite to individual production blocs is given by


 pwfj ¼ pwfref þ gwb Zj  Zref

where gwb is the wellbore fluid gravity, Zj is the depth of the jth grid blocks, Zref is the reference grid block depth, pwfref is the wellbore pressure at reference depth. Thus the production rate at the jth grid block can be calculated by

(24) qg;scj

where Np is the number of production grid blocks, and (krg)j. However, this approach cannot satisfy our requirement in this study since it cannot take into account the hydrostatic pressure difference between layers. Furthermore, it doesn’t consider the influence of near wellbore damaged permeabilities on early-time transient layer flow rates. Therefore, we implement the try-and-

(25)

Gw ¼  mg Bg

!




pj  pwfj



j


 2pDzj 1  pj  pwfj
¼  alg  1 alg lnalg m B g g ln ln j lnalg alg  1 þ n n b b k r1;j k r2;j

Fig. 2. Solution diagram of the forward model.

(26)

64

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

Fig. 3. Flowing bottomhole pressure calculated by Saphir software and our model. n where Gwj is the geometry factor and b k r1;j is the non-Darcy permeability of the grid block next to the wellbore. Additionally, the sum of individual layer flow rates must be equal to the total flow rates,

qg;sc ¼

Np X

qg;scj

(27)

j¼1

The wellbore pressure at reference depth can be determined by combining Eqs. (26) and (27),

PNp j¼1

pwfref ¼

(

Gw mg Bg

!

h




pj  gwb Zj  Zref

j

PNp j¼1

Gw mg Bg

i

) þ qg;sc

!

(28)

j

Therefore, Eq. (28) is used to estimate pwfref using the pressure value at old time level n, and then calculate qg,scj at the time level of n þ 1, iteration level of v ¼ 0. If the calculated and given qg,sc values are not consistent with each other, the pj is renewed and the same procedure is repeated for the iteration level of v ¼ 1 until the discrepancy between calculated and given qg,sc within 1% accuracy requirement. 2.2. Reservoir thermal model The reservoir thermal model is derived from the energy balance equation, and has been described in great details in our previous work (Sui et al., 2008), the governing equation is given as

r Cb p






vT vp b p VT þ bT T  1 vg $Vp þ K V2 T;  fbg T ¼ rg vg $ C T g g vt vt (29)

Fig. 4. Layer flow rates calculated by Saphir software and our model.

Fig. 5. Gas JeT coefficient variations with well depth.

b p ¼ fr C b b where r C g pg þ ð1  fÞrr C pr is the average formation property (density and specific heat capacity) of gas and formation T bp rock, bg is the gas thermal compressibility, rg is the gas density, C g is the gas specific heat capacity, KT is the heat conductivity coefficient of the formation rock saturated with gas. The initial reservoir temperature is set to be geothermal temperature. The upper, lower, and outer boundaries of the computation region are set to be constant geothermal temperature. At the inner boundary, radiation boundary condition is used to describe the heat exchange between the wellbore and the formation. The reservoir thermal model is also discretized using radial logarithmic grid blocks consistent with the reservoir flow model. Due to the strong compressibility of gas, the gas properties in the reservoir thermal model usually cause strong nonlinearities. To solve this problem, the gas properties in the coefficients are linearized using simple iteration method, and the reservoir thermal model is solved implicitly. 2.3. Wellbore flow model Wellbore flow model consists of mass balance equation and momentum balance equation (Yoshioka, 2007). The mass balance equation is written as

 2grI;g vI;g d
r vg ¼ ; dz g R

(30)

where R is the pipe inner radius, rI,g and vI,g are the density and volumetric velocity of gas entering into wellbore from formation, g Table 1 Input reservoir and fluid properties. Initial reservoir pressure at reference depth, MPa Initial reservoir temperature at reference depth,  C Surface temperature,  C

25.0/40.8

Wellbore radius, m

0.08

60.9/78.9

Porosity, fraction

0.15

15

1675

Geothermal gradient,  C/m Reservoir pressure gradient, kPa/m Reference depth, m

0.018 10.18

Drainage radius, m

1500

Gas specific heat capacity at standard Condition, J/kg- C Gas heat conductivity coefficient, W/m- C Gas viscosity at standard condition, mPas Gas gravity at standard condition

JouleeThomson coefficient of gas at reference depth,  C/MPa

1.7/0.8

2550/3550

0.035 1.08e-2

0.6

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

65

Table 2 Layer properties for studying non-Darcy effects.

Layer 1 Layer 2

Depth (m)

k (103 mm2)

ks (103 mm2)

rs (m)

s

3510e3525 3535e3550

3 3

0.41 e

0.41 e

10 0

is the open ratio of the pipe. The momentum balance equation is given by


rg v2g f d rg v2g dp ¼   rg g;  dz dx R

(31)

where f is the friction factor. Fig. 7. Reduction percentage of bottomhole pressure and temperature (at 3530 m) for different surface flow rates (5  104 m3/d w 50  104 m3/d).

2.4. Wellbore thermal model Wellbore thermal model is derived from the energy balance equation (Sui, 2009),


  vT 2g rI;g vI;g
2ð1  gÞ Tr jr¼rw  T þ UT Tr jr¼rw  T ¼ b p;g vz R rg vg Rrg vg C vp g þ KJT;g ;  b p;g vz C

(32)

where Trjr ¼ rw is the formation temperature at right outside of wellbore, KJT,g is the JouleeThomson coefficient of gas. It should be noted that the wellbore thermal model used here is a transient model, and the wellbore flow model is a sequential steady-state model, which will be updated at every each time step after solving the reservoir counterpart layer flow rates. 2.5. Solution procedure of the forward model Due to the different transfer mechanisms of pressure and temperature, the wellbore/reservoir flow model is solved decoupled from the wellbore/reservoir thermal model, which could avoid solving a giant sparse matrix. The solution diagram of the forward model is shown in Fig. 2. First, the wellbore/reservoir flow model is solved after initialization for obtaining fluid velocity and pressure distribution inside wellbore and reservoir, which will be provided to the wellbore/reservoir thermal model. In the coupled flow model, there is only one set of grid block system used for wellbore and reservoir, the inner boundary points of the reservoir are used to

represent the grid blocks of wellbore. It should be noted that such discretization methodology is only applicable for 2D cylindrical single-well model (Abou-Kassem et al., 2006). From Fig. 2, we can see that there are three major iterative loops for solving the forward model, which is required by the following situations, a) individual layer flow rates cannot be assigned explicitly. Iteration procedure is used to obtain accurate inflow profiles; b) ‘effective permeability’ b k g is a function of fluid velocity, which must be solved iteratively with the reservoir flow model; c) Due to the strong nonlinearity of all the governing equations in forward model, fluid properties must be updated using the newest pressure and temperature fields. The iterative loops are controlled by the defined upper bounds of relative errors. The relative error of surface flow rate is defined as

Eerr;q ¼

    qg;sc  q+g;sc  q+g;sc

and 3 q ¼ 0:01;

(33)

where q+g;sc is the real surface flow rate. The relative errors of b k g and wellbore fluid temperature are defined by

Eerr;k ¼

Eerr;T ¼

  b  k g;new   k g;old  b b k g;new

and 3 k ¼ 0:01;

(34)

jTold  Tnew j and 3 T ¼ 0:01: Tnew

(35)

where b k g;old and Told are the effective permeability and wellbore fluid temperature at last iteration step, b k g;new and Tnew are the effective permeability and wellbore fluid temperature at current iteration step. In the end, the coupled wellbore/reservoir flow model has been validated using the gas reservoir multilayer testing model in

Table 3 Layer properties for studying the effect of damaged radius.

Fig. 6. Reduction ratio of permeability in reservoir (Layer 2) for different surface flow rates (5  104 m3/d w 50  104 m3/d).

Case 1

Depth (m)

k (103 mm2)

ks (103 mm2)

ds (m)

s

Layer 1

2510e2525

Layer 2

2535e2550

3 3 3 3 3

0.09 0.22 0.39 0.52 e

0.03 0.10 0.30 0.60 e

10 10 10 10 0

66

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

Fig. 8. Bottomhole pressure change and pressure change derivatives for different damage radii in Case 1.

a commercial well testing software Saphir (KAPPA, 2010). The transient bottomhole pressure behavior and layer flow rate variations from our model are consistent with the simulation results given by the software. The flowing bottomhole pressure and layer flow rate simulation results are shown in Figs. 3 and 4. The wellbore/reservoir thermal model cannot be validated at this time due to the lack of published model and commercial software with this function. 3. Result and discussion With the developed mathematical model, we studied the characteristics of transient temperature behavior during the transient test with varying formation property values in gas wells with commingled zones. The study results show that transient temperature behavior is very sensitive to different formation property values although near-wellbore non-Darcy effects add more complications into interpretation. In the end, such testing and analysis techniques are applied to two hypothetical examples to illustrate the model performance in gas reservoirs.

Fig. 10. Wellbore fluid temperature change derivatives at depth of 2505 m for different damage radii in Case 1.

studying transient temperature behavior in gas wells, it should be clearly realized that gas JouleeThomson effect has some different features for different production times and well production depth. In recent published research works, the JouleeThomson effect is considered as the major reason of producing fluid temperature deviation from original geothermal temperature. Since it is common to observe a large temperature drop opposite to the production zone from the wellbore temperature distribution, gas JouleeThomson (JeT) coefficient magnitude is believed to be much higher than oil and water case. However, such observations usually come from long-time production wells, where the formation has been cooled down by previous gas flow, and JeT effects can show maximum cooling. Under such condition, Steffensen and Smith presented the general range of JeT coefficient for natural gas is 2.42e4.83  C/MPa (0.03e0.06  F/psi). On the other hand, it should be noted that gas JeT coefficient is strongly dependent on pressure and temperature. Thus it is a function of well depth. The gas JeT coefficient magnitude can be very different for a shallow well and a deep well. From the thermodynamic definition of JeT coefficient (Bird et al., 2002),

3.1. Gas JouleeThomson coefficient Before we delve into the specific transient temperature behavior in gas wells, we would provide some background information and considerations for Gas JouleeThomson coefficient, which has a great influence on the investigated topic. Since this paper is

Fig. 9. Wellbore fluid temperature change at depth of 2505 m for different damage radii in Case 1.

KJT;g ¼

bg ðp; TÞ$T  1 rg ðp; TÞ$Cp;g ðp; TÞ

Fig. 11. Permeability fields in Layer 1 (ds ¼ 0.30 m, s ¼ 10).

(36)

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

67

Table 4 Layer properties for studying the effect of skin location.

Case 2a Layer 1 Layer 2 Case 2b Layer 1 Layer 2

Fig. 12. Permeability fields in Layer 2 (s ¼ 0).

where the gas thermal expansion coefficient, gas density and gas heat capacity are considered as the function of pressure and temperature, which can be calculated by published equations or correlations (John and Wattenbarger, 1996). We generated the gas JeT coefficient (CH4 component) variations with well depth, which is shown in Fig. 5. From Fig. 5, we can see that gas JeT coefficient at depth of 1000 m can be almost 10 times higher than that at depth of 4000 m, where higher pressure and temperature make gas properties close to liquid. In this paper, the gas well producing intervals are between 2510 m and 2550 m. Similar considerations for JouleeThomson coefficient in gas wells were also presented by Wang et al. (2008) who showed the smaller or even reverse temperature anomalies at producing zones in deeper gas wells. 3.2. Influence of non-Darcy flow on transient pressure/temperature behavior Before investigating characteristics of transient temperature behavior, we would study the influence of non-Darcy flow during the test first. The input reservoir and fluid properties are shown in Table 1 and 2, which will be used for all the case studies in this work unless otherwise specified. Considering the surface gas flow rate varies in the range of 5  104 m3/d w50  104 m3/d, and the gas production well is

Fig. 13. Temperature change derivatives at depth of 2530 m for different damage radii in Case 1.

Depth (m)

k (103 mm2)

ks (103 mm2)

ds (m)

s

2510e2525 2535e2550

3 3

e 0.39

e 0.30

0 10

2510e2525 2535e2550

3 3

0.39 e

0.30 e

10 0

operated at such constant rates for 1000 h, the reduction ratio of effective permeability in the non-damage layer is calculated using our model by considering and without considering non-Darcy effects. By observing Fig. 6, we can see that the results indicate great effective permeability deductions with high flow rate nonDarcy effects, which is consistent with previous study results (Fligelman et al., 1989). The reduction percentage of bottomhole pressure and temperature (at 3530 m) are calculated and shown in Fig. 7. The reduction percentages are defined as the proportion between the extra bottomhole pressure and temperature changes caused by non-Darcy effects and the total bottomhole pressure and temperature changes. Fig. 7 shows that such reduction percentage increases with surface production rate, and it also 18.7% and 32.7% reduction of pressure and temperature can happen for a high-rate gas well. Thus, non-Darcy effects should be taken into account in this study. 3.3. Characteristics of transient temperature behavior in gas reservoirs From our previous research work, we found that transient temperature behavior showed strong sensitivities to different damage radius, skin factor, and permeability in oil reservoirs. Such sensitivities made it possible to inverse the true values of formation properties from transient temperature and pressure measurement. We believe similar characteristics also exist in gas reservoirs. In this part, we will investigate transient temperature behavior during the transient flow test with varying formation properties. 3.3.1. Case1: sensitivity of the transient temperature to the damage radius A hypothetical two-layer gas reservoir is used for studying the sensitivity of the transient temperature to the damage depth. The detailed reservoir information is listed in Table 3. Four different damage radii are used to represent different formation damage

Fig. 14. Temperature behavior during test in Case 2a (z ¼ 2505 m and 2530 m).

68

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

Fig. 15. Temperature behavior during test in Case 2b (z ¼ 2505 m and 2530 m).

scenarios from shallow to deep. The parameter ds is the damage depth. Assume the gas well has been producing at a rate of 20  104 m3/d for a while and then the surface production rate is cut back to be 10  104 m3/d. According to the data acquisition configuration we showed in Fig. 1, during the transient flow test, downhole pressure at 2505 m and 2530 m should be recorded for analysis. For this hypothetical case, the ‘temperature and pressure measurement data’ are generated using the forward model by simulating the transient flow test. First, the bottomhole pressure changes and their derivatives for four different scenarios are shown in Fig. 8, we can see that the pressure behavior do not tell differences between different scenarios. Such dilemma are alleviated by exploring temperature changes shown in Fig. 9, which shows transient temperature behavior at 2505 m are sensitive to different damage radii. The transient temperature sensitivity is also magnified by generating temperature change derivative curves ðjdTjÞ=ðdlnðDTÞÞ. From Fig. 10, we can see that temperature change derivative curves clearly show similar humps but with different peak time for different damage scenarios, which is a significant fingerprint for transient temperature diagnostics. The temperature derivative humps are caused by the varying speed of temperature changes, which is the outcome of permeability variations in near-wellbore region. The peak time of temperature derivative humps mainly depends on the damaged depth, i.e. the smaller the damaged depth,

Fig. 17. Temperature change derivative behavior during test in Case 2b (z ¼ 2505 m and 2530 m).

the earlier the peak appears. However, non-Darcy effects in gas reservoirs do impact the transient temperature behaviors. In gas reservoirs, besides the permeability reduction due to the formation damage, permeability is also reduced by the non-Darcy flow in near well region. Therefore, there is ‘non-Darcy skin’ in both damaged layer and undamaged layer, which makes transient temperature behavior in gas reservoirs more complex than oil cases. For the case of ds ¼ 0.30 m, the reservoir permeability fields in both layers are shown in Figs. 11 and 12. From these figures, we can see clearly that there is a permeability altered zone caused by non-Darcy flow in both layers. Due to zero skin and a higher layer rate, the non-Darcy permeability field even goes deeper in Layer 2. The temperature change derivative curves at 2530 m is generated and presented in Fig. 13, from which shows that the temperature change derivative curves for four scenarios show similar humps with almost same peak time, because Layer 2 produces same flow rates in those damaged scenarios which yields same non-Darcy permeability fields. 3.3.2. Case 2: sensitivity of the transient temperature to skin locations In Case 1, we have studied the sensitivity of transient temperature to the damage depth, and the skin factor is located in Layer 1. In this case, we will study how the skin location affects transient temperature performance. The damage skin location in multilayer reservoirs is valuable information for stimulation work design and usually cannot be interpreted from conventional transient pressure analysis workflow. The same reservoir diagram and production history as Case 1 are used for this study. The detailed formation properties are listed in Table 4. The formation damage is located in the bottom layer and the upper layer in Case 2a and 2b respectively. For this case, we calculated the temperature changes and their derivatives at different measured locations (2505 m and 2530 m) for Case 2a and 2b, which are shown in Figs. 14 through 17. The key point underlying in these figures is that the temperature change

Table 5 Layer properties for studying the effect of permeability.

Fig. 16. Temperature change derivative behavior during test in Case 2a (z ¼ 2505 m and 2530 m).

Depth (m)

k (103 mm2)

s

Layer 1 Layer 2

3490e3510 3520e3524

Layer 3

3530e3550

6 (a) 0.1 (b) 10 (c) 50 3

0 0 0 0 0

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

69

Fig. 20. Wellbore fluid temperature history during test in Case 3b

3.3.3. Case 3: sensitivity of transient temperature to layer permeabilities The preceding cases addressed the sensitivity of transient temperature to the damage radius and skin locations. Actually,

transient wellbore fluid temperature behavior is also sensitive to contrasting layer permeabilities. In this part, we will use a hypothetical three-layer gas reservoir to study the sensitivity of transient temperature to layer permeabilities. In this reservoir, the upper and bottom layer are set to be both 20 m thick with permeability of 6 md and 3 md respectively. The middle layer is a 4 m-thickness zone with either relatively very high or very low permeabilities. The formation damage skin factor is set to be zero for all layers to eliminate the skin effects on this study. The detailed layer properties are listed in Table 5. Assume the gas well has been producing at a rate of 40  104 m3/d for a while and then the surface production rate is cut back to be 20  104 m3/d. To study the transient temperature variation, the wellbore fluid temperature changes above each producing layer are recorded during the transient test. The specific depths of measurement are 3486 m, 3516 m, and 3528 m. From Case 1 and 2, we know that the near well permeability altered zone usually causes a greater pressure gradient and further lead to larger temperature changes. Since we do not have any damage skin effect in this case, the final wellbore mixing fluid temperature will mainly depend on the layer flow rate distribution. Using different permeability values in Layer 2, we calculate the flow profile at the end of test and show the results in Fig. 18. With permeability ranging between 0.1 md and 50 md, the layer flow rate percentage of Layer 2 varies from less than 1% to almost 50%. Such contrasting flow profiles yield completely different transient

Fig. 19. Wellbore fluid temperature history during test in Case 3a.

Fig. 21. Wellbore fluid temperature history during test in Case 3c

Fig. 18. Flow profiles at the end of test in Case 3a, 3b, and 3c.

derivative curve right above the damaged layer showing a higher hump compared to the derivative curves at other locations, for example the ‘2530 m’ curve and the ‘2505 m’ curve in Figs. 16 and 17 respectively. This is a very helpful hint for us to determine the damage skin location in practical interpretation. The derivative humps at other locations are caused by different reasons. When skin factor is located in Layer 2 (Fig. 16), the temperature change derivative curves above both producing layers both can see the damage skin effect by showing humps; when skin factor is located in Layer 1 (Fig. 17), only the temperature change derivative curve above Layer 1 can see the damage skin effect. The hump at ‘2530 m’ curve is caused by the non-Darcy permeability alterations in Layer 2. Another characteristic shown in the temperature change derivative curves is the final curve drop, which means the temperature is achieving steady-state status gradually.

70

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

Table 6 Regression results for a two-layer case. True values

Layer 1 Layer 2

Initial guesses

Regression results

k (103 mm2)

ks (103 mm2)

rs (m)

s

k (103 mm2)

ks (103 mm2)

rs (m)

s

k (103 mm2)

ks (103 mm2)

rs (m)

s

6.0 3.0

e 0.37

e 0.40

0 10

5.0 5.0

0.83 0.83

0.2 0.2

5.0 5.0

5.98 2.99

5.98 0.40

0.64 0.4

0.0001 9.98

Table 7 Regression results for a three-layer case. True values k (10 Layer 1 Layer 2 Layer 3

1.0 1.0 10.0

3

mm ) 2

Initial guesses ks (10 e e 0.87

3

mm ) 2

3

rs (m)

s

k (10

e e 0.2

0 0 10

4.0 4.0 4.0

mm ) 2

Regression results 3

ks (10 0.67 0.67 0.67

wellbore fluid temperature behavior during test. The temperature histories for three cases are shown in Figs. 19 through 21. If we take Case 3a as an example, we can see that the temperature curve above Layer 2 (3516 m) is almost identical to the temperature curve above Layer 3 (3528 m), which is because there is only little formation fluids coming from Layer 2, and such small quantity of fluids does not affect wellbore fluid temperature prominently. On the contrary, the big flow rate contribution from Layer 2 in Case 3c make the temperature curve at 3516 m much different from the temperature curve at 3528 m. While the situation of Case 3b lies between Case 3a and 3c. 3.4. Regression results of the inverse problem After investigating most characteristics of transient temperature behavior, the inverse problem will be solved by using nonlinear regression method. The objective function f(x) in Eq. (27) is defined to describe the discrepancy between observed data d and simulated data g(x).

fðxÞ ¼

1 1 kd  gðxÞk22 ¼ ðd  gðxÞÞT C1 n ðd  gðxÞÞ 2 2

(37)

By minimizing the objective function, the true values of formation properties represented by vector x can be obtained. The regression vector x is defined as

x ¼ ½k1 ; k2 ; /; kN ; ks1 ; ks2 ; /; ksN ; rs1 ; rs2 ; /; rsN T3N1

Fig. 22. Bottomhole pressure data with noise in the 3-layer inverse case.

(38)

mm ) 2

rs (m)

s

k (103 mm2)

ks (103 mm2)

rs (m)

s

0.15 0.15 0.15

3 3 3

0.98 0.72 8.12

0.96 0.26 0.97

0.4 0.09 (¼rw) 0.4

0.02 0 10.3

where N is the number of producing layers. During the regression procedure, the LevenburgeMarqurdt algorithm (C.T. Kelly, 1999; Weibo Sui, 2009) is applied to renew the iteration parameters. A hypothetical two-layer and a three-layer example are presented here to illustrate the performance of the inverse method. The objective of the interpretation is mainly to find out the specific location and magnitude of the formation damage skin factor. First, the proposed forward model is used to generate the ‘observed temperature and pressure data’. Using the generated observation data and the initial guesses of the regression parameters (k, ks, and rs), we could calculate the true values of formation properties by doing nonlinear regression. The detailed interpretation procedure can be referred to our previous work (Sui et al., 2008). The true values, initial guesses, and final regression results for both examples are shown in Tables 6 and 7 respectively. From Table 6, we can see the regression results of the two-layer case match the true values of the formation properties very well. The magnitude and location of the damage skin factor are interpreted accurately. While from Table 7, we can see that the regression results of the threelayer case do not perfectly match the true values of regression parameters. The imperfectness is mainly caused by the layer permeability alterations caused by non-Darcy effects. Considering the practical complications caused by data noise, we regenerated the ‘observed temperature and pressure data’ for the 3layer case by adding random data noises. The temperature data noise amplitude is 0.1  C, and the pressure data noise amplitude is set to be 50Pa. The temperature and pressure data with noise are shown in Figs. 22 and 23. With the same initial guesses, the inverse

Fig. 23. Downhole Temperature data with noise in the 3-layer inverse case.

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

71

Table 8 Regression results for a three-layer case with noise data. True values

Layer 1 Layer 2 Layer 3

Initial guesses

Regression results

k (103 mm2)

ks (103 mm2)

rs (m)

s

k (103 mm2)

ks (103 mm2)

rs (m)

s

k (103 mm2)

ks (103 mm2)

rs (m)

s

1.0 1.0 10.0

e e 0.87

e e 0.2

0 0 10

4.0 4.0 4.0

0.67 0.67 0.67

0.15 0.15 0.15

3 3 3

1.17 0.88 12.15

0.59 0.38 1.32

0.12 0.33 0.36

0.37 1.81 12.10

annular space with constant gas production Q0 ¼ 470  103 m3/d or 5.52 m3/s at standard condition. Observation of the temperature change was during 2 h. During this time pressure decreased by 14.7 MPa, temperature decreased by 25  C. With provided limited formation and well production information, we simulated the transient testing procedure and the matching results are shown in Fig. 24. The reservoir permeability was determined in this procedure to be 3.7 md. At the same time, the transient temperature change derivative plot is generated and shown in Fig. 25. From Fig. 25, we can see the distinct hump characteristic for near-wellbore permeability alteration. With considering the high-rate non-Darcy effect, we used our inverse simulator to find out the formation damage scenario. The regression result is s ¼ 4 and rs ¼ 0.32 m. 5. Conclusions Fig. 24. Downhole pressure and temperature data with simulation results.

results are shown in Table 8 where we can see that the regression accuracy decreases a little bit, but still yield very valuable information. More detailed works on sensor resolution and data noise impacts can be found in our previous work (Sui et al., 2008). 4. Field example Although downhole permanent temperature and pressure gauges are applied widely in recent years, due to the lack of effective interpretation techniques, transient temperature field data were rarely published especially for multilayer reservoirs. Due to the failure of obtaining such field data, here what we only have is the transient bottomhole temperature data in a single-layer gas well on Shebelinskoe field (Chekalyuk, 1965). We would like to use this field example to prove our findings and interpretation approach in some extent. The completion section in this gas well is 1476 me1499 m, wellbore diameter on bit is 25 cm, 2.500 flowing tubes was pulled down to the depth of 1495 m. The gas well operated through the

In this paper, a 2D numerical wellbore/reservoir coupled model is established, which extends the previous proposed multilayer testing method for single-phase oil reservoir to gas reservoirs successfully. The near well non-Darcy flow caused by high flow rate is taken into account by using Forchheimer equation in reservoir flow model, and the non-Darcy effects on transient temperature and pressure have been studied. The study results show that the non-Darcy effects become greater in high rate gas reservoirs and it mainly happens within the near well region. In this study, the nonDarcy effects are considered as the permeability alteration. With the developed model, some interesting characteristics of transient temperature behavior in gas reservoirs have been investigated. Research results show that the transient temperature behavior in gas reservoirs has sensitivities to variations of the damage radius, skin location, and layer permeability. However, due to the non-Darcy effects in gas reservoirs, the temperature performance in gas reservoirs does show more complexities than that in oil reservoirs. For synthetic illustration, a two-layer and a threelayer case are presented to show the performance of the inversion method with considering data noise impacts. From the calculation results, we can see that the skin location and the magnitude can be interpreted successfully for both cases. In the end, a singlelayer field example was shown to prove the characteristic of transient temperature behavior and calculate the damage scenarios. Nomenclature B bp C

Fig. 25. Simulated temperature change derivatives in field case.

Cn d Eerr,k Eerr,q Eerr,T f f(x) g g(x)

formation volume factor, m3/m3 specific heat capacity, J/(kgK) covariance matrix observation data relative error of effective permeability relative error of surface flow rate relative error of temperature friction factor objective function gravity acceleration, m/s2 simulation data

72

K KJT k b k ks p pri qsc R r re rs rwb s T TI t Vb v x z Zmin Zmax

a b a b bT 3

g f m r

W. Sui, D. Zhu / Journal of Natural Gas Science and Engineering 9 (2012) 60e72

thermal conductivity, W/m2K JouleeThomson coefficient, K/Pa permeability, md, m2 effective permeability due to non-Darcy flow, md, m2 damaged permeability, md, m2 pressure, Pa initial reservoir pressure, Pa production rate, m3/d wellbore radius, m radial coordinate outer boundary of the computation region, m damaged radius, m wellbore radius, m skin factor temperature, K; transmissibility coefficient inflow temperature, K time, s volume of the grid block, m3 velocity vector regression vector vertical coordinate the depth of the upper boundary of the computation region, m the depth of the lower boundary of the computation region, m intermediate parameter intermediate parameter non-Darcy coefficient, m1 thermal expansion factor, K1 defined upper bound of the relative error pipe open ratio porosity viscosity, Pas density, kg/m3

Superscripts v the vth iteration step at current time step n the nth time step

Subscripts g gas I inflow i position index j layer index N number of producing layers ND non-Darcy r radial direction z vertical direction References Abou-Kassem, J.H., Farouq Ali, S.M., Islam, M.R., 2006. Petroleum Reservoir Simulation-A Basic Approach. Gulf Publishing Company, Houston, TX. Achnivu, O.I., Zhu, D., 2008. Field Application of an Interpretation Method of Downhole Temperature and Pressure Data for Detecting Water Entry in Inclined

Gas Wells. In: Paper SPE 115753 Presented at SPE Annual Technical Conference and Exhibition, Denver, Colorado, 21e24 September. http://dx.doi.org/10.2118/ 115753-MS. Almutairi, F.H., Davies, D.R., 2008. Modification of Temperature Prediction Model to Accommodate I-Well Complexities. In: Paper SPE 113594 Presented at Europec/ EAGE Conference and Exhibition, Rome, Italy, 9e12 June. http://dx.doi.org/10. 2118/113594-MS. Bird, R.B., Stewart, W.E., Lightfoot, E.N., 2002. Transport Phenomena. John Wiley and Sons, New York. Chekalyuk, E.B., 1965. Thermodynamics of Oil Formation. Nedra, Moscow. Fligelman, H., et al., 1989. Pressure Drawdown Test Analysis of a Gas Well-Application of New Correlations, vol. 4. SPEFE (3) pp. 406e412. Forchheimer, P., 1901. Wasserbewegung durch Boden. Z. Ver. Deutch 45. 1972e1788. Glasbergen, G., Gualtieri, D., Van Domelen, M., Sierra, J., 2009. Real-time Fluid Distribution Determination in Matrix Treatments Using DTS. In: SPEPO, vol. 24 (1), pp. 135e146. SPE-107775-PA. http://dx.doi.org/10.2118/107775-PA. Huckabee, P., 2009. Optic Fiber Distributed Temperature for Fracture Stimulation Diagnostics and Well Performance Evaluation. In: Paper SPE 118831 Presented at SPE Hydraulic Fracturing Technology Conference, the Woodlands, Texas, 19e21 January. http://dx.doi.org/10.2118/118831-MS. Huebsch, H., Moss, M., Trilsbeck, T., Brown, G., et al., 2008. Monitoring Inflow Distribution in Multi-zone, Velocity String Gas Wells Using Slickline Deployed Fiber Optic Distributed Temperature Measurements. In: Paper SPE 115816 Presented at SPE Annual Technical Conference and Exhibition, Denver, Colorado, 21e24 September. http://dx.doi.org/10.2118/115816-MS. John, L., Wattenbarger, A., 1996. Gas Reservoir Engineering. Society of Petroleum Engineers, Richardson, TX. Jones, S.C., 1987. Using the Inertial Coefficient, B, to Characterize Heterogeneity in Reservoir Rock. Paper SPE 16949 presented at SPE Annual Technical Conference and Exhibition, New Orleans, October 5e8. Kabir, C.S., Izgec, B., Hasan, A.R., Wang, X., et al., 2008. Real-time Estimation of Total Flow Rate and Flow Profiling in DTS-Instrumented Wells. In: Paper IPTC 12343 Presented at International Petroleum Technology Conference, Kuala Lumpur, Malaysia, 3e5 December. http://dx.doi.org/10.2523/12343-MS. KAPPA Technical Reference and User’s Guide. 2010. KAPPA. Kelly, C.T., 1999. Iterative Methods for Optimization. Society for Industrial and Applied Mathematics, Philadelphia. Li, Z., Zhu, D., 2009. Predicting Flow Profile of Horizontal Well by Downhole Pressure and DTS Data for Water-Drive Reservoir. In: Paper 124873 Presented at SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 4e7 October. http://dx.doi.org/10.2118/124873-MS. Muradov, K.M., Davies, D.R., 2009. Temperature Modeling and Analysis of Wells with Advanced Completion. In: Paper SPE 121054 Presented at EUROPEC/EAGE Conference and Exhibition, Amsterdam, the Netherlands, 8e11 June. http://dx. doi.org/10.2118/121054-MS. Pedrosa Jr., O.A., Aziz, K., 1986. Use of hybrid grid in reservoir simulation. SPERE 1 (6), 611e621. Pereira, C.A., Kazemi, H., Ozkan, E., 2006. Combined effect of non-Darcy flow and formation damage on gas-well performance of dual-porosity and dual-permeability reservoirs. SPEREE 9 (5), 543e552. http://dx.doi.org/10.2118/ 90623-PA. Smith, R.V., 1961. Unsteady-State gas flow into gas wells. J. Pet. Technol. 13 (11), 1151e1159. Swift, G.W., Kiel, O.G., 1962. The prediction of gas-well performance including the effect of non-Darcy flow. J. Pet. Technol. 14 (7), 791e798. Sui, W., Zhu, D., Hill, A.D., Ehlig-Economides, C., 2008. Determining Multilayer Formation Properties from Transient Temperature and Pressure Measurements. In: Paper SPE 116270 Presented at SPE Annual Technical Conference and Exhibition, Denver, Colorado, 21e24 September. http://dx.doi.org/10.2118/ 116270-MS. Sui, Weibo. 2009. Determining Multilayer Formation Properties from Transient Temperature and Pressure Measurements. PhD dissertation, Texas A&M U., College Station, TX. Wang, H., 2008. Application of Temperature Observation Wells during SAGD Operations in a Medium Deep Burial Extra Heavy Oil Reservoir. In: Paper PETSOC 2008-118 Presented at Canadian International Petroleum Conference, Calgary, Alberta, Canada, Jun 17e19. http://dx.doi.org/10.2118/2008-118. Wang, X., Lee, J., Vachon, G., 2008. Distributed Temperature Sensor (DTS) System Modeling and Application. Paper SPE 120805 Presented SPE Saudi Arabia Section Technical Symposium, Al-Khobar, Saudi Arabia, 10e12 May. Yoshioka, K. 2007. Detection of Water or Gas Entry into Horizontal Wells by Using Permanent Downhole Monitoring System. PhD dissertation, Texas A&M U., College Station, TX.