A mathematical model, numerical algorithm and program code for simulation and long-term forecasting of changes in permafrost as a result of operation of a multiple well pad of 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, thermal stabilization of soil, possible insulation of the objects, seasonal uctuations in air temperature, and freezing and thawing of the upper soil layer. A parallel algorithm of decomposition with splitting by spatial variables is presented.
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 elds 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 e ect 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 e ects in permafrost, and due to thermal stabilization of soil, to increase operational safety of oil and gas elds, 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 signi cantly reduces the list of initial parameters. This \anchor" is contained in the nonlinear boundary conditions and allows to simulate natural thermal elds 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 eld and showed a good agreement (di erence 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 signi cantly 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 elds for a are system [12] takes up to 10 hours of computing. 2
Problem statement and mathematical model 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. α q + b(Tair − T z=0) = εσ (T 4 z=0 − Ta4ir ) + λ ∂T
∂z z=0 z solar radiation x ∂∂Ty = 0 y x u lfr a e n li ∂∂Tx = 0 Ωi,i = 1,n insulating shells y iiit ssv m ∂∂Tx = 0 e Ω
∂∂Ty = 0 ∂∂Tz = 0
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 and possible devices.
Simulation of processes of heat distribution is reduced to solution of threedimensional 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 c (T ) + k (T
T ) with initial condition
T (0; x; y; z) = T0(x; y; z):
Here
is density [kg=m3], T is temperature of phase transition [K],
c (T ) =
cc21((xx;; yy;; zz));; TT >< TT ;; is speci c heat [J/kg K];
(
21((xx;; yy;; zz));; TT >< TT ;; is thermal conductivity coe cient [W/m K ], k = k(x; y; z) is speci c heat of phase transition, is Dirac delta function.
Balance of heat uxes at the surface z = 0 de nes the corresponding nonlinear boundary conditions q + b(Tair
T (x; y; 0; t)) = " (T 4(x; y; 0; t)
Ta4ir) +
:
(
To determine the parameters in boundary condition (
In condition (
The others parameters in condition (
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 g. 1. These surfaces are inner boundaries with conditions
T
= Ti(t); i = 1; : : : ; n: i
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 de ned by positive numbers Lx, Ly, Lz: Lx x Lx, Ly y Ly, Lz z 0.
At the boundaries of the computational domain the boundary conditions are given =
In (
Among the mathematical models, which are closer to the considered (
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 (
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,
(
2.
y 3.
z the method of front straightening, the method of smoothing coe cients 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 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 [5] 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.
I II
processors 1 2
N 0
With using these ideas [5, 6], to solve problem (
Approaches to parallelization
On the base of ideas in [6] a nite di erence method is used with splitting by the
spatial variables in three-dimensional domain to solve the problem (
In Fig.3 the stencil of the scheme is presented. The scheme is divided into 3 steps: sweeping by x-variable with xed 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 di erent 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. 4
Numerical results
In Fig.5 thermal elds 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 in uence of wells is enhanced. 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 producing wells. Numerical calculations based on the model for the arrangement of well pads can improve safety and e ciency of northern oil elds due to optimal location of wells and other technical systems in the area and provides signi cant economic e ect 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.