Introduction

According to the papers [1,2] a mathematical model is suggested for long-time forecasting of impacts of development and exploitation of oil and gas fields located in areas of permafrost, as well as on Arctic shelf. New methods for simulation and studying of permafrost degradation in well pads areas which related with permafrost heating by various technical system are developed, taking into account climate changes, solar radiation, as well as taking into account a combined effect of all engineering facilities and technical systems, which are located on well pads. An optimal arrangement of objects on a well pad allows to minimize the temperature effects in permafrost, and due to thermal stabilization of soil, to increase operational safety of oil and gas fields, and to considerably reduce the costs and risks. Simulated problems take into account a number of climatic, natural and man-made factors forming long-term prognosis of degradation of permafrost.

The computer realization of this original methods also uses an numerical algorithms to "anchor" geographical coordinates of the area, as well as a climatic database developed on the base of open databases NASA, and significantly reduces the list of initial parameters. This "anchor" is contained in the nonlinear boundary conditions and allows to simulate natural thermal fields related with seasonal changes in the upper layer of soil. The novelty of parameters adaptation allowed to compare of numerical and experimental data obtained for "Russkoye" oil field and showed a good agreement (difference is about 5%) between the considered model and the practice.

The program code has to be oriented to carry out high-performance computations because of long-time period to be simulated. The computational system is a hybrid computer cluster "Uran" with MPI. Note that complete simulation of all technical systems that located in a well pad makes it necessary to solve such problems in a significantly larger area with three-dimensional computational grid, resulting an essential increasing time of computations (up to 100 hours). For example, a detailed simulation of thermal fields for a flare system [12] takes up to 10 hours of computing.

Problem statement and mathematical model

Simulation of unsteady three-dimensional thermal fields, such as oil and gas fields (the well pads) located in the area of permafrost, is required to take into account the different climatic, physical and technological factors. The first 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, flare systems, pipelines, foundations of buildings, etc. In addition, it is necessary to take into account parameters of used thermal insulation and possible devices. Simulation of processes of heat distribution is reduced to solution of threedimensional diffusivity equation with non-uniform coefficients including localized heat of phase transition -an approach to solve the problem of Stefan type, without the explicit separation of the phase transition in Ω (fig. 1). The equation has the form

ρ c ν (T ) + kδ(T − T * ) ∂T ∂t = ∇ (λ(T )∇T ),(1)

with initial condition

T (0, x, y, z) = T 0 (x, y, z).(2)

Here

ρ is density [kg/m 3 ], T * is temperature of phase transition [K], c ν (T ) = c 1 (x, y, z), T < T * , c 2 (x, y, z), T > T * , is specific heat [J/kg K], λ(T ) = λ 1 (x, y, z), T < T * , λ 2 (x, y, z), T > T * , is thermal conductivity coefficient [W/m K ], k = k(x, y, z) is specific heat of phase transition, δ is Dirac delta function.

Balance of heat fluxes at the surface z = 0 defines the corresponding nonlinear boundary conditions

γq + b(T air − T (x, y, 0, t)) = εσ(T 4 (x, y, 0, t) − T 4 air ) + λ ∂T (x, y, 0, t) ∂z .(3)

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 NASA climate data. Fig. 2a shows the data for the considered field.

The others parameters in condition (3) are determined as a result of geophysical research of oil and gas field. Fig. 2b shows temperature distribution in an exploratory well. Applying the developed iterative algorithm [3,4] to define 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 different climate and natural features of the considered geographical location.

Let n objects(technical systems) be included in Ω which are heat (foundations, producing insulated wells, pipelines) or cold (SCDs) sources. The surfaces of these objects are Ω i (x, y, z), i = 1, . . . , n in fig. 1. These surfaces are inner boundaries with conditions

T Ωi = T i (t), i = 1, . . . , n.(4)

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. We assume that the size of the box Ω is defined by positive numbers L x , L y , L z :

−L x ≤ x ≤ L x , −L y ≤ y ≤ L y , −L z ≤ z ≤ 0.

At the boundaries of the computational domain the boundary conditions are given

∂T ∂x x=±Lx = ∂T ∂y y=±Ly = 0, ∂T ∂z z=−Lz = γ.(5)

In ( 5) γ is a positive number, corresponding to a geothermal flux value. As a rule γ is a small number and it is possible to be set zero in calculations.

