*Article* **Probabilistic Analysis as a Method for Ground Freezing Depth Estimation**

**Tomasz Godlewski 1,\*, Łukasz Wodzy ´nski <sup>2</sup> and Małgorzata Wsz ˛edyrówny-Nast <sup>1</sup>**

**Abstract:** Accurate frost depth prediction is an important aspect in different engineering designs such as for pavements, buildings, bridge foundations, and utility lines. This paper presents a probabilistic method of assessment of the depth of soil freezing. Annual (winter) maxima of the position of the zero centigrade temperature measured in the soil were approximated by Gumbel probability distribution. Its parameters were estimated using maximum likelihood method. The results received on the basis of data from 36 meteorological stations in Poland and 50 years of observations, as characteristic values with 50-year return period, reflect the influence of the climatic conditions on the freezing depth. On the other hand, the soil structure and its conditions also play an important role in freezing. Nowadays they may be taken into account using correction coefficients. It is concluded that this method is more precise than a method using the air freezing index because through the use of direct measurements it takes into account additional factors affecting the actual depth of freezing. The obtained results are not the same as those given in the older Polish Standard which was based on the simplified and limited data. The results confirm the impact of climate change on ground freezing depth.

**Keywords:** soil freezing; probability; return period; reliability; correction coefficient; depth of ground freezing

#### **1. Introduction**

Freezing of the soil in wintertime is one of the climatic actions which must be considered when designing structures founded on soils. The actual frost depth is affected by the material type, soil thermal properties, soil water content, and climatic conditions such as temperature, wind speed, precipitation, and solar radiation. Frost depth can be estimated using numerical or analytical modeling technique (Gontaszewska [1], Konrad & Shen [2], Rajaei & Baladi [3], Venäläinen et al. [4], ). In the past few years different numerical techniques have been used for modeling transient heat flow in pavement layers. These methods require a large number of input data [4,5]. In most cases the input data are either not available or just too expensive to collect. In the case of missing data, higher or similar accuracy can be obtained from simpler analytical and semi-empirical models to estimate frost depth prediction—for e.g., Stefan model, modified Berggren model or Chisholm and Phang empirical model [3,5]. In those methods, the depth of frost penetration has usually been defined using the so-called air freezing index. However, it is obvious that the frost penetration into the soil is a complex random process in which many factors affect the final results, which are often divergent from the observations [6]. Those are not only air temperature but also the depth and the thermal properties of snow cover, which is changing during winter, heat contained in the soil [7,8] and the soil properties, its structure and water content. All the climatic conditions have a random character and should be analyzed using probabilistic methods. Nowadays, the simplest approach is the probabilistic analysis of the annual (winter) maximum depth of the zero centigrade temperature. However, it

**Citation:** Godlewski, T.; Wodzy ´nski, Ł.; Wsz ˛edyrówny-Nast, M. Probabilistic Analysis as a Method for Ground Freezing Depth Estimation. *Appl. Sci.* **2021**, *11*, 8194. https:// doi.org/10.3390/app11178194

Academic Editor: Daniel Dias

Received: 26 July 2021 Accepted: 30 August 2021 Published: 3 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> Building Research Institute, 02-656 Warsaw, Poland; m.wszedyrowny@itb.pl

<sup>2</sup> Institute of Civil Engineering, Warsaw University of Life Science, 02-776 Warsaw, Poland; lukasz\_wodzynski@sggw.pl

**<sup>\*</sup>** Correspondence: t.godlewski@itb.pl; Tel.: +48-22-56-64-163

should be emphasized that it is not the same as the depth of the soil freezing. The purpose of the paper is to present the probabilistic analysis of the maximum annual depth of the zero centigrade temperature of the soils, the implemented method and some examples of the obtained results.

Since 1955, the maps of the soil freezing depths have been given in Polish Standards. The values were calculated in centimeters according to the formula given in Soviet recommendations [9]:

$$h\_z = 23 \cdot \sqrt{\sum m + 2} \tag{1}$$

where: *m*—an absolute value of the winter monthly mean negative air temperature; only negative values are summarized.

The numerical factor 23 means the type of soil. The value 23 is for loams and clays. Numbers for other types of soil were not given in the Polish Standards. According to Zavarina [10], this formula was derived on the basis of the experiments made at six Soviet meteorological stations. A similar formula but with more detail was elaborated at the US Army Cold Regions Research and Engineering Laboratory [6].

