At present, much attention is paid to the development of the Arctic regions. Prediction of stability of various engineering structures in permafrost propagation zones under the impact of various heat sources on environment, which can be located both inside the ground and on its surface, are in the top. Such heat sources are generated, for example, in operation of northern oil and gas fields. In this case, the correct setting of the boundary conditions on the surface of the soil during computer simulation of thermal field propagation in the ground plays an important role. The main attention is paid to the nonlinear boundary condition on the ground surface and its role in thermal field propagation in the near-surface layer of the soil. The results of numerical calculations showing the possibility of numerical algorithms adaptation to specific geographic coordinates with taking into account soil properties at the location are presented.
The main characteristics of permafrost, which are usually taken into account in the construction of production
wells and various engineering structures under these conditions, are the type of distribution (continuous,
intermittent, island), the type of cryogenic structure (massive, layered, mesh), and the degree of iciness. The frozen
soils have different physicochemical properties, which can vary in all directions. In summer, due to positive
air temperatures and solar radiation, seasonal thawing of the upper layer of the soil takes place, in winter the
reverse freezing process is observed. It was noted that in the northern high latitudes the average air temperature
increased faster than the average global temperature, that affects on the state of permafrost. Permafrosts that
have a negative temperature below the seasonal thawing zone occupy about 25% of the globe and the study of
the dynamics of its boundaries is important for various structures constructing in these territories and is also
associated with the climate changes [
Note that more than 75% of all Russian buildings and structures in the permafrost zone are constructed
and operated on the base of principle of conservation of frozen soil foundation [
To solve the problem, a number of mathematical models are developed [
Let T = T (t; x; y; z) be soil temperature at the point (x; y; z) at the time moment t. 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 technological (figure 1a) and climatic (figure 1b) factors. (1) (2) THERMAL EXCHANGE
SOLAR RADIATION
EMISSIVITY THERMAL CONDUCTIVITY (b) 6,00 4,00 2,00 0,00 -2,00 -4,00 (c) -6,00 m 1 m 2 m 3 m 4 m 5 m 6 m 7 sep okt nov Aver(sep-nov) MesData (a) with initial condition is density [kg=m3], T is temperature of phase transition [K], c (T ) =
{ cc21((xx;; yy;; zz));; TT >< TT ;; is specific heat [J/kg K]; (T ) = {
12((xx;; yy;; zz));; TT <> TT ;; is thermal conductivity coefficient [W/m K ], k = k(x; y; z) is specific heat of phase transition, is Dirac delta function. 2.1
Boundary Conditions The main heat flow associated with climatic factors on the surface z = 0 is shown in figure 1a.
Let consider in detail the boundary condition obtaining on the surface of the ground, since a condition of this type is rarely used in solving the problems under consideration. As the boundary condition on the surface of the ground (the main zone of natural thermal fields generating), the equation of the balance of flows that bring and take away energy is used, which takes into account the main climatic factors: the average monthly air temperature and solar radiation power (mainly on the spring and summer months). Let Tair = 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. Let compose the heat balance on the surface of the soil z = 0. We’ve got: q1 = q is the part solar radiation to heat the ground; q2 = b(Tair − T (t; x; y; 0)) is the heat exchange with air on the surface of the ground z = 0; q3 = " (T 4(x; y; 0; t) − Ta4ir) is emissivity of the heated soil surface; q4 = − @T (x@; yz; 0; t) is the heat transfer due to possible internal heat flow on the surface of the soil.
Here = 5; 67 · 10 8Wt=(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 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 atmospheric conditions, angle of incidence of solar radiation, i.e. latitude and time. The balance of flows on surface z = 0 has the following form: q1 + q2 = q3 + q4 and determines the corresponding nonlinear boundary condition @T (x; y; 0; t) q + b(Tair − T (x; y; 0; t)) = " (T 4(x; y; 0; t) − Ta4ir) + : (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 [
Taking into account the absorbed and reflected solar energy in condition (3) is a complex problem, since there are many undetermination in the parameters and ". Ideally, for a complete quantitative description of solar radiation absorbed and reflected by soil surface in a specific geographical area, many factors have be taken into account, which, as a rule, are not known. In particular, the specific parameter p is a part of long-wave radiation reflected by atmosphere toward the Earth’s surface. If, generally, p = 0:84 for the Earth, then for a specific geographic point on the earth’s surface this value will be very different from this value. In order to obtain more accurate estimates of these parameters, it is necessary to take into account, with other things, the number of sunny and rainy days during a year, as well as the other climatic parameters that are often probabilistic, and their monthly averages are required.
Moreover, the value of " parameter does not always accurately may be determined for a specific geographic
point. The value of this parameter can significantly affect on the boundary condition (3). The same arguments
are valid and for other parameters contained in boundary condition (3). Also, all these parameters change during
the time in accordance with the climate cycle (in our model such a change is considered). For example, it is
necessary to consider the layers of snow cover, the structure of the snow, the presence on the surface of a layer of
humus, vegetation, etc. The study of the influence of snow cover is available in many works (for example, [
Thus, an important task is to determine the parameters in the nonlinear boundary condition (3). We will assume that for a given geographic point we know the values represented in figure 1c. The initial temperature distribution is usually determined by the data from an exploration well at the time moment t0 at the point (x0; y0; 0) on the surface of the computational domain Ω. Let denote this initial temperature distribution in the exploratory well by T (t0; x0; y0; z). For definiteness, we assume that t0 is the number of days since the beginning of the year in which measurements are taken. Our next task is to refine (select) the parameters in the boundary T (t0; x0; y0; z) ≈ T (t0 + 365j; x0; y0; z); j = 1; 2; : : : ; J: That is the solution of problem (1)–(3) is periodically repeated (in accordance with the annual temperature cycle) with a certain approximation of a given initial distribution T (t0; x0; y0; z). And the more precisely this relation is satisfied during a longer time interval (for a larger value of J ), the more accurately it is possible to select the parameters in (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 figure 1c the calculated temperature distribution in the exploratory well is compared with the measured temperature distribution at a given time point t0.
At the boundaries of the computational domain the boundary conditions are given =
In (5) is a positive number, corresponding to a geothermal flux value. As a rule is possible to be set zero in calculations. is a small number and it 3
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. (4) (5)
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 [
With using these ideas [
After the algorithm has adapted to a specific geographic location in one of the northern oil and gas fields, the detailed numerical calculations were carried out, related not only to predict the permafrost thawing from the producing wells, but also for the optimum choice of pad riprap.
In the figures 2–5, the numbers indicate the month’s number, the curves show the temperature distribution for the month with the z depth.
The temperature distribution in spring (curves 3, 4, 5) and summer (curves 6, 7, 8) in the soil without and with a riprap are shown in figure 2.
The temperature distribution in autumn (curves 9, 10, 11) and winter (curves 12, 1, 2) in the soil without and with a riprap are shown in figure 3.
Figures 4–5 show the temperature changes depending on the used riprap. 4
Thus, the developed mathematical model and software product allow to carry out detailed numerical calculations on long-term forecasting of temperature field changes in the near-surface layer of soil. The simulations take into account the most significant climatic and physical factors, which in general are difficult to be described in detail and contain many undetermined parameters. The computations allow to choos the optimal version of ripraps for the development of an oil and gas field with using available materials. On the other hand, the developed approach adequately describes long-term dynamics of changes in the active layer (ALT), taking into account various scenarios of climate change, which allows a long-term forecast related to changes in the permafrost boundaries.
This work was supported by Russian Foundation for Basic Research projects (16–01–00401) and program of scientific research UrB RAS 18–1–1–8.