THM-coupled numerical analysis of temperature and groundwater level in-situ measurements in artificial ground freezing

of for water-saturated equations of the model on balance laws energy and momentum for a fully saturated porous media. Clausius-Clayperon equation and poroelastic constitutive relations to coupled processes in water and ice pore pressure, porosity and a stress-strain state of freezing soil. The proposed model used to predict equivalent water content measured in Mizoguchi’s test and frost heave in a one-sided freezing test. Numerical simulation of ground freezing in the Petrikov mining complex located in that the model is able field measurements of inside forming the mismatch between hydro- and thermo-monitoring data obtained during the artificial freezing is analyzed. et al. [39] and Lai et al. [40], the mechanical behavior of the freezing soil is described by a change in porosity which is defined by the mass balance equation. Soil porosity, equivalent pore pressure, and a stress-strain state of the freezing soil are evaluated using constitutive relations of poroelasticity proposed by Coussy [35-36] and effective stress conception developed by Bishop [41]. Clausius – Clapeyron equation is used for estimation of pore ice pressure and cryogenic suction. Phase transition of water into ice is incorporated in the model by a soil freezing characteristics curve. Coupled set of nonlinear equations of the model were implemented in Comsol Multiphysics software. The effectiveness of the proposed model was demonstrated by numerical simulation of two laboratory tests. The first one is a well-known Mizoguchi’s test in which an evolution of equivalent water content in sandy loam samples was measured during freezing in a closed system. The second test is one-sided freezing of silty sand samples in an open system to measure frost heave displacement. The validated model was applied to numerical simulation of pore pressure evolution in unfrozen soil inside a closed cylindrical frozen wall. Numerical predictions were compared to field measurements of groundwater level recorded by two hydro-observation wells of different depths during AGF in the Petrikov mining complex (Republic of Belarus). A mismatch between field measurements of temperature and groundwater level collected in different layers is analyzed. In addition, an effect of water migration to the inner boundary of the frozen wall on the pore pressure inside the frozen wall is discussed. Comparison of the numerical results and experimental data has shown that the proposed model is able qualitatively to describe both a porosity reduction due to water migration and thermal shrinkage of the solid skeleton as well as porosity rise in the frozen zone induced by the frost heave. The second and third stages show better agreement between calculated and experimental results than the first one. In the first stage, maximum deviation between the plots is observed for 24 hours of the freezing and is equal to 0.04. As regards the second stage, the maximum deviation is 0.016 and also corresponds to the freezing for 24 hours. The obtained deviation is less than 10 %, which indicates an acceptable agreement between numerical results and experimental measurements. laboratory test


