Numerical simulation of vertical shaft sinking using artificial ground freezing

. Artificial ground freezing (AGF) is used worldwide for vertical shaft sinking in difficult hydrogeological conditions. The modern tendency is to determine the design parameters of the freezing technique based on numerical simulation. This work is devoted to the numerical simulation of the formation of an ice-soil wall in the soil stratum due to the AGF and shaft sinking under the protection of the wall. For this purpose, a fully coupled thermo-hydro-mechanical model of soil freezing has been developed on the basis of the theory of poromechanics. The developed model considers important features of the freezing process, such as the phase change, pore water migration due to cryogenic suction, frost heave, and consolidation of the soil. The results have shown that the model allows to predict the distribution of ice content, assess stress and strain in the ice-soil wall, and estimate displacement of the excavation wall.


INTRODUCTION
In the modern mining industry, construction of vertical mine shafts for elaboration of potash deposits is frequently accompanied by complex geotechnical problems related to groundwater seepage and unstable soft soil stratums. One of the most effective techniques applied worldwide in geotechnical engineering to shaft sinking under hard hydrogeological conditions is artificial ground freezing (AGF) [1]. The aim of the method is the formation of a temporary ice-soil retaining structure around the intended excavation. Freezing of water-saturated soft soils leads to the reduction of its permeability and improves its strength and stiffness properties. As the result, the ice-soil wall serves as a temporary excavation support, which eliminates groundwater flow into the excavation and prevents squeezing of the excavation wall by the rock pressure.
Determination of the design parameters of an ice-soil wall and cooling regimes requests taking into account the main features of the mechanical behavior of freezing saturated soils. Mechanical behavior of frozen soils is characterized by the rheological properties [2,3]. Therefore, the rock pressure could induce a significant creep deformation of the ice-wall during shaft sinking. In geotechnical engineering to predict the creep deformation of the wall, Vyalov's constitutive relations are widely used. In particular, the well-known analytical expression for assessing the thickness of the ice-soil wall by the critical deformation criterion was derived on the basis of the relations [4,5]. Experimental validations of the relations were performed by uniaxial and multiaxial creep tests of frozen soils at various negative temperatures and ice content [6].
Another important feature of the mechanical behavior of the soil during artificial freezing is a change of the stress and strain fields due to the frost heave and consolidation [7]. The transformation of pore water into ice and the migration of water from unfrozen soil into the freezing front leads to a volumetric expansion of the pore space in the frozen zone of the soil. At the same time, the water outflow causes a reduction of the the pore space volume in the unfrozen zone. During the freezing, the change of the volume can induce a significant rise in stress. As a result, hazardous geotechnical problems, such as a fracture of the ice-soil wall and soil inrush into the excavation can arise.
To take into account the main features of the freezing process of soils, thermo-hydromechanical models have been developed. Mu &Ladanyi [8] proposed one of the first models to describe coupling between heat transfer, water migration, and evolution of stress and strain fields during the freezing process. Additional volumetric strain was used in the model to describe the frost heave of the freezing soil. To estimate creep deformation of the frozen soil, the Prandtl-Reuss law was used. In a model of Liu & Yu [9], volumetric strain due to the phase change was determined in an analogous manner, but volume increment induced by water migration was estimated by the change of the matric potential. As a result, not only the frost heave deformation in the frozen zone of soil was described, but also the consolidation deformation in the unfrozen zone. The model was used to numerically study of the influence of seasoning temperature variations on the pavement. Another way to derive the thermo-hydro-mechanical models of soil freezing is based on the effective stress principle. Neaupane&Yamabe [10] applied the principle to develop a coupled model considering nonlinear elasto-plastic deformation of the freezing soil. To describe the inelastic strain, the Mohr-Coulmb's criterion was adopted. The applicability of the model was demonstrated by a numerical simulation of the deformation process of an underground cavern. In a thermo-hydro-mechanical model developed by Nishimura et al. [11], the mechanical behavior of frozen and unfrozen soil was described based on the extended Barcelona basic model. The results of a numerical simulation of frost heave deformation of soil placed around a buried gas pipeline showed good agreement with field measurements. Bekele et al. [12] proposed a coupled model of soil freezing which included the Clausius-Clapeyron equation and an expression for estimating the volumetric expansion and contraction strains. The model was calibrated on field measurements of a soil-pipeline interaction. Panteleev et al. [13] developed a model of artificial freezing of a rock mass taking into account heat transfer, water filtration, and the evolution of an elastic stressstrain state. Conducted numerical simulations allowed to determine the influence of thermogravitational convection and seepage flow on the formation of an ice-soil wall. Zhou &Meschke [14] formulated the governing equations for the soil freezing model in the framework of thermo-poromechanics theory by Coussy [15,16]. This approach provided for deriving an equation for estimating porosity depending on volumetric strain, temperature, and pore pressure. The applicability of the model was demonstrated by numerical simulations of the use of AGF for construction of tunnels. Another numerical model based on the Coussy's theory was proposed by Tounsi et al. [17]. The model was adapted to large-scale simulations of response of rock masses to the AGF and excavation works. Numerical predictions with the model have shown good agreement with field measurements of temperature and frost heave displacements. A promising approach to modelling soil freezing was developed by Zhou & Li [18]. In this approach, the pore pressure was estimated by the void ratio of soil and mechanical stress, and the void ratio was computed from the mass balance equation. Also, a concept of separating void ratio as a criterion of the formation of ice lens was proposed. In the thermo-hydro-mechanical model proposed by Lay et al., the void ratio of soil in the governing equations was changed on the porosity of the soil. The model validation has been carried out on the basis of the data of one-side freezing tests of soil samples. It was shown that the numerical results are in good agreement with measurements of the evolution of frost heave deformation and temperature, and the distribution of the equivalent water content.
A specific feature of the models by Zhou & Li, Lay et al. is accurate simulation of the evolution of porosity, distribution of the equivalent water content, deformation of soil, and changes in the pore pressure during the process of freezing. However, these models can be applied only to study one-dimensional problems. In the present study, an attempt to extend the models to solving three-dimensional problems of artificial freezing of water-saturated soil was made. The mechanical behavior of saturated soil under freezing conditions is described within the framework of the Coussy theory of porous media using the effective stress concept. The Clausius-Clapeyron equation is applied to determine the relation between temperature, ice pressure, and water pressure in the frozen zone. Temperaturedependent ice saturation is used for modeling the phase change. Heat transfer by conduction and convection is taken into account in the energy conservation equation. Also

