Influence of air temperature and snow cover accumulation regimes on ground freezing depth variations in Moscow region

According to developed algorithm and calculating scheme the calculations of ground freezing depth variations for observation sites of meteorological stations of Moscow region (Mozhaysk, Kolomna) were performed on basis of meteorological data on air temperature and snow cover thickness for winter periods of 1988/89-2018/19. The comparison of calculated and available in open access observation data on ground freezing depth for these winter periods was also conducted and indicated good correspondence. The calculating scheme for ground freezing is constructed on basis of three layer media heat conductivity problem (snow cover, frozen and thawed ground) with phase transition on the boundary of frozen and unfrozen ground. Heat balance equation includes phase transition energy, inflow of heat from unfrozen ground and outflow to frozen ground, snow cover and atmosphere. The heat flux is calculated on basis of Fourier law as a product of heat conductivity and temperature gradient. It is supposed, that temperature changes in each media linearly.


Introduction
Thermal and snow cover accumulation regimes of winter periods mostly determine the thermal regime of underlying grounds and its freezing depth. According to the known data on long-term averaged regime of snow accumulation and seasonal variations of air temperature for the particular region and on the base of construction norms the depth of freezing and placement of underground pipelining are determined. However variations in the process of intraseasonal snowfall deposition, accumulation of snow cover and seasonal variations of air temperature in relation to mean values lead to variations of ground temperature, variations of ground freezing depth and hazards for underground pipelining. Also variations of absolute values and the dynamics on the ground freezing depth are important for development of microbiota and root systems of plants. And also, as known, accumulation of snow cover may lead to thawing of permanently frozen ground and increase of seasonal thawed layer and disappearance of permafrost, but the absence of snow cover may lead to occurrence of permafrost in the zones of its sporadic presence or absence.
That is why for determination of influence of air-temperature and snow cover thickness absolute values variations on the dynamics on the ground freezing and thawing depth, the number of models and calculating schemes were developed. For example, V.A. Kudriavcev [1] characterized warming and cooling action of snow cover on the ground depending on snow accumulation regime and on its duration and suggested an equation for estimation of ground freezing depth including snow cover thickness, its thermal properties and amplitude of yearly air temperature oscillations. In the works of A.V. Pavlov [2] snow cover on the ground surface is considered as additional layer with thermal resistance and this is used in the formula for ground freeing depth calculations.

Material and methods
In our case calculating scheme for estimation of ground freezing depth under covered with the snow cover ground surface on basis of air temperature and snow cover thickness data [3][4] is constructed as in [5][6][7][8] and the examples of calculations are performed for the meteorological stations of Moscow region (Mozhaysk, Kolomna) for winter periods of 1988/89-2018/19. The comparison of results of calculating scheme and available data of observations at meteorological stations is performed and indicated good correspondence. The calculating scheme for estimation of ground freezing process is constructed on basis of three layer media heat conductivity problem (snow cover, frozen and thawed ground) with phase transition on the boundary of frozen and unfrozen ground. Heat balance equation includes phase transition energy, inflow of heat from unfrozen ground and outflow to frozen ground, snow cover and atmosphere. The heat flux is calculated on basis of Fourier law as a product of heat conductivity and temperature gradient. It is supposed, that temperature changes in each media linearly like in [9].
The calculations of freezing for covered with snow cover ground surface in winter period on basis of daily data on air temperature and snow thickness allows the estimation of the movement rate of ground freezing interface during the winter period. The rate of movement of ground freezing interface can be expressed with the formulas: or F1 -is heat outflow through frozen ground (and snow cover) from ground freezing interface (W/m 2 ) into atmosphere; cLV = c L dhfg /dτ -heat value for phase transition in the ground, c -ground moisture content (1-4 kg/cm*m 2 ), (last value correspond to full filling of porous by water for light clay with density 2000 kg/m 3 and porosity coefficient 0.617 [10]) L -energy of H2O phase transition (335 kJ/kg), V -rate of movement of ground freezing interface (cm/s); F2 -heat outflow for cooling of thawed ground in front of ground freezing interface (W/m 2 ).
Heat flux is expressed according to Fourier law by means of temperature gradient and heat conductivity as F=-λ (grad T). Heat conductivity and heat flux through combination of two media (snow and frozen ground) according to [11] can be expressed as: Here Tair -air temperature, hs и hfg -snow cover thickness and ground freezing depth, and λs and λfg -heat conductivity of snow and frozen ground and this expression valid also for hs =0.
It was supposed, that on the depth of 10 m in ground there is a point of zero annual temperature oscillation with temperature value T0 about 7°C. That is why Here λthg -heat conductivity of thawed ground. According to [10] averaged heat conductivity of the thawed and frozen ground λthg and λfg -was assumed to be equal 1.4 and 1.8 W/m °C correspondingly.
The differential scheme hfg(tn+1)= hfg(tn)+ Δt V(tn) was constructed by Euler approximation for the equation for the rate of ground freezing depth. By obtained differential scheme the calculations of ground freezing depth were done for the winter periods of 1988/19-2018/19 and the comparisons of obtained results with the available observation data for covered with snow cover site surface of the meteorological stations of Moscow region (Mozhaysk, Kolomna) were performed. And this way the calculations were done with the step-size of one day. For initial conditions, it was supposed that frozen ground thickness hfg was equal 0.5 cm. For each time step (for each day) the rate of movement of freezing interface V and the value frozen ground thickness hfg for the next day (time-step) were calculated.

Results and discussion
The results of calculations of maximal ground freezing depth for the sites of the meteorological observatory of Moscow region (Mozhaysk, Kolomna) for the winter periods of 1988/89-2018/19 and theirs comparison with the observed data are displayed on the Fig.  1 and indicate general consistency. The example of results of calculations for ground freezing depth for the covered with snow cover site of meteorological station of Moscow region (Mozhaysk) for the winter period 2009/10 and theirs comparison with the observed data is displayed on the Fig. 2.

Conclusions
The calculating scheme was constructed by means of approximation of derived ordinary differential equation for ground freezing depth variation with numerical scheme. On basis of obtained numerical scheme the calculations were performed and the comparison of obtained results with observed values were done.
The applied method is well physically consistent. The obtained by method solution indicates good correspondence with the process of ground freezing variation during the winter period. The important point for successful work of method is good setting up of initial conditions. Considered in this work method of linear gradients (like in [9]) differs from considered before in the classical book of A.N. Tikhonov and A.A. Samarskii [12] (first edition from 1951) or in the works of A.V. Pavlov [2] or presented in the other works method, where for calculation of seasonal dynamics of ground freezing depth the heat conductivity partial differential equation of second order for space and first order for time is used. In this work there is only reduced ordinary differential equation of first order for time is used. Solving of this reduced first order ODE is simpler and could be done within Excel program.
The possible way of improvement of the method's accuracy is taking into consideration of density and heat conductivity of each layer of snow cover and determination of bulk heat conductivity according to the formula of multi-layer media.
The work was performed in a frame and under support of state topic «Mapping, modelling and assessment of risk of hazardous natural processes» AAAA-A16-116032810093-2.