3D simulation of the thermal regime of a group of gas wells in the Srednetyungskoe field Текст научной статьи по специальности «Физика»

Математические заметки СВФУ Июль—сентябрь, 2019. Том 26, № 3

UDC 622.324.5

3D SIMULATION OF THE THERMAL REGIME OF A GROUP OF GAS WELLS IN THE SREDNETYUNGSKOE FIELD V. A. Ivanov, P. V. Sivtsev, and I. I. Rozhin

Abstract. The thermal interaction between a group of gas wells and the permafrost ground was simulated using the finite element method for the conditions of the Srednety-ungskoe field. For the spacings between the wellheads of 10 m, 15 m, and 20 m, thermal regimes for the gas and surrounding ground massif are forecasted for 30 years of operation. The feature of the present work is that atmospheric conditions and the effect of the wells are modeled simultaneously in 3D, allowing for accurate characterization of the wellhead area. Field data on the ground temperature, the bottomhole pressure and temperature, the production rate, the gas composition, thermophysical properties of the ground, and weather conditions are used as input parameters. Temperature profiles of the ground, the position of the thawing front, and timing of merger of the thawing halos around the wells are estimated. The results of the survey can be used in development planning of the field.

DOI: 10.25587/SVFU.2019.26.34.009 Keywords: permafrost, gas well, Stefan problem, numerical simulation, Srednetyung-skoe field, zonal isolation

Introduction

Srednetyungskoe is a large gas and condensate field located in a remote area of the Sakha Republic, far from transport routes. Currently, small amount of gas is produced there for local use and the field is considered as a perspective object for future development.

The region belongs to the subarctic zone where the ground is frozen from the surface down to the depth of about 500-600 meters, forming permafrost. Thus, oil and gas wells there have to pass through the frozen soil in order to reach the reservoir. When the production starts, the fluid from the reservoir flows up through the permafrost layer causing its thaw. The thawing is in fact inevitable because the fluid comes from deeper horizons and therefore has higher temperature. The consequences of this process could be destabilization of the upper section of the well, breaking of the casing integrity, erosion of the zonal isolation leading to gas seepage [1-3]. Most probably, the named issues will emerge after a long stretch of time at an unexpected moment. And the cost of the remedial works as well as the

The research was supported by the Government of the Russian Federation Megagrant No. 14.Y26.31.0013.

© 2019 V. A. Ivanov, P. V. Sivtsev, I. I. Rozhin

environmental impact would be rather high. Therefore, we need to carefully predict the thermal processes in and out of the petroleum wells through their lifecycle. This can be done by means of numerical simulation based on a mathematical model describing heat transfer in frozen rocks.

In [4, 5] a problem of permafrost thawing around a single well with constant temperature is solved analytically. While this kind of approach is convenient for rapid estimates, it cannot be used for precise quantitative calculations. The main reason being that the gas temperature in the well can vary with time and depth, thus affecting the thermal regime of the surrounding rocks considerably. In [3, 6, 7] a conjugate numerical problem in an axisymmetric formulation is solved for single oil and gas wells. Heat transfer in the producing petroleum well and in the surrounding rocks is simulated taking into account the key factors influencing it. The numerical approach is quite simple and effective. However, in the modern oil and gas projects the wells are often located very closely so that their zones of thermal influence impose on each other and the axisymmetric formulation cannot be used. In [8] a 3D Stefan problem, describing the dynamics of ground thawing around two adjacent producing wells, is solved using the finite difference method. The authors divided the thawing process into three stages: 1) the wells do not affect each other and act as if they were isolated and single; 2) there is a merger of thawing halos around wells; 3) transition to a quasi-stationary regime with slowing down the movement of the phase boundary. The research doesn't account for the cooling of the gas during its flow upwards through the well. Meanwhile the throttling effect and heat losses to the surrounding rocks can cause an intensive cooling. In [9] similar problem solved using finite element method [10-13].

The objective of the present work is to build a 3D conjugate numerical model of thermal regime of a well cluster and surrounding permafrost rocks using the finite element method. The developed model can be used to make a quantitative estimate of the destabilized thawed zone around group of wells. At the same time, the temperature of the produced gas will be calculated. Therefore, a gas hydrate formation can be predicted and prevented as well as certain measures against destructive warming processes in the permafrost can be taken.

1. Method

Consider a group of 5 slant wells whose heads are arranged in a row (Fig. 1). The wells produce natural gas from great depths of several kilometers. The upper sections located in permafrost zone are vertical, while lower sections in the thawed zone are slant with different azimuths. Thus, the wells in their lower sections can be modeled as single, but in their upper sections they are located in close proximity and interact with each other thermally. Therefore, in permafrost the group of wells must be treated as whole.

As the gas is mainly used for heating buildings, the production rates vary with season of the year. In winter all the wells are under production and in summer many of them shut down.

group of gas wells

permafrost zone

-500 m

producing layer

thawed zone

Fig. 1. Computational domain.

The studied physical system consists of three processes: 1) heat transfer with the gas moving upwards inside the wells; 2) heat transfer in the surrounding rocks; 3) heat transfer between the gas and the rocks.

The first process is described using the differential equations from [13], given initial conditions at the bottom hole. The equations are stationary as the transition processes in the flowing gas have much higher rates than in the rock massif. In this case, the state of the gas will be immediately reacting to slow change in the outside temperature of the rocks Tg.

The equations of the mathematical model are:

dp dx

= -pg sin y

4 pS2*

dT dp dx dx

^Dot rri\ 9 •

p(0)= po, T (0) = To.

(1)

(2) (3)

Here x is the coordinate along the axis of the well, p is the gas pressure, p is the gas density, g is the gravity, p is the well inclination measured from horizontal, ^ is the hydraulic resistance coefficient, M is the mass flow rate, S is the cross-section of the well, T is the gas temperature, e is the throttling coefficient, D is the inner diameter of the well, at is the total heat transfer coefficient through the wellbore wall, cp is the specific heat of gas at constant pressure, p0 is the pressure at the bottom hole, and T0 is the temperature at the bottom hole.

The real gas law is defined according to Latonov, Gurevich [14]:

p=

p

ZRT'

Z= ( 0,17376In ( Y^j + 0,73^ PC +0,1^-.

(4)

(5)

Here Tc is the critical temperature and pc is the critical pressure which depend on the gas composition.

The throttling coefficient e is dependent on the imperfection coefficient Z and is defined as follows:

RT2 iaz\ <«<

To describe the heat transfer in the rocks, both frozen and thawed, we neglect mass transfer processes. Thus, we can use the heat equation as follows [15]:

~ dT

C(Tg)^=V(\(Tg)VTg)

(7)

Here C(Tg) is the piecewise constant volumetric heat capacity of rocks, depending on the temperature Tg, which contains implicit heat of the phase transition by means of the Dirac delta function; A(Tg) is the piecewise constant coefficient of thermal conductivity of rocks, also dependent on temperature Tg; t is time. We use the following approximations:

As

at Tg <Tf - AT,

KTg) = { As - (As - A;)[Tg'2TATAT] at Tf - XT < Tg < Tf + AT, (8)

Ai

C(Tg )

Cs

Cf + (Cs - Cf )

Cf + (Ci - Cf ) Ci

-Jl

at

zIL

at

at Tg > Tf + AT, at Tg < Tf - AT,

2

at Tf - AT < Tg < Tf,

at Tf < Tg < Tf + AT, , at Tg > Tf + AT.

(9)

Here Tf is the temperature of the ice-water phase transition in the ground; AT is the temperature half-interval for smoothing the delta function; As and Ai are the thermal conductivities of the ground in the solid and liquid phases, respectively; Cs and Ci are the volumetric heat capacities of the soil in the solid and liquid phases, respectively; and Cf is the volumetric heat capacity at the phase transition temperature, which is found by the following formula:

(10)

Here W is the implicit heat of phase transition of a unit volume of wet ground.

In the rock massif from the bottom hole of the well to the bottom of permafrost, thermophysical coefficients in the heat equation (7) are piecewise constant and the equation can be solved using standard methods. Because the wells in the lower thawed zone are considered as single, there we use axisymmetric formulation of the heat equation as in the work [3].

To make up a boundary problem, the equation (7) is completed with three boundary conditions. First, on the far side of the computational domain heat flow is set equal to zero:

dTg dn

0.

(11)

g

2

g

Second, at the wellbore wall we have Newton's law of cooling with the gas:

dT

\{Tg)-5±=at{Tg-T). (12)

Third, at the upper surface we state the 3rd type boundary condition to account for the weather influence:

dT

KTg)-^ = as(Tg - Tair) + qs (13)

Here at is the total coefficient of heat transfer from gas to rocks, as is the convection coefficient of the ground surface, Ta;r is the air temperature, and qs is the total heat flow on the ground surface, n is the normal to the surface.

We use the finite element method to solve the heat transfer problem (7)-(13) and state it in a variational form. Let us define Q as a computational domain in the ground consisting of both thawed and frozen zones and H(Q) as a Sobolev space. Then the goal is to find T G H(Q), complying with the following:

a(Tn+1,v) = L(v), v G H(Q), (14)

a{Tn+\v) = i [ C(Tg)nTn+1vdx+ f (A(Tff)n gradTn+1, gradf) cte

T

n

+ y atTn+1vds + ^ asTn+1vds, (15)

IW

L(v) = - f C{Tg)nTnvdx+ f atTgasvds+ f asTn+1vds. (16)

T

Here H(Q) = {v G H(Q) : v(x) = 0, x G rout}, rout is the outer boundary of Q, rw is the boundary between the wellbore and the ground, rair is the upper surface boundary, Tn+1 is the ground temperature at the next time step tn+1, and Tn is the ground temperature at the current time step tn. The linear basis functions are used.

When solving the time-dependent equations (1)-(16) numerically, each time step consists of calculating the temperature and pressure of the gas and then of calculating the temperature distribution in the rock massif. Numerical model based on these equations has been tested on synthetic case against analytical solution.

The computational experiment is divided into two stages:

1st stage. In the lower thawed zone, the thermal boundary problem is solved in the axisymmetric formulation for single wells. Here, the computational domain in the rock massif is a cylinder with radius of 70 m and with height of the entire length of the inclined section of the well from the bottom to the permafrost base, 3420 m. The well is located in the center of the cylinder, making it hollow. Boundary condition (11) is set on outer sides of the cylinder, and boundary condition (12) is set on the borehole wall. The initial temperature of the rocks is set according to the

logging data. The problem is solved using the finite difference method. Distribution of temperature and pressure of the gas in the well is found by solving the first-order Cauchy problem for equations (1)-(3).

2nd stage. In the upper frozen zone, a Stefan problem is solved in 3D formulation for the entire well cluster. The computational domain is a 3D box with dimensions of 70 x 50 m in horizontal plane and of height from the permafrost base to the surface, 670 m (Fig. 1). The wells are located in a row on the side with dimensions 70 x 670 m. Boundary condition (11) is set on all outer sides of the 3D box, except for the upper side: on the far sides, due to the distance from the wells; on the sides where the wells are located, due to the symmetry of the problem formulation; on the bottom, due to the vertical heat flux being neglected. On the upper side we have the climatic condition (13). Boundary condition (12) is set on the borehole walls. Initial temperature distribution of the rocks is set according to the logging data. The problem is solved using the finite element method. Distributions of temperature and pressure of the gas in the wells are found by solving the firstorder Cauchy problem for equations (1)-(3), where the results of the 1st stage are taken as initial conditions.

The temperature regime of the rocks and the group of wells is projected for 30 years of operation.

Information on composition of the rocks, their thermal characteristics, thermal logging data, and the gas properties are used as input parameters for the calculations. The density, water content, thermal conductivity, and thermal capacity of sedimentary rocks are examined from a core of the Lena-Vilyuy oil and gas province (Table 1) [16], where the reviewed Srednetyungskoe field is located. The climatic data, which is presented in Table 2, is taken from [17].

The depth of the productive formation is 3100 m, the depth of the permafrost base is 670 m, the geothermal gradient in the thawed zone is 36.9°C/1000 m, the geothermal gradient in the frozen zone is 4.1°C/1000m, the formation temperature is 73°C, the formation pressure is 320 atm, and the production rate is 8.5kg/s. The total production rate of the 5 wells varies depending on the season of the year, reaching the maximum in January and the minimum in July. This is caught up by shutting down the wells one by one in spring and then putting them back into production in autumn. The schedule of the wells' operation is arranged so that the accumulated productions of the wells are leveled.

The gas properties are as follows: the hydraulic resistance coefficient of the gas is 0.02, the specific heat at constant pressure is 2300 J/(kg • K), the critical temperature is 217.44°K, the critical pressure is 46.13 • 105 Pa, the gas law constant is 398.47 J/(kg • K), and the gas density at bottom hole is 0.15 kg/m3. The well design data is the following: the radius of the production tubing is 57 mm, the radius of the inner casing is 89 mm, the radius of the outer casing is 245 mm, the thermal conductivity of the cement is 0.169 W/(m • K), and the thermal conductivity of the gas is 0.52 W/(m • K).

Table 1. Properties of the ground

Depth interval, m Density, kg/rrfi Water content Thermal conductivity of thawed ground, W/(m • K) Thermal conductivity of frozen ground, W/(m • K) Volumetric heat capacity of thawed ground, 106 J/(m3 • K) Volumetric heat capacity of frozen ground, 106 J/(m3 • K)

0-86 2000 0.12 1.69 1.93 2.57 2.31

86-488 2000 0.12 1.62 1.86 2.68 2.42

488-670 2120 0.1 1.44 1.67 2.72 2.51

Table 2. Climatic data

Month I II III IV V VI VII VIII IX X XI XII

Thermal radiation balance of an earth surface, W/m2 -9.65 -4.63 -5.02 45.14 113.04 126.16 119.6 79.48 35.49 0 -14.27 -9.65

Monthly average atmospherical temperature (1991-2003), °C -37.5 -32.7 -21.4 -5 7.3 16.6 19.8 15.4 5.8 -7.4 -27.6 -38

Wind velocity, rn/s 1.4 1.4 2 2.8 3.4 3.3 3 2.8 2.6 2.6 2 1.3

2. Results

In Fig. 2 the gas temperature for spacing of 10 m between wellheads is presented for four time points compared to the initial temperature of the ground. The gas cools while flowing up the well due to interaction with the surrounding ground and the throttling effect. The simulation shows that the flow rate is high enough to prevent leveling of the gas temperature with that of the ground. One can also see that the gas temperature changes insignificantly during 30 years of operation.

In Fig. 3 the vertical temperature profiles of the ground surrounding five gas wells near the surface are presented. Upper 1-3 meters of the ground experience thawing and freezing according to season of the year. Whereas the deeper ground continuously thaws due to the thermal effect of the group of wells. Also we can see that 3 wells are shut down in September and 2 wells are shut down in March. These wells can be found out by colder color of the surrounding ground: red is for the thawed zone, blue is for the frozen zone, and the colors in between are for the phase transition zone.

In Fig. 4 the radius of thawing is plotted against the depth for the spacing of 20 m between wellheads at four time points. A constant growth of the thawing halo is clearly seen. Because the gas has much higher temperature than the frozen

Fig. 2. Initial temperature of the ground and the gas temperature for spacing 10 m.

Fig. 3. Vertical temperature profile of the ground surrounding five gas wells near the surface at four time points with spacing of 10 m between the wellheads: (a) in March of the 1st year, (b) in September of the 1st year, (c) in March of the 3rd year, and (d) in September of the 3rd year.

Fig. 4. Permafrost thawing radius around the gas producing well at four time points of the field operation.

Fig. 5. Time of merger of the thawing halos around adjacent wells for three spacings between wellheads.

ground we always have positive heat flux as soon as the well is in active production. The shape of the curves is determined mostly by the initial temperature of the ground. The only area where the radius is not increasing is the upper part of the computational domain subject to climatic interference.

In Fig. 5 the time of merger of the thawing halos around the adjacent wells is plotted against the depth for three spacings between wellheads. For the spacing of 10 m, the thawing halos merge by 4-6 years after the start of the field development. For the spacing of 15 m, the timing is 12-15 years and for the spacing of 20 m, it is 22-25 years. Note that in all of the considered cases the thawing halos experience a merger, sooner or later. As with the radii of thawing, the shape of the curves is determined mostly by the initial ground temperature.

Conclusion

Under the given conditions, the thawing halos around the production wells continuously grow and eventually merge into one. Time of the merger depends on the spacings between the wells: for 10 m, 4-6 years after the start of the field development, for 15 m, 12-15 years, and for 20 m, 22-25 years. This circumstance must be taken into account in the future development plan of the Srednetyungskoe field as it can cause certain environmental issues.

REFERENCES

1. Permyakov P. P., "Vliyanie kriolitozony v osnovanii podvodnogo perehoda gazoprovoda cherez r. Lenu [in Russian]," Gaz. Promyshl., No. 2, 59-61 (2013).

2. Sheveleva D. V., Dinamika Slozhnogo Teplovogo Vzaimodeystviya Neftyanyh i Gazovyh Skva-zhin s Mnogoletnemyorzlymi Porodami [in Russian], Diss. ... Kand. Fiz.-Mat. Nauk, Tyumen (2008).

3. Argunova K. K., Bondarev E. A., and Rozhin I. I., "Teplovoe vzaimodeystvie neftedobyvayu-shchikh skvazhin s mnogoletnemerzlymi gornymi porodami [in Russian]," Nauka i Obraz., No. 4, 78-83 (2008).

4. Polozkov K. A., "Vybor rasstoyaniy mezhdu kustovymi dobyvayushhimi skvazhinami v zonakh MMP [in Russian]," Stroitel'stvo Neftyanyh i Gazovykh Skvazhin na Sushe i na More, No. 3, 21-29 (2008).

5. PAO Gazprom, VRD 39-1.9-015-2000, Rukovodstvo po Termometricheskim Metodam Kontro-lya Kachestva Stroitel'stva, Krepleniya Skvazhin v Mnogoletnemerzlykh i Nizkotemperatur-nykh Porodakh [in Russian], OOO VNIIGAZ, Moscow (2001).

6. Bondarev E. A., Rozhin I. I., and Argunova K. K., "Modeling the formation of hydrates in gas wells in their thermal interaction with rocks," J. Eng. Phys. Thermophys., 87, No. 4, 900-907 (2014).

7. Bykov I. Yu., Pushkin V. N., and Pushkin V. V., "Dynamics of axis-symmetric frosting of wellhead area in conditions of ventilated shaft directions with account of cement ring influence [in Russian]," Stroitel'stvo Neftyanyh i Gazovyh Skvazhin na Sushe i na More, No. 3, 15-19 (2011).

8. Gorelik J. B., Shabarov A. B., and Sysoyev Yu. S., "The dynamics of frozen ground melting in the influence zone of two wells [in Russian]," Kriosfera Zemli, 12, No. 1, 59-65 (2008).

9. Afanaseva N. M. and Kolesov A. E., "Numerical solution of the thermal influence of oil well cluster on permafrost," in: AIP Conf. Proc., 1773, No. 1, 110001 (2016).

10. Vasilyeva M. V., Zakharov P. E., Sivtsev P. V., and Spiridonov D. A., "Numerical modeling of thermoelasticity problems for constructions with inner heat source [in Russian]," Mat. Zamet. SVFU, 24, No. 3, 52-64 (2017).

11. Vasilyeva M. V. and Prokopiev G. A., "Numerical solution to the problem of two-phase filtration with heterogeneous coefficients by the finite element method [in Russian]," Mat. Zamet. SVFU, 24, No. 2, 46-62 (2017).

12. Grigoriev V. V. and Zakharov P. E., "Numerical modeling of the two-dimensional Rayleigh-Benard convection [in Russian]," Mat. Zamet. SVFU, 24, No. 1, 87-98 (2017).

13. Bondarev E. A., Vasil'ev V. I., Voevodin A. F., Pavlov N. N., and Shadrina A. P., Termogidro-dinamika Sistem Dobychi i Transporta Gaza [in Russian], Nauka, Novosibirsk (1988).

14. Latonov V. V. and Gurevich G. R., "Raschet koeffitsienta szhimaemosti prirodnykh gazov [in Russian]," Gaz. Promyshl., No. 2, 7-9 (1969).

15. Samarskii A. A. and Moiseenko B. D., "An economic continuous calculation scheme for the Stefan multidimensional problem," U. S. S. R. Comput. Math. Math. Phys., 5, No. 5, 43-58 (1965).

16. Semyonov V. P., Geotemperaturnoe Pole i Kriolitozona Vilyuyskoy Sineklizy [in Russian], Diss. ... Kand. Geol.-Min. Nauk, FGBUN Inst. Merzlotovedeniya im. P. I. Mel'nikova SO RAN, Yakutsk (2018).

17. Gavrilova M. K., Klimat Tsentral'noy Yakutii [in Russian], Yakut. Knizh. Izdat., Yakutsk (1973).

Submitted July 19, 2019 Revised August 29, 2019 Accepted September 3, 2019

Victor A. Ivanov

Institute of Oil and Gas Problems SB RAS, 1, Oktyabrskaya Street, Yakutsk 677000, Russia [email protected] Petr V. Sivtsev

International Research Laboratory "Multiscale Model Reduction", Ammosov North-Eastern Federal University, 42 Kulakovsky Street, Yakutsk 677000, Russia [email protected] Igor I. Rozhin

Institute of Oil and Gas Problems SB RAS, Oktyabrskaya str., 1, Yakutsk 677000, Russia [email protected]