УДК 517.9
A Numerical Model of the Seasonal Thawing of Permafrost in the Swamp-lake Landscapes
Victor M. Belolipetskii*
Institute of Computational Modelling SB RAS Akademgorodok, 50/44, Krasnoyarsk, 660036 Institute of Mathematics and Computer Science Siberian Federal University Svobodny, 79, Krasnoyarsk, 660041
Russia
Svetlana N. Genova^
Institute of Computational Modelling SB RAS Akademgorodok, 50/44, Krasnoyarsk, 660036
Russia
Received 25.11.2015, received in revised form 30.01.2016, accepted 20.02.2016 The theoretical description of the temperature field in the soils during freezing or thawing is carried out using solutions of Stefan's problem. A mathematical model based on the equations of thermal conductivity for frozen and thawed layers. We consider the areas in which there are lakes or bogs. We distinguished the following layers in the vertical structure of the zone of permafrost: thawed soil, frozen soil, water, ice, snow. We offer a simplified numerical algorithm for solving of one-dimensional (in the vertical direction) heat conduction problems with moving boundaries of phase transition with the formation of new and cancellation of existing layers. A simplified numerical algorithm for solving one-dimensional (in the vertical direction) heat conduction problems with moving boundaries of phase transition with the formation of new and cancellation of existing layers is offering.
Keywords: permafrost, Stefan's problem, thawed and frozen soil, small dimensional numerical model. DOI: 10.17516/1997-1397-2016-9-2-158-165.
Introduction
In connection with the change in global temperature is of interest assess the response of permafrost to climate change. We consider the areas in which there is a lake or swamp. Since the vertical temperature gradients are usually larger than the horizontal one, so all physical process assumes one-dimensional in the vertical direction in the description of the heat transfer. We distinguish the following layers in the vertical direction: the snow, ice, water, thawed soil , frozen soil. The theoretical description of the temperature field in the water and soils during their freezing and thawing is carried out using solutions of Stefan problem [1]. A mathematical model based on the equations of thermal conductivity for frozen and thawed areas. At the borders of phase transition (freezing-thawing) the conditions of equality of temperatures to the phase
* [email protected] [email protected] © Siberian Federal University. All rights reserved
transition temperature and the Stefan condition are posed. There is an extensive literature on the mathematical modeling of permafrost (see, Eg, [2-8]).
In this paper the numerical model of small dimension is proposed for discription vertical temperature distribution in thawed and frozen layers taking into account the formation of new and revocation of the existing layers ( [9]). One can allocate two types of water bodies: a) reservoirs, which freeze to the bottom in the winter, and the ice melt throughout the depth and the top layer of the bottom defrost in the summer ("shallow" water bodies); b) reservoirs that do not freeze to the bottom in the winter ("deep" water). There are various options for the location of frozen and thawed layers (Tab. 1). For "deep" waters there are three variants. Variant 1: water layer, layer of thawed soil, layer of frozen soil, layer of thawed soil. Variant 2: snow-ice layer, water layer, layer of thawed soil, layer of frozen soil, layer of thawed soil.
Variant 3: water layer, ice layer, water layer, layer of thawed soil, layer of frozen soil, layer of thawed soil.
For "shallow" water bodies seven variants are considered. When switching from one variant to another the layers are added or deleted.
Table 1. Variants of the location for frozen and thawed layers in the swamp-lake landscapes
variant water ice water thawed frozen thawed frozen thawed
N (snow) soil soil soil soil soil
1 + + + +
2 + + + + +
3 + + + + +
4 + + +
5 + + + +
6 + + + + + +
7 + + + + + +
8 + + + + + +
Fig. 1 illustrates a vertical structure for variants 6 and 7, Fig. 2 shows a scenario of switching from one version to another for the "shallow" reservoirs.
1. Mathematical model of dynamics of freezing and thawing of permafrost
The vertical temperature distribution in every layer are determined by solving the heat equation, satisfying the appropriate boundary conditions:
dT d2T■ f=* -è- «
Here T(t, z) is a temperature of i-th layer (h-1 < z < h), * is a thermal diffusivity, t is time, z is a vertical coordinate (downward).
Fig. 1. Variants 6 (a) and 7 (b)
Fig. 2. The scheme of switching from one variant to another
Boundary conditions. The condition on the boundary of the atmosphere - water (z = 0)
K^ 1 =- ,
m
cw pw
the condition on the boundary of the atmosphere - snow+ice (z = 0)
T = Tice -
(2)
(3)
The condition at the bottom of the pond (boundary of i-th hi-1 < z < hi and (i + 1)-th
hi < z < hi+1 layers, z = hi = hbt )
t=T+1 m,=(€)
(4)
i+l
On moving boundaries of phase transition (z = hj )
j ^ = ). - (^ )
Tj = T.
j+i
Tp
ph 7
p.
dt
dz J
(5)
j+i
Here pw is water density, cw is the specific heat capacity of water, Fn is the total heat flow on the boundary of the atmosphere-water, Tice is temperature on ice surface, z = hi is coordinate of the boundary between i-th and (i + 1)-th layers, A is coefficient of heat conductivity, Tph is phase transition temperature, p. is density of j-th layer, L.
Q.
Lw ■ Wj is volumetric latent heat
of melting environment in j-th layer, Wj = -Q— is soil humidity in j-th layer, Qjw is the volume
of water in the soil, Qj is the volume of the soil. The total heat flow (Fn) and temperature on the boundary snow-ice (Tice) is defining according to the formulas, described in [10].
On the lower boundary of permafrost layer the conditions (5) are supplemented with setting
dT
of a geothermal gradient in the layer on thawed soil G = -7—, (G = 0.02 — 0.05° C/m [11]).
dz
The initial conditions:
Ti(0,z) = TO, Si = S°. (6)
Consider an arbitrary i-th layer: hi-i ^ zi ^ hi, Si = hi — hi-i, where Si is the thickness of i-th layer. We introduce new independent variables t, Z:
zi — hi-1
Since
dZi
wi = ~m
written as
1
Si
d =
dzi
(Zi — 1)
t = 17
1 d
Si dZi' dhi-i dt
Zi =
dz2
51
1 dL
52 dz27
0 < Zi < 1.
d
d = d
dt = dt + Wta&
dhi — Zi d
dTi
-dt + Wi dz
7 then equation (1) and boundary conditions (2)-(5) dTi Ki d2 Ti
K m
Si dZi
Fn
ti=°
cw pw
S2 dz2
Ti\ii=° = Tice7
Ti\ii = 1 = Ti+iki+l =
(Ai OTA f Ai+i dTi+A T| T . T
U -dTi)i=i " IS+1 az+i)5i+i=°7 Tj=1 = Tj+ll'j+1=° = Tph7
dhj _ f A^-LJ
Pj Ljd = I
(Ai j —( Aj+i dTj+A V s. dzj J j=i V Sj+i dzj+i) j+i=°.
(7)
(8) (9)
(10)
2. Numerical algorithm
Let us consider the solution of the formulated problems on coarse (three-point) grid for each of layers: Zii = 0, Zi2 = 0.5, Zi3 = 1.
We approximate the equation (7) using the implicit difference scheme and directional differences for the convective terms [12]. Grid equations corresponding to the differential equation (7) for i-th layer are of the form:
Ti,2 — Ti At
2 + w Ti,3 — 2w0Ti,2 + w+Tii
Ki(Tii — 2Ti2 + Ti3)
(0.5Sn)2
(11)
= -
(dhi + dhi-i \ \dt dt J
/W),
\wi^\, w = Wi,2 -\wifl\, w+ = wi,2 + here At is
time step, t„+i = tn + At, T^ = Ti(tn,&,k), Ti,k = Ti(tn+i,&,k).
Boundary conditions (8)-(10) for the difference grid are represented as: on the border of the atmosphere-water
1 +
8KiAt
TW.
8KiAt ±i'1 (gn )2 '2
nn . 4AtFn
Tn +
cw pU
(12)
on the border of the atmosphere- snow+ice
Ti,1 Tic
(13)
on the water bottom z = hi
Ti,3 = Ti
i+i,i,
Ai /rp rp \ _ Ai+1
Oi Oi+1
(Ti
i+1,2
T
Ti
i+1,1)
At the borders of phase transition z = hj
(14)
(15)
T,
j,3
Tj
j+1,1
T,
Pj Li
hn+1 - hn
At
Xj
Tj,3 - t,
j,2
0.50n
ph,
- Xj+1
Tj + 1,2 - Tj+1,1 0-50n+1
(16)
(17)
If Sj becomes less than a predetermined small value e (Oj < e), the corresponding layer will be cancelled. If as a result of melting or freezing a new k-th layer is formed, the initial thickness Sk = e and temperatures is put for this layer.
The following algorithm is used for the obtained tasks. Let is known all parameters (temperature distributions in layers under study and the positions of phase transition) on n-th time step. Then a finding the unknown parameters at time tn+1 is performed in two stages. The first step is to define the temperature distributions in the selected layers (taking into account the relations (11)-(16)) by solution of systems of linear algebraic equations of small dimension. In the second phase it is clarified the position of the phase boundary by the numerical solution of equation (17).
0
w
3. Calculations samples
Model calculations were performed for a swamp depth of 25 cm, 50 cm, 75 cm and 100 cm by weather data 2010-2011 years for the weather station Dudinka. Year 2010 can be considered as "cold" (about 240 days with mean daily negative temperature, the average temperature for the period from 1 October to 30 April was —20.72oC); year 2011 is "warm" (190 days with mean daily negative temperature and the average temperature for October to April was —15.22oC). There were determined the depth of seasonal thawing of permafrost soils, the temperature of the melt layer, the time intervals of the existence for layer of thawed soil. In Tabs. 2, 3 are given the results of calculations shown the impact of the swamp depths and weather data on the main characteristics of the permafrost.
From the freezing-thawing of soil calculation results it follows that with increasing depth of "shallow" water it is increased the period of existence of the thawed soil layer. For sufficiently
Table 2. The thickness of the layer of thawing frozen soil (cm)
swamp depth (cm) 25 50 75 100
2010 122.9 118.5 114.2 108.5
2011 143.2 139.38 135.04 144.6
Table 3. The count of days existence of the thawed soil layer
swamp depth (cm) 25 50 75 100
2010 212 228 243 264
2011 242 292 311 365
great depths of the lake a layer of thawed soil does not freeze completely during the year. With further increase of the depth of the reservoir it does not freeze to the bottom ("deep" water) and thawed layer exists throughout the year.
In Fig. 3 the calculated vertical temperature distributions are given. In Fig. 3a is given temperature distribution for variant 7, in Fig. 3b is given for variant 2.
10 0 10 20 -30 -20 -10 0 10
Fig. 3. The calculated temperature in layers for different weather data: a) for summer (variant 7), b) for winter (variant 2). Water depth is 1 m
The proposed numerical algorithm allows to describe the dynamics of annual freezing-thawing
permafrost in the swamp-lake district landscapes and to assess the impact of climate change.
The work is executed at partial financial support of Basic Research Program of Presidium of
RAS "Searching basic research for the development of Russian Arctic" (project 12).
References
[1] A.M.Meirmanov, The Stefan Problem, Novosibirsk, Nauka, 1986.
[2] O.A.Anisimov, F.E.Nelson, Forecast of changes in permafrost conditions in the Northern hemisphere: application of results of balance and transitive calculations on models of the General circulation of the atmosphere, Kriosfera Zemli, (1998), no. 2, 53-57 (in Russian).
[3] S.P.Malevsky-Malevich, E.K.Molkentin, E.D.Nadezhina, Model estimates of the evolution of permafrost and the distribution of the layer of seasonal thaw, depending on the climatic conditions in the Northern regions of West Siberia, Kriosfera Zemli, 4(2000), no. 4, 49-57 (in Russian).
[4] M.M.Arzhanov, A.V.Eliseev, P.F.Demchenko, I.I.Mokhov, Modeling of changes in temperature and hydrological regimes of subsurface permafrost, using the climate data (reanalysis), Kriosfera Zemli, (2007), no.4, 65-69 (in Russian).
[5] T.S.Sazonova, V.E.Romanovsky, A model for regionalscale estimation of temporal and spatial variability of activelayer thickness and mean annual ground temperatures, Permafrost and Periglacial Processes, (2003), no. 2, 125-140.
[6] L.E.Goodrich, The influence of snow cover on the ground thermal regime, Canadian Geotech-nical Journal, 19(1982), no. 4, 421-432.
[7] D.J.Nicolsky, V.E.Romanovsky, V.A.Alekseev, D.M.Lawrence, Improved modeling of permafrost dynamics in a GCM land-surface scheme, Geophysical Research Letters, 34(2007), no. 8, L08501.
[8] T.J.Zhang, O.W.Frauenfeld, M.C.Serreze, A.Etringer, C.Oelke, J.McCreight, R.G.Barry, D.Gilichinsky, D.Q.Yang, H.C.Ye, F.Ling, S.Chudinova, Spatial and temporal variability in active layer thickness over the Russian Arctic drainage basin, Journal of Geophysical Research-Atmospheres, 110(2005), Art. no. D16101.
[9] V.M.Belolipetskii, S.N.Genova, Numerical model of dynamics Permafrost in the bog-lake landscapes, Fundamental Problems of Water Resources: Proceedings of IV Russian Scientific Conference, Moscow, Water Problems Institute RAS, 2015, 95-97.
[10] V.M.Belolipetskii, S.N.Genova, V.B.Tugovikov, Y.I.Shokin, Numerical modelling of problems of hydro-icethermics of currents, SB RAS, Novosibirsk, 1994 (in Russian).
[11] M.S.Krass, V.G.Merzlikin, Radiation physics of snow and ice, Leningrad, Gidrometeoizdat, 1990 (in Russian).
[12] V.I.Polezhaev, A.V.Bune, N.A.Verezub, etc., Mathematical modeling of convective heat and mass transfer based on the Navier-Stokes equations, Moscow, Nauka, 1987 (in Russian).
Численная модель сезонного оттаивания вечной мерзлоты в болотно-озерных ландшафтах
Виктор М. Белолипецкий
Институт вычислительного моделирования СО РАН Академгородок, 50/44, Красноярск, 660036 Институт математики и фундаментальной информатики Сибирский федеральный университет Свободный, 79, Красноярск, 660041
Россия
Светлана Н. Генова
Институт вычислительного моделирования СО РАН Академгородок, 50/44, Красноярск, 660036
Россия
Теоретическое описание температурного поля в почвах при их промерзании или оттаивании осуществляется с помощью решений задач Стефана. Математическая модель основывается на уравнениях теплопроводности для мерзлых и талых слоев. Рассматриваются территории, на которых имеются озера или болота. Выделяются следующие слои в вертикальной структуре зоны вечной мерзлоты: талый грунт, мерзлый грунт, вода, лед, снег. Предлагается упрощенный численный алгоритм решения одномерных (в вертикальном направлении) задач теплопроводности с подвижными границами фазового перехода с образованием новых и аннулированием существующих слоев.
Ключевые слова: вечная мерзлота, задачи Стефана, мерзлые и талые слои, малоразмерная численная модель.