=Paper=
{{Paper
|id=Vol-2543/spaper06
|storemode=property
|title=Supercomputer Technologies for Long-term Modeling of Permafrost Changes
|pdfUrl=https://ceur-ws.org/Vol-2543/spaper06.pdf
|volume=Vol-2543
|authors=Mikhail Filimonov,Nataliya Vaganova,Elena Akimova,Vladimir Misilov
|dblpUrl=https://dblp.org/rec/conf/ssi/FilimonovVAM19
}}
==Supercomputer Technologies for Long-term Modeling of Permafrost Changes==
https://ceur-ws.org/Vol-2543/spaper06.pdf
Supercomputer Technologies for Long-term Modeling
of Permafrost Changes
M. Yu. Filimonov1,2[0000-1111-2222-3333], N. A. Vaganova1,2[0000-0001-6966-9050],
E. N. Akimova1,2[0000-0002-4462-5817] and V. E. Misilov1,2[0000-0002-5565-0583]
1 Ural Federal University, 19 Mira street, Yekaterinburg, Russia
2 Krasovskii Institute of Mathematics and Mechanics of UrB RAS, 16 S.Kovalevskaya
Str., Yekaterinburg, Russia
fmy@imm.uran.ru
Abstract. A model of propagation of thermal fields in permafrost from various
engineering objects operating in Arctic regions is considered. The proposed
model includes the most significant technical and climatic parameters affecting
the formation of thermal fields in the surface layer of the soil. The main objec-
tive of the study is a long-term forecasting of changes in the dynamics of per-
mafrost boundaries during operation of cluster sites of northern oil and gas
fields. Such a forecast is obtained by simulation of complex system consisting
of heat or cold sources and frozen soil, thawing of which can lead to the loss of
the bearing capacity and possible technogenic and environmental accidents. For
example, the sources of heat can be production wells, and the sources of cold
can be seasonal cooling devices that are used to stabilize the soil. To minimize
the impact of heat sources on permafrost, various options for thermal insulation
are used, and to preserve the original temperature regime of the top layer of
soil, riprap materials consisting of sand, concrete, foam concrete, or other heat
insulating material are used. The developed set of programs was used in the de-
sign of 12 northern oil and gas fields. To solve the described problem in a com-
plex three-dimensional area, substantial computational resources are required.
The computing time of one variant can often exceed 10–20 hours of machine
time on a supercomputer. To speed up the numerical calculations, multicore
processors are used. Numerical calculations illustrate the possibility of a devel-
oped set of programs for making long-term forecasts for determining changes in
the boundaries of the permafrost zones, and show that on multicore processors
it is possible to achieve acceleration close to the theoretical one.
Keywords: Heat and Mass Transfer, Cryolithozone, Simulation, Parallel Com-
puting, OpenMP.
1 Introduction
The term “permafrost” was introduced into the English literature by S.W. Muller [1]
as an abbreviation of the original Russian term “Permanently frozen ground (soil)”
suggested by M.I. Sumgin [2]. The permafrost zone occupies about 25% of all the
land of the globe [3], [4] and is located in the northern or mountainous regions under
Copyright © 2020 for this paper by its authors.
Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
362
the zone of seasonal thawing of the soil, which is determined by geographical coordi-
nates, the intensity of solar radiation, and other climatic factors. The climate changes
related with, for example, the warming may lead to the essential permafrost degrada-
tion especially in high-latitudes regions, which may influence on the environment of
the globe [4], in particular, due to the release of methane.
The development of the permafrost regions is especially important for Russia,
since about 93% of Russian natural gas and 80% of oil are produced in these areas. At
the same time, the development of these regions leads to negative consequences asso-
ciated with the thawing of permafrost, which leads to the formation of dangerous
geological phenomena called thermokarst. The average thickness of the permafrost
varies from 10 to 800 meters, and the components of the permafrost soils have various
physicochemical properties that can vary in all directions. In summer, as a result of
positive temperatures and solar radiation, seasonal thawing of the upper soil layer
occurs; in winter, the reverse process is observed. Computer simulation of such sea-
sonal processes is described in [5], and, in particular, with using multicore processors
and OpenMP technology in [6]. A significant effect on the formation of thermal fields
in the soil is also caused by anthropogenic conditions, often leading to more extensive
changes in the boundaries of the permafrost [7], [8]. In well pads of northern oil and
gas fields the following technical systems may have the same effect of a massive heat
source: production and injection wells [9], [10], flare systems [11] having some peri-
odical operation conditions [12] and other engineering facilities. For thermal stabiliza-
tion of the soil (including additional cooling or refreezing), the cooling devices
(SCDs) are used that can reduce thermal effects from heat sources on the surrounding
soil and can be used to prepare a construction site in the permafrost distribution zone
by freezing the upper part of the soil [13].
To solve the described problems, mathematical models, numerical methods of so-
lution, and codes are developed for modeling non-stationary thermal fields in the
near-surface soil layer in a complex three-dimensional area. The computational time
of the problems solution, for which a large number of variants with different parame-
ters are required, could exceed tens of hours of computer time on a supercomputer.
The parallel approaches are considered [14], and cloud technologies were also used,
which allow to carry out remote computing to solve certain problems arising from the
placement and operation of northern oil and gas fields [15], [16]. In this work, we
studied the use of multicore processors and OpenMP technology to solve the prob-
lems of long-term forecasting of changes in the thawing boundaries of permafrost
around a production well. The results of numerical calculations and evaluation of the
effectiveness of multicore processors are presented.
2 Problem Statement and Mathematical Model
Let 𝑇 = 𝑇(𝑡, 𝑥, 𝑦, 𝑧) be the soil temperature at the time moment 𝑡, where (𝑥, 𝑦, 𝑧) is
a point of the computational area Ω = {(𝑥, 𝑦, 𝑧): 0 ≤ 𝑥 ≤ 𝐿𝑥 , 0 ≤ 𝑦 ≤ 𝐿𝑦 , −𝐿𝑧 ≤
𝑧 ≤ 0}. The area is a 3D “box” (see Fig. 1), in which the cylinder-shaped well Ω1 and
the cold sources (SDCs) Ω2 are excluded. The axes x and y are parallel to the soil
363
surface, and the axe z is directed vertically down into Ω. In Fig. 1 there presented two
variants of the unit combinations: a single well and a well surrounded by SCDs.
thermal exchange
emissivity
solar radiation
x insulation
y
z
Ω
Ω2
Ω1
Fig. 1. The domain Ω, well Ω1 , and seasonal cooling devices Ω2 .
Thermal conductivity processes are described by the following equation:
T
с (T ) k (T T * ) div (T ) grad T (1)
t
with the initial condition
T (0, x, y, z) T0 ( x, y, z), (2)
where ρ=ρ(x,y,z) is density [kg/m3], T*=T*(x,y,z) is temperature of phase transition,
с1 ( x, y, z ), for T T *
сν(T) = is specific heat [J/(kg K)],
с2 ( x, y, z ), for T T
*
364
1 (x, y, z ), for T T *
λ(T)= is thermal conductivity [W/(m K)], δ is Dirac δ-
2 (x, y, z ), for T T
*
function, k=k(x,y,z) is specific heat of phase transition. The justification of the ap-
plicability of this equation for solving problems of the Stefan type is presented in
[17], [18]. In [19] as a result of the heat balance on the surface, the following bounda-
ry condition is proposed:
T
q(t ) b(Tair (t ) T z 0 ) (T 4 Tair4 (t )) (3)
z z 0
and an algorithm for adaption of the mathematical model to a specific geographic
location is described. In condition (3) 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, σ = 5,67∙10-8W/(m2K4) is Stefan-Boltzmann constant, b=b(t,x,y) is heat
transfer coefficient, ε=ε(t,x,y) is the coefficient of emissivity. The coefficients of heat
transfer and emissivity depend on the type and condition of the soil surface. Total
solar radiation q(t) is the sum of direct solar radiation and diffuse radiation. Soil is
absorbed only a part of the total radiation which equal to αq(t), where α=α(t,x,y) is the
part of energy that is formed to heat the soil, which in general depends on atmospher-
ic conditions, angle of incidence of solar radiation, i.e. latitude and time. The tech-
nical objects Ω1 and Ω2 are additional sources of heat or cold in the permafrost soil.
At these inner boundaries we set the following conditions:
T Ti (t ), i 1, 2. (4)
i
To use the numerical methods it is necessary to set the boundary conditions at the
lateral boundaries of Ω. Let suppose
T T T
0, 0, 0. (5)
x x Lx y y L z z L
z
y
In this case, the computational domain should be chosen large enough to avoid the
influence of the boundary conditions (5) on the thermal fields in the computational
domain Ω created by the objects Ωi.
3 Numerical Calculations and Performance Study
To apply a numerical method of solution of the problem (1)–(5) it is necessary to
evaluate the parameters included into boundary condition (3). In [19], [20] an iterative
algorithm is described to determine these parameters. The choice of these parameters
allows to indirectly take into account the influence of snow cover, climatic, and envi-
ronmental conditions associated with the geographical coordinates of the considered
well pad. In the numerical implementation of the solution of problem (1)–(5), the
finite-difference method is used, which allows to apply the method of splitting into
365
spatial variables for better organization of numerical calculations. Following [17] and
[18] equation (1), for each of the spatial directions, the equation is approximated by
an implicit central-difference three-point scheme, and the system of difference linear
algebraic equations having the tri-diagonal form is solved by the sweep method. On
the ground surface, in view of condition (3), a fourth-degree algebraic equation arises,
for the solution of which the Newton method is used [21]. In the computations, an
orthogonal grid is used, which is uniform or condensing near the soil surface, or the
surfaces Ωi.
0 0
T T
40 40
-2 20 -2 20
10 10
5 5
-4 -4
1 1
0 0
-6 -1 -6 -1
-5 -5
-8 -8
-10 -10
z
z
-12 -12
-14 -14
-16 -16
-18 -18
-20 -20
0 5 10 15 20 0 5 10 15 20
x x
(a) (b)
Fig. 2. The thermal fields after 5 years of operation for the non-insulated well (a) and for the
well with a complex shell (b).
0 0
T T
40 40
-2 20 -2 20
10 10
5 5
-4 -4
1 1
0 0
-6 -1 -6 -1
-5 -5
-8 -8
-10 -10
z
z
-12 -12
-14 -14
-16 -16
-18 -18
-20 -20
0 5 10 15 20 0 5 10 15 20
x x
(a) (b)
Fig. 3. The thermal fields after 5 years of operation for the non-insulated well (a) and for the
well with a complex shell and a system of SCDs (b).
The described algorithm is implemented in the “Wellfrost” certified software pack-
age and in the various modifications, which was tested at 12 northern oil and gas
366
fields. Also the experimental data and the results of numerical calculations was com-
pared in view to determining the boundary of the location of the zero isotherm, which
determines the boundary of the area of soil thawing around the producing well. In
2012 the accuracy of numerical calculations was verified for the Russkoye oil field
(Yamalo-Nenets Autonomous Okrug, Russia), for which the obtained numerical re-
sults are in accordance with the experimental results, and the accuracy reached 5%
after 3 years from the start of the field operation.
In Figures 2 and 3 the calculated thermal fields are shown for 5 years of operation
around the well without heat-insulating shells (a) and with a combined heat-insulating
shell (b), as well as with a system of 8 SCDs.
The permafrost soil temperature is -0.7°С, the temperature of the fluid in the well
is 50°С. Calculations allow to evaluate the influence of heat or cold sources, the zone
of seasonal thawing and freezing of the upper soil layer. To assess the long-term im-
pact, calculations are carried out for a time period of up to 50 years, the time step does
not exceed 24 hours. The calculation time depends both on the size of the computa-
tional grid and on the complexity of the objects introduced into the computational
domain.
Table 1 shows the considered variants of the calculations. Fig. 4 shows a graph of
the change of calculation time for each of the model variants. The calculations were
carried out on a grid with 91x91x51 nodes.
Table 1. Descriptions of the model variants.
No. Description
1 Freezing and thawing of soil without inserted objects.
2 Thawing of the soil around the well without thermal insulation.
3 Thawing of the soil around the well with a simple insulation.
4 Thawing of the soil around the well with a complex insulation.
5 Freezing of the soil with 9 SCDs.
Freezing and thawing of the soil around the well without thermal insulation and 8
6
SCDs.
Freezing and thawing of the soil around the well with a simple insulation and 8
7
SCDs.
Freezing and thawing of the soil around the well with a complex insulation and 8
8
SCDs.
The computing time increases nonlinearly when the number of the included object
grows, as well as, the complexity of the objects, even if the grid size remains the
same.
367
160
Time of computation, s
110
60 T1, s
T6, s
10
1 2 3 4 5 6 7 8
-40
The model variants
Fig. 4. Computational times of 1 year for different variants: sequential program execution time
(red), execution time on 6 cores (blue).
Table 2 contains the computing times of the parallel program for simulating 1 year
on the 6-core AMD Ryzen 5 1600X processor. 𝑇1 is the computing time of the serial
program, 𝑇6 is the computing time of the parallel program running on 6 cores.
Table 2. Computing time of the serial and parallel programs for various models.
No. 𝑇1 [s] 𝑇2 [s] Speedup 𝑆6 Efficiency 𝐸6
1 59 29.7 1.98 0.33
2 56.2 28.8 1.95 0.32
3 59.6 29.5 2.01 0.34
4 64.3 33.6 1.91 0.32
5 80.8 38.1 2.12 0.35
6 81.8 40.5 2.01 0.34
7 103.7 52.4 1.98 0.33
8 153.8 84.9 1.81 0.30
For studying the parallel algorithm running on 𝑛 cores, the speedup Sn T1 and
Tn
efficiency En Sn coefficients were used. The resulted values are consistent with
n
1
the theoretical ones, obtained by the Amdahl’s law Sn 2.2 , where
1
n
0.35 is the serial code proportion.
368
4 Conclusion
The developed models and algorithms make it possible to simulate the propagation of
unsteady thermal fields in frozen ground from producing wells at well pads, taking
into account the possible placement of cooling devices around the wells. The numeri-
cal experiments show the possibility of obtaining a long-term forecast in the dynamics
of permafrost boundaries for various options for well operation using additional tech-
nical systems and heat-insulating materials. The complications and additional con-
structions inserted into the model by taking into account various additional sources of
cold, for example the SCDs, and the use of thermal insulation increase the computa-
tional time. Parallel computations approaches significantly reduced the computational
time in solving such problems. A parallel software package for multicore processors
using OpenMP technology has been developed.
Acknowledgments
This work was supported by the RFBR (project No. 19-07-00435).
References
1. Muller, S.W.: Permafrost or Permanently Frozen Ground and Related Engineering Prob-
lems. Ann Arbor: J.W. Edwards (1947).
2. Sumgin, M. I.: Permafrost in the USSR Boundary. Vladivostok: Far-East Geophysical Ob-
servatory (1927) [in Russian].
3. Nelson, F.E., Anisimov, O.A., Shiklomanov, N.I.: Subsidence risk from thawing perma-
frost. Nature 410, pp. 889–890 (2001).
4. Nelson, F.E., Anisimov, O.A., Shiklomanov, N.I.: Climate Change and Hazard Zonation in
the Circum-Arctic Permafrost Regions. Natural Hazards 26, pp. 203–225 (2002).
5. Akimova, E.N, Filimonov, M.Y., Misilov, V.E., Vaganova, N.A.: Supercomputer model-
ling of thermal stabilization processes of permafrost soils. In: 18th Intern. Conf. Geoin-
formatics: Theoretical and Applied Aspects, Geoinformatics 2019, 15482. Kyiv (2019).
6. Akimova, E. N., Filimonov, M.Yu., Misilov, V.E., Vaganova, N.A.: Simulation of thermal
processes in permafrost: parallel implementation on multicore CPU. In: CEUR Workshop
Proceeding, vol. 2274, pp. 1–9 (2018).
7. Vaganova, N., Filimonov, M.: Simulation of freezing and thawing of soil in Arctic regions.
Conference Series: Earth and Environmental Science 72, 012005 (2017).
8. Vaganova, N.A., Filimonov, M.Y.: Computer simulation of nonstationary thermal fields in
design and operation of northern oil and gas fields. AIP Conference Proceedings 1690
020016 (2015).
9. Vaganova, N.A., Filimonov, M.Yu.: Simulation of Engineering Systems in Permafrost.
Vestnik Novosibirskogo Gosudarstvennogo Universiteta. Seriya Matematika, Mekhanika,
Informatika 13(4), pp. 37–42 (2013) [in Russian].
10. Filimonov, M.Y., Vaganova, N.A.: Simulation of technogenic and climatic influences in
permafrost for northern oil fields exploitation. Lecture Notes in Computer Science (includ-
ing subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)
9045, pp. 185–192 (2015).
369
11. Filimonov, M., Vaganova, N.: Short and Long Scale Regimes of Horizontal Flare System
Exploitation in Permafrost. In: CEUR Workshop Proceedings, vol. 1662, pp. 253–260
(2016).
12. Filimonov, M. Yu., Vaganova, N. A.: Simulation of Influence of Special Regimes of Hori-
zontal Flare Systems on Permafrost. Lecture Notes in Computer Science (including subse-
ries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 11386,
pp. 233–240 (2019).
13. Vaganova, N. A., Filimonov, M. Y.: Simulation of Cooling Devices and Effect for Ther-
mal 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, pp. 580–587 (2019).
14. Vaganova, N., Filimonov, M.: Parallel splitting and decomposition method for computa-
tions of heat distribution in permafrost. In: CEUR Workshop Proceedings, vol. 1513,
pp. 42–49 (2015).
15. Bersenev, A. Iu., Vaganova, N. A., Vasev, P. A., Igumnov, A. S., Filimonov, M.:
Klasternye vychisleniia kak servis na primere zadachi modelirovaniia teplovykh polei ot
skvazhin na severnykh neftegazovykh mestorozhdeniiakh. In.: Nauchnyi servis v seti In-
ternet: mnogoobrazie superkompiuternykh mirov: trudy Mezhdunarodnoi superkompiuter-
noi konferentsii (22-27 sentiabria 2014 g., g. Novorossiisk), pp. 147–151. Izd-vo MGU,
Moscow, Russia (2014) [in Russian].
16. Vaganova, N. A., Vasev, P. A., Gusarova, V. V., Igumnov, S. T., Filimonov, M.: Ispol-
zovanie oblachnykh tekhnologii pri modelirovanii ekspluatatsii severnykh neftegazovykh
mestorozhdenii. In: Trudy IMekh UrO RAN «Problemy mekhaniki i materialovedeniia».
Materialy konferentsii «Aktualnye problemy matematiki, mekhaniki, informatiki», pp. 23–
28. IM UrO RAN, Izhevsk, Russia (2014) [in Russian].
17. Samarsky, A. A., Vabishchevich, P. N.: Computational Heat Transfer. Vol. 2. The Finite
Difference Methodology.Wiley, Chichester (1995).
18. Samarskii, A. A., Moiseenko B. D.: Ekonomicheskaia skhema skvoznogo scheta dlia
mnogomernoi zadachi Stefana. ZhVMiMF 5(5), 816–827 (1965) [in Russian].
19. Filimonov, M. Yu., Vaganova, N. A.: On Boundary Conditions Setting for Numerical
Simulation of Thermal Fields Propagation in Permafrost Soils. In: CEUR Workshop Pro-
ceedings, vol. 2109, pp. 115–122 (2018).
20. Vaganova, N. A., Filimonov, M. Yu.: Dolgosrochnoe prognozirovanie dinamiki zon ot-
taivaniia mnogoletnemerzlykh porod v uste kusta dobyvaiushchikh skvazhin. In: XXXI
Sibirskii teplofizicheskii seminar, posviashchennyi 100-letiiu akademika S.S. Kutateladze:
sb. tr. Vserossiiskoi konferentsii, pp. 42–48. IT SO RAN, Novosibirsk, Russia (2014) [in
Russian].
21. Bashurov, Vl. V., Vaganova, N. A., Filimonov, M. Yu.: Chislennoe modelirovanie
protsessov teploobmena v grunte s uchetom filtratsii zhidkosti. Vychislitelnye tekhnologii
16(4), pp. 3–18 (2011) [in Russian].