INTRODUCTION
verlying stratums of potash salt deposits in Belarus consist of aquifers and unstable soils. Therefore, artificial ground freezing (AGF) technology is required to provide groundwater control and enhance the mechanical properties of such soils. AGF is performed before and during the shaft sinking by a series of freezing wells drilled around the projected contour of the shaft. A circulation of refrigerant in the pipes placed into the wells induces heat extraction from the soils and the formation of a frozen wall. For a safe shaft sink operation, it is required to achieve a solid impermeable frozen wall of the designed thickness. Nowadays control of a frozen wall state is carried out by field measurements of ground temperature and groundwater level. However, analysis of the measurements is frequently complicated by a lack of data on initial hydrogeological conditions in the excavation site and freezing conditions. A promising way to study the freezing process and analyze the monitoring data is through mathematical modeling of AGF. Research of soil freezing was begun about one hundred years ago. The pioneering works of [1][2][3] were related to studies of unfrozen water content in soils at negative temperatures. Taber [4] and Beskow [5] observed that the frost heave in freezing soil was induced by ice lens formation associated with water migration from the unfrozen zone to the freezing front. Further efforts were intended to detailed investigation and a quantitative description of such processes as cryogenic suction in soils at negative temperatures [6][7], the dependence of unfrozen water content on temperature [8][9], frost heave of freezing soils [10][11][12] and mechanical behavior of frozen soils [13][14]. As a result, coupled model of heat-fluid transport in freezing soil [15][16], the ice rigid model [17], and a conception of the segregation potential [18][19] were developed. These models provide a theoretical foundation for evaluation of temperature, water velocity and frost heave strain due to ice segregation in freezing soils. It is assumed that water and ice coexist in phase equilibrium. The unfrozen water content is expressed by soil-freezing characteristic curves. The studies related to the soil-freezing characteristic curves can be found in many works [20][21][22][23]. Modern thermo-hydraulic models are able to predict temperature evolution and distribution of water, and ice content in laboratory tests [24][25][26] and in engineering applications of the AGF [27][28]. However, analysis of frost heave effect on surrounding unfrozen soil and structures requires stress and strain fields evaluation. Currently, fully coupled thermo-hydromechanical models for geotechnical application are derived using the theory of porous media. A model proposed by Nishimura et al. [29] was applied to a study of a frost heave effect on buried pipelines. In the model, the mechanical behaviour of freezing soil is described by extended Barcelona Basic Model with two-stress variable constitutive relations. In Liu and Yu [30] performance of a thermo-hydro-mechanical model was tested on field data collected in pavement of seasonally frozen soils. A volumetric expansion of freezing soil is incorporated through additional inelastic strains. Fully coupled thermo-hydro-mechanical models of AGF are presented in [31][32][33][34]. According to these models, the stress-strain state of the freezing ground is simulated using constitutive relations of poroelasticity proposed by Coussy [35][36]. In the studies of Zhou and Meschke [31] and Tounsi et al. [34] numerical simulation of AGF is performed for horizontal excavations. Panteleev et al. [32][33] consider AGF for a vertical shaft sinking in the Petrikov potash deposit. However, in the model, the mass balance law is considered only for unfrozen soil. This work is a continuation of [37][38] where the pure mechanical behavior of ice wall at the Petrikov mining complex (Republic of Belarus) has been investigated. In this paper, a thermo-hydro-mechanical model of freezing for saturated soil is proposed. The model was applied to the analysis of field measurements of ground temperature and groundwater level obtained by hydro-observation wells. The governing equations of the model include the mass balance equation, the energy conservation equation, and the momentum balance equation. According to thermo-hydro-mechanical models of frost heave developed by Zhou et al. [39] and Lai et al. [40], the mechanical behavior of the freezing soil is described by a change in porosity which is defined by the mass balance equation. Soil porosity, equivalent pore pressure, and a stress-strain state of the freezing soil are evaluated using constitutive relations of poroelasticity proposed by Coussy [35][36] and effective stress conception developed by Bishop [41]. Clausius -Clapeyron equation is used for estimation of pore ice pressure and cryogenic suction. Phase transition of water into ice is incorporated in the model by a soil freezing characteristics curve. Coupled set of nonlinear equations of the model were implemented in Comsol Multiphysics software. The effectiveness of the proposed model was demonstrated by numerical simulation of two laboratory tests. The first one is a well-known Mizoguchi's test in which an evolution of equivalent water content in sandy loam samples was measured during freezing in a closed system. The second test is one-sided freezing of silty sand samples in an open system to measure frost heave displacement. The validated model was applied to numerical simulation of pore pressure evolution in unfrozen soil inside a closed cylindrical frozen wall. Numerical predictions were compared to field measurements of groundwater level recorded by two hydro-observation wells of different depths during AGF in the Petrikov mining complex (Republic of Belarus). A mismatch between field measurements of temperature and groundwater level collected in different layers is analyzed. In addition, an effect of water migration to the inner boundary of the frozen wall on the pore pressure inside the frozen wall is discussed.

