Analysis of Different Freezing/Thawing Parameterizations using the UTOPIA Model

Abstract

Soil moisture changes are generally due to external factors (precipitation, evaporation, etc.) and internal forces (gravitational force, capillarity, transpiration, etc.). When soil temperatures remain below 0 °C for a long time (hours or even entire consecutive days), part of the liquid water content of the soil can freeze, thus freezing/thawing effects must be taken into account in those conditions. The present work is devoted to the numerical modeling of the water phase change in the soil. The model used in this study for the land surface processes is UTOPIA (University of TOrino land Process Interaction in Atmosphere) model, which is the updated version of LSPM (Land Surface Process Model). Scientific literature proposes some formulations to account for freezing/thawing processes. Three different parameterizations have been compared using a synthetic dataset in order to assess which one performs best from a physical point of view. Parameterizing freezing/thawing processes creates numerical instability and water overproduction in the UTOPIA model. These problems have been solved and described in the paper by means of synthetic data created to test the new parameterizations. The results show that UTOPIA is able to capture the freezing/thawing physical processes.

1. Introduction

It is well known that the fluxes between soil and atmosphere can trigger precipitation and convection, allow fog or dew development in the atmospheric surface layer, enable frost on vegetation and soil, and strongly modify the soil water content and temperature [1].

A wrong or an imperfect soil moisture estimation may lead to errors in the boundary layer estimation [2]. Moreover soil property changes (hydraulic conductivity, thermal capacity etc.) during the freezing/thawing transient stage may lead to variations in the hydrological balance components [3,4].

The goal of this study is to improve the SVAT (Soil Vegetation Atmosphere Transfer) scheme UTOPIA (University of TOrino land Process Interaction in Atmosphere) parameterizing the freezing/thawing in the soil. This model is an evolution of the LSPM (Land Surface Process Model), which was developed at the University of Torino [5], and is now used operationally at the Regional Agency for Environmental Protection—Arpa Piemonte, and also for climatological studies [6].

The improvement of the performances of such kind of models has practical importance in several disciplines. In meteorology, a better estimation of the boundary layer development is linked to the prediction of the atmospheric stability, and leads to better performances of GCMs (General Circulation Models) or LAMs (Limited Area Models) in predicting two meters temperature, relative humidity, wind speed and precipitation [7]. In agriculture, the knowledge of the temperatures in the root zone and in the canopy layer may eventually allow activation of practices in order to avoid plant freezing. In hydrology, it could allow a better estimation of the runoff associated with rainfall and snowmelt on frozen soil, or to determine the portion of soil frozen (for instance, related to the permafrost).

2. The UTOPIA Model

The UTOPIA (University of TOrino land Process Interaction in Atmosphere) model was evolved from another model called LSPM (Land Surface Process Model). It is a diagnostic one-dimensional model developed at the University of Torino [5] since 1989, and can be categorized in the class of the SVAT schemes. UTOPIA can be used both as a stand-alone model or coupled with an (general or mesoscale) atmospheric circulation model, behaving as its lower boundary condition.

The UTOPIA domain is vertically subdivided into three main zones: the soil, the vegetation and the atmospheric layer within and above the vegetation canopy layer. Variables are mainly diagnosed in the soil and in the vegetation layers. The canopy itself is represented as a single uniform layer (big leaf approximation), whose properties are described by vegetation cover and height, leaf area index, albedo, minimum stomatal resistance, leaf dimension, emissivity and root depth. The soil state is described by its temperature and moisture content, and its main parameters are: thermal and hydraulic conductivities, soil porosity, permanent wilting point, dry soil heat capacity, surface albedo and emissivity. These variables are calculated by integrating the heat diffusion (Fourier) and the water mass conservation equations. The soil is discretized using a multi-layer structure and solving the equations using the Crank-Nicholson tridiagonal method. The possible presence of snow, considered as a single layer and characterized by density, albedo, water equivalent and ice content, is taken in account with specific routines.

Momentum, heat and water vapor fluxes between the soil (bare and vegetated, and eventually covered by snow) and the reference height, above the canopy, are computed using an electric analog formulation, for which they are directly proportional to the gradients of the related scalars and inversely proportional to the adequate resistance.

Since the UTOPIA is a diagnostic model, observations in the atmospheric layer at its upper border (i.e., the surface layer above the canopy) are needed as boundary conditions. These are: air temperature, humidity, pressure, wind speed, cloud cover, long- and short-wave incoming radiation, and precipitation rate. Usually these data are measured, eventually reconstructing the missing data with adequate interpolation techniques, or furnished by the atmospheric model in case of coupling.

All specific details about UTOPIA model physical parameterizations, its use and its features are fully described in [8].

3. The Freezing Parameterization in UTOPIA

Three different parameterizations have been implemented in UTOPIA to describe freezing processes in the soil. The first one, developed by Schrödin et al. [9], is used in the frame of the COSMO model. The second one, developed by Viterbo et al. [2], is used in the frame of the ECMWF general circulation model. The third one, developed by Cassardo [8], is based on the previous one but contains different parameterizations. They are briefly summarized below.

3.1. Energy vs. Temperature

This parameterization, hereafter called “S” and based on the paper of Schrödin et al. [9], imposes the energy balance conservation in the soil. The energy available for the phase change is:

$$\Delta E={(\rho c)}_s\Delta z(T-T_0)$$

(1)

where (ρc)_s is the soil volumetric thermal capacity, Δz is the soil layer depth and T₀ is the freezing point of water (0 °C). The maximum possible change of liquid water (ice) Δη_w,max (Δη_i,max) is:

$$\Delta {\eta }_{w,max}=-\Delta {\eta }_{i,max}=\frac{\Delta E}{L_f{\rho }_w}$$

(2)

where L_f is the latent heat of fusion and ρ_w is the water density of water (999.8 kg/m³). The ice content η_i, depends, as shown by the formula (3), on the soil energy available for the phase change, also on the soil water content η_w (if freezing is taking place) or on the soil ice content η_i (if thawing is taking place):

$$\begin{array}{l}\Delta {\eta }_i=-MIN(-\Delta {\eta }_{i,\max },{\eta }_i)\begin{array}{c}if\end{array}\begin{array}{c}\Delta {\eta }_{i,\max }<0\end{array}\\ \Delta {\eta }_i=MIN(\Delta {\eta }_{i,\max },{\eta }_w)\begin{array}{c}if\begin{array}{c}\Delta {\eta }_{i,\max }>0\end{array}\end{array}\end{array}$$

(3)

Finally, the temperature at the current time-step is

$$T=T_0+(\Delta {\eta }_i-\Delta {\eta }_{i,max})\frac{L_f{\rho }_w}{{(\rho c)}_s\Delta z}$$

(4)

where Δη_i is the actual change of ice content.

3.2. Thermal Capacity vs. Temperature

This parameterization, hereafter called “V” and based on the paper of Viterbo et al. [2], modifies the soil thermal capacity according with the quantity of water and ice in the soil. The soil energy equation in the presence of soil water phase change is evaluated as:

$${\left(\rho C\right)}_s\frac{\partial T}{\partial t}=\frac{\partial }{\partial z}\left(k_T\frac{\partial T}{\partial z}\right)+L_f{\rho }_w\frac{\partial {\eta }_i}{\partial t}$$

(5)

where k_T is the soil thermal conductivity. The volumetric ice content is evaluated by:

$${\eta }_i=f\left(T\right){\sigma }_v{\eta }_{fc}$$

(6)

where σ_v is the vegetation cover, η_fc is the soil water content at the field capacity, and f(T) is a function of temperature given by:

$$f\left(T\right)=\{\begin{array}{ll}0& T>T_1\\ 0.5[1-sin(\frac{\pi \left(T-0.5T_1-0.5T_2\right)}{T_1-T_2}\left)\right]& T_2\le T\le T_1\\ 1& T<T_2\end{array}$$

(7)

where T = T₀ + 1K e T₂ = T₀ - 3K, and T₀ is the freezing point of water. The behavior of f(T) is shown in Figure 1.

Figure 1. Behavior of the function f(T) in the Viterbo parameterization.

Figure 1. Behavior of the function f(T) in the Viterbo parameterization.

3.3. Thermal Capacity vs. Temperature

This parameterization, hereafter called “C” and based on the paper of Cassardo [8], is based on the same physical considerations of the one of Viterbo et al. [2], but the parameterization of the ice water content η_i is given by:

$${\eta }_i=f(T)\eta -{\eta }_{min}$$

(8)

where η_w,min is the minimum quantity of water in soil (set to 0.01 m³/m³ in order to avoid undesired numerical problem when soil water content tends to zero) and η = η_w + η_i. Thus, in this case, the ice water content does not depend on the vegetation cover, and is limited by the actual liquid water content.

4. Sensitivity Experiments Using Synthetic Data

A set of data, hereafter called “synthetic data”, has been created as a workbench in order to verify the correct estimation of the physical mechanisms taking place during the water phase change.

This data consist in an idealized year (360 days) subdivided into three well distinct periods. The initial period, lasting 27 days, is characterized by a daily mean temperature permanently above 0 °C (about 4.5 °C, Figure 2) and a continuous precipitation (daily mean precipitation rate of about 0.001 mm/s corresponding to 86.4 mm/day, Figure 3), in order to saturate the upper soil layer and achieve the numerical stability. The intermediate period, lasting 280 days, is characterized by the absence of precipitation (in order to avoid snow accumulation on the ground which might affect the ice behavior in the soil) and by a mean air temperature well below 0 °C (about −4.5 °C), in order to provoke the freezing in the upper soil layers. The final period, lasting 50 days, is characterized by the same conditions of the first period (average air temperature above 0 °C and continuous precipitation), in order to force the thawing of the soil ice formed during the intermediate period.

Cloud cover has been set to overcast conditions when precipitation occurs, clear sky conditions elsewhere.

The remaining meteorological fields necessary to run UTOPIA model have been chosen in order to repeat the same day throughout the year. In other words, the same diurnal cycle is repeated during all the simulations not considering an annual cycle typical of the real world. The reason for such a choice is to investigate the main mechanisms taking part of the water phase change avoiding those strictly related to the annual cycle. More in-depth, the minimum and maximum values for the main meteorological fields during their repeated diurnal cycle are:

• Relative humidity (%): min = 61.1; max = 91.5.
• Wind velocity (m/s): min = 2.0; max = 5.0.
• Overcast conditions (W/m²): max radiation = 74.65; clear sky radiation max = 248.85 .
• Pressure: min = 1022.17; max = 1023.73.

Figure 2. Daily mean 2 m temperatures used in the synthetic data.

Figure 2. Daily mean 2 m temperatures used in the synthetic data.

Figure 3. Daily mean precipitations used in the synthetic data.

Figure 3. Daily mean precipitations used in the synthetic data.

Table 1. Model setting for the simulation with the synthetic data.

SYNTHETIC DATA – UTOPIA setting
Depth of the six soil layers (cm, top to bottom): 20, 20, 40, 80, 160, 320
t = 60 s	Initial w = 0.4 m3/m3	Initial T = 1 °C
Soil type: silty clay
Soil porosity s [m3/m3]	Wilting point wi [m3/m3]	Field capacity fc [m3/m3]
0.492	0.283	0.4
Vegetation type: short grass Vegetation cover v = 0.8

Table 1. Model setting for the simulation with the synthetic data.

SYNTHETIC DATA – UTOPIA setting
Depth of the six soil layers (cm, top to bottom): 20, 20, 40, 80, 160, 320
t = 60 s	Initial w = 0.4 m3/m3	Initial T = 1 °C
Soil type: silty clay
Soil porosity s [m3/m3]	Wilting point wi [m3/m3]	Field capacity fc [m3/m3]
0.492	0.283	0.4
Vegetation type: short grass Vegetation cover v = 0.8

5. Results

The results (Figure 4 and Figure 5) show that the model is able to reproduce the physical processes of both freezing and thawing using all three parameterizations above mentioned. In particular, Figure 4 shows the warming of the layer affected by freezing with respect to ‘no freezing’ (“NF”) case (i.e., when freezing is not activated), due to the latent heat release, while Figure 5 shows the decrease of the liquid water content η_w during the freezing period. The warming of the layer in which ice forms is proportional to the quantity of ice present in the soil: this fact is not surprising, as the heating is linked to the latent heat release. In this sense, V and C parameterizations simulate the warmest temperatures and the biggest ice quantity.

Figure 4. Simulated soil temperature in the 10 cm soil layer.

Figure 4. Simulated soil temperature in the 10 cm soil layer.

Figure 5. Simulated soil moisture in the 10 cm soil layer. The cyan line is the porosity η_s, the pink line is the wilting point η_wi.

Figure 5. Simulated soil moisture in the 10 cm soil layer. The cyan line is the porosity η_s, the pink line is the wilting point η_wi.

5..1 Numerical Instability during the Thawing Phase

Figure 4 and Figure 5 show soil temperature and soil moisture oscillations during the thawing phase. This is due to two main factors: the first is that precipitation begins during the thawing period; the second is related to numerical instability, especially for the parameterization C.

Inspection of this problem revealed that the numerical part of this instability is caused by a single term.

In fact, the soil moisture is calculated from the conservation equation:

$$\frac{\partial {\eta }_w}{\partial t}=-\frac{\partial q_{wz}}{\partial z}$$

(9)

where q_wz is the soil water flux.

In this equation, the soil water flux, responsible for the above-mentioned numerical instability, is expressed by the following equation implemented in the UTOPIA model [10]:

$$q_{wz}=-(D_{w\eta }+D_{v\eta })\frac{\partial {\eta }_w}{\partial z}-K_{\eta }-D_{vT}\frac{\partial T}{\partial z}$$

(10)

The first term on the right-hand-side of the equation is the soil water flux, related to the humidity gradient. This flux depends on the hydraulic diffusivity of liquid water D_wη and from the hydraulic diffusivity of water vapor D_vη. The second term is the water flux due to gravitational drainage, while the last one is the water vapor flux due to the temperature gradient.

In particular, the term responsible of numeric instability of soil moisture and temperature is the water vapor flux due to humidity gradient ($-D_{v\eta }\frac{\partial {\eta }_w}{\partial z}$). Niu et al. [11] showed that this flux can be considered negligible for values of volumetric water content larger than 0.06 m³/m³. For this reason, the water vapor hydraulic diffusivity in the equation of soil water vertical flux has been set to zero for η_w > 0.06 m³/m³, leading to a noticeable reduction of the instability (Figure 6 and Figure 7).

Figure 6. Soil temperature in the 10 cm soil layer (as for Figure 4) for the simulation performed when the water vapor hydraulic diffusivity is set to zero.

Figure 6. Soil temperature in the 10 cm soil layer (as for Figure 4) for the simulation performed when the water vapor hydraulic diffusivity is set to zero.

Figure 7. Soil moisture in the 10 cm soil layer (as for Figure 5) for the simulation performed when the water vapor hydraulic diffusivity is set to zero

Figure 7. Soil moisture in the 10 cm soil layer (as for Figure 5) for the simulation performed when the water vapor hydraulic diffusivity is set to zero

5.2. The Problem of Water Overproduction during Soil Freezing

A second relevant problem encountered during the sensitivity tests concerning the phase change of water in the soil is the lack of conservation of the total soil water mass. This fact appears evident by considering the total volumetric water content (water and ice) in all model soil layers considered in a certain simulation.

$${\eta }_{tot}={\displaystyle \sum _{i=1}^N{\eta }_w(i)+{\eta }_i(i)}$$

(11)

In the η_TOT calculation, we considered the water density equal to 999.8 kg/m³ and the ice density equal to 917.0 kg/m³.

Figure 8 shows η_TOT values corresponding to the simulation referring to Figure 6 and Figure 7. It is evident that, during the period of the soil freezing, the increase of η_TOT is not justified on the basis of considerations related to hydrological soil balance. In fact, in the freezing period, there is no precipitation, and the transpiration stops because the temperatures drop below 0 °C, thus minimizing also the evaporation. Considering also the drainage from the bottom layer, a decrease, and not an increase, of total water content may be considered realistic.

Several analyses have shown that the term responsible of the observed water overproduction is the liquid water flux due to the humidity gradient ($-D_{l\eta }\frac{\partial {\eta }_w}{\partial z}$). This last term describes capillarity. In fact, when the soil starts to freeze, capillarity takes place because the ice formed in the soil occupies part of the soil pores volume.

To demonstrate the relation between the water overproduction and the capillarity, we set all terms of the hydrological balance equation to zero, except this flux. If only surface and intermediate runoff are allowed (Figure 9), an unreasonable increase of η_tot is still present.

Figure 8. η_tot trend for the simulation with synthetic data.

Figure 8. η_tot trend for the simulation with synthetic data.

Figure 9. Representation of the moisture fluxes between the different soil layers and the external environment. During the sensitivity test discussed in Section 5.2, transpiration evaporation, infiltration and gravitational drainage have been set to zero.

Figure 9. Representation of the moisture fluxes between the different soil layers and the external environment. During the sensitivity test discussed in Section 5.2, transpiration evaporation, infiltration and gravitational drainage have been set to zero.

This behavior has been identified as a problem of numeric instability due to very fast variation of volumetric water content in the soil layers during the water-ice phase change. A reduction of the integration time-step in the UTOPIA model has resulted in a decrease of the water overproduction, as shown in Figure 10. However, this approach for the solution of this problem cannot be adopted because it requires a very long computing time for the simulations.

Figure 10. η _tot trends evaluated in the sensitivity test summarized in Figure 9. The three simulations have been executed using three different time step (DTSEC) values. The numerical convergence toward a constant value is just slightly shown. The complete numerical convergence would have been reached shorter values of DTSEC, however this is not shown because it requires too large a computation time.

Figure 10. η _tot trends evaluated in the sensitivity test summarized in Figure 9. The three simulations have been executed using three different time step (DTSEC) values. The numerical convergence toward a constant value is just slightly shown. The complete numerical convergence would have been reached shorter values of DTSEC, however this is not shown because it requires too large a computation time.

An alternative approach has been adopted considering for each time-step the calculation of two variables for η_tot (in units of volumetric water content). The first one, called η_tot;correct, conserves the water mass during the freezing, while the second variable, called η_{tot;effective}, does not conserve the water mass.

Simulation results show (Figure 11 and Figure 12) that the physically and hydrologically unreasonable increment of the total VSWC occur in presence of a water flux, due to the gradient humidity, towards soil layers with a significant ice production (η_i > 0.01 m³/m³). Thus η_tot begin to increase when the total VSWC of the second layer, where the volumetric ice content is greater than 0.01 m³/m³, is considered.

Figure 11. η_tot trend relative to the sensitivity test summarized in Figure 9, considering: the three deepest soil layers (top left), the four deepest soil layers (top right), the five deepest soil layers (bottom left), and all soil layers (bottom right).

Figure 11. η_tot trend relative to the sensitivity test summarized in Figure 9, considering: the three deepest soil layers (top left), the four deepest soil layers (top right), the five deepest soil layers (bottom left), and all soil layers (bottom right).

Figure 12. η_i (volumetric ice content) trends relative to the sensitivity test summarized in Figure 9, for the fourth layer (top left), the third layer (top right), second layer (bottom left) and first layer (bottom right).

Figure 12. η_i (volumetric ice content) trends relative to the sensitivity test summarized in Figure 9, for the fourth layer (top left), the third layer (top right), second layer (bottom left) and first layer (bottom right).

At the beginning of each time step, the quantity η_i,tot, i.e., the total volumetric ice content in all soil layers, is calculated. If η_i,tot is greater than the minimum value 0.01 m³/m³ (i.e., if there is a significant quantity of ice in soil), the variable η_tot,correct is computed by setting to zero the hydraulic diffusivity (D_wη= 0), in order to neglect the water fluxes among the soil layers due to the humidity gradient. In this way, the total VSWC (η’_w+η’_i) of each layer contributes to conserve the total VSWC. Thus, the variable η_tot;correct is calculated as:

$${\eta }_{tot,correct}={\displaystyle \sum _{j=1}^N\eta {\text{'}}_w(j)+\eta {\text{'}}_i(j)}$$

(12)

at the same time, the total VSWC of all soil layers η_{tot,effective} is calculated by accounting the true hydraulic diffusivity (D_wη≠ 0) and it does not conserve the total VSWC:

$${\eta }_{tot,effective=}{\displaystyle \sum _{j=1}^N{\eta }_w(j)+{\eta }_i(j)}$$

(13)

At this point, as stated before, among all soil layers, those who do not conserve the total VSWC possess a significant ice production. Thus, the following correction on η_w is performed only in those soil layers in which the volumetric ice content is significant (with η_i > 0.01 m³/m³):

$${\eta }_{w,new=}{\eta }_{w,old}\frac{{\eta }_{tot,correct}}{{\eta }_{tot,effective}}$$

(14)

This last correction is inserted in a cycle of 30 iterations per time step. For every iteration the variable η_{tot,effective} is recalculated from the η_w,new values in the previous iteration based on the following scheme:

$${\eta }_{w,new=}{\eta }_{w,old}\frac{{\eta }_{tot,correct}}{{\eta }_{tot,effective}}{\eta }_{w,old}={\eta }_{w,new}\times 30\text{ times}$$

(15)

This cycle was introduced to ensure the convergence of η_{tot,effective}, as shown in Figure 13.

Figure 13. Trends of η_tot;correct for parameterization V without correction (black line) and of η_{tot;effective} with corrections (other lines) and with different iterations numbers.

Figure 13. Trends of η_tot;correct for parameterization V without correction (black line) and of η_{tot;effective} with corrections (other lines) and with different iterations numbers.

Figure 14, Figure 15 and Figure 16 show the total volumetric water content for all soil layers, without (blue line) and with (red line) the above-mentioned correction. In the second case, the trend of η_tot is justified on the basis of consideration related to hydrological soil balance. The decrease of η_tot is mainly due to drainage from the bottom layer. In fact, in the period with temperature below zero degrees, the transpiration is stopped and evaporation is reduced.

Figure 14. Trends of η_tot;correct for simulations without correction (blue line) and of η_{tot;effective} with correction (red line). In these simulations each soil layer has been initialized with a value of η_w = 0.4 m³/m³.

Figure 14. Trends of η_tot;correct for simulations without correction (blue line) and of η_{tot;effective} with correction (red line). In these simulations each soil layer has been initialized with a value of η_w = 0.4 m³/m³.

Figure 15. Same as Figure 14 but for the parameterization V.

Figure 15. Same as Figure 14 but for the parameterization V.

Figure 16. Same as Figure 14 but for the parameterization C.

Figure 16. Same as Figure 14 but for the parameterization C.

Another diagnostic to analyze the water overproduction problem is to look at the hydrological balance given by:

$$\frac{\Delta \eta }{\Delta t}-(I+E+R+D)=0$$

(16)

In the previous formula $\Delta \eta ={\eta }_{tot}^i-{\eta }_{tot}^{i-1}$, where i and i-1 are the current and the previous timesteps, respectively, Δt is the model timestep, I is the precipitation, E the evapotranspiration, R the sum of intermediate and surface runoff and D is the drainage from the bottom soil layer.

Figure 17 shows the hydrological balance and its components for simulation with and without the correction of the water overproduction correction for the parameterization V. Without the correction there is an unbalance of the hydrological balance during the transition phases (freezing and thawing). Implementing the correction in UTOPIA model, the hydrological balance and hence the mass conservation are respected.

Figure 17. Hydrological balance and its components with (a, b) and without (c, d) the correction for the water overproduction and relative to the parameterization V. Notice the different scales of the plots (10^-6 for the hydrological balance, 10^-7 for the components).

Figure 17. Hydrological balance and its components with (a, b) and without (c, d) the correction for the water overproduction and relative to the parameterization V. Notice the different scales of the plots (10^-6 for the hydrological balance, 10^-7 for the components).

6. Conclusions

Water phase changes in the soil are very important for several reasons. In this work, some freezing parameterizations have been implemented and tested in the UTOPIA land surface model, searching the most performing one, if any. To this end, an idealized experiment has been designed and executed in order to initialize and guide the model.

The introduction of the freezing/thawing parameterizations in UTOPIA has revealed that they work properly, but two numerical problems became evident: the intrinsic numerical instability of the scheme, and the failure of the total soil water content conservation during the freezing/thawing phases.

The correction designed to eliminate the numerical instability during the thawing phase, and the algorithm designed to avoid the water overproduction in the soil freezing processes have allowed a more precise estimate of the soil moisture during freezing and thawing periods. This factor is very important in order to get more accurate estimations of hydrological resources in the winter and springs periods.

The final version of UTOPIA is able to catch the transient phase related to the soil water phase changes. All parameterizations show the impact of the latent heat release/absorption during the freezing/thawing phases on the soil temperature [12].

Comparisons with observations are useful and necessary in order to estimate the importance of the soil freezing/thawing for a correct estimation of the energy and hydrologic budgets in the surface layer.

In the future developments, another aspect to be taken into account is the inclusion of more appropriate parameterizations for the soil property changes induced by the freezing/thawing phenomena, such as the thermal and hydrological changes of the soil due to the presence of ice [13,14].