A mathematical model, numerical algorithm and program code for simulation thermal changes in permafrost as a result of operation of a horizontal °are system in northern oil and gas ¯eld are presented. In the model the most signi¯cant climatic and physical factors are taken into account such as solar radiation, determined by speci¯c geographical location, heterogeneous structure of frozen soil, seasonal °uctuations in air temperature, and freezing and thawing of the upper soil layer. Results of computing are presented.
the surface heating is stronger, and on the upper boundary we may assume a ¯xed temperature.
Optimization of exploitation of productive wells in permafrost area are considered [5, 6]. Operating modes of production wells are selected in such a way to reduce the thermal e®ect of hot oil on the permafrost, and in case of shut-in reverse to slow freezing of engineering constructions. In this paper short and long continued operation of the °are system and simulated the e®ect of di®erent thermal cycles on the ground are considered. 2
Simulation of unsteady three-dimensional thermal ¯elds, such as oil and gas ¯elds (the well pads) located in the area of permafrost, is required to take into account the di®erent climatic, physical and technological factors.
The ¯rst group of factors is related with solar radiation, seasonal changes in air temperature, resulting a periodic thawing (freezing) of soil, and possible snow layer. The second group factors includes parameters of soil: thermal, dependent with humidity, structure and temperature. The third group of factors are the possible source of heat as production and injection wells, °are systems, pipelines, foundations of buildings, etc. In addition, it is necessary to take into account parameters of used thermal insulation [7] and possible devices used for thermal stabilization (cooling) of the soil such as seasonal cooling devices, as in [5].
We consider a horizontal °are, which is simulated by a heat source on the surface of the permafrost soil (¯g. 1). The initial temperature distribution in the soil is presented in ¯g. 2. The basic thermal parameters of soil are presented in Table 1. In ¯g. 3 monthly averages air temperature and solar radiation through a year used for simulation the annual temperature cycle in the soil are shown.
y z x thermal exchange flare riprap permafrost soil (a) 20С 70С 400 (b) 27х18 m 12х18 m 6х9 m 2х3m 2.1
Simulation of processes of heat distribution is reduced to solution of three-dimensional di®usivity equation with non-uniform coe±cients including localized heat of phase transition | an approach to solve the problem of Stefan type, without the explicit separation of the phase transition in (¯g. 1). The equation has the form with initial condition
Here ½ is density [kg=m3], T ¤ is temperature of phase transition [K], ½¡cº (T ) + k±(T ¡ T ¤)¢ @T = r (¸(T )rT );
@t
T (0; x; y; z) = T0(x; y; z): k = k(x; y; z) is speci¯c heat of phase transition, ± is Dirac delta function.
cº (T ) = ½ c1(x; y; z); T < T ¤; c2(x; y; z); T > T ¤;
is speci¯c heat [J/kg K]; ¸(T ) = ½ ¸1(x; y; z); T < T ¤; ¸2(x; y; z); T > T ¤;
is thermal conductivity coe±cient [W/m K ], 1. 2. 3. 4. 5.
Layers of riprap and soil concrete slab 0-0.15 m. breakstone 0.15-0.25 m. sand 0.25-2.35 m. soil 2.35-10.35 . soil 10.35-40.0 .
Thermal conductivity, W/(m K) frozen melted 1.93 1.57
kJ/(m3 K) frozen melted 2150.0 3490.0
transition, kJ/(m3 K) 141504.0
Temperature
of phaze transition, C 0.0 1.75 1.83 1.89 2.21 1.63 1.39 1.68 1.75 2160.0 1536.0 2200.0 2350.0 2630.0 1983.0 2780.0 2750.0 The computational domain is a three-dimensional box , where x and y axes are parallel to the ground surface and the z axis is directed downward (¯g 1a). We assume that the size of the box is de¯ned by positive numbers Lx, Ly, Lz: ¡Lx · x · Lx, ¡Ly · y · Ly, ¡Lz · z · 0.
The °are platform is simulated by a region of ¯xed temperature in a upper boundary (¯g 1b). Balance of heat °uxes at the surface z = 0 outside of the °are platform is de¯ned by corresponding nonlinear boundary conditions @T (x; y; 0; t) °q + b(Tair ¡ T (x; y; 0; t)) = "¾(T 4(x; y; 0; t) ¡ Ta4ir) + ¸ : (3) @z
Nonlinear boundary conditions of fourth degree is often used for simulations of process where there is a heat exchange as solar radiation or other type of heatet surfaces interaction, for example, in [8].
To determine the parameters in boundary condition (3), an iterative algorithm is developed that takes into account the geographic coordinates of considered area, lithology of soil and other features of the selected location.
In condition (3) values of intensity of solar radiation and seasonal changes in air temperature are obtained by weather stations or on the base of an open climate data for a °are system to be simulated location. Fig. 2 shows the data for the considered ¯eld.
The others parameters in condition (3) are determined as a result of geophysical research of oil and gas ¯eld. Fig. 2 shows temperature distribution in an exploratory well. Applying the developed iterative algorithm [5, 7] to de¯ne some of the parameters in nonlinear boundary condition (3) it is possible to identify them so that the temperature distribution in the soil found as a solution of equation (1){(3) to be periodically repeated over the next few years, that allows to implicitly take into account di®erent climate and natural features of the considered geographical location.
At the boundaries of the computational domain the boundary conditions are given = (4)
In (4) ° is a positive number, corresponding to a geothermal °ux value. As a rule ° is a small number and it is possible to be set zero in calculations. 2.3
Numerical methods of solving problems are the most e®ective and universal method of research for models considered in this paper. A large number of works is devoted to development of di®erence methods for solving boundary value problems for the heat equation To solve (1){(4) a ¯nite{di®erence method is used.
At present there are the following di®erence methods for solving Stefan type problems: the method of front localization by the di®erence grid node, the method of front straightening, the method of smoothing coe±cients and schemas of through computation [9]. The method of front localization in the mesh node is used only for one-dimensional single-front problems and method of front straightening for the multi-front problems. A basic feature of these methods is that the di®erence schemes are constructed with explicit separation of the front of phase transformation. It should be noted that the methods with explicit separation of unknown boundary of the phase transformation for the case of cyclic temperature changes on the boundary are not suitable, because the number of non-monotonically moving fronts may be more than one, and some of them may merge with each other or disappear.
In [10] an e®ective scheme of through computations is developed with smoothing of discontinuous coe±cients in the equation of thermal conductivity by temperature in the neighborhood of the phase transformation. Through calculation scheme is characterized by that the boundary of phase separation is explicitly not allocated, and the homogeneous di®erence schemes may be used. The heat of phase transformation is introduced with using the Dirac ±-function as a concentrated heat of phase transition in the speci¯c heat ratio. Thus obtained discontinuous function then \shared" with respect to temperature, and does not depend on the number of measurements and phases. Collocation and least residuals method is also used for such equations [11].
With using these ideas [9, 10], to solve problem (1){(4) in three-dimensional box a ¯nite di®erence method is used with splitting by the spatial variables and taking into account the inner boundaries i. Solvability of the same di®erence problems approximating (1){(4) is proved in [12, 13]. 1 0 -1 -2 -3 -4 z-5 -6 -7 -8 -9 -10 -210 -3-1 -3 -2 -1 0 1 2 3 103070 0 -2-3 -3
-3 40 x march 1 0 3 1 2 -1 0 -2 -3 -2 -4 z-5 -6 -7 -8 -9 -10 30 -1 40 x june 1 0 -1
The developed mathematical model allows to take into account the most signi¯cant physical and climatic factors
in°uencing on formation of temperature ¯elds in permafrost during operation of °are system [14]. The algorithm
used to determination of annual changing of upper boundary parameters allows to decrease amount of initial
data and to expect the program code to be presented in remote and clouds simulations [
Let consider a 3D computational domain with Lx = 87m, Ly = 87m, Lz = 40m. The °are system in in center of the soil surface (¯g 1). The °are platform has 3 layers: 0.15 m. of concrete slab, 0.10 m. of breakstone, and 2.10 m. of sand. The permafrost is combined of two di®erent frozen soils. Parameters of soils and riprap layers are presented in table 1. We will compare 4 types of °are system operation mode (table 2). Let assume that the month of starting operation is may.
In ¯gures 4{7 temperature ¯eld are shown in section by y axis near the °are platform for di®erent operating modes. In ¯gures 8 and 9 annual temperature changing are presented at the depth of 1m and 3m. In ¯gures 10 and 11 the detalization of temperatures under °are platform are shown for summer and winter. The \ground" line in the ¯gures denotes the natural temerature disribution without in°uence of °are.
E®ect of heating for I operating mode is not too intense, but is permanent, and its in°uence is presented even at the depth of 3m (¯g 8). Operating mode II is obviously has a signi¯cant e®ect on the temperature distribution under the °are platform, the deepest soil thawing is observed in this mode (¯g 5). When the °are heated for a month (operating mode III), despite the fact that high peak temperature is observed (¯g 9), but signi¯cant period of rest allows ground to be cooled up to near background values of temperature (¯g 10 and On the base of a mathematical model of thermal ¯eld in a nonuniform frozen soil from a heat source located on the surface of the ground, a horizontal °are system used to gas °aring in northern oil and gas ¯elds is simulated. Forecasts of permafrost degradation around this system is obtained for di®erent operating modes. The purpose of simulations is making recommendations on the optimal choice of the structure of heat-insulating riprap layers to reduce the thermal e®ects on the permafrost. By computer simulation an acceptable depth is determined, which may be used for heat insulating materials with limited operating temperature range (e.g., for penoplex the temperature has to be no greater than 70oC). Also the features of thermal ¯eld propagation in the soil are considered in detail, taking into account di®erent physical and climatic factors. Operating modes, when the °are platform time could be possible to cool down without heat accumulation, are preferable.
This work was supported by Russian Foundation for Basic Research projects (16-01-00401 and 14-01-20141) and Program of UB RAS (project 15-7-1-13). [1] Y. Zhang, W. Chen, J. Cihlar. A process-based model for quantifying the impact of climate change on permafrost thermal regimes. Journal of Geophysical Research, 108(D22): ACL5-1{ACL5-16, 2003. [2] M. Filimonov, N. Vaganova. Prediction of changes in permafrost as a result technogenic e®ects and climate.
Academic Journal of Science, 3(1): 121{128, 2014. [3] M. Y. Filimonov, N. A. Vaganova. Simulation of technogenic and climatic in°uences in permafrost for northern oil ¯elds exploitation. Lecture Notes in Computer Science, 9045:185{192, 2015. march [4] M. Yu. Filimonov, N. A. Vaganova. Simulation of thermal e®ect of a vertical °are system on permafrost.
Ecology. Economy. Informatics. Systems analysis and mathematical modeling of ecological and economic systems. Collection of articles, 1:292{299, 2015. [5] M. Yu. Filimonov, N. A. Vaganova. Simulation of thermal stabilization of soil around various technical systems operating in permafrost. Applied Mathematical Sciences, 7(141-144): 7151{7160, 2013. [6] N. A. Vaganova, M. Yu. Filimonov. Computer simulation of nonstationary thermal ¯elds in design and operation of northern oil and gas ¯elds. AIP Conf. Proc., 1690:020016, 2015. [7] N. A. Vaganova, M. Yu. Filimonov. Simulation of Engineering Systems in Permafrost. Vestnik Novosibirskogo
Gosudarstvennogo Universiteta. Seriya Matematika, Mekhanika, Informatika, 13(4): 37{42, 2013. [8] A. N. Cherepanov, V. P. Shapeev, V. I. Isaev. Simulation of Heat Transfer Processes in Laser Welding of
Dissimilar Metals with an Insert. High Temperature, 53(6):841{846, 2015. [9] A. A. Samarsky, P. N. Vabishchevich. Computational Heat Transfer, Volume 2, The Finite Di®erence
Methodology. Wiley, October 1995. [10] A. A. Samarsky, B. D. Moiseenko. Economic scheme of the level set method for the multidimensional Stefan problem. Computational Mathematics and Mathematical Physics, 5(5): 816{827, 1965. [11] V. Shapeev. Collocation and Least Residuals Method and Its Applications. EPJ Web of Conferences, 108(01009):1{12, 2016. equation with a nonlinear boundary condition. Proceedings of the Steklov Institute of Mathematics, S1: to describing natural and technogenic processes in soils in northern ¯elds well pads.
Proceedings of XI All-Russian Congress on Fundamental Problems of Theoretical and Applied Mechanics. Izhevsk: 3917{3919, 2015.