THERMO-HYDRO-MECHANICAL MODEL OF SOIL FREEZING
reezing soil is modelled as a fully saturated three-phase porous medium consisting of solid grains (index s), liquid water (index l) and ice (index i). In the initial state, the pore space is filled only with water. The air phase is ignored. According to [31], [34] the following hypotheses are applied: 1. All phases of the porous media are in the local thermal equilibrium, so the temperature of all phases is the same. 2. Phase densities are assumed to be constant during freezing. 3. Effects of pore water salinity and external load on the freezing temperature are not considered. 4. Ice moves together with solid grains, so a velocity of ice relative to solid grains is zero. 5. Ice lens formation is negligible. 6. The soil is elastic and isotropic. Deformation of the soil is estimated using the small strain formalism. Balance equations for a three-phase porous media can be formulated on the base of these hypotheses. The mass balance equation is written as where jSjn is the mass content per unit of volume of water (j = l) and ice (j = i) at time t; j is the density of the phase j; Sj is saturation of the phase j; n is the porosity; v l is the velocity of water relative to the solid skeleton. The energy conservation law has the following form where T is the temperature of the porous media, C and λ are the volumetric heat capacity and the thermal conductivity of the porous media; С l is the volumetric heat capacity of water; Q ph is the heat source related to the latent heat due the water phase change. The momentum balance equation for the porous media is given as: σ is the total stress tensor; γ is the unit weight of the porous media. The ice saturation S i is determined by a soil freezing characteristic curve which is expressed as a power function of the temperature T [42]: where Tph is the freezing temperature of pore water and α is the experimental parameter. The water saturation is found from the condition of the full saturation Sl = 1 -Si.
The relative velocity of water v l can be described by the Darcy law as where k is the hydraulic conductivity and ψ is the soil-water potential. Following Nixon [43], the hydraulic conductivity k is given as F where k 0 is the hydraulic conductivity of the unfrozen soil, β is the experimental parameter. The soil water potential ψ driving pore water migration is defined according to Bernoulli's equation: where p l is the pressure of the pore water; g is the gravitational acceleration; z is the vertical coordinate. The volumetric heat sources in the energy conservation Eqn.
(2) are written as where L is the latent heat. The volumetric heat capacity C and the thermal conductivity λ are determined as [44]    where cj, λj (j = s,l,i) are the specific heat capacities and the thermal conductivities of the phase j. Unit weight of the saturated porous media γ is written as The total stress σ of the porous media is defined as where σ′ is the effective stress; p is the equivalent pore pressure; b is the effective Biot coefficient; I is the identity tensor. The effective stress σ′ is given by the Hooke's law for isotropic linear-elastic media where K is the effective bulk modulus; G is the effective shear modulus;  e is the elastic strain tensor;  e vol is the volumetric part of the tensor. According to the additive decomposition of the total strain , the elastic strain  e can be expressed as where  th is the thermal strain. Total strain  is defined by the geometric relation for an infinitesimal deformation: where u is the displacement vector. Thermal strain is written as where a T is the volumetric thermal dilation coefficient and T 0 is the initial temperature of the unfrozen soil.
The equivalent pore pressure p is assumed to be the weighted sum of the pore water pressure p l and the pore ice pressure p i [17], [41]: where χ is the pore pressure parameter defined as The pore ice pressure pi is expressed by the Clausius-Clapeyron equation as follows where p hydr is the reference pressure in soil before the freezing process. Following [39][40], expression for the pore water pressure pl can be obtained from (18), (20) as According to theory of poroelasticity [31], [35], [36] the equivalent pore pressure p can be expressed by the porosity and the volumetric strain where n 0 is the initial porosity; N is the effective Biot tangent modulus. The effective mechanical properties are estimated as [31]   where X is the effective value, X fr and X un are the values for the frozen and the unfrozen states.