Among the mathematical models, which are closer to the considered (1)-( 5), we have to mention the works of researches from USA and Canada in which there is used one-dimensional heat equation, and take into account various factors: snow cover, vegetation, etc. (see review in [9]). It is assumed that there are no engineering systems located in the permafrost zone. Taking into account solar energy it was shown that short-wave part of the radiation can penetrate into the thick snow into a considerable depth, ranging in depth by the Bouguer-Lambert law. In the proposed three-dimensional model snow cover, vegetation and other factors are taken into account by a special iterative algorithm variating some coefficients in nonlinear boundary conditions on the soil surface. This approach allows to user to simplify the task of initial data setting, for example, it is not necessary to know thickness of snow cover, changes in the thermal properties of snow, depending on the solar radiation, etc.

Numerical methods of solving problems are the most effective and universal method of research for models considered in this paper. A large number of works is devoted to development of difference methods for solving boundary value problems for the heat equation To solve (1)-( 5) a finite-difference method is used.

At present there are the following difference methods for solving Stefan type problems: the method of front localization by the difference grid node, the method of front straightening, the method of smoothing coefficients and schemas of through computation [6]. 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 difference 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 [5] an effective scheme of through computations is developed with smoothing of discontinuous coefficients 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 difference 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 specific heat ratio. Thus obtained discontinuous function then "shared" with respect to temperature, and does not depend on the number of measurements and phases. With using these ideas [5,6], to solve problem (1)-( 5) in three-dimensional box a finite difference method is used. Solvability of the same difference problems approximating (1)-( 5) is proved in [8,10,11] in the case of thermal traces of of underground pipelines without phase transition in soil [7].

Approaches to parallelization

On the base of ideas in [6] a finite difference method is used with splitting by the spatial variables in three-dimensional domain to solve the problem (1)-( 5). We construct an orthogonal grid, uniform, or condensing near the ground surface or to the surfaces of Ω i . The original equation is approximated by an additive onedimensional implicit central-difference scheme and a three-point sweep method to solve a system of linear differential algebraic equations is used.

In Fig. 3 the stencil of the scheme is presented. The scheme is divided into 3 steps: sweeping by x-variable with fixed y and z, sweeping by y, and sweeping by z. These three steps are successively carried out, but it is possible to compute it in different grid lines simultaneously so to perform a decomposition of Ω with no overlapping.

Fig. 4 shows two basic steps of the computational algorithm. Ist step is parallel sweeping by x and y on N processors, the "zero" processor works to read initial and upper boundary parameters and to compute sweeping by z. The processors are exchanged by the values of temperatures and use these to compute the next sweeping step. In Fig. 5 thermal fields from two heated wells are shown for 3, 5, and 10 years of exploitation. The temperature of the wells re assumed to be 45 • C, permafrost temperature is −1 • C. In upper layers there are seasonal melting of frozen soil. The melted zones around the wells merge and raise so the influence of wells is enhanced. Fig. 6 shows the computational time of an 1 year of two wells exploitation in a domain by using a series of grids with 201x151x151, 401x301x301, and 401x301x601 nodes. The used computational system is a hybrid computer cluster "Uran" with MPI. The computations are carried out on different numbers of processors. The time step of numerical algorithm is 24 hours. To compute 1 year with using single processor in 201x151x151 nodes more than 2 and a half hours are required (the diamond point in the figure). The squares denote the computational time of the same task by using the presented parallel algorithm. There are some increasing of computational time due to non-balanced loading of the processors. When the numbers of the grid nodes is enlarge, then the computational times is stabilised. The acceleration is not so big as expected and it is related with the domain splitting by Z. The investigations and the presented algorithm optimization will be continued.

Numerical results

Conclusion

The developed mathematical model allows to take into account the most significant physical and climatic factors influencing on formation of temperature fields in permafrost during operation producing wells. Numerical calculations based on the model for the arrangement of well pads can improve safety and efficiency of northern oil fields due to optimal location of wells and other technical systems in the area and provides significant economic effect already at the design stage. The suggested approach of splitting and decomposition allows to use distributed and parallel computations and, as a result, essentially increase complexity and detailed elaboration of the objects to be simulated. This work was supported by Russian Foundation for Basic Research 14-01-00155, by the Program of UB RAS "Mathematical models, algorithms, highperformance computational and information technologies and applications" (prj. 15-7-1-13), and by Act 211 Government of the Russian Federation, contract N 02.A03.21.0006.