CC BY
6
DOI: 10.17516/1997-1397-2022-15-6-753-762 УДК 517.96
Modelling of the Ice Cover Dynamics of a Freshwater Reservoir
Evgeniy N. Vasil'ev*
Institute of Computational Modelling of SB RAS Krasnoyarsk, Russian Federation
Received 22.04.2022, received in revised form 28.08.2022, accepted 20.10.2022
Abstract. A non-stationary computational model for determining the temperature regime of a reservoir and the dynamics of the ice cover in the cold season is considered in the paper. The model is based on the numerical solution of a generalized Stefan problem. Average daily air temperature, wind speed and solar radiation absorbed by the surface of the reservoir were taken into account. Computational modelling of the temperature dynamics of the reservoir was carried out for climatic conditions in the region of Krasnoyarsk. The influence of heat fluxes and snow cover thickness on the time dependence of the temperature in the reservoir and dynamics of the ice growth was analysed. The results of the computational experiment are compared with the measurement data of the ice cover thickness.
Keywords: Stefan problem, numerical simulation, thermal conductivity, ice cover dynamics, temperature regime of the reservoir.
Citation: E.N. Vasil'ev, Modelling of the Ice Cover Dynamics of a Freshwater Reservoir, J. Sib. Fed. Univ. Math. Phys., 2022, 15(6), 753-762. DOI: 10.17516/1997-1397-2022-15-6-753-762.
Introduction
Modelling of temperature regimes and ice cover dynamics on lakes and reservoirs is of great practical importance. The simulation results can be used to predict and evaluate various ice phenomena, such as the onset of freeze-up, an increase in ice thickness, a decrease in the thickness and strength of the ice cover with the onset of warming. Such data are important for determining the timing of navigation on water bodies to ensure the safety of the operation of ice crossings and winter roads.
To predict the state of water bodies in the cold period various methods are used. One of the approaches for assessing ice growth is based on the use of the sum of negative air temperatures for the cold period of the year [1,2]. Another approach to determine the thermal evolution of the ice cover is based on the heat balance equation for a two-layer "ice-snow" medium [3,4]. Non-stationary mathematical models based on the solution of the Stefan problem provide more accurate results. They take into account detailed climatic data which allow for the main factors affecting the dynamics of ice covers [5]. This paper presents a mathematical model and the results of computational modelling of the temperature regime and the ice cover dynamics of a reservoir for the climatic conditions of the city of Krasnoyarsk.
* ven@icm.krasn.ru https://orcid.org/0000-0003-0689-2962 © Siberian Federal University. All rights reserved
1. Computational model of the ice cover dynamics of the reservoir
With the onset of the cold period a reverse temperature stratification is established in lakes and reservoirs. It means that there is an inhomogeneity of temperature distribution only in the vertical direction [6]. At the same time the water temperature is close to 0 °C in the surface layer while in the deep bottom layer it can rise to 3-4 C. In this mode, the water density at the depth of the reservoir has the maximum value. When approaching the surface the water density decreases. It is assumed that there is no convection and the only mechanism of heat transfer between water layers is thermal conductivity.
Modelling of the ice cover dynamics of the reservoir is carried out on the basis of the numerical solution of the Stefan problem. The problem is presented in a generalized formulation when the heat of the "ice-water" phase transition is taken into account in the value of the effective heat capacity of the medium [7]
t s t sdT d / VT) ce(z)p(z) — = — I \{z) — ) + qv (z), (1)
here ce, p are specific effective heat capacity and density of the medium; T — temperature; A — thermal conductivity coefficient; t — time; z is the vertical coordinate directed downward from the surface of the reservoir, qv is the power of the volumetric heat source associated with the absorption of the solar radiation flux in the depth of the reservoir. The inhomogeneity of ce, p and A is due to the presence of a phase transition, the difference in the thermophysical properties of water and ice and also the relationship between these properties and temperature.
The heat transfer equation (1) is supplemented by the boundary condition of the third kind at the upper boundary of the reservoir
dT
dz
= q, (2)
=0
here q is the total value of the heat flux on the surface of reservoir. The boundary condition at the lower boundary is dictated by the depth of the reservoir. For deep reservoirs it is assumed that temperature is constant at the bottom of reservoir (it is 4 C or another value from the measurement data). For shallow reservoirs the condition of heat exchange with the bottom of reservoir is specified [6].
The heat capacity of a medium located at the phase boundary at temperature T* has a jump associated with the presence of the latent heat of the phase transition q*. To take into account q* in the effective heat capacity the ¿-function is used and
Ce(T) = c(T)+ q*5(T - T*). (3)
In the numerical solution of the problem the smoothing method [8] was applied, and discontinuous ¿-function was replaced by a smooth function f (T). The function is represented as one half-cycle of the sine function in the temperature interval (T*-AT, T*+AT). It also satisfies the normalization condition [9]
T *+AT
/ f (T) dT =1. (4)
JT *-AT
This method assumes that liquid and solid phases are simultaneously present in the medium in a given temperature range. This intermediate state of the medium has certain physical interpretation and corresponds to wet ice. Such dynamic state is especially characteristic during the
z
melting of ice when at near zero temperature the medium is no longer a solid body, but not yet liquid water. Polycrystalline ice consists of individual single crystals. Ice melting begins along the boundaries of ice crystals. As a result, a volumetric network of water inter-layers is formed in the medium between the crystals [10]. With the growth of ice on the lower boundary of the cover, under certain conditions a layer with an intermediate state can also form that includes lamellar ice crystals with water enclosed between them [11]. The values of the thermal conductivity coefficient A and the density p in this temperature range is set according to a linear law from the values corresponding to water at temperature T*+AT to the values for ice at temperature T* — AT. The numerical solution of equations (1)-(4) is carried out using the implicit finite difference algorithm [8,9].
The temperature regime of the reservoir and the change in the ice cover thickness are primarily determined by the value of the total heat flux q at the upper boundary of the reservoir. This is represented by boundary condition (2). The terms in (5) describe the main mechanisms of heat transfer on the surface of the reservoir. Solar radiation incident on the surface is partially reflected into the atmosphere, and the rest of the radiation power qr is absorbed by the reservoir. The thermal radiation of the atmosphere qa is absorbed by the surface of the reservoir. The radiation of the reservoir qs is emitted from this surface. The intensity of convective heat exchange between the reservoir surface and the surrounding air depends on the heat transfer coefficient a and the temperature difference between the atmospheric air Ta and the reservoir surface Ts. Then
Here A is the reflection coefficient of solar radiation (albedo), e is the emissivity of the reservoir surface. The value of the heat transfer coefficient a depends on the state of the reservoir surface and the wind velocity u. For the water surface of an unfrozen reservoir the relation a=5.8(u + 0.3)0'5 is used [6]. For a frozen reservoir the heat transfer coefficient is calculated using the expression for the snow cover a=3.4 + 2.2u [13].
The power of the radiation of the reservoir is determined as follows qs=ea(Ts + 273)4, where a is the Stefan-Boltzmann constant. Thermal radiation of the atmosphere qa is [14]
To determine qr data from long-term measurements of the total energy of solar radiation incident on a horizontal surface under actual cloudiness conditions in Krasnoyarsk was used. The measurement results are available as monthly averages of total (direct and diffuse) solar radiation. However, values of daily average solar radiation flux density qr are required for computational modelling. They were determined from the monthly average values of total solar radiation using spline interpolation. Approximation accuracy was checked by numerical integration of the obtained qr. As a result, monthly average values of the total solar energy were obtained, and they differ from the initial values by less than 1%. The solid curve in Fig. 1 shows the calculated qr (t). The horizontal segments for the corresponding months show the monthly average values of the solar flux density.
2. The effect of snow cover on the temperature regime of a reservoir
The presence of a layer of snow on the ice cover significantly affects the thermal balance of the reservoir. Firstly, in the presence of snow the albedo and emissivity values of the reservoir surface
q = (1 - A)qr + eqa - qs + a(Ta - Ts).
(5)
qa = 0.925a(Ta + 273)4 - 0.03.
(6)
Fig. 1. Time dependence of the solar radiation flux density in the cold season for the city of Krasnoyarsk
are radically changed. Secondly, snow has low thermal conductivity, and this leads to additional thermal insulation of ice. A temperature difference is established on the snow layer while the temperature of the lower layer of snow cover is equal to the temperature of the upper surface of ice. The temperature of the upper snow layer is determined at each time step of problem (1)-(4) from the numerical solution of the non-linear heat balance equation (7). This equation includes incoming qa and outgoing qsn=ea(Tsn + 273)4 radiation fluxes on the snow surface and vertical heat fluxes towards the atmosphere and ice:
(1 - A)qr + eqa - qsn + a(Ta - Tsn ) + (T - Tsn)/Rsn = 0,
(7)
where Tsn is the temperature of the upper surface of the snow cover, Rsn is the thermal resistance of the snow layer. The value of Rsn is calculated from the values of the snow cover thickness Ssn and the thermal conductivity of snow Asn using the relation Rsn=Ssn/Asn, where the value of Asn was determined by the empirical formula [15]
Xsn = 0.3824 x 10-3psn +0.1362.
(8)
Temperature is in degrees Celsius in relations (1)-(8). The SI system of units is used for other parameters.
3. Temperature regime of the reservoir and dynamics of ice cover growth
To determine the temperature regimes of the reservoir and the dynamics of the ice cover various factors should be taken into account. Climatic conditions have obvious impact on the value of heat fluxes. The optical and thermophysical parameters A, e, A and p strongly depend on the state of ice and snow. For example, the characteristic albedo of pure ice, freshly fallen snow, melting or polluted snow, and water are 0.3, 0.8, 0.5, and 0.08 [5], respectively. Moreover, the values of these thermophysical parameters vary with time during the cold period. With the onset
of negative temperatures the water surface is covered with the layer of ice. When precipitation falls a snow cover appears on the layer of ice. The thermophysical properties and thickness of the layer of snow can vary with time. Some bodies of water may not have permanent snow cover as the snow is blown off the ice surface by the wind.
To simulate the dynamics of the ice cover in the absence of snow the following values of optical and thermophysical parameters were set: for water A=0.08, e=0.98, A=0.6 W/(m-K), for ice A=0.3, e=0.98, A=2.2 W/(m-K), q*=335 kJ/(kg-K). Fig. 2 shows the time dependences of temperature during the cold season (October-March). The average daily climatic data are based on the results of long-term measurements. The figure shows the time dependences of temperature of atmospheric air (line 1), temperature on the surface of reservoir (line 2), and temperature at the depth of 0.1 m (line 3).
t, months
Fig. 2. Time dependences of air temperature (1), temperature on the surface of reservoir (2) and temperature at the depth of 0.1 m (3)
The temperature of the ice surface is several degrees higher than the air temperature because of the radiation flux and heat coming from the depth of the reservoir. Minimum ice temperatures are set for the coldest climate conditions. Regions of positive temperature (lines 2 and 3 in Fig. 2) correspond to the liquid state. At negative temperature water is in frozen state. The intermediate state at T«0 corresponds to loose ice saturated with water. In mid-March the average temperature of the surface layer of ice reaches zero degrees while average daily air temperature is negative. This is due to heating by solar radiation. When temperature becomes positive on the surface of the reservoir the layer of water begins to form. At the depth of 0.1 m the ice temperature during the second half of March is about 0 °C that corresponds to the state of wet ice. At the same time, the surface temperature of the ice cover quickly passes the zero level and takes on a positive value.
The presence of a layer of snow on the ice cover significantly affects the thermal balance of the reservoir due to changes in albedo and initiation of additional thermal resistance between air and ice. To simulate the dynamics of the ice cover with a layer of snow the following values of the optical and thermophysical parameters of snow were set: A=0.3, e=0.98, p=300 kg/m3, A=0.25 W/(m-K). At the same time, the values of average daily climatic data remain the same. Under the conditions of a real reservoir, the thickness of the snow cover varies with time during
the cold season. To assess the effect of snow thickness a model profile of the time dependence of snow thickness on the surface of the reservoir was set. From the moment the ice cover appeared the thickness of the snow layer linearly increases for 30 days to a predetermined value and then it remains at a constant level. Fig. 3 shows the simulation results for snow thickness of 0.1 m. The time dependences of the temperature of the atmospheric air (line 1), of the temperature of the snow surface (line 2), of the temperature of the ice surface (line 3) and of the temperature at the depth of 0.1 m (line 4) are shown.
t, months
Fig. 3. Time dependences of the temperature of air (1), of the temperature of the snow surface (2), of the temperature of the surface of the reservoir (3) and of the temperature at the depth of 0.1 m (4)
Comparison Fig. 2 and 3 allows us to estimate the effect of snow cover on the temperature regime of the reservoir. The thermal resistance of the snow layer causes the temperature of the ice surface to be several degrees higher than the temperature of the body of water not covered with snow for the period from the beginning of December to the end of February. In March, on the contrary, the temperature of the ice under the snow is lower than the temperature of the open ice. At the same time, the temperature of the snow surface slightly differs from that of the air. This is due to minor influence of solar radiation which is mainly reflected from the surface because of the high value of the snow albedo.
The temperature regime of the reservoir and the dynamics of ice cover growth are determined by the balance of heat flows. The balance depends to a large extent on the presence of snow and thickness of the snow layer Ssn. The calculated time dependences of the ice thickness d(t) are shown in Fig. 4. Numbers on the graphs indicate thickness of the snow cover in centimetres. It should be noted that at small thickness of snow cover the removal of heat from the reservoir increases and it causes faster increase in the thickness of ice layer than in the case of complete absence of snow (6sn=0). This is because thermal resistance a thin layer of snow is small while the reflection of solar radiation increases sharply. As the snow thickness increases the thermal resistance increases and the associated decrease in the heat flux leads to a slowdown in the growth of the ice layer. A layer of ice of maximum thickness approximately 1 m is observed at ¿sn=0.05 m. The thickness of the ice layer does not exceed 0.6 m at Ssn=0.2 m. In the second half of March when ¿sn=0 the surface layer of ice begins to melt and significant change in the
slope of the graph is observed. It is associated with a decrease in albedo from 0.3 to 0.08 when a layer of water appears on the surface of the ice cover, and solar radiation absorbed by the surface increases.
1.21---t-i-,-T-
t, months
Fig. 4. Thickness of the ice layer versus time, numbers on the curves indicate the thickness of the snow layer
The adequacy of the computational model was verified by comparing the thickness of the ice cover of a fresh lake and a river channel located within the city of Krasnoyarsk with experimental data [16]. The determination of ice thickness was carried out both using remote sensing technology based on satellite navigation systems GLONASS and GPS and by direct measurement using drilling. The ice thickness measured using drilling on February 3 2015 in the lake was 0.66 m, and measured on March 20 2015 in the channel was 0.7 m. The measured values of snow thickness and density were approximately 0.07-0.08 m and 300 kg/m3, respectively. The value of the ice cover thickness that corresponds to the end of winter or the beginning of spring is an important integral indicator. It reflects the impact of the entire set of climatic conditions and other thermophysical factors during the cold period on the temperature regime of the reservoir. Climatic data from the Minino weather station (latitude 56.07, longitude 92.73) were used in simulation. The station is located in the immediate vicinity of the lake where measurements were taken. The time dependence of temperature for the cold season (October 2014 - March 2015) is shown in Fig. 5 (line 1). The coordinates of the measurement site on the river channel are not indicated in the article [16]. The thickness of the snow layer was taken equal to 0.075 m in simulation.
The time dependences of temperature and the thickness of the ice cover were determined. Time dependence of temperature at the depth of 0.1 m is shown in Fig. 5. The time dependence of the ice thickness is shown in Fig. 6, where circles present the measurement data [16]. The results of modelling and measurements differ by approximately 0.04-0.05 m. Various factors can cause the discrepancy between measurement data and results of modelling. As can be seen from Fig. 4, the thickness of the snow layer significantly affects the growth of the ice cover. Therefore, it is the change in the thickness of the snow cover during the cold period that can have the greatest influence. Such changes in the thickness of the snow cover can be caused by both snowfall and wind drift. In addition, the values of the optical and thermophysical parameters of snow can
t, months
Fig. 5. Time dependences of the temperature of atmospheric air (1) and the temperature at the depth of 0.1 m (2)
vary in the case of surface impurity, changes in snow density or some other reasons.
0 Okt. 1 Nov. 2 Dec. 3 Jan. 4 Feb. 5 Mar. 6 t, months
Fig. 6. The ice thickness versus time, circles present the measurement data
Conclusion
Computational modelling of the temperature regimes of the reservoir and the dynamics of the ice cover was carried out on the basis of the numerical solution of the generalized Stefan problem. To describe thermal processes average daily climatic data that correspond to the cold season for the city of Krasnoyarsk were used. Calculations of the temperature regime of the reservoir were made both for the ice cover not covered with snow and for the presence of snow cover of various thicknesses. The influence of the snow cover thickness on the balance of heat
fluxes and the dynamics of ice growth was analysed. It was shown that in the absence of snow the process of ice melting on the reservoir begins earlier. In addition, it was found that with an increase in the snow layer the thickness of the ice cover first increases, and then begins to monotonously decrease. A satisfactory agreement between the results of the computational experiment and measurement data for the thickness of the ice cover on water bodies within the city of Krasnoyarsk was obtained.
References
[1] Z.K.Naurozbayeva , V.A.Lobanov, Methods of short-term forecasting of ice thickness growth in the north-eastern part of the Caspian sea, Geographical bulletin, 54(2020), no. 3, 81-97 (in Russian). DOI 10.17072/2079-7877-2020-3-81-97
[2] I.O.Dumanskaya, Long-term forecast methodology for the ice conditions on the European seas of Russia, Meteorology and hydrology, (2011), no. 11, 64-77 (in Russian).
[3] S.V.Klyachkin, R.B.Guzenko, R.I.May, Numerical model of the ice cover evolution in Arctic Seas for the operational forecasting, Ice and Snow, 55(2015), no. 3, 83-96 (in Russian). D0I.org/10.15356/2076-6734-2015-3-83-96
[4] Yu.M.Georgievskiy, Short-term and long-term forecasts of ice phenomena on rivers, lakes and reservoirs, LGMI, Leningrad, 1986 (in Russian).
[5] M.S.Krass, V.G.Merzlikin, Radiative thermal physics of snow and ice, Gidrometeoizdat, Leningrad, 1990 (in Russian).
[6] L.A.Behovyh, S.V.Makarychev, I.V.Shorina, Fundamentals hydrophysics, Publishing house of the Altai state agrarian University, Barnaul, 2008 (in Russian).
[7] A.A.Samarskij, P.N.Vabishhevich, Computational Heat Transfer, Editorial URSS, Moscow, 2003 (in Russian).
[8] A.A.Samarskii, The theory of difference schemes, Nauka, Moscow, 1989 (in Russian).
[9] E.N.Vasil'ev, V.A.Derevyanko, The dynamics of phase changes in a heat storage of thermal control system for onboard radio-electronic equipment, Thermophysics and Aeromechanics, 25(2018), no. 3, 461-467. DOI: 10.1134/S0869864318030125
[10] G.S.Bordonskii, Electromagnetic properties of ice near the water-ice phase transition temperature, Physics of the Solid State, 47(2005), 715-719.
[11] V.A.Borodkin, S.M.Kovalev, A.I.Shushlebin, Change of structure and some physical properties of level fast ice during the spring and summer period of 2014 in the vicinity the research station "Ice base Cape Baranov", Arctic and Antarctic Research, 65(2019), no. 1, 63-76 (in Russian). DOI:10.30758/0555-2648-2019-65-1-63-76
[12] V.M.Kotlyakov, A.V.Sosnovsky, Estimation of the thermal resistance of snow cover based on the ground temperature, Ice and Snow, 61(2021), no. 2, 195-205 (in Russian). DOI:10.31857/S2076673421020081
[13] N.I.Osokin, A.V.Sosnovsky, Impact of dynamics of air temperature and snow cover thickness on the ground freezing, Earth's Cryosphere., 19(2015), no. 1, 99-105 (in Russian).
[14] T.V.Kirillova, Radiation regime of lakes and reservoirs, Gidrometeoizdat, Leningrad, 1970 (in Russian).
[15] V.M.Kotlyakov, A.V.Sosnovsky, N.I.Osokin, Estimation of thermal conductivity of snow by its density and hardness in Svalbard, Ice and Snow, 58(2018), no. 3, 343-352 (in Russian). DOI: 10.15356/2076-6734-2018-3-343-352
[16] M.I.Mikhailov, K.V.Muzalevskiy, V.L.Mironov, Ice thickness measurements at freshwater lake and river using GLONASS and GPS signals, Current Problems in Remote Sensing of the Earth from Space, 14(2017), no. 2, 167-174 (in Russian).
DOI: 10.21046/2070-7401-2017-14-2-167-174
Моделирование динамики ледового покрова пресноводного водоема
Евгений Н. Васильев
Институт вычислительного моделирования СО РАН Красноярск, Российская Федерация
Аннотация. В работе рассмотрена нестационарная вычислительная модель для определения температурного режима водоема и динамики ледового покрова в холодное время года. Модель основана на численном решении задачи Стефана в обобщенной формулировке. Климатические условия учитывались в виде зависимостей среднесуточных значений температуры воздуха, скорости ветра и солнечной радиации, поглощенной поверхностью водоема. Вычислительное моделирование температурной динамики водоема проводилось для климатических условий в районе г. Красноярска. Выполнен анализ влияния тепловых потоков и толщины снежного покрова на временные зависимости температуры в водоеме и динамику нарастания льда. Проведено сравнение результатов вычислительного эксперимента с данными измерений толщины ледового покрова.
Ключевые слова: задача Стефана, вычислительное моделирование, теплопроводность, динамика ледового покрова, температурный режим водоема.