The results of the calculations according to the Formula (1) were implemented into the Polish Standard in 1955 (PN-55/B-03020) and its next editions with small modification in 1974. For the western half of Poland, the depth of freezing was 0.8 m, for the eastern part it was 1.0 m and for small area around Suwałki region it was 1.4 m.

The approach, based on the air freezing index, is still recommended [1]. However, now it is possible to define the depth of the zero centigrade temperature on the basis of direct measurements of the soil temperature.

#### **2. Material and Methodology**

#### *2.1. Measurements of Soil Temperature*

The Earth's surface is heated by the absorption of solar radiation. The variety of the substrate results in different heating mechanism. Snow and ice covering the ground reflects most of the radiation [11]. Water, due to its partial transparency, absorbs radiation throughout the layer to which it reaches, and the ground absorbs radiation in a thin surface layer [12]. Changes in ground temperature depend on its vegetation coverage, chemical composition, and to a large extent on water and air content. The source of heat in the ground can also be heat exchange with the atmosphere; thawing and freezing processes of water contained in the ground; precipitation heat and biological, chemical, and physical processes in the soil [13]. Soil, which is not covered with vegetation, is most strongly heated during the day in clear weather; at night in the same conditions, it intensely radiates heat. The size of the heat flux absorbed or given away by the soil depends on the temperature difference and the thermal capacity of the soil [14,15]. The thermal capacity of the soil is influenced by its mineral composition, water and air content. There is a very large impact of the snow cover on the depth of the zero isotherm position. The maximum positions of zero isotherms in the ground without snow cover are even more than twice as deep as in snow-covered soils [16,17].

Soil freezing (frost) depth can be determined using a number of direct approaches, for example, by direct measurement by thermometers installed in the ground. Another approach is to use the concept of 'degree days'—freezing degree days (FDDs) and thawing degree days (TDDs)—as a surrogate for frost depth. The advantage of the 'degree day' approach is that meteorological data is easier to obtain and less intrusive than the other described approaches. A limitation of the degree-day approach is the representativeness of temperature data taken off-site [18].

In Poland, the soil temperature is measured every day at the meteorological stations of the Institute for Meteorology and Water Management—State Research Institute (IMGW-PIB). The measurements were taken in the last 40 years using mercury-in-glass thermometers, at 6:00 and 18:00 UTC at the depths of 5 cm, 10 cm, 20 cm, and 50 cm and at the depth of 1.0 m but only at 12:00 UTC [19]. Measurements at this last depth have been made only since 1982. Currently, the measurements are automated, performed every

minute with an accuracy of 0.1 ◦C. Measurements of the minimum temperature above the surface of the exposed soil and the ground temperature are made on the testing site. This measuring field is a part of the ground with area of 2 × 4 m; the longer sides of the plot are directed in the direction E-W. Experimental plot are free of the grass. The measurement technique for the soil temperature and the method of calculation of the zero centigrade depth were described in the instruction for meteorological measurements [20,21].

Freezing the ground is a physical process consisting of the freezing of water present in the ground and its molecules into a hard permafrost. The water contained in the soil is a solution of different salts of varying concentrations, the higher their concentration, the lower the ground freezing temperature. This temperature may drop to −2 ◦C and below. Freezing of the ground depends on the air temperature, the duration of the negative temperatures, the thickness and type of snow cover, the terrain profile, and the structure and humidity of the ground. The low air temperature that persists for a long time lowers the temperature of deeper and deeper layers of the ground, while strong but short-term temperature drops have less impact on the freezing depth of the ground. An important role is played by the snow cover: a thicker, less dense and even layer of snow inhibits the process of freezing the ground. Different types of soils are characterized by different ease of freezing. The process of ground freezing consists of solar radiation, the advection of warm air masses, and heat penetrating from deeper layers of soil. The timing of freezing of the first layers of soil, the length of permafrost persistence and the maximum freezing depth are of great importance in construction, agriculture, and water management. At weather stations measuring ground temperature at depths of 5, 10, 20, 50, and 100 cm, isotherm "0" in the ground is calculated. Isotherm (ET) "0" is the plane connecting the points at 0 ◦C in the ground [19]. Depending on the course of the weather, there may be one or more "0" isotherms in the ground. On the basis of measured everyday values of soil temperature personnel of weather stations calculate the position of the zero centigrade temperature. Following formula is used:

$$Z = A + \frac{a \cdot \Delta T\_A}{\Delta T} \tag{2}$$

where: *A*—the depth of the lowest place of thermometer which measured negative temperature; *a*—the distance between depths of thermometers measuring negative and positive temperature; Δ*TA*—temperature at the level *A*; Δ*T*—the difference between temperature values at the levels with zero centigrade temperature between them.

An example is presented in Figure 1 below: *A* = 0.2 m, *a* = 0.5 − 0.2 = 0.3 m, Δ*TA* = −0.4 ◦C, ΔT = −0.8 ◦C. Zero-degree temperature is at the depth of 0.35 m.

**Figure 1.** Example of temperature of the soil at the meteorological station in Białystok for one day (15 March 2005) adapted from [19].

The problem may occur when the negative soil temperature exceeds 1.0 m depth. In such a case it is recorded as below 1.0 m. The true value can be defined using extrapolation of values measured at smaller depths. An example is presented in Figure 2. However, it must be pointed out that such cases were not frequent.

**Figure 2.** Example of extrapolation of the data on the soil temperature made in order to define maximal winter value adapted from [19].

When determining the position of a zero isotherm by extrapolation, it should be taken into account that as the temperature rises in the upper layers of the ground at an unchanged value at the depth of the largest measurement, the slope of the straight regression decreases. The use of such a set may lead to an error in the determination of the extrapolated value.

#### *2.2. Methods of Probabilistic Analysis*

The aim of the analysis is to find such values of the position of zero centigrade temperature which would be exceeded with the accepted probability in appropriate reference period. Nowadays, the commonly accepted annual probability of exceedance of climatic actions is 0.02 which means 50-year return period. Such an assumption may also be used for the probabilistic calculation of freezing depths. The probability that a value *Zk* will be exceeded once in *t* years is:

$$P(Z > Z\_k) = \frac{1}{t} \tag{3}$$

and the probability of not exceeding it is:

$$P(Z \le Z\_k) = 1 - \frac{1}{t} = F(Z\_k) \tag{4}$$

In order to find *Zk* one has to use the annual (winter) maxima of the position of the zero centigrade temperature and the appropriate probability distribution. The Gumbel distribution [22] has been preferred here. The Gumbel extreme—value distribution was used to analyze the data on locations of the zero isotherm. The purpose of this study was to fit an extreme-value distribution to a given set of winter maxima of soil freezing depth denoted by *Zi*. The cumulative distribution function (CDF) is given as:

$$F(Z) = \exp\{-\exp\left[-\mathfrak{a}(Z - lI)\right]\}\tag{5}$$

where: *α* and *U* are parameters of the distribution estimated using the maximum likelihood method.

This distribution and procedure have been widely used to forecast wind velocity and snow load, e.g., [23]. It was also used to forecast soil freezing [24]. After being doubly logarithmized, it takes the form of a linear function with respect to *α*, *U*. The value *Zk* of the depth of zero centigrade temperature and 50-year return period is calculated using the following equation:

$$Z = \mathcal{U} - \frac{1}{\mathfrak{a}} \ln[-\ln F(Z)] \tag{6}$$

The unknown Gumbel distribution parameters (*α*, *U)* were estimated using the least squares method, the maximum likelihood method, the method of moments, and the method described by Julius Lieblein in his work [25].

#### *2.3. Estimation of Gumbel Distribution Parameters*

In order to determine the parameters, the measurement results are depicted on the Gumbel probability grid depending on the position each observation occupies in the statistics order. The measurement results are treated as the values of the function which arguments have the form, where the ordinal number is in the ascending order of the n-element measurement data set [26].

#### 2.3.1. Least Square Estimation

This method is relatively easy to use as double logarithm of Gumbel distribution results in the equation of a straight line. The parameters of the regression function are approximations of the Gumbel distribution parameters. Good approximation is obtained when noticeable linear relationship is present. The estimation method of Gumbel distribution parameters based on the least square method gives the best results as the coefficient of correlation R ≥ 0.95. However, it is very sensitive to the cases when some points are located far outside the straight line. Such cases are often occurring, for example, in snow load on the ground. Such points are increasing the steepness of the regression line, which leads to the overestimation of the predicted values. In extreme cases, the linear regression do not fit the approximated empirical distribution at all, which then disqualifies this method.

