Parallel Splitting and Decomposition Method for
Computations of Heat Distribution in
Permafrost
Nataliia Vaganova1 and Mikhail Filimonov1,2
1
IMM UB RAS, Yekaterinburg, Russia
2
Ural Federal University, Yekaterinburg, Russia
{vna,fmy}@imm.uran.ru
Abstract. 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
field are presented. In the model the most significant climatic and phys-
ical factors are taken into account such as solar radiation, determined
by specific geographical location, heterogeneous structure of frozen soil,
thermal stabilization of soil, possible insulation of the objects, seasonal
fluctuations in air temperature, and freezing and thawing of the upper
soil layer. A parallel algorithm of decomposition with splitting by spatial
variables is presented.
Keywords: parallel computations · splitting by spatial variables · do-
main decomposition · simulation of heat distribution
1 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 lo-
cated in areas of permafrost, as well as on Arctic shelf. New methods for simu-
lation 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 com-
bined 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 al-
gorithms to “anchor” geographical coordinates of the area, as well as a climatic
database developed on the base of open databases NASA, and significantly re-
duces the list of initial parameters. This “anchor” is contained in the nonlinear
Computations of Heat Distribution in Permafrost 43
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 con-
sidered model and the practice.
The program code has to be oriented to carry out high-performance compu-
tations 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 computa-
tional 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.
2 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.
∂T
α q + b(Tair − T z =0 ) = εσ (T 4 z =0 − Tair4 ) + λ
∂z z =0
emissivity
linear flux
solar radiation
∂T
=0
∂x
∂T
=0 ∂T
∂y Ω =0
x y ∂y
z ∂T
=0
∂T ∂z
=0 insulating shells
∂x Ω i , i = 1, n
Fig. 1. A computational area with boundary conditions
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.
44 Nataliia Vaganova and Mikhail Filimonov
(a) (b)
Fig. 2. Intensity of solar radiation (solid) q and average air temperature (dashed) Tair
for the considered location (a). Soil temperature (b).
Simulation of processes of heat distribution is reduced to solution of three-
dimensional 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
∂T
ρ cν (T ) + kδ(T − T ∗ ) = ∇ (λ(T )∇T ), (1)
∂t
with initial condition
T (0, x, y, z) = T0 (x, y, z). (2)
Here ρ is density [kg/m3 ], T ∗ is temperature of phase transition [K],
c1 (x, y, z), T < T ∗ ,
cν (T ) = is specific heat [J/kg K],
c2 (x, y, z), T > T ∗ ,
λ1 (x, y, z), T < T ∗ ,
λ(T ) = is thermal conductivity coefficient [W/m K ],
λ2 (x, y, z), T > T ∗ ,
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
∂T (x, y, 0, t)
γq + b(Tair − T (x, y, 0, t)) = εσ(T 4 (x, y, 0, t) − Tair
4
)+λ . (3)
∂z
To determine the parameters in boundary condition (3), an iterative algo-
rithm is developed that takes into account the geographic coordinates of consid-
ered 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.
Computations of Heat Distribution in Permafrost 45
The others parameters in condition (3) are determined as a result of geo-
physical 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 solu-
tion 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 (founda-
tions, 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 = Ti (t), i = 1, . . . , n. (4)
Ω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 defined 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
∂T ∂T ∂T
= = 0, = γ. (5)
∂x x=±Lx ∂y y=±Ly ∂z z=−Lz
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,
46 Nataliia Vaganova and Mikhail Filimonov
z
y
splitting
x by spatial
x z
variables
1.
y 2. 3.
Fig. 3. Stencil splitting by spatial variables: 1 — sweeping by x, 2 — by y, 3 — by z
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 meth-
ods is that the difference schemes are constructed with explicit separation of
the front of phase transformation. It should be noted that the methods with ex-
plicit 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 ho-
mogeneous 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.
processors
I 1 2 N
II
0
Fig. 4. Two basic parts of the parallel algorithm: I. Sweeping by X and Y . II. Sweeping
by Z.
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
Computations of Heat Distribution in Permafrost 47
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].
3 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 one-
dimensional 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.
4 Numerical results
5 5 5
T: -5-3-1 1 3 5 T: -5-3-1 1 3 5 T: -5-3-1 1 3 5
0 0 0
3.0
-1.0 0.0 2.0 2.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 3.0 0.0
1.0 -1.0 -1.0 2.0 -1.0 -1.0 2.0 -1.0
0.0
0.0
-5 -5 -5
0.0
1.0
1.0
0.0
0.0
1.0
-1.0
0.0
-1.0 -1.0
0.0
-1.0 -1.0 -1.0
0.0
-10 0.0 -10 -10
0.0
0.0
0.0
1.0
0.0
-15 -15 0.0 -15 2.0
1.0
3.0
1.0
1.0
z
z
z
3.0 2.0 4.0
5.0
5.0
3.0
-20 0.0 -20 4.0 -20
4.0
1.0
1.0
0.0
5.0
1.0
4.0
0.0
1.0
5.0
3.0
1.0
0.0
3.0 5.0
5.0
2.0 3.0
-25 -25 -25
3.0
0.0
1.0
2.0
0.0
2.0
4.0
5.0
-30 -30 -30
5.0
5.0
3.0
2.0
4.0
-35 -35 -35
1.0
0.0
5.0
4.0
0.0
5.0
0.0
0.0
3.0
5.0
-40 -40 -40
0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50
x x x
(a) (b) (c)
Fig. 5. Two wells and computed thermal fields for 3 (a), 5 (b), and 10 (c) years of
operating
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.
48 Nataliia Vaganova and Mikhail Filimonov
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.
Fig. 6. Computational time of 1 year of wells operation for different numbers of pro-
cessors and grid nodes
5 Conclusion
The developed mathematical model allows to take into account the most signifi-
cant 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, high-
performance computational and information technologies and applications” (prj.
Computations of Heat Distribution in Permafrost 49
15-7-1-13), and by Act 211 Government of the Russian Federation, contract
N 02.A03.21.0006.
References
1. Filimonov, M.Yu., Vaganova, N.A.: Simulation of thermal fields in the permafrost
with seasonal cooling devices, 9th International Pipeline Conference(IPC 2012).
ASME Conference Proceedings. 4, 133-141 (1999)
2. Filimonov, M., Vaganova, N.: Prediction of changes in permafrost as a result
technogenic effects and climate. Academic Journal of Science. 3(1), 121–128 (2014)
3. Filimonov, M.Yu, Vaganova, N.A.: Simulation of thermal stabilization of soil
around various technical systems operating in permafrost. Applied Mathematical
Sciences. 7(141-144), 7151–7160 (2013)
4. Vaganova, N.A., Filimonov, M.Yu.: Simulation of Engineering Systems in Per-
mafrost. Vestnik Novosibirskogo Gosudarstvennogo Universiteta. Seriya Matem-
atika, Mekhanika, Informatika. 13(4), 37–42 (2013)
5. Samarsky, A.A., Moiseenko, B.D.: Economic scheme of the level set method for the
multidimensional Stefan problem. Computational Mathematics and Mathematical
Physics. 5(5), 816–827 (1965)
6. Samarsky, A.A., Vabishchevich, P.N.: Computational Heat Transfer, Volume 2,
The Finite Difference Methodology (1995)
7. Vaganova, N.: Mathematical model of testing of pipeline integrity by thermal fields.
AIP Conference Proceedings (AMEE’14). 1631(1), 37–41 (2014)
8. Bashurov, V.V., Vaganova, N.A., Filimonov, M.Y.: Numerical simulation of ther-
mal conductivity processes with fluid filtration in soil. Journal of Computational
technologies. 16 (4), 3–18 (2011)
9. Zhang, Y., Chen, W., Cihlar, J.: 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)
10. Vaganova, N.A.: Existence of a solution of an initial-boundary value difference
problem for a linear heat equation with a nonlinear boundary condition. Proceed-
ings of the Steklov Institute of Mathematics. S1, S260–S272 (2008)
11. Filimonov, M.Yu., Vaganova, N.A.: Simulation and Numerical Investigation of
Temperature Fields in an Open Geothermal System. Lecture Notes in Computer
Science: Finite Difference Methods, Theory and Applications. 9045, 393–399 (2015)
12. Filimonov, M.Yu., Vaganova, N.A. Simulation of thermal effect of a vertical flare
system on permafrost. Ecology. Economy. Informatics. Systems analysis and math-
ematical modeling of ecological and economic systems. Collection of articles: 1,
292–299 (2015)