Short and long scale regimes of horizontal flare system
exploitation in permafrost
Mikhail Filimonov1,2 Nataliia Vaganova1,2
fmy@imm.uran.ru vna@imm.uran.ru
1 – Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia)
2 – Ural Federal University (Yekaterinburg, Russia)
Abstract
A mathematical model, numerical algorithm and program code for sim-
ulation thermal changes in permafrost as a result of operation of a hor-
izontal flare system in northern oil and gas field are presented. In the
model the most significant climatic and physical factors are taken into
account such as solar radiation, determined by specific geographical
location, heterogeneous structure of frozen soil, seasonal fluctuations
in air temperature, and freezing and thawing of the upper soil layer.
Results of computing are presented.
1 Processes in Permafrost
Oil and gas production has a significant effect on permafrost due to heat from tubes with oil in the wells leads
to thawing. Permafrost degradation leads to significant difficulties in the construction and operation of various
engineering structures, some of which are already in emergency condition. This problem is compounded by
processes of global warming [1] because of request to ensure reliability of the conservation status of frozen soils
during the construction phase. In this paper, for simulation of heat distribution from wells and other facilities
in permafrost a mathematical model is proposed that takes into account not only climatic(seasonal changes
in temperature and intensity of solar radiation due to geographical location of the area) and physical factors
(different thermal characteristics of inhomogeneous soil changing over time), but also engineering features of
technical systems.
One type of thermal influence is gas flaring in the normal and emergency modes. To reduce the thermal
effect and the permafrost distruction there used different insulation matherials and flare pads. Additionally, it
is necessary take into account and regulate the mode and capacity of the heat flow from the gas torch.
To solve this problem, it is possible to use a mathematical simulation. Choosing an adequate mathematical
model of heat propagation in nonuniform frozen soil, taking into account most significant factors, and carying out
neerical experiments is practicably analize different long and short scale modes of a flare system exploitation and
choose the most safe and reliable. This is an urgent task for oil and gas industry and construction in permafrost
areas.
According to the papers [2, 3] a mathematical model is suggested for long-time forecasting of impacts of
development and exploitation of a flare pad located in areas of permafrost.
In [4] is considered a vertical flare system. To simulate influence of a vertical flare as a rule a radiation heat
transfer model is used on the soil surface. The height of the flame in a horizontal flare system is lower, so that
c by the paper’s authors. Copying permitted for private and academic purposes.
Copyright °
In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the 47th International Youth School-conference “Modern Problems in
Mathematics and its Applications”, Yekaterinburg, Russia, 02-Feb-2016, published at http://ceur-ws.org
253
the surface heating is stronger, and on the upper boundary we may assume a fixed 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 effect 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 flare system and simulated the effect of different thermal cycles on the ground are considered.
2 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 [7] and possible devices used for thermal
stabilization (cooling) of the soil such as seasonal cooling devices, as in [5].
We consider a horizontal flare, which is simulated by a heat source on the surface of the permafrost soil (fig. 1).
The initial temperature distribution in the soil is presented in fig. 2. The basic thermal parameters of soil are
presented in Table 1. In fig. 3 monthly averages air temperature and solar radiation through a year used for
simulation the annual temperature cycle in the soil are shown.
thermal exchange
20С 27х18 m
70С
12х18 m
400
flare 6х9 m
x riprap 2х3m
y
z permafrost
soil
(a) (b)
Figure 1: (a) — considered area, (b) — temperature distribution in flare platform.
Figure 3: Average air temperature and solar
Figure 2: Initial temperature in a frozen soil.
radiation by months.
254
2.1 Basic Equations
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.
Table 1: Basic thermal parameters of riprap layers and soil
Layers of Thermal conduc- Volumetric heat, Heat of phaze Temperature
riprap tivity, W/(m K) kJ/(m3 K) transition, of phaze
and soil frozen melted frozen melted kJ/(m3 K) transition, C
1. concrete slab 1.93 1.57 2150.0 3490.0 141504.0 0.0
0-0.15 m.
2. breakstone 1.75 1.63 2160.0 2630.0 83616.0 -0.8
0.15-0.25 m.
3. sand 1.83 1.39 1536.0 1983.0 69505.8 -0.1
0.25-2.35 m.
4. soil 1.89 1.68 2200.0 2780.0 90983.7 -0.8
2.35-10.35 .
5. soil 2.21 1.75 2350.0 2750.0 64461.7 -0.8
10.35-40.0 .
2.2 Initial and Boundary Conditions
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 (fig 1a). 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.
The flare platform is simulated by a region of fixed temperature in a upper boundary (fig 1b). Balance of
heat fluxes at the surface z = 0 outside of the flare platform is defined by 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
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 flare system to be simulated location. Fig. 2 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. 2 shows temperature distribution in an exploratory well. Applying the developed iterative algorithm [5, 7]
255
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.
At the boundaries of the computational domain the boundary conditions are given
¯ ¯ ¯
∂T ¯¯ ∂T ¯¯ ∂T ¯¯
= = 0, = γ. (4)
∂x ¯x=±Lx ∂y ¯y=±Ly ∂z ¯z=−Lz
In (4) γ 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.
2.3 Methods of Solutions
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)–(4) 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 [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 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 [10] 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. 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 finite difference method is
used with splitting by the spatial variables and taking into account the inner boundaries Ωi . Solvability of the
same difference problems approximating (1)–(4) is proved in [12, 13].
1 1 1
-3 -2 -1 0 1 2 3 10 30 70 -3 -2 -1 0 1 2 3 10 30 70 -3 -2 -1 0 1 2 3 10 30 70
0 3 3 0 -1 -3 -2 0 2 10 103 -2 -1 -3 0
1-2
0
1 2 1 1
0 -3
-1 0 10 2 -1 1 -1 -1
-3-2
3 1
1 0 -1
3
-2 -1 -2 -2
0 2 1 0
-3 -3 -3
-2 -1 -1
-4 -2 -4 -4
z
z
z
-1
-5 -5 -5 -3 -3
-6 -6 -6
-2 -2
-7 -7 -7
-8 -8 -8
-9 -9 -9
-10 -10 -10
30 40 50 60 30 40 50 60 30 40 50 60
x x x
june october march
Figure 4: I operating mode. Temperature field under flare platform on june, october, and march.
256
1 1 1
-3 -2 -1 0 1 2 3 10 30 70 -3 -2 -1 0 1 2 3 10 30 70 -3 -2 -1 0 1 2 3 10 30 70
0 3 0 -3 -103 10 -2 -3 0
70
10
2 -1 10 10
2
30
30 1
70
10 1
0 1 1 70 -2
-2
-1 -1 -1
-1
-3
2
30
3 1
-3
2
1
10
0 30
-2 -2 0 -2
2
-1 0
13
-1
3
2
1
-3 -3 -3
3 2
-2 0 -1 10
-2 10 -1
0
23
-4 -4 -4
3
1
1
-1
0
z
z
z
-5 -5 2 -5
-2
-3
-2
0
-2 -3
-1
1
-6 -6 -6
0
-2 -1 -2
-7 -7 -7
-1
-8 -8 -8
-9 -9 -9
-2
-2
-2
-10 -10 -10
30 40 50 60 30 40 50 60 30 40 50 60
x x x
june october march
Figure 5: II operating mode. Temperature field under flare platform on june, october, and march.
3 Numerical Results
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 of flare 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 [15]. Also the used
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 [16]. Implicit method of
solution allows to use different scales of time steps in numerical simulation.
Let consider a 3D computational domain with Lx = 87m, Ly = 87m, Lz = 40m. The flare system in in center
of the soil surface (fig 1). The flare 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 different frozen soils. Parameters of soils and riprap layers
are presented in table 1. We will compare 4 types of flare system operation mode (table 2). Let assume that the
month of starting operation is may.
In figures 4–7 temperature field are shown in section by y axis near the flare platform for different operating
modes. In figures 8 and 9 annual temperature changing are presented at the depth of 1m and 3m. In figures 10
and 11 the detalization of temperatures under flare platform are shown for summer and winter. The “ground”
line in the figures denotes the natural temerature disribution without influence of flare.
Effect of heating for I operating mode is not too intense, but is permanent, and its influence is presented
even at the depth of 3m (fig 8). Operating mode II is obviously has a significant effect on the temperature
distribution under the flare platform, the deepest soil thawing is observed in this mode (fig 5). When the flare
heated for a month (operating mode III), despite the fact that high peak temperature is observed (fig 9), but
significant period of rest allows ground to be cooled up to near background values of temperature (fig 10 and
1 1 1
-3 -2 -1 0 1 2 3 10 30 70 -3 -2 -1 0 1 2 3 10 30 70 -3 -2 -1 0 1 2 3 10 30 70
0 0 -2 -3 -1 -2 -3 0
1 1 0 10 -1 0
1 2 1
23
1
2
-1 -1 -1
10
1
3 -1
2
0
-2 70 -2 0 0 -2
0
-3 30 -3 1 -3
10 1
-1 -1
1
-4 32 -4 -1 -4
z
z
z
-1
0
-5 0 -5 -5 -3
-2 -2
-3 -3
-6 -6 -1 -6
-2
-7 -7 -2 -7
-8 -8 -8
-9 -9 -9
-2
-10 -10 -10
30 40 50 60 30 40 50 60 30 40 50 60
x x x
june october march
Figure 6: III operating mode. Temperature field under flare platform on june, october, and march.
257
Table 2: Types of flare system operatng mode
time in action periodicity
I operating mode 1 hour 6 days
II operating mode 24 hours 6 days
III operating mode 1 month 12 years
IV operating mode 1 month 6 months
fig 11). Moreover, heating in the same type cyclically (every 6 months, operating mode IV) has no a significant
thermal trace in the soil temperature (fig 7).
4 Conclusion
On the base of a mathematical model of thermal field in a nonuniform frozen soil from a heat source located on
the surface of the ground, a horizontal flare system used to gas flaring in northern oil and gas fields is simulated.
Forecasts of permafrost degradation around this system is obtained for different 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 effects 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 70o C). Also the features of thermal field propagation in the soil are
considered in detail, taking into account different physical and climatic factors. Operating modes, when the flare
platform time could be possible to cool down without heat accumulation, are preferable.
Acknowledgements
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).
References
[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 effects and climate.
Academic Journal of Science, 3(1): 121–128, 2014.
[3] M. Y. Filimonov, N. A. Vaganova. Simulation of technogenic and climatic influences in permafrost for
northern oil fields exploitation. Lecture Notes in Computer Science, 9045:185–192, 2015.
1 1 1
-3 -2 -1 0 1 2 3 10 30 70 -3 -2 -1 0 1 2 3 10 30 70 -3 -2 -1 0 1 2 3 10 30 70
0 70 0 -3
-2 -3
-2 0
10
1 1 0 -1 10 2 -1 0
2
1 1
0
10
-1 -1 -1
70
3
3
2
-2 -2 1 -2
1
2
-3
-1
-1
30
-3 -3 -3
-2 10 0
10 -1
1 0
0
-2
-4 2 -4 -1 -4
3 -3
z
z
z
-2
-5 -5 -5
-1 -3
-2
-6 -2 -6 -6
-1
-7 -7 -7
-8 -8 -2 -8
-2
-2
-9 -9 -9
-2
-10 -10 -10
30 40 50 60 30 40 50 60 30 40 50 60
x x x
june october march
Figure 7: IV operating mode. Temperature field under flare platform on june, october, and march.
258
% #
$ "
#! !
"
!
(depth is 1 m) (depth is 3 m)
Figure 8: I and II operating mode. Temperature changes under flare platform.
% %
$ $
#! #!
" "
! !
(depth is 1 m) (depth is 3 m)
Figure 9: III and IV operating mode. Temperature changes under flare platform.
[4] M. Yu. Filimonov, N. A. Vaganova. Simulation of thermal effect of a vertical flare 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 fields in design and
operation of northern oil and gas fields. 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 Difference
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.
259
' '
& &
%#$ %#$
# #
!" !"
(depth is 1 m) (depth is 3 m)
Figure 10: III and IV operating mode. Temperature changes under flare platform on summer.
% &
$ %
#! $"
" #
! "
!
(depth is 1 m) (depth is 3 m)
Figure 11: III and IV operating mode. Temperature changes under flare platform on winter.
[12] N. A. Vaganova. Existence of a solution of an initial-boundary value difference problem for a linear heat
equation with a nonlinear boundary condition. Proceedings of the Steklov Institute of Mathematics, S1:
S260–S272, 2008.
[13] V. V. Bashurov, N. A. Vaganova, M. Y. Filimonov. Numerical simulation of thermal conductivity processes
with fluid filtration in soil. Journal of Computational technologies, 16(4):3–18, 2011.
[14] M. Yu. Filimonov, N. A. Vaganova, V. V. Gusarova. Elaboration of mathematical models and software
to describing natural and technogenic processes in soils in northern fields well pads. Proceedings of XI
All-Russian Congress on Fundamental Problems of Theoretical and Applied Mechanics. Izhevsk: 3917–3919,
2015.
[15] N. A. Vaganova, P. A. Vasev, V. V. Gusarova, A. S. Igumnov, M. Yu. Filimonov. Application of cloud
technologies for simulations of northern oil and gas fields exploitation. Proc. of IMech UrO RAN. Conf.
“Actual problems of mathematics, mechanics, computer science”: 23–26, 2014.
[16] N. A. Vaganova, M. Yu. Filimonov. Parallel Splitting and Decomposition Method for Computations of Heat
Distribution in Permafrost. Proceedings of the Ural-PDC 2015, Vol-1513: 42–49, 2015.
260