The straight line was fitted to the points placed on the probability grid so that it runs as close as possible to all points, i.e., the sum of the squares of vertical distance of points from this straight line is as small as possible. The inverse of the slope coefficient of this line and its Z-intercept are the parameters of the Gumbel distribution—being the results of the least squares estimation (LSM).

#### 2.3.2. Maximum Likelihood Method

Maximum likelihood estimation is recommended due to its many advantages. In this method, the assessment of the estimated parameters is based on values maximizing the credibility of the sample. Due to the lack of analytical solution for assessment of the distribution parameters, numerical methods are used. Estimators obtained with this method are consistent, asymptotically most efficient and they have normal asymptotical distribution. In practice, this method is rarely used due to calculation complexity, as the values of parameter distribution are assessed iteratively. The implementation of the maximum likelihood method (Max LM) consisted of seeking the maximum of the likelihood function given by the formula:

$$l = \prod\_{i=1}^{n} a \cdot \exp\{-a(Z\_i - lI) - \exp[-a(Z\_i - lI)]\}\tag{7}$$

The values *α*, *U* for which the Function (7) reaches its global maximum are the parameters of the Gumbel distribution estimated by the maximum likelihood method.

#### 2.3.3. Method of Moments

The moment method is simpler. In order to determine the parameters (*α*, *U)* using the method of moments, certain sample moments referring to measured values *Zi* were compared to the corresponding moments of the Gumbel distribution [27,28]. The arithmetic mean *Zm* and variance *σ*<sup>2</sup> as properties of the Gumbel distribution are expressed by the following formulas:

$$Z\_{m} = \mathcal{U} + \frac{1}{a} \cdot \gamma\_{\prime} \cdot \sigma^{2} = \frac{\pi^{2}}{6 \cdot a^{2}} \tag{8}$$

where: *γ* = 0.5772 is the Euler constant.

Comparing the Formulas (8) to corresponding quantities related to the n-element sample of winter maxima of freezing depth, we obtain the values *α*, *U* of the Gumbel distribution parameters estimated by the method of moments (MM).

#### 2.3.4. Lieblein Method—BLUE

Julius Lieblein [25] proposed a method of estimating Gumbel distribution parameters based on order statistics, requiring the sample elements *Zi* to be placed in ascending order. The estimators of the investigated distribution parameters take the form of a linear function with respect to positional statistics:

$$
\hat{\mathfrak{d}}\_{(n)} = \frac{1}{\hat{\mathfrak{d}}\_{(n)}} = \sum\_{i=1}^{n} b\_i Z\_{i\prime} \text{ \"if } = \sum\_{i=1}^{n} a\_i Z\_i \tag{9}
$$

where: Best Linear Unbiased Estimators (BLUE) are with minimal variance. The coefficients *bi*, *ai* were collated in the paper [25] for sample size of *n* = 1, 2, ... , 16. The aforementioned work also describes a method of estimating parameters *b*, *U* for arbitrarily large sample with a size of *n* having given coefficients for an *m*-element sample if *m<n*. For an *n*-element sample, the coefficients *b i* , *a <sup>i</sup>* can be calculated using the following formulas:

$$\begin{aligned} b'\_i &= \sum\_{j=1}^m b\_j \cdot (j/i) p(n\_\prime m\_\prime i\_\prime j)\_\prime \\ &\quad i = 1, 2, \dots, n \\ a'\_i &= \sum\_{j=1}^m a\_j \cdot (j/i) p(n\_\prime m\_\prime i\_\prime j)\_\prime \end{aligned} \tag{10}$$

where:

$$p(n,m,i,j) = \binom{i}{j} \binom{n-i}{m-j} / \binom{n}{m} \tag{11}$$

is the hypergeometric probability function. The formulas below allow the calculation of "good" estimators ˆ *b* (*n*) , *U*ˆ (*n*) for an *n*-element sample:

$$
\hat{\sigma}'\_{(n)} = \frac{1}{\hat{a}'\_{(n)}} = \sum\_{i=1}^{n} b'\_i \cdot Z\_{i\prime} \,\hat{\mathcal{U}}'\_{(n)} = \sum\_{i=1}^{n} a'\_i \cdot Z\_i \tag{12}
$$

The samples analyzed in this work contain significantly more than 16 measurement results; therefore, a set of BLUE coefficients *bj*, *aj* for *m =* 16 was used in the calculations.

