CC BY
8
DOI: 10.17516/1997-1397-2021-14-6-795-804 УДК 519.6
Thawing of Permafrost During the Operation of Wells of North-Mukerkamyl Oil and Gas Field
Mikhail Yu. Filimonov* Nataliia A. Vaganova^
Krasovskii Institute of Mathematics and Mechanics Yekaterinburg, Russian Federation Ural Federal University Yekaterinburg, Russian Federation
Received 05.05.2021, received in revised form 31.06.2021, accepted 10.08.2021 Abstract. Thawing of ice-saturated rocks due to climate change or various technological impacts will be accompanied by subsidence of the earth's surface and development of dangerous permafrost geological processes called thermokarst, leading to accidents, which may destruct the wells. Currently, the investment programs of the development of new northern oil and gas fields are restricted. In this regard, reducing the cost of developing the oil and gas fields is an urgent problem. For example, diminishing the area of well pads and maintaining efficiency in the northern oil and gas fields can significantly reduce the costs, in particular, during the design stage. A model of unsteady thermal fields propagation in frozen soil from new well construction for the North Mukerkamyl oil and gas field is developed, taking into account the construction features, the annulus, and the complex lithology of the soil surrounding the well. It is planned to take into account climatic and technological factors, in particular, an annual rest period of well operation, which held from several hours to two weeks. The paper discusses the computational features of the thermal fields calculating in frozen ground from wells and explores the influence of various parameters, which in the computations may lead to a significant increasing of thawing area in the well pads.
Keywords: mathematical modelling, heat transfer, permafrost.
Citation: M.Yu. Filimonov, N.A. Vaganova, Thawing of Permafrost During the Operation of Wells of North-Mukerkamyl Oil and Gas Field, J. Sib. Fed. Univ. Math. Phys., 2021, 14(6), 795-804. DOI: 10.17516/1997-1397-2021-14-6-795-804.
Introduction
In Russia, the reserves of underground ice of the permafrost zone with amount up to about 19000 km3, which gives the right to call permafrost underground glaciations. These areas cover 10 million km2. Permafrost thawing due to climate change and various technogenic heat sources can lead to development of thermokarst and major accidents, which are especially dangerous to permafrost [1-3]. Modeling of thermal fields in frozen soil is an important problem [4-6], since the stability of various engineering structures depends on the preservation of the negative temperatures in the soil. There are various technical systems that affect both the rate of thawing of permafrost [7-9] and the thermal stabilization of the soil [10,11].
The most of the oil and gas fields of Russia are located in the permafrost zone. The operation has the specifics and restriction. In accordance with Russian Constructive Standards, before starting the operation and development of an oil and gas field, it is necessary to determine the
* fmy@imm.uran.ru t vna@imm.uran.ru © Siberian Federal University. All rights reserved
minimum distance between the producing wells located on well pads. The following restriction is used: two wells cannot be drilled at the distance lower than two radii of thawing from each other (i.e., lower than the distance which will be covered by zero isotherm for 30 years of operation of a single well). This means that a reduction in the radius of thawing will lead to significant savings of expensive materials used to build the well pads, and therefore to decreasing the cost of produced oil, which, subject to a possible reduction in oil prices, can be a determining factor.
Currently, the investment programs of development of new northern oil and gas fields are restricted. In this regard, reducing the cost of designing and optimal operation of new and existing oil and gas fields is an urgent problem. For example, diminishing the area of well pads and maintaining efficiency in the northern oil and gas fields can significantly reduce the costs, in particular, during the design stage. When the existing northern fields operate, there are various options to carry out some planned technological operations related, in particular, with wells rest period, or flare systems using for the gas utilization, etc. Optimal strategies of operation are important for minimization of exploitation costs.
In this paper a mathematical model and software for calculating the propagation of unsteady thermal fields in frozen soil from new well designs for the North-Mukerkamyl oil and gas field are proposed. In numerical calculations, the design features of the well are taken into account, which makes it possible to more accurately estimate the thawing radius from the producing well and is recommended to use heat-insulating materials for certain modes of its operation. The results of numerical calculations and recommendations related to the accuracy of obtaining numerical solutions for such problems are presented.
1. Statement of the problem and mathematical model for two wells
Following [3], we consider a mathematical model of non-stationary thermal fields propagation in frozen ground from two producing wells. To simulate non-stationary thermal fields in the permafrost zone various climatic, physical and technological factors should be taken into account when different technical systems operate.
Let the computational domain is a three-dimensional box in which the vertical wells are considered as inner boundaries (Fig. 1). The considered soil is possible to be non-homogenous and may include a number of elements, such as ice-rock lenses, different layers, engineering constructions (riprap), and insulating layers around the wells (for example, Fig. 2).
A mathematical model with including localized heat of phase transition is considered. This approach allows to solve the problem of Stefan type without the explicit separation of the phase transition [12]. The heat of phase transformation is introduced with using Dirac ¿-function as a concentrated heat of phase transition in the specific heat ratio.
First we consider heat exchange on a flat ground surface illuminated by the sun. The computational domain is a three-dimensional box Q, where x and y axes are parallel to the ground surface and the z axis is directed downward. We assume that the size of the box Q is defined by positive numbers Lx, Ly, Lz: -Lx ^ x ^ Lx, -Ly ^ y ^ Ly, -Lz ^ z ^ 0. Let T = T(t, x, y, z) be soil temperature at the point (x,y, z) at the time moment t. The main heat flow associated with climatic factors on the surface z = 0 is shown in Fig. 1.
The ground surface z = 0 is a main zone of formation of natural thermal fields. On this surface the equation of balance of flows is used as a boundary condition, with taking into account the main climate factors: air temperature and solar radiation. Tair = Tair (t) denotes the temperature in the surface layer of air, which varies from time to time in accordance with the annual cycle of temperature; a = 5.67 • 10~8 W/(m2K4) is Stefan-Boltzmann constant; b = b(t,x,y) is heat transfer coefficient; A = X(T) is thermal conductivity coefficient; e = e(t,x,y) is the coefficient of emissivity. The coefficients of heat transfer and emissivity depend on the type and condition of
Fig. 1. The main heat flows and boundary conditions
the soil surface. 0 can include a number of engineering structures. Suppose that in 0 there are 2 objects that are heat sources (for example producing insulated wells). We denote the surface of these objects by 0j = 0i(x,y,z), i = 1,2 (Fig. 1).
Thus, the modeling of thawing in the soil leads to the solution of the following heat equation [12] in 0:
dT
p(cv(T) + kS(T - T*)) — = V (X(T)AT), (1)
where p is density [kg/m3], T* is temperature of phase transition [K],
cv(T) = { ^XyZi T > T** is specific heat [J/(kg • K)],
f A (x z) T < T* X(T) a Х У z) T > T* is thermal conductivity coefficient [W/(m • K) ],
k = k(x,y, z) is specific heat of phase transition, S is Dirac delta function.
Thus, it is necessary to solve equation (1) in the area Q with initial condition
T(0, x, y, z) = To(x,y,z) (2)
and boundary conditions [13]
Yq + b(Tair - T (x, y, 0, t)) = ea(T 4(x, y, 0, t) - Ta4ir ) + \9T (x^ M), (3)
T
= Ti(i), T
Qio
= T2 (t), (4)
^20
Fig. 2. The composition of the shell of the producing well
dT dx
x=±Lx
dT dy
V = ±Lv
dT
dz
0.
(5)
-—Lz
Initial condition (2) describes an initial distribution of soil temperature at the time moment from which we plan to start the numerical calculation. The method of determining the parameters included in the nonlinear boundary condition (3) is described in [3,13].
Thus, to solve the problem by computer modeling it is necessary to consider the following factors:
— solar radiation and emissivity, which lead to nonlinear boundary condition (3) on the upper
surface;
— various thermophysical parameters of soils that form permafrost that vary not only in the
vertical plane;
— seasonal variations in air temperatures in the area under investigation (seasonal thawing of
the upper layer of the soil in summer, the reverse freezing process in winter);
— possible phase transitions in the soil;
— the structure of the base thermal insulation of the well pad;
— technological features of wells.
2
2. Numerical results
A large number of works are devoted to development of numerical methods for solving boundary value problems of heat conduction. Basics of finite difference methods are detailed in the works [14,15].
To solve the problem (1)-(5), a finite-difference method is developed for solving problems of Stefan type [12,16], which allows using the splitting method with respect to the space variables. At the same time, the algorithm is "adapted" to the geographic location of a particular oil and gas field [3,13].
On the surface z = 0 an algebraic equation of the fourth degree is generated, to be solved by Newton method. The solvability of the corresponding difference equations approximating the original problem was proved in the paper [17].
Calculations based on the proposed models and developed programs were used for 13 northern oil and gas fields. When the source of heat in the frozen ground was a well, a comparison was made between numerical data on the distribution of the boundary of the melting of frozen ground (zero isotherm) and experimental data. The comparison showed that the numerical data coincide with the experimental data with accuracy of 5%.
Numerical calculations were carried out in accordance with the technical characteristics of wells in Fig. 2 for the North-Mukerkamyl oil and gas field and the climate data at the surface 2 = 0 in accordance with Fig. 3. For this oil and gas field, the temperature of the phase transition T* in the soil is equal to -0.2oC.
Fig. 3. Climate data
In numerical calculations, an orthogonal condensing grid is used. In the {x, y}-plane, the computational grid is condensed in the shell of the well. In the vertical plane (along the 2 coordinate), the mesh is refined near the boundary of various soil layers. Fragments of a condensing grid in the {x, y}-plane are shown in Fig. 4. In the calculations, the following mesh sizes are used: 97 x 97 x 21 = 197589 nodes, 197 x 197 x 21 = 814989 nodes, 397 x 397 x 21 = 3309789 nodes, 697 x 697 x 21 = 10201989 nodes.
Numerical calculations showed that the thawing radius depends not only on the total number of calculated nodes, but also on the number of nodes that concentrated in a neighborhood of a radius of 1.5m around the well. In Fig. 5 a graph of the thawing radius depending on the time of operation of the well is shown.
A special role in the accuracy of finding the radius of thawing also belongs to the number of nodes near the well. In Fig. 5 the graphs of radius of thawing for a different number of nodes in the neighborhood of the well (from 51 nodes to 151 nodes) are preented. Calculations show that we can use 101 nodes in the neighborhood of the well. Therefore, we take a grid consisting of 397 x 397 x 21 = 3309789, as a computational grid. The calculation time of one year of production simulation for a single production well on a supercomputer ranged from 30 minutes to 90 minutes of machine time.
In Fig. 6 thawing radius are presented for the following well insulation options: (1) corresponds to the case of the all thermal insulations consist only of ordinary cement; (2) is for the simple thermal insulation of the well with using cement shells; (3) is the special case of the combined insulation of the well for the North-Mukerkamyl oil and gas field. Graph (3) shows
Fig. 4. Computational grid in {x,y}- plane
40m_197_51
40m_297_51
40m_297_101
40m_397_101
40m_397_151
40m_697_101
Fig. 5. Thawing radius
that taking into account the design features of the wells can significantly refine the evaluation of the thawing radius and the distance between the wells.
To estimate the distance between the neighboring wells and the moment of the taliks junction, it is also necessary to take into account the mutual thermal effect of these wells. In [3], it was investigated and shown that when designing an oil and gas field, taking into account only the
Fig. 6. Thawing radius for three well designs
Fig. 7. Continuous well operation for 30 years
thawing radius for a single well may be insufficient. Consider a simulation of the operation of two neighboring wells located at the distance of 7.6 meters from each other. Numerical calculations were carried out on the following grids: 657 x 347 x 45 = 8822475 nodes and 651 x 397 x 41 = 10596327 nodes. The calculation time of one year of modeling the operation of two neighboring wells ranged from 70 to 140 minutes. Fig. 7 shows the isotherms (— 0.2oC) in the vertical plane {x, z} for two continuously operating neighboring wells for 30 years in October, which intersect before the end of the operation. Fig. 8 shows the temperature fields for the case corresponding to Fig. 7, but with annual preventive rest period for 14 days. In this case, the taliks are not joined.
N
-15
-10
-5
0
-20
0 5 10 15 20
25 30 35 +0 +5
0 5 10 15 20
25 30 35 40 45
X
X
Fig. 8. Well operation for 30 years with annual preventive rest period for 14 days
Conclusion
The numerical calculations for a specific oil and gas field in the permafrost zone result the following:
— when thermal fields are simulated in a frozen ground, it is necessary to choose a sufficiently
large computational domain to avoid the influence of boundary conditions of the lateral boundaries of the computational region;
— it is necessary to use a condensed grid so that a sufficient number of grid nodes get into the
vicinity of the well. It was found experimentally that in our case it is enough that at least 101 nodes of the computational grid in a horizontal plane are in the vicinity of 1.5 m near the well;
— computer simulation of the operation of two wells over 30 years showed that when the wells
are located at the distance of 7.6 meters from each other and with the annual prevention rest period for 14 days, the thawing radius of these wells are not joined.
The work was supported by Russian Foundation for Basic Research 19-07-00435.
References
[1] F.E.Nelson, O.A.Anisimov, N.I.Shiklomanov, Subsidence risk from thawing permafrost, Nature, 410(2001), 889-890. DOI: 10.1038/35073746
[2] R.R.Muskett, L-Band InSAR Penetration Depth Experiment, North Slope Alaska, Journal of Geoscience and Environment Protection, 5(2017), no. 3, 14-30.
DOI: 10.4236/gep.2018.69016
[3] M.Y.Filimonov, N.A.Vaganova, Simulation of Technogenic and Climatic Influences in Permafrost for Northern Oil Fields Exploitation, Lecture Notes in Computer Science, 9045(2015), 185-192. DOI: 10.1007/978-3-319-20239-6_18
[4] N.Vaganova, M.Filimonov, Simulation of freezing and thawing of soil in Arctic regions, IOP Conference Series: Earth and Environmental Science, 72(2017), 012005.
DOI: 10.1088/1755-1315/72/1/012005
[5] M.S.Krass, V.G.Merzlikin, The Models of Permafrost for Arid and Alpine Regions, IOP Conference Series: Earth and Environmental Science, 272(2019), 022157.
DOI: 10.1088/1755-1315/272/2/022157
[6] N.A.Vaganova, M.Y.Filimonov, Computer Simulation of Nonstationary Thermal Fields in Design and Operation of Northern Oil and Gas Fields, AIP Conference Proceedings, 1690(2015), 020016. DOI: 10.1063/1.4936694
[7] N.Vaganova, M.Y.Filimonov, AIP Conference Proceedings, 1789(2016), 020019. DOI: 10.1063/1.4968440
[8] M.Y.Zemenkova, A.Shabarov, U.Shastunova, A.Kislitsyn, A.Shuvaev, Physical and mathematical modeling of process of frozen ground thawing under hot tank, IOP Conference Series: Materials Science and Engineering, 357(2018), 012007.
DOI: 10.1088/1757-899X/357/1/012007
[9] M.Yu.Filimonov, N.A.Vaganova, Flare Systems Exploitation and Impact on Permafrost, Journal of Physics: Conference Series, 899(2017), 092004.
DOI: 10.1088/1742-6596/899/9/092004
[10] F.Kiselyov, F.Sergeyev, Prediction of construction bases frozen soil temperature development under intense heating, Journal of Physics: Conference Series, 1425(2020), 012208. DOI: 10.1088/1742-6596/1425/1/012208
[11] N.A.Vaganova, M.Yu.Filimonov, Simulation of Cooling Devices and Effect for Thermal Stabilization of Soil in a Cryolithozone with Anthropogenic Impact, Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 11386(2019), 580-587. DOI: 10.1007/978-3-030-11539-5_68
[12] A.A.Samarsky, P.N.Vabishchevich, Computational Heat Transfer, Volume 2, The Finite Difference Methodology, New York, Chichester, Wiley, 1995.
[13] M.Yu.Filimonov, N.A.Vaganova, On Boundary Conditions Setting for Numerical Simulation of Thermal Fields Propagation in Permafrost Soils, CEUR-WS Proceedings, 2109(2018), 18-24.
[14] N.N.Yanenko, The method of fractional steps (The solution of problems of mathematical physics in several variables), Berlin, Springer-Verl., 1971.
[15] S.V.Patankar, Numerical Heat Transfer and Fluid Flow, New York , Hemisphere, 1980.
[16] A.A.Samarskii, B.D.Moiseyenko, An economic continuous calculation scheme for the Stefan multidimensional problem, it USSR Computational Mathematics and Mathematical Physics(1965), 5(1965), 43-58 (in Russian).
[17] N.A.Vaganova, Existence of a solution of an initial-boundary value difference problem for a linear heat equation with a nonlinear boundary condition, Proceedings of the Steklov Institute of Mathematics, 261(2008), no. 1, 260-71. DOI: 10.1134/S0081543808050209
Таяние вечной мерзлоты при эксплуатации скважин Северо-Мукеркамыльского нефтегазового месторождения
Михаил Ю. Филимонов Наталия А. Ваганова
Институт математики и механики им. Н. Н. Красовского УрО РАН
Екатеринбург, Российская Федерация Уральский федеральный университет Екатеринбург, Российская Федерация
Аннотация. Таяние многолетнемерзлых пород в результате изменения климата или различных технологических воздействий будет сопровождаться оседанием земной поверхности и развитием опасных геологических процессов вечной мерзлоты, называемых термокарстом, что приведет к авариям, которые могут вызвать разрушение скважин. В настоящее время ограничены инвестиционные программы освоения новых северных нефтегазовых месторождений. В связи с этим удешевление разработки нефтегазовых месторождений является актуальной проблемой. Например, уменьшение площади кустовых площадок и использование оптимальной эксплуатации северных нефтегазовых месторождениях могут значительно снизить затраты, в частности, на стадии проектирования. Разработана модель распространения нестационарных тепловых полей в мерзлых грунтах от скважин для Северо-Мукеркамыльского нефтегазового месторождения с учетом особенностей конструкции скважин, затрубного пространства и сложной литологии грунта, окружающего скважину. Планируется учесть климатические и технологические факторы, в частности ежегодный перерыв в работе скважин, который составляет от нескольких часов до двух недель. В статье обсуждаются вычислительные особенности расчета тепловых полей в мерзлых грунтах из скважин и исследуется влияние различных параметров, которые в расчетах могут привести к значительному увеличению площади оттаивания на кустовых площадках.
Ключевые слова: математическое моделирование, тепломассоперенос, многолетнемерзлые породы.