COMPUTER IMPLEMENTATION OF THE MODEL
he Eqns. (1)- (22) were solved in finite-element software Comsol Multiphysics® using Weak Form PDE Interface, Heat Transfer Module and Solid Mechanics Module. The porosity n, the temperature T and the displacement u were chosen as primary variables. In [39][40] have been shown that solution of the mass balance equation relative to the porosity n enables accurate prediction of water and ice distributions. To incorporate the mass balance Eqn. (1) in Comsol a weak form of the equation was obtained. By taking partial derivative with respect to time, using condition of fully saturation of the porous media and the Darcy's law (5) we can rewrite the Eqn.
(1) as Eqn. (23) was multiplied by a test function F and integrated over the domain Ω using integration by parts to obtain the weak form. The final form of the equation is where m=(m x ;m y ) is the outward unit normal to the boundary Γ of the domain Ω.
Eqn. (24) was solved in Comsol using Weak Form PDE Interface. The energy conservation Eqn. (2) was solved using Heat Transfer in Solids Interface where the heat source Qph and the convective term cl vl ·gradT were included by Heat Source nodes. The momentum balance Eqn. (3) together with constitutive relations (13)- (16) were solved by Solid Mechanics Interface. Influence of the pore pressure p on a stress-strain state was taking into account by External Stress node. Thermal strain  th was added to the total strain  by Thermal Expansion node. For the spatial discretization of the porosity Eqn. (24), quadratic Lagrange shape functions were used whereas linear Lagrange shape functions were applied for the heat transfer equation. Displacement vector u was approximated using quadratic serendipity shape functions. Optimal element sizes were achieved by successive refinement of the mesh in each calculation. Time-dependent problem (1)- (22) was integrated according to the Backward Differentiation Formula. Segregated solver was used for the computation, where the temperature field T was defined in the first step while the porosity n and displacement u were obtained in the second step. System of non-linear algebraic equations in every step was solved by a damped Newton method with a constant damping factor. Direct PARDISO solver was employed for the linearized system of equations.

VALIDATION OF THE PROPOSED MODEL
he validation of the proposed model has been carried out by simulation of two different laboratory tests. The first is the benchmark test carried out by Mizoguchi [24][25], [30]. This test gives us information about water and ice content during one-sided freezing of four identical cylindrical samples packed with sandy loam in a closed system. One cylinder was taken as a reference, the rest were frozen for 12, 24, and 50 hours. The test was conducted in a closed system under the following conditions. The cylindrical sample had a height of 20 cm and a radius of 4 cm. In the initial time moment, the soil was fully saturated with water. Volumetric water content which coincided in this case with the porosity was uniformly distributed along the sample height and was equal to 0.35. Initial temperature of the cylinder was 6.7 O C. The top surface of the cylinder was subjected to a constant temperature of -6 O C. The other surfaces were thermally insulated. Due to the radial symmetry of the problem, half of the cylinder's cross-section was simulated. Therefore, a computational domain had a rectangular shape with a width of 4 cm and a height of 20 cm. According to Zhou et al. (2012) it is assumed that the top surface of the soil was subjected to instantaneous freezing without water migration. The difference in the densities of the water and ice induces increase in the volume of the pore water by 9%. Therefore, the value of the porosity at the top boundary of the sample is assumed to be 1.09·n 0 =0.3815, where n 0 =0.35. Zero flux boundary condition was applied at the other sides of the domain. The displacement vector at the bottom boundary was constrained in all directions. Symmetry boundary condition was given at the axis of the symmetry and only vertical displacement was permitted at the lateral boundary. The parameters of the soil used in the simulation are listed in Tab. 1. The considered area was divided into quadrilateral elements. The optimal mesh size was determined by series of calculations with various element sizes. Reference numerical solution was obtained on a computational mesh with 30000 elements which ensures sufficiently fine partition of the computational domain. The relative tolerance tol has been determined by a formula: where subscript r denotes the reference solution, L represents the middle line of the computational domain along which the porosity n has been defined.    Fig. 2 presents porosity distribution along the sample height measured in the test and obtained by numerical simulation for 12, 24, and 50 hours of the freezing. The presented plots can be divided into three stages. In the first stage, the porosity rises and the temperature of the soil is less than the freezing temperature Tph. Therefore, at this stage frost heave of the soil is occurred due to the phase transition of water into ice and water inflow. In the second stage, the temperature is also less than T ph , but the porosity significantly reduces to a minimum value. It indicates that there is an intensive water migration into the freezing front induced by cryogenic suction. The water migration together with the mechanical impact of the frozen zone leads to a volumetric shrinkage of the soil. Since a temperature gradient in the sample decreases in the process of freezing, the length of this stage increases. The third stage corresponds to unfrozen soil. At this stage the porosity is also reduced by water migration to the freezing front. Comparison of the numerical results and experimental data has shown that the proposed model is able qualitatively to describe both a porosity reduction due to water migration and thermal shrinkage of the solid skeleton as well as porosity rise in the frozen zone induced by the frost heave. The second and third stages show better agreement between calculated and experimental results than the first one. In the first stage, maximum deviation between the plots is observed for 24 hours of the freezing and is equal to 0.04. As regards the second stage, the maximum deviation is 0.016 and also corresponds to the freezing for 24 hours. The obtained deviation is less than 10 %, which indicates an acceptable agreement between numerical results and experimental measurements.    Fig. 3 presents the distribution of the vertical displacement along the sample after 2 hours of the freezing (Fig. 3(a)) and variation of the vertical displacement of the top surface with time ( Fig. 3(b)). It can be seen that shrinkage of the soil arises near the freezing front. On the other hand, a frost heave induces a vertical uplift of the soil in the frozen zone. The mechanical behavior indicates a rise in the equivalent pore pressure of the frozen zone and a decrease in the pressure of the unfrozen zone. Similar distributions were obtained by [39]. Fig. 3 (b) shows that upward displacement induced by frost heave achieved 4.6 mm as it was measured at the end of the freezing test. Calculated evolution of the vertical displacement qualitatively close to experimental curves recorded in one-sided freezing experiments carried out by Lai et al. [40]. At the beginning of the freezing, a small frozen shrinkage of the soil occurs. After that, the frost heave strain increases due to the freezing of the initial and migrated water.