Mathematical formulation of the model
Upon freezing, soil is supposed to be a three-phase porous medium consisting of solid grains (index s), liquid water (index l), and crystals of ice (index i). To develop a thermohydro-mechanical model, the following assumptions were made according to Mu  The effect of salinity of pore water on the freezing process is not considered.

5.
Ice movement relative to solid grains is not taken into account, i.e., ice and solid grains are moving as a whole.
6. The segregation of ice during freezing is ignored. 7.
It is assumed that saturated soil is an isotropic and homogeneous porous material that undergoes small deformation.
8. Rheological properties of saturated soil are described by viscoelastic strain. On the basis of the above assumptions, the governing equations of the model can be written as follows: The mass balance equation: The equilibrium equation: The energy conservation equation: div grad grad where ρ j S j n is the mass content of water (j = l) and ice (j =i) at time t; ρ j is the density and S j is saturation of the phase j; n is the porosity; v l is the velocity of water relative to the solid skeleton; σ is the total stress tensor; γ is the unit weight of the porous medium; T is the temperature; C is the volumetric heat capacity and λ is the thermal conductivity of the porous media; Q ph is a heat source related to the latent heat due the water phase change; Q sk is a heat source caused by the soil deformation.
In the mass balance equation S i is expressed by a power function of the temperature T [20] Where T ph is the freezing temperature of pore water and αis an experimental parameter. The condition of the fully saturation implies that S l = 1 -S i . The relative velocity of water v l can be described by the Darcy law as Where k is the hydraulic conductivity and ψ is the water head. The hydraulic conductivity k is supposed to be dependent on the temperature T by the following way [21]: Where k 0 is the hydraulic conductivity of the unfrozen soil, β is an experimental parameter. The water head ψ driving the pore water migration is defined as Where p l is the pressure of the pore water, g is the gravitational acceleration, z is the vertical coordinate. The unit weight of the saturated porous media γcan be written as According to the effective stress principle, the total stress σ in the equilibrium equation is defined by the following way where  σ is the effective stress, p por is the pore pressure, b is the effective Biot coefficient, I is the identity tensor. The effective stress  σ is given by the relation of isotropic linear elasticity as Where K is the effective bulk modulus, G is the effective shear modulus, e ε is the elastic strain and e vol  is the volumetric part of the tensor e ε . According to the principle of the additive decomposition of the total strain ε , the elastic strain e ε can be expressed as Where ve ε is the viscoelastic strain due to rheological properties of the saturated soil and th ε is the thermal strain. The total strain ε is defined by the geometric relation for a case of infinitesimal deformation: Where u is the displacement vector. The thermal strain is written as Where a T is the volumetric thermal dilation coefficient and T 0 is an initial temperature of the unfrozen soil. The viscoelastic strain ve ε is given by the Vyalov's relations, which for multiaxial stress conditions can be expressed as Where ξ, m, ω are material parameters, σ eq is the equivalent stress. The equivalent stress σ eq is defined by where dev  σ is a deviator part of the tensor  σ .
The pore pressure p por is assumed to be weighted sum of the pore water pressure p l and the pore ice pressure p i [22] (1 ) Where χ is the pore pressure parameter, such that (Lay et al. 2014) 1.5 (1 ) Thepore ice pressure p i is expressed by the Clausius-Clapeyron equation as follows Wher ep hydr is the hydrostatic pressure in soil before the freezing process. The expression for the pore water pressure p l can be obtained from (16) An estimation of the pore pressure p por is conducted on the expression derived in the theory of the porous media by Coussy (2005Coussy ( , 2010 Where n 0 is the initial porosity; N is the effective Biot tangent modulus; B is a material parameter which determines an influence of the viscoelastic strain  (14) is zero, in this case it does not influence on the porosity n and the pore pressure p por .
The effective mechanical properties are computed as Where X is the effective value, X fr and X un are values for the frozen and the unfrozen states.
The volumetric heat sources in the energy conservation equation is written as follows The volumetric heat capacity C and the thermal conductivity λ are determined as [8] (1 ) s s Where c j , λ j (j = s,l,i) are the specific heat capacities and the thermal conductivities of the phase j.

Computer implementation
The governing equations were implemented using the ComsolMultiphysics® software. Porosity n, displacement vector u and temperature Twere considered as primary field variables. The spatial discretization was performed using the finite element method. The field variables were approximated by linear Lagrange shape functions. For temporal discretization, the backward Euler scheme was used.

NUMERICAL SIMULATION
The proposed model was used for numerical simulation of the formation the ice-soil wall by the method of artificial freezing and the stress-strain state of the excavation wall inside the formed ice-soil wall. The simulation was carried out for a silty sand stratum lying at a depth of 74.5-75.97 m in a potash deposit in the Republic of Belarus. The thermal regimes of artificial freezing and design parameters of the ice-soil wall were specified on the basis of the technical and design documentation of the Belaruskali Company for the construction of vertical mine shafts.
The layout of forty-one freeze boreholes used for the AGF is presented in Figure 1a. The boreholes are located in a circle with a radius of 8.25 m. The radius of the freezing borehole is 7.3•10 -2 m. The distance between two boreholes is 1.11 m. It is assumed that angles of inclination of the boreholes from the vertical direction can be neglected. Therefore, the artificial freezing process can be studied in a domain bounded by two symmetry planes, projections of which are shown in Figure 1a with red lines. The first plane passes through the center of the freezing borehole. The second plane passes through the middle of the distance between the chosen borehole and the neighboring one. In what follows, we will call this plane the middle plane.
The geometry of the computational domain and a layout of the boundary conditions are shown in Figure 1b. The distance from the center of the projected excavation to the outer boundary is 16.5 m. On the outer boundary of the lateral pressure P lat is given. Temperature T and porosity n on the surface are supposed to be constant and equal to the initial values T 0 and n 0 . The overburden pressure P ob acts on the top boundary. Vertical displacement u z on the bottom surface is fixed. Displacement u on the borehole wall is constrained in the horizontal direction. The temperature of freezing T well is set constant on the borehole wall. Porosity n on the wall is assumed equal to 1.09n 0 . Due to the symmetry condition, displacement u at the inner edge is constrained in the horizontal direction. Material parameters used in the numerical simulation were determined on the basis of laboratory tests performed by the Belaruskali Company. The parameters are listed in Table  1 and Table 2. The initial and the boundary conditions are listed in Table 3. Simulation of artificial freezing was conducted for a period of 85 days.
The distributions of temperature T and porosity nare shown in Figure 2, 3 after 7 days, 20 days, and 85 days of the freezing. It can be seen that the formation of an ice-soil wall is accompanied by intensive frost heave and soil consolidation. The porosity in the frozen zone is increased by 40%. At the same time, in the unfrozen zone, the porosity decreases by 17.5%. Since water migrates from the unfrozen zone to the frozen one, the propagation of the freezing front to the center of the projected excavation increases the shrinkage of the soil. As a result, during freezing, the porosity at the inner side of the ice-soil wall decreases by more than 7.5 % than that at the outer one. Figure 4a shows the distribution of the porosity throughout the ice-soil wall by the 85th day of the freezing. Also, the field of water velocity v l is presented. The distribution of the porosity is non-homogeneous. The water flow is directed to the freezing front. Due to cryogenic suction, the flow stops only when the temperature is below the freezing point T ph . This fact agrees with experimental observations. As for the outer side, there is free water supply from the surrounding soil, which induces a more intensive frost heave.     As a result, the porosity at the outer side is 16% higher than at the inner one. Also, a domain with reduced porosity is observed near the middle plane. As the ice-soil wall moves to the middle plane, cryogenic suction causes water to flow out of the domain towards the phase transition front. Therefore, shrinkage of the soil occurs in the domain. Figure 4b shows the profile of the volumetric ice content along a line segment from the wellbore to the middle plane. It can be seen that the volumetric ice content in the domain is 58% less than that of the wellbore. Figure 5 presents the distributions of volumetric strain ε vol and mechanical pressure P after 85 days of the freezing. The distribution of strain qualitatively coincides with the distribution of porosity, which emphasizes the relationship between heat transfer, mass transfer, and stress-strain state. The distributions of mechanical pressure and volumetric strain reflect the processes of frost heave and soil consoldation. It can be seen that artificial freezing causes a significant change in the stress-strain state of the soil stratum. Due to the volumetric expansion of the soil in the frozen zone, the pressure takes a negative value, and the strain is positive. At the same time, the volumetric shrinkage of the soil in the unfrozen zone leads to the opposite situation. In the domain with reduced porosity near the middle plane, the pressure is positive. Therefore, it can be concluded that the soil in the domain remains in the consolidation state at a negative temperature.
The pressure profile along the line segment from the center of the projected excavation to the middle point on the outer boundary is shown in Figure 6a. The pressure near the icesoil wall increases. Near the inner side, the pressure is 32% higher than on the outside. During shaft sinking, high pressure could lead to significant deformation of the excavation wall.
From Figures 5 (a) and 6 it can be seen that the mechanical pressure and the equivalent stress reach high values. This may be related to an intensive volumetric expansion of the silty sand due to frost heave and high stiffness of the soil in the frozen state. However, it is known from the experimental data that the frost heaving pressure is not so significant. Therefore, some type of inelastic strain arises in the freezing soil due to frost heave, and the constitutive relations of poroelasticitymust be expanded to account for this strain. The  Further study was devoted to the influence of excavation works on the stress-strain state of the retaining ice-soil structure. For this, only the mechanical part of the developed model was used.
The retaining ice-soil structure is taken into account on the 85th day of freezing. The excavation radius is 5.25 m. The depth coincides to the thickness of the soil stratum. The excavation wall is free of any constraints on the displacement and loadings. The simulation is performed within 24 hours, which requires the installation of a mine shaft lining.
In Figure 7a, the distributions of the radial displacement and the vector filed of the displacement u over the soil stratum are shown at the start of subsidence. Figure 7b shows the displacement distributions 24 hours after the start of the deformation process of the excavation wall. It can be seen that on the outside of the ice-soil wall, the mechanical pressure is so significant that the frozen soil moves towards the surrounding soil, despite the lateral pressure. Thus, it can be concluded that the thickness of the ice wall is sufficient to limit the lateral pressure. In the middle of the wall, the radial displacement is minimal. On the inner side of the wall the displacement directs to center of the cut. The performed shaft sinking has a significant impact on the ice-soil wall. The displacement on the inner side of the wall is increased by 20 % and on the outside it is reduced by 7%. Due to the high mechanical pressure inside the ice-soil wall, the soil moves into the excavation. The radial displacement reaches a maximum value of 9.3 cm on the excavation wall. It should be noted that the movement of the soil is caused not by the lateral pressure, but frost heave caused by freezing.

CONCLUSION
The paper proposes a three-dimensional thermo-hydro-mechanical model of freezing of water saturated soil. The governing equations of the model include the mass balance equation, the equilibrium equation, and the energy conservation equation. The phase transition of water into ice is modeled through the temperature-dependent ice saturation. To describe the interaction between the transfer of pore water and the evolution of the stressstrain state of the soil, the Coussy's theory of porous media and the principle of the effective stress are applied. Along with the Clausius-Clapeyron equation, this approach makes it possible to establish the relationship between pore pressure, temperature, porosity and volumetric strain. The rheological behavior of saturated soil is described by viscoelastic strain based on the Vyalov's constitutive relations. The nonlinear equations of the model are implemented in the ComsolMultiphysics® software and solved numerically using the finite element method.
The developed model was applied to numerically simulate the formation of an ice-soil wall by the AGF and shaft sinking into a silty sand stratum at a depth of 74.5-75.97 m in a potash deposit. The results of the simulation show that artificial freezing causes intensive frost heave of the soil in the frozen zone and consolidation of the soil in the unfrozen zone. Due to frost heave, the porosity increases by 34 %. The distribution of volumetric ice content along the ice-soil wall is nonhomogeneous. The free flow of water from the surrounding soil to the outside of the wall increases the ice content. As a result, the ice content on the side is 16% higher than on the inside. Also, the migration of water to the freezing front during the formation of the ice-wall induces the consolidation of soil in the area near the middle plane. The ice content in the area is 58% less in comparison to that near the wellbore. A numerical study of the mechanical behavior of the soil allows us to conclude that freezing causes a significant change in the stress-strain state of the soil stratum. An increase in mechanical pressure and equivalent stress is observed near the sides of the ice-soil wall. The pressure on the inside is 32 % higher than on the outside. That is due to the intensive consolidation of unfrozen soil inside the wall. The results of the simulation of shaft sinking show that the excavation in the soil stratum has a significant impact on the ice-soil wall. Under the impact of the pressure, unfrozen soil moves into the excavation. Radial displacements on the excavation wall reach 9.3 cm. At the same time, significant pressure on the outer side of the ice-soil wall opposes the influence of the lateral pressure, therefore, the frozen soil moves towards the unfrozen zone.
The pressure and stress obtained in the numerical simulation reach high values, which may be associated with the intensive volumetric expansion of the soil due to frost heave, high stiffness of the soil in the frozen state, and the use of the constitutive relations of poroelasticity in the model. Expansion of the relationships that take into account inelastic strain could allow a more correct description of the stress state of the freezing soil.
The reported study was funded by RFBR, project number 19-31-90107.