The evaluation of the results obtained with four methods proved to be more difficult than the estimation of the Gumbel distribution parameters. As we have four Gumbel distributions, we face the question which of them best approximates the empirical distribution and can be used to determine the characteristic value *Zk* of the location of the zero isotherm in the ground [26]. The test statistics of chi-squared, Kolmogorov–Smirnov and Cramer von Mises were adopted as a measure of adjustment of the theoretical to the empirical distribution. Low values of these statistics indicate that the given theoretical distribution truly reflects the empirical distribution.

Amid the known compliance tests, the chi squared test is mentioned first. It requires the decomposition of the *Z* axis into *r* + 1 intervals *I*1, *I*2, ··· , *Ir*+<sup>1</sup> by the numbers −∞ = *g*<sup>0</sup> < *g*<sup>1</sup> < ··· < *gr* < *gr*+<sup>1</sup> = ∞. The difference in the value of the cumulative distribution function *F* at the ends of the interval is equal to the probability *p* that the random variable *Z* will take a value belonging to that interval:

$$p\_j = F(\mathcal{g}\_j) - F(\mathcal{g}\_{j-1}), \; j = 1, 2, \dots, r+1 \tag{13}$$

The chi squared test statistic for a measurement sequence consisting of *n* maximum freezing depths is given by the formula:

$$\chi^2 = \sum\_{j=1}^{r+1} \frac{\left(n\_j - np\_j\right)^2}{np\_j} \tag{14}$$

where: *nj*—the number of measured values of *Z* within *Ij*, while *npj*—is the expected number of measurement results which, according to the considered distribution, should be in *Ij* [29]. An examination of distribution compliance by means of statistics (14) leaves the choice of the number of intervals *Ij* into which the axis *Z* is divided and the points *gj* constituting the limits of these intervals. In this work because of the ambiguity mentioned chi-squared test was implemented in two variants. The first one assumes the division of the *Z* variable axis into one-element intervals, the limits of which are arithmetic means of two values adjacent in the order statistics. The only exception here are points *g*0, *gr+*<sup>1</sup> equal to minus and plus infinity, respectively. In the second variant, the division into five intervals with equal probability *pj* was used. Point *g1* was set as the arithmetic mean of the fifth and sixth elements in a series of ascending order data: *g*<sup>1</sup> = (*Z*<sup>5</sup> + *Z*6)/2. The value of the cumulative distribution function at point *g1* is equal to the probability of the first interval: *F*(*g*1) *= p*. The same probability of all intervals is ensured by the next boundary points *gj* being defined as:

$$\mathcal{G}\_{\vec{j}} = \mathcal{U} - \frac{1}{\mathfrak{a}} \ln \left[ -\ln \left( p + F(\mathcal{g}\_{j-1}) \right) \right], \dots, j = 1, 2, \dots, r \tag{15}$$

except the last: *gr*+<sup>1</sup> = ∞.

Other tests refer to hypothetical cumulative distribution function and verify its compliance with empirical cumulative distribution function. Kolmogorov–Smirnov statistics *Dn*<sup>1</sup> verifies the quality of the adjustment on the basis of the largest absolute value of the hypothetical *F*(*Z*) and empirical *Fn*(*Z*) cumulative distribution function difference:

$$D\_{n1} = \sup\_{z} |F\_n(Z) - F(Z)|\tag{16}$$

In Formula (16), unlike in [30], the <sup>√</sup>*<sup>n</sup>* factor is omitted, because only equally numerous sets of values *Fn*(*Z*) and *F*(*Z*) are compared with each other and without any reference to samples of a different size *n*. Another variation of Kolmogorov–Smirnov statistics found in the literature was also used, namely:

$$D\_n^+ = \max\_{1 \le i \le n} \left| \frac{i}{n} - F(Z\_{i:n}) \right|, \\ D\_n^- = \max\_{1 \le i \le n} \left| F(Z\_{i:n}) - \frac{i-1}{n} \right| \tag{17}$$

from where:

$$D\_{n2} = \max\{D\_{n}^{+}, D\_{n}^{-}\} \tag{18}$$

The above Formulas (16)–(18) define the maximum discrepancy between the postulated Gumbel distribution and the measured depths of the zero isotherm location. So, the idea arose to do a similar check for medium differences of hypothetical and geometrical cumulative distribution function. The subsequent tests are equivalent to statistics (16)–(18) and check not the maximum, but the average distance of the hypothetical and empirical cumulative distribution function. According to this reasoning, statistics (16) representing the average spacing of the cumulative distribution function is given by the formula:

$$D\_{n1s} = \sum\_{i=1}^{n} \frac{|F\_n(Z\_i) - F(Z\_i)|}{n} \tag{19}$$

Whereas the equivalent of statistics (18) in relation to the average value is expressed by the formulas:

$$D\_n^+ = \frac{1}{n} \sum\_{i=1}^n \left| \frac{i}{n} - F(Z\_{i:n}) \right|, \\ D\_n^- = \frac{1}{n} \sum\_{i=1}^n \left| F(Z\_{i:n}) - \frac{i-1}{n} \right| \tag{20}$$

and next:

$$D\_{\rm n2s} = \frac{1}{2} \left( D\_{\rm n}^{+} + D\_{\rm n}^{-} \right) \tag{21}$$

Another test used in this work, comparing the cumulative distribution function with order statistics of measurement results, is Cramer von Mises's statistics [31] in the form:

$$\mathcal{W}\_{\rm CM} = \frac{1}{12n} + \sum\_{i=1}^{n} \left( F(Z\_i) - \frac{2i - 1}{2n} \right)^2 \tag{22}$$

The best estimation method was considered to be the one indicated by most of the mentioned tests. For example, if the chi squared statistic is the smallest for the Gumbel distribution parameters *α*, *U* determined by the maximum likelihood method, it means that the chi squared test indicates the maximum likelihood method as the best of four applied methods.

#### **3. Results**

#### *3.1. Estimation Models*

From all recorded values of the position of zero centigrade temperature from measurements during every winter, one can extract the maximum depth and use in probabilistic calculations. Using the above mentioned Gumbel distribution, its parameters were estimated and the values of the predicted position of the zero isotherm in the ground were calculated (Figures 3–6). Calculations were made for data from years 1980–2011.

**Figure 3.** Estimation models for results from Białystok (**left**) and Elbl ˛ag (**right**) station—predicted position of zero isotherm for return period.

**Figure 4.** Estimation models for results from Jelenia Góra (**left**) and Kalisz (**right**) station—predicted position of zero isotherm for return period.

**Figure 5.** Estimation models for results from K ˛etrzyn (**left**) and Kłodzko (**right**) station—predicted position of zero isotherm for return period.

**Figure 6.** Estimation models for results from Łeba (**left**) and Warszawa-Bielany (**right**) station—predicted position of zero isotherm for return period.

This results in the location of the zero isotherm for the ground conditions of the 36 weather stations (first approach) where the measurements were made—Table 1.

#### *3.2. Validation of Chosen Model of Estimation (Lieblein)*

It was assumed that, as in the case of climatic interactions, the characteristic values of the zero isotherm position should have a return period of 50 years—as required reference level to ensure the reliability in design.

For the next step of calculation, only one model of Gumbel's method of analysis was used, with Lieblein's best linear unbiased estimators (BLUE). These estimators are recommended in European Standard for climate actions, e.g., derivation of wind speeds from measurements at meteorological stations. To assess the correctness of the predictions indicated in the first approach, a recalculation was performed for the enlarged dataset. The new data related to freezing depth measurements, obtained from 2011 to 2020 (black points on graphs), were added. On this basis, the previous forecasts have been compared with the values based on current observations (Figures 7–10). A new trendline has been set for the new data (1980–2020). With regard to the previous prediction, the change in the empirically determined freezing depth for the assumed return period (50 years) can thus be assessed.


**Table 1.** Results of calculation for all localizations.

**Figure 7.** Comparison of predicted and measurement value for chosen location: Białystok (**left**) and Elbl ˛ag (**right**) station.

**Figure 8.** Comparison of predicted and measurement value for chosen location: Jelenia Góra (**left**) and Kalisz (**right**) station.

**Figure 9.** Comparison of predicted and measurement value for chosen location: K ˛etrzyn (**left**) and Kłodzko (**right**) station.

**Figure 10.** Comparison of predicted and measurement value for chosen location: Łeba (**left**) and Warszawa-Bielany (**right**) station.

#### **4. Discussion and Conclusions**

The use of the described method allows for a safe and reliable prediction of the depth of freezing of the ground. Real measurements over many years and an analysis of trends show that these values are not constant and subject to change. This is mainly due to the increasingly intense climate change associated with global warming. For the chosen locations, new predictions show that the empirical depth of freezing is general descending. This is one symptom of global warming.

The presented method of analyzing the depth of soil freezing takes into account actual soil temperature as a function of all external climatic actions with their coincidence and duration: air temperature, snow cover and precipitations. For this reason, it is more precise than the method using the so called air freezing index [16,26].

A comparison of calculation results related to the location of zero isotherm, conducted based on ground temperature analysis and freezing index (according to Gontaszewska [1]) leads to a conclusion that the second method has a lower accuracy. It does not account for snow cover, leading to much deeper freezing depths. Only in two cases (Leszno, Koszalin), almost the same results (difference of <2%) were obtained. The largest differences occurred in the mountains, where snowfall is the biggest (Jelenia Góra and Zakopane; values based on AFI are larger by 86%) and in North-East Poland (Białystok, 48%). Significant differences were also obtained in Warsaw (difference of 49%) and Kraków (77%).

Moreover, the probabilistic approach permits us to analyze the soil freezing in the terms of the reliability of structures. As the influence of climate can be easily analyzed using probabilistic approach, the influence of soil structure on the depth of freezing is much more complicated. It is due to the heterogeneity of soil and its structure. The depth of freezing is influenced by the type of the soil, its state of saturation, porosity, and mineral composition, and the arrangement of layers (Xu et al. [32], Zurawski & Godlewski [ ˙ 33], Yu, et al. [30]). This problem can be solved now only on the basis of the available literature using correction coefficients. One of the sources is available to authors now: former Soviet correction coefficients to the Formula (1) (Ickiewicz, Pogorzelski [34]): 23 for clays and silts, 28 for clay-sands, dusts and fine-grained sands, 30 for coarse-grained sands, and 34 for gravels. On the base of these values, it is possible to assume that gravel is a reference soil and the ratio of above mentioned values to 34 will be correction coefficients. In order to reduce the values based on different types of soil the inverse values of these coefficients will be used. An example is given below.

At Suwałki meteorological station the soil can be classified as silts and sands. In this case the actual depth of freezing reduced to the gravels as a standard reference soil can be received as a product of the calculated value 1.05 m for the 50-year return period and a factor of 34/28 = 1.21. If the soil in Suwałki region is assumed as a reference gravel, then the freezing depths will be 1.05 × 1.21 = 1.27 m. On the contrary, if we want to calculate the depth for clays and foams, we will receive 1.27 × (23/34) = 0.86 m.

Also, the lack of snow cover under conditions of reduced temperature can significantly increase the freezing range. Snow cover greatly inhibits freezing depths in soils, such that freezing depths can vary from year-to-year or spatially within a given time period, solely dependent on the presence or absence of a snow cover. The implications of a future climate, without or an increasingly intermittent snow cover, are greater freezing depths or varying depths, respectively. The currently observed snowless winters (in Poland and central part of Europe), with low rainfall (hydrological drought) and possible sudden waves of frost can lead to an increase in the range of freezing of the ground, which is already observed in the results of measurements and presented predictions.

At present, new analysis is currently being done in order to update the old provisions presented in Polish Standards. The influence of the soil structure will be taken into account using correction coefficients, as listed above. It should be noted here that the new prediction ranges and the depths of the zero isotherm are to be considered as design values representative for reliability analysis for the return period of 50 years.

Further research on the influence of the soil structure on its freezing is necessary. Such investigations should be carried out for different types of soil. They should also examine the difference between the depth of the position of zero centigrade temperature and actual frost penetration into the soil. Meanwhile, the correction coefficients will be used as presented above.

**Author Contributions:** Conceptualization, T.G. and Ł.W.; methodology, T.G., Ł.W. and M.W.-N.; software, T.G. and Ł.W.; validation, Ł.W., and M.W.-N.; formal analysis, T.G. and M.W.-N.; investigation, Ł.W. and M.W.-N.; resources, T.G. and M.W.-N.; data curation, Ł.W.; writing—original draft

preparation, Ł.W. and M.W.-N.; writing—review and editing, T.G.; visualization, Ł.W. and M.W.-N.; supervision, T.G.; project administration, T.G.; funding acquisition, T.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We appreciate the soil temperature data provided by the Institute for Meteorology and Water Management—State Research Institute.

**Conflicts of Interest:** The authors declare no conflict of interest. The authors do not have any competing financial, professional, or personal interests from other parties.

#### **References**