ARTIFICIAL GROUND FREEZING FOR SHAFT SINKING IN PETRIKOV MINING COMPLEX
he validated thermo-hydro-mechanical model of soil freezing was applied for the study of pore pressure inside a cylindrical frozen wall formed by AGF at the Petrikov mining complex.

T
The image of the freezing system is given in Fig. 4(a). The system consisted of vertical freezing wells drilled along the projected contour of the shaft. Freezing pipes were located circumferentially with a radius r of the freezing contour equal to 8.25 m. The radius of each well was 7.3 cm. The distance between the neighboring wells was 1.1 m. The depth of the freezing wells was 250 meters. Steel pipes were installed in the wells. The temperature of refrigerant pumped in the pipes by the cooling station was -27 ºС. A detailed description of artificial freezing at the Petrikov potash deposit is presented in [46]. The location of freezing and hydro-observation wells is shown in layout diagram given by Fig. 4(b). The wells were drilled inside the freezing contour at distances of 3.2 m and 3.5 m from the center. The wells are marked as HW1, HW2 and have depths of 82 m and 200 m respectively. Feeding zones of the wells were located in two different fine-grained soil layers of glauconite sand and slightly clayey quartz-feldspar sand. Fig. 4 (c) presents data on groundwater levels measured by HW1 and HW2 during 30 days from 03.04.2016 to 02.05.2016. An increase in a frozen wall thickness induces shrinkage of unfrozen soil inside the circle of the freezing wells due to the frost heave. When the frozen wall is not closed, pore water can flow freely to the unfrozen soil through gaps in the wall. If the frozen wall is solid (impermeable to the water flow), the amount of water in the unfrozen soil is constant. As a result, the unfrozen soil compression leads to an increase in the pore pressure and, as a consequence, in groundwater level in a hydro-observation well which can be seen in Fig. 4 (c). A difference in depths of hydro-observation wells enables the control of the frozen wall integrity along the vertical direction. In the Petrikov mining complex, HW2 was equipped with fiber-optic sensors. Fig. 4 (d) shows the ground temperature profile along the depth of the HW2 measured on the first day of the monitoring. Data given by Fig. 4 (c) shows that temperature reduces with the depth. Therefore, the frozen wall is expected firstly to close at the depth of 200 m rather than at the depth of 82 m. However, the groundwater level rises up to a maximum value after 17 days in HW1 and only after 22 days in HW2. Hence, it can be concluded that freezing at the depth of 82 m is more intensive than at the depth of 200 m. The resulting mismatch between the field measurements is analyzed in the following section.

NUMERICAL ANALYSIS OF PORE PRESSURE VARIATION INSIDE A CYLINDRICAL FROZEN WALL
n order to study an evolution of pore pressure in unfrozen soil inside a frozen wall and analyze possible reasons for the inconsistency between the temperature monitoring and groundwatel level data obtained by hydro-observation wells at the Petrikov mining complex, numerical simulation of freezing of saturated soil layers was carried out in Comsol Multiphysics® software. To perform the simulation some assumptions were made. The feeding zones of hydro-observation wells HW1 and HW2 are located at two different soil layers (glauconite sand and slightly clayey quartz-feldspar sand respectively). Results of three-dimensional numerical simulation of AGF at the Petrikov mining complex performed by Panteleev et al. [32][33] have shown that a frozen wall has a cylindrical form within separate soil stratums. Heat transfer along the vertical direction and water flow due to the gravity potential can be neglected. Also, as the distance between neighboring freezing well is small in comparison with the radius of the frozen wall, the frozen wall before the closure has a small influence on the unfrozen soil inside the circle of the freezing wells. Therefore, we may consider the frozen wall formation in a horizontal section of the soil stratums at the depths corresponding to the feeding zones of HW1 (the first simulation case) and HW2 (the second simulation case) and perform the simulation for 2D domain. The material parameters used for the simulation are listed in Tab. 3. Thermal, filtration, and elastic properties of the soils were obtained by the INM of NAS. In both simulation cases, the computational domain was a circle with a radius of 8.25 m. A reference day of the calculation corresponds to the day when the freezing system had already begun to operate but the closure of the frozen wall has not yet come. According to the data presented in Fig. 4 (d), the initial temperature corresponding to the depth of HW1 is 5.2 O C and the initial temperature corresponding to the depth of HW2 is 1 O C. The freezing temperature was applied at the boundary of the computational domain. Assuming a linear variation of the coolant temperature with a depth we obtain that the freezing temperature at the depth of 82 m is 3 O C higher than at the depth of 200 m and is equal to -24 O C. A constant value of the porosity was given at the boundary of the circle. In the first simulation case, it was equal to 1.09·0.66=0.72 and in the second simulation case, it was equal to 1.09·0.24=0.26. In both cases, the displacements were constrained at the boundary of the circle. Triangular elements were applied to discretize simulation domain. A mesh convergence study was performed by a series of calculations with varying numbers of finite elements. Reference numerical solution was obtained on a computational mesh with 25970 elements. The relative tolerance tol has been calculated with regards to the pore pressure as: where subscript r denotes the reference solution. It should be noted that value of pore pressure p in (26) Table 3: Parameters of soils and initial conditions in the feeding zones of the hydro-observation wells.
Results presented in Fig. 5 show that relative tolerance reaches plateau when total number of elements exceeds 6550 elements. This finite-element mesh corresponds to the maximum size of the element equal to 0.33 m and minimum size equal to 0.00124 m. Temperature and pore pressure distributions in the computational domain for the first and the second simulation case are presented in Fig. 6,7. The results are shown for the 17 th day and 22 nd day of the calculation which correspond to the closure of the frozen wall according to the data provided by hydro-observation wells (Fig. 4(c)). It can be seen that the thickness of the frozen wall in the second case is higher than in the first case, which is confirmed by the temperature measurements recorded by HW2 (Fig. 4(d)). However, hydro-observation wells predict the earlier rise in the groundwater level and, as a consequence, in pore pressure value for HW1 (Fig. 4 (c)).  Fig. 8 shows a variation of the pore pressure normalized by its maximum value at the points of the computational domain, which correspond to locations of HW1 and HW2. Experimental pressures corresponding to hydrostatic pressure were obtained from the data provided by Fig. 4 (c). The trend predicted by the model is similar to the experimental data. In the first simulation case, the maximum value of the pore pressure is reached by the 17th day of the freezing. In the second case, the maximum value of the pore pressure is attained by the 22 nd day of the freezing. This delay can be explained by a difference in material properties of the considered soils as well as by discrepancy in the freezing conditions. The effect of these parameters on pore pressure variation will be analyzed below. Disagreement between experimental and computational results for HW2 from 15th to 19th day can be explained by measurement errors, which lead to the oscillation of the pore pressure and its decline in this time while the model predicts constant rise. The model predicts a uniform formation of the frozen wall, so the pore pressure monotonically rises.  Fig. 4(c)).
In order to analyze possible reasons for the mismatch between the temperature and groundwater level measurements in HW2, the effect of some model parameters on the pore pressure at the point corresponding to the location of HW2 has been studied. Fig. 9 presents the effect of the coolant temperature on the pore pressure variation. It is known that the temperature of the refrigerant is distributed non-uniformly along the freezing well. The shallower is the depth, the higher is the freezing temperature. Results of simulation show that a significant difference in the pore pressure values occurs only for the high discrepancy between coolant temperatures. Calculation shows that difference in coolant temperatures at the depth of 82 m and 200 m is 3 O C. According to Fig. 9, this mismatch can provide only a slight difference in pore pressure variation. Figure 9: Effect of the coolant temperature on pore pressure variation. Fig. 10 shows the effect of  and  on the variation of pore pressure with time. Parameter  is responsible for the shape of the soil freezing characteristic curve defined by Eqn. (4). Parameter  determines the rate of decrease in hydraulic conductivity with temperature. It can be seen that an increase in the magnitude of  leads to a rise in pore pressure value ( Fig. 10 (a)). For higher values of  a fraction of water converted into ice rises, so frost heave in the frozen zone evolves more intensively. The opposite effect can be observed for  (Fig. 10 (b)). Relatively high values of magnitude provide lower values of the pore pressure and delay its abrupt rise. Low values of  provide more intensive water migration into the frozen zone due to cryogenic suction. As a result, the frost heave intensifies increasing compressive stress acting on the unfrozen soil. This leads to a faster increase in the pore pressure value compared to the higher values of . However,  in the first simulation case was lower than in the second case. Therefore, higher pore pressure in the first simulation case could be provided by . Except for , Eqn. (6) includes one more parameter k 0 . This parameter describes the hydraulic conductivity of the unfrozen soil. The effect of k0 on pore pressure variation is presented in Fig. 11. Higher values of k0 lead to the increase in pore pressure value compared to the lower ones. Therefore, the less permeable soils can provide a significant delay in pore pressure rise. The first reason for this delay is similar to the delay induced by high magnitudes of . Higher values of k0 lead to a rise of amount water migrated into the frozen zone. As a result, a more amount of water transforms into ice, so the compressive loadings acting on the unfrozen soil increases and, as a consequence, the pore pressure rises. The second reason is that in low permeable soils water flows from the compressive zone of the unfrozen soil to a feeding area of a hydroobservation well with lower velocity leading to a delay in pore pressure rise. From Figs.8, 11 we can conclude that even though k 0 in the second simulation case was higher it could not provide a significant increase in pore pressure. Fig. 11 (b) shows the effect of initial porosity on pore pressure evolution. We can observe that higher initial porosity induces higher pore pressure at the hydro-observation well. High-porous soils contain large amount of water, so during freezing, an intensive frost heave occurs. This situation is similar to the first simulation case.
(a) (b) Figure 11: Effect of k 0 (a) and n 0 (b) on pore pressure variation. Fig. 12 demonstrates influence of the Biot tangent modulus N un of the unfrozen soil on pore pressure at the hydroobservation well. According to Eqn. (21) the pore pressure is directly proportional to this parameter. Therefore, higher values of Biot tangent modulus of soil lead to higher pore pressure value as it was obtained in numerical simulation.
Presented results have shown that an abrupt rise in pore pressure inside the frozen wall depends on the effect of frost heave in the frozen zone on the unfrozen soil inside the wall. The intensity of the effect rises with the thickness of the frozen wall. Therefore, the pressure rise can occur with a perceptible delay relative to the closure of the frozen wall. Also, it should be noted, that water migration into the frozen zone contributes to the frost heave of the soil. However, significant water migration can also lead to pore pressure reduction. Simulation results show that in the presence of a strong cryogenic suction, which defines by Nun =200 MPa , k0 =1.5·10 -9 m/s and = -2.53 parameters we can obtain the situation when pore pressure decreases with an increase in thickness of the frozen wall as illustrated in Fig. 13.  Aside from the material properties, the pore pressure inside of a frozen wall can also depend on a radius r of the freezing wells location. Fig. 12 illustrates the effect of the freezing сontour radius on pore pressure variation. It can be seen that a decrease in the radius leads to an increase in the pore pressure inside the frozen wall. Therefore, freezing wells located at the larger radius provide fewer values of the pore pressure. Commonly, the freezing wells are drilled with a substantial vertical deviation, which can reduce or enlarge the radius of the freezing contour. Thus, drillhole deviation also can play a significant role in the pore pressure variation. Figure 13: Effect of the freezing contour radius r on pore pressure variation.
The conducted study has illustrated the strong dependence of pore pressure evolution on material properties of the soil as well as freezing conditions. On the one hand, the obtained results have shown that earlier rise in pore pressure level in HW1 can be induced by the higher value of initial porosity of glauconite sand, a larger magnitude of  parameter, or higher value of Biot tangent modulus of the unfrozen soil. At the same time, an abrupt rise in HW1 cannot be explained by  parameter and k0. On the other hand, vertical deviation of freezing wells can also explain the mismatch between temperature and groundwater level data while coolant temperature has a minor effect on this discrepancy.

CONCLUSION
GF for a vertical shaft sinking is a complex engineering activity that induces a series of large-scale thermal, hydraulic and mechanical processes inside the ground at depths up to several hundred meters. To study AGF we have developed a coupled thermo-hydro-mechanical model of saturated freezing soil, which takes into account phase transition of water into ice with the associated release of the latent heat, water migration into the frozen zone due to cryogenic suction as well as nonhomogeneous distribution of water and ice content. The model was verified by two different laboratory experiments. The results have shown that the model is able to describe porosity evolution related to frost heave and consolidation of the freezing soil as well as to predict upward displacement induced by frost heave. Application of the model to numerical simulation of AGF in the Petrikov mining complex enables us to describe pore pressure variation in hydro-observation wells located inside the frozen wall and also propose an explanation for the mismatch between measurements of the temperature and groundwater level in fine-grained sand stratum laid at depth of 82 m and 200 m. According to the obtained results, variation of the pore pressure at the hydro-observation wells is strongly dependent on the material properties of soils and freezing conditions. Analysis of the results has shown that delay in a pore pressure increase in the hydro-observation well HW2 (200 m) in comparison with HW1 (82 m) can be related to less intensive frost heave in the frozen zone due to lower values of the initial porosity,  parameter responsible for a rate of water freezing and Biot tangent modulus of the unfrozen soil. The effect of the coolant temperature is less significant. Moreover, a decrease in the radius of the freezing well contour induces an increase in the pore pressure. Therefore, drillhole deviation could also lead to a mismatch between the field measurements.