3D SIMULATION OF THE THERMAL REGIME OF A GROUP OF GAS WELLS IN THE SREDNETYUNGSKOE FIELD

. The thermal interaction between a group of gas wells and the permafrost ground was simulated using the ﬁnite element method for the conditions of the Srednety-ungskoe ﬁeld. 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 eﬀect 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 proﬁles 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 ﬁeld.


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][2][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 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][11][12][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.

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.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 T g .
The equations of the mathematical model are: Here x is the coordinate along the axis of the well, p is the gas pressure, ρ is the gas density, g is the gravity, ϕ 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, ε is the throttling coefficient, D is the inner diameter of the well, α t is the total heat transfer coefficient through the wellbore wall, c p is the specific heat of gas at constant pressure, p 0 is the pressure at the bottom hole, and T 0 is the temperature at the bottom hole.The real gas law is defined according to Latonov, Gurevich [14]: Here T c is the critical temperature and p c is the critical pressure which depend on the gas composition.
The throttling coefficient ε is dependent on the imperfection coefficient Z and is defined as follows: 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]: Here C(T g ) is the piecewise constant volumetric heat capacity of rocks, depending on the temperature T g , which contains implicit heat of the phase transition by means of the Dirac delta function; λ(T g ) is the piecewise constant coefficient of thermal conductivity of rocks, also dependent on temperature T g ; t is time.We use the following approximations: Here T f is the temperature of the ice-water phase transition in the ground; T is the temperature half-interval for smoothing the delta function; λ s and λ l are the thermal conductivities of the ground in the solid and liquid phases, respectively; C s and C l are the volumetric heat capacities of the soil in the solid and liquid phases, respectively; and C f is the volumetric heat capacity at the phase transition temperature, which is found by the following formula: 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: Second, at the wellbore wall we have Newton's law of cooling with the gas: Third, at the upper surface we state the 3rd type boundary condition to account for the weather influence: Here α t is the total coefficient of heat transfer from gas to rocks, α s is the convection coefficient of the ground surface, T air is the air temperature, and q s 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 as a computational domain in the ground consisting of both thawed and frozen zones and H( ) as a Sobolev space.Then the goal is to find T ∈ H( ), complying with the following: Here H( ) = {v ∈ H( ) : v(x) = 0, x ∈ out }, out is the outer boundary of , w is the boundary between the wellbore and the ground, air is the upper surface boundary, T n+1 is the ground temperature at the next time step t n+1 , and T n is the ground temperature at the current time step t n .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 × 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 × 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/1000 m, the formation temperature is 73 • C, the formation pressure is 320 atm, and the production rate is 8.5 kg/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 • 10 5 Pa, the gas law constant is 398.47 J/(kg • K), and the gas density at bottom hole is 0.15 kg/m 3 .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).

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   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.

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

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

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