2. Mining Institute, Ural Branch of Russian Academy of Sciences, Perm 614007, Russia
One of the characteristic features of the modern underground construction is an increase in the depth of deposits and, as a consequence, complication of geotechnical conditions. This circumstance is directly connected with an increase in opening depth of unstable fluidsaturated rocks requiring special building methods. One such method is artificial ground freezing which is used for vertical shaft sinking. Effectiveness of this technique is confirmed by longterm experience of its application in Russia and other countries and is determined by the reliability of the ice wall with calculated thickness. The Verkhnekamsk potash deposit in the Urals, Starobinsk potash deposit in Belorussia and the Gremyachinsk potash deposit in the Volga region were developed with the use of artificial ground freezing method.
Modern development of artificial ground freezing technique is directly related to an increase in freezing depth. In turn, this is connected with economic efficiency issues and problems of guaranteed closure of individual ice bodies to a single wall along with achieving the given strength as soon as possible. Application of this technique was developed in response to the 1990's, with recommendations and analytical techniques for the definition of freezing time, the required thickness of the ice wall, and diameter of the circle of the freezing columns location. The current recommendations and analytical techniques often lead to the two limiting cases: creation of nontight ice wall and formation of an ice wall with excessive thickness which is not effective from an economic point of view. One of the possible solutions to both problems is direct numerical simulation of the ice wall formation in fluidsaturated rocks under specific hydrogeological conditions.
Solution to the problem of artificial ground freezing with the use of direct numerical simulation depends on two factors: initial data on the thermal and mechanical properties of rocks, and adequacy of the used physical model. Artificial ground freezing is a complex multiphasic phenomenon including closely related thermal, mechanical and hydrodynamic processes. Modern physical models which are used in civil engineering, extractive industry, soil science and agricultural engineering can be divided into the following groups: rigidice (Miller, 1978; O'Neill and Miller, 1985; Sheng et al., 1995 ), thermodynamic (Harlan, 1973; Guymon and Luthin, 1974; Jame and Norum, 1980; Konrad, 1994; Noborio et al., 1996a ,b; Hansson et al., 2004 ; Nishimura et al., 2009 ), semiempirical (Konrad and Morgenstern, 1981; Nixon, 1992) and poromechanical models (Coussy, 2004; Coussy and Monteiro, 2007, 2008). Thermodynamic models can be divided into two classes: freezing models of fully saturated porous media (Mikkola and Hartikainen, 2001; Kruschwitz and Bluhm, 2005) and models which take into account unfrozen water after phasetransition (Rempel et al., 2004 ; Zhou and Meschke, 2013). Despite the fact that later models describes physical processes of individual grains and pores more accurately in order to describe phase transition at large spatial scales, separate investigations are still required of unfrozen water on the phase transition rate and coupled filtrationmechanical processes within porous media upon freezing.
This study is devoted to the development of a coupled thermohydromechanical model of fluidsaturated poroelastic media under an assumption of the full transition of liquid into ice. This model is used for numerical simulation of artificial ground freezing Starobinsk deposit taking into account available geophysical and physical data.
2 Mathematical formulation of thermohydromechanical boundary value problemThe object of simulation is layered porous media. In general, this media is a threephased material consisting of a dry skeleton, fluid and ice that completely fill the pore space. In the initial configuration porous media is completely filled with water. All liquid transforms into ice and dry skeleton remains unchanged during the phase transition process. Developed model is restricted to isotropic linear theory of thermoporoelasticity within which it is assumed that the media undergoes small deformations, and changes in density and porosity in each phase are small compared with the initial configuration. It is also assumed that the liquid is nearly incompressible.
The complete threedimensional formulation of the model under the aforementioned assumptions includes the heat Equation (1), Fourier law (2), equilibrium Equation (3), constitutive Equations for the mechanical behavior description (4), geometric relation for linear strain tensor (5), continuity Equation (6), and the Darcy law (7) which can be written in the following form:
$\rho {c_p}\frac{{\partial T}}{{\partial t}} + {\rho _f}{c_{p,f}}\bar v \cdot \bar \nabla T + \bar \nabla \cdot \bar q = 0,$

(1) 
$\bar q =  k\bar \nabla T,$

(2) 
$\bar \nabla \cdot \tilde \sigma = \rho \bar g,$

(3) 
$\tilde{\sigma }=\tilde{\tilde{C}}:\left( \tilde{\varepsilon }{{{\tilde{\varepsilon }}}_{T}} \right){{\alpha }_{B}}{{p}_{f}}\tilde{E},$

(4) 
$\tilde \varepsilon = \frac{1}{2}\left[ {\bar \nabla \bar u + {{\left( {\bar \nabla \bar u} \right)}^T}} \right],$

(5) 
${\rho _f}S\frac{{\partial {p_f}}}{{\partial t}} + \bar \nabla \cdot \left( {{\rho _f}\bar v} \right) =  {\rho _f}{\alpha _B}\frac{{\partial {\varepsilon _{vol}}}}{{\partial t}},$

(6) 
$\bar v =  \frac{{k'}}{\mu }\left( {\bar \nabla {p_f} + {\rho _f}g\bar \nabla z} \right).$

(7) 
where ρ is effective density of the system "dry skeletonliquidice" (kg/m^{3}), c_{p} is effective heat capacity of the system "dry skeletonliquidice" at constant pressure (J/(kg·K)), k is effective thermal conductivity coefficient of the system "dry skeletonliquidice" (W/(m·K)), T is absolute temperature (K),
Expressions for density of the considered system taking into account phase transition from state 1 (fluid) to state 2 (ice) and influence of the pore pressure p_{f} and temperature T have the form:
$\rho = \theta {\rho _{sf}} + (1  \theta ){\rho _{si}},$

(8) 
${\rho _{sf}} = {\rho _f}\left[ {1 + {\alpha _f}(T  {T_0}) + {\chi _f}{p_f}} \right]n + {\rho _s}(1  n),$

(9) 
${\rho _{si}} = {\rho _i}\left[ {1 + {\alpha _i}(T  {T_0}) + {\chi _i}{p_i}} \right]n + {\rho _s}(1  n),$

(10) 
where θ characterizes fraction of fluid phase in the material and (1−θ) is fraction of ice, ρ_{s} is density of the dry skeleton (kg/m^{3}), χ_{i} is ice compressibility coefficient (Pa^{−1}), ρ_{i} is ice density (kg/m^{3}), α_{f} is thermal expansion coefficient of the fluid (K^{−1}), α_{i} is thermal expansion coefficient of ice (K^{−1}), χ_{i} is compressibility coefficient of ice (Pa^{−1}), and T_{0} is initial temperature (K).
The following relation was used in order to define effective specific heat at constant pressure:
$\begin{array}{l}{c_p} = \displaystyle \frac{1}{\rho }\left[ {\theta {\rho _{sf}}\left\{ {n{c_{p,f}} + \left( {1  n} \right){c_{p,s}}} \right\}} \right.\\\quad \quad \left. { + (1  \theta ){\rho _{si}}\left\{ {n{c_{p,i}} + (1  n){c_{p,s}}} \right\}} \right] +\displaystyle L\frac{{\partial {\alpha _m}}}{{\partial T}},\end{array}$

(11) 
where
Effective thermal conductivity was defined as:
$k = \theta \left\{ {n{k_f} + (1  n){k_s}} \right\} + \left( {1  \theta } \right)\left\{ {n{k_i} + (1  n){k_s}} \right\},$

(12) 
where k_{f} = thermal conductivity of the fluid (W/(m·K)), k_{s} = thermal conductivity of the dry skeleton (W/(m·K)), and k_{i} = thermal conductivity of ice (W/(m·K)).
Influence of the θ on the permeability coefficient k' was taken into account with the use of Heaviside function:
$k'\left( \theta \right) = \left\{ {_{ \displaystyle 0,\theta = 0}^{ \displaystyle k',\theta = 1}} \right.$

(13) 
Permeability coefficient k' was estimated by the values of the filtration coefficient k_{f}' of every layer, obtained during hydrogeological investigations of this rock mass:
$k' = \frac{{{k_f}'\mu }}{{\rho g}},$

(14) 
where μ = 10^{−3} = dynamic viscosity of fluid (Pa·s).
Effect of the volumetric strain, pressure and temperature was taken into account as Equation (4):
$n = {n_0} + {\alpha _B}{\varepsilon _{vol}} + \frac{1}{N}{p_f}  3{\alpha _s}\left( {{\alpha _B}  {n_0}} \right)\left( {T  {T_0}} \right),$

(15) 
where N = Biot tangent modulus (Pa), α_{s} = thermal expansion coefficient of the dry skeleton (K^{−1}), and n_{0} = initial porosity.
As it has been mentioned above, every layer is isotropic. Results of laboratory studies of the elastic properties K_{si} and G_{si} of watersaturated rock samples under negative temperature were used to define bulk modulus K and shear modulus G of the dry skeleton. Assuming that effective elastic modulus of the frozen rock mass is determined by the mixture rule, elastic modulus of the dry skeleton K and G can be found as
$\begin{array}{l}K = \displaystyle \frac{{{K_{si}}  n{K_i}}}{{1  n}},\$12pt]G = \displaystyle \frac{{{G_{si}}  n{G_i}}}{{1  n}},\end{array}$

where K_{i} = 1.8 × 10^{9} Pa, G_{i} = 1.4 × 10^{9} Pa, which are bulk modulus and shear modulus of ice. Values of K and G for every considered layer are presented in Table 1.
The Biot coefficient and the Biot tangent modulus were defined as:
$\begin{array}{l}{\alpha _B} = \displaystyle1  \frac{{{K_d}}}{K}\$12pt] \displaystyle \frac{1}{N} = \frac{{{\alpha _B}  n}}{K}\end{array}$

where K_{d} is bulk modulus of the dry porous skeleton calculated by the formula:
${K_d} = \frac{K}{{\displaystyle 1 + \frac{{3\left( {1  \nu } \right)n}}{{2\left( {1  2\nu } \right)\left( {1  n} \right)}}}}$

where
Thermal strain
${\tilde \varepsilon _T} = {\alpha _s}\left( {T  {T_0}} \right)\tilde E$

(16) 
Equations (1)~(15) are supplemented with the following initial and boundary conditions:
${p_f}\left( {t = 0} \right) = 0$

(17) 
$\bar v\left( {t = 0} \right) = \bar 0$

(18) 
$\bar u\left( {t = 0} \right) = \bar 0$

(19) 
${T_i}\left( {t = 0} \right) = {T_{{o_i}}},\;\;i = 1,..., n$

(20) 
${\left. {{u_x}} \right_{{\Gamma _j},{\Gamma _j}}} = 0,\;\;j = 1,..., m$

(21) 
${\left. {{u_y}} \right_{{\Gamma _j},{\Gamma _l}}} = 0,\;\;j = 1,..., m$

(22) 
${\left. {\bar u} \right_{{\Gamma _b}}} = \bar 0$

(23) 
${\left. {\tilde \sigma \cdot \bar n} \right_{{\Gamma _b}}} = \rho \bar gz$

(24) 
${\left. {\tilde \sigma \cdot \bar n} \right_{{\Gamma _u}}} = \rho \bar gz$

(25) 
$ \bar n \cdot \bar q\left {_{{\Gamma _l},{\Gamma _u},{\Gamma _b}}} \right. = 0$

(26) 
$ \bar n \cdot \rho \bar v\left {_\Gamma = 0} \right.$

(27) 
$\begin{aligned}&{\left. T \right_{{\Gamma _j}}} = {T_1}\left( {t,x} \right) \\&{T_1}\left( {t,x} \right) =  ({T_{pod}}(t)  {T_{obr}}(t))x/{h_1} + {T_{pod}}(t),\;j = 1,..., m\end{aligned}$

(28) 
where
The initial data were generated close to real observed situation. Artificial freezing zone includes 41 freezing well. Centers of these wells are located on the circle which has a diameter of 10.5 m and are located at a distance of 1,608 mm from each other. The diameter of the freezing well is 146 mm and its depth is 225 m. To simulate influence of the real geometry of the wells after drilling, artificial inclinometer data were generated. Deviation of the bottom of each well from the vertical was taken into account during the simulation. Inclinometer data are presented in Table 1.
Data on the coolant temperature supply and return presented in Table 2 were used to define boundary conditions.
According to geological data, sedimentary cover consists of 30 layers up to a depth of 225 meters. All 30 layers were grouped into 13 layers with similar thermophysical properties and permeability. Averaged properties of each layer were calculated as the arithmetic mean of the physical properties combined in one group. The averaged thermophysical properties of the layers are presented in Table 3. The following abbreviations are used in the Table 3: ρ_{p} is density of soil particles (kg/m^{3}), T_{cr} is crystallization temperature of water in the layer (K).
Averaged mechanical properties of the layers are presented in Table 4.
The considered area is a cylinder with a depth of 225 meters and diameter of 26.5 meters. This area has been divided by finite elements which have a form of rectangular prisms. The size of the finite element was 9~11 centimeters in the vicinity of the well and 6~7 meters in the periphery. The height of the elements was no more than 6 meters. Finite element model of one of the layers is demonstrated in Figure 1. The total number of elements was about 1 million.
Numerical solution of the coupled boundary value problem of artificial ground freezing for the time interval of 300 days was obtained. Figure 2 and 3 show evolution of the phase transition front.
Average thickness of the ice wall after 282 days was calculated in order to estimate thickness variation. The results are presented in Figure 4. Average thickness of the ice wall is 3.32 meters (standard deviation is 0.86 m).
An important factor affecting the process of ice wall formation is seepage flow of the fluid filling pore space. Moving fluid is a convective source of heat transfer. A change in the fluid pressure or density due to the phase transition leads to change in the stressstrain state of the solid skeleton.
Figure 5 shows distribution of the vertical and horizontal component of the Darcy's velocity in the xOz section in various time moments after the beginning of the freezing process. Z axis is directed downward.
Distribution of the horizontal component of the Darcy's velocity shows characteristic "mosaic" structure of the alternating areas with a multidirectional fluid motion. The velocity amplitude varies from layer to layer within ±1.6 millimeters per day.
In the vertical direction the fluid rises up along the wells and boundaries of the ice wall. After that, it changes direction. Amplitude of the Darcy's velocity exceeds the amplitude of the horizontal component and varies within ±2.5 millimeters per day. The presented figures of the Darcy's velocity field indicate the presence of convective flows of a certain topology near the wells. Structure of these flows can be seen more clearly in Figure 6.
Analysis of the obtained stream lines shows that after the beginning of the phase transition process formation, two toroidal convective cells are observed. One is formed inside of the ice wall and the other is formed outside of the ice wall. Fluid filling pore space rises up along the cold walls of the freezing arch inside of each cell. Furthermore, the fluid moves along the layer boundary and goes down in the middle of the inner surface of the freezing arch in the second cell and along the lateral boundaries in the first cell. Thus, after the beginning of the phase transition process, seepage flows induced by thermogravitational convection are observed in each layer with the nonzero permeability.
4.4 Effect of constant seepage flow in the layer on the ice wall formationThe hypothesis of the immobility of groundwater in the aquifers (relative to time if ice wall formation) is quite rough because the presence of the constant seepage flow can be induced by the peculiarities of the geological structure of the region as well as its geodynamic behavior. Series of numerical experiments on the ice wall formation in the most heatconducting layer (layer 10) with a sufficiently high permeability were carried out in order to investigate the effect of these flows on the phase transition rate. Moreover, the 10^{th} layer has a pronounced spreading of the freezing wells location in the vertical direction relative to location of the wellheads.
Numerical simulation was carried out in the aforementioned coupled mathematical formulation. Constant horizontal seepage flow in the layer was simulated by applying Darcy's velocity to the left and right boundaries of the layer. The calculations were made for the three values of the seepage flow velocities: 17.3, 65 and 173 millimeters per day. In case of the zero seepage flow the ice wall closure is observed on the 70 th day. Temperature distribution after 50 and 150 days from the beginning of the freezing process is presented in Figure 7. All results are presented for the central section of the computational area by the xOy plane.
In case of the presence of horizontal seepage flow with a velocity of 17.3 millimeters per day (minimum of the considered velocities) ice wall closure is observed 38 days later than in the case of the absence of such flow (108 days after the beginning of the freezing process). Figures 8 and 9 represent evolution of temperature and second phase (ice) distributions.
These distributions show that the fluid flow induces insignificant transfer of cold water from the zone of the ice wall formation. In case of the presence of several nonclosed zones, flow inside of the ice wall shifts for the outflow of the fluid through these zones, thereby making the closure even more difficult. Figure 10 presents distribution of the ycomponent of the Darcy's velocity. It can be noted that fluid velocity in these zones increases significantly. After 50 days from the beginning of the freezing process fluid velocity in the nonclosed areas exceeds prescribed flow on the boundaries 8.5 times (fluid velocity is about 140 mm per day).
Abutting into the ice wall fluid flow induces exceeding of the pore pressure which is leveled in the region outside of the ice wall (Figure 11). The value of the pressure jump does not exceed 5% of the hydrostatic pressure at this depth. In case of constant flows with velocities of 65 and 173 mm per day, complete closure of the ice wall are not observed for the considered time (300 days) and influence of the flow on the ice wall formation is the most pronounced (Figure 12).
According to the obtained results, transfer of the cooled fluid by the flow is very significant. Therefore, there are closure and growth of the two disconnected ice zones on both sides of the flow and creation of the two separate ice clusters between them, integration of which is difficult. Fluid velocity between separate clusters increases from 6 (prescribed velocity is 173 mm per day) to 8 (prescribed velocity is 65 mm per day) times as well as in the case of the minimum prescribed velocity (Figure 13). In addition to the increase in the value of the horizontal fluid velocity component there is observed insignificant vertical component which does not exceed 0.2% of the y Darcy's velocity component (Figure 14).
Due to the nonuniform thickness of the ice wall (because of the vertical temperature gradient of the refrigerant) before the phase transition front vertical component of the Darcy's velocity field has a positive direction (Figure 14) which does not depend on the velocity of the seepage flow in the layer. There is observed local increase in pore pressure ahead of the ice wall at velocities of 65 and 173 mm per day as well as in the case of the minimal seepage flow (Figure 15). Pore pressure variations do not exceed units of percent of the hydrostatic pressure corresponding to this depth.
Thereby, results of the simulation show that presence of the horizontal seepage flow significantly affects ice wall formation by slowing its growth.
5 SummaryThis study is devoted to the development of a multiphysics model for the artificial ground freezing process. The proposed model includes solution to several problems: the problem of nonstationary thermal conductivity with phase transition, the problem of linear filtration, and the problem of thermoelasticity. This model is based on the following hypotheses: isotropy of physical properties in the rock layers, the absence of the unfrozen water after phasetransition, small deformation of the porous media, and weak compressibility of the fluid.
The proposed model was used for the numerical simulation of the ice wall formation in waterinfiltrated soil designed for the construction of a vertical mine at Starobinsk deposits taking into account accompanying filtration and mechanical processes. Numerical simulation was carried out by the finiteelement method. This paper has detailed information on issues concerning determination of material parameters.
Results of numerical simulation allowed us to define time of the ice wall closure, thickness of the ice wall in each layer, and average thickness of the whole ice wall. It has been shown that ice wall forms nonuniformly. This can be connected with deviations of the wells from the vertical as well as with the difference in thermal and filtration properties of the layers. Analysis of the stream lines distribution shows that after the beginning of the phase transition, formation of two toroidal convective cells was observed. One is created inside the ice wall while the other is created outside of the wall. The fluid filling the pore space rises up along the cold walls of the ice wall, moves along the layer boundary and goes down in the middle of the inner region of the ice wall in the second cell and along the lateral boundaries in the first cell. Thereby, there are observed seepage flows induced by the thermogravitational convection after the beginning of the phase transition in each layer with the nonzero permeability.
The influence of constant horizontal seepage flow on the ice wall formation was investigated. Numerical simulation of the ice wall formation in the layer with a seepage flow with several prescribed velocities (17.3, 65 and 173 millimeters per day) was carried out. Results of the simulation show that in the case of the absence of the flow, the closure of the ice wall is observed after 70 days from the beginning of the freezing process. In case of the presence of the flow with the minimal considered velocity (17.3 millimeters per day), the closure is observed 38 days later (after 108 days). In case of the fluid flow with the velocities of 65 and 173 millimeters per day, complete closure of the ice wall was not observed for the considered time (300 days). Thereby, results of the simulation show that presence of the horizontal seepage flow significantly affects ice wall formation.
Acknowledgments:This work was supported by the Russian Science Foundation (Grant No. 171101204).
Coussy O, 2004. Poromechanics. Chichester: Wiley, pp. 298.

Coussy O, Monteiro P. 2007. Unsaturated poroelasticity for crystallization in pores. Computers and Geotechnics, 34(4): 279290. DOI:10.1016/j.compgeo.2007.02.007 
Coussy O, Monteiro P. 2008. Poroelastic model for concrete exposed to freezing temperature. Cement and Concrete Research, 38(1): 4048. DOI:10.1016/j.cemconres.2007.06.006 
Guymon GL, Luthin JN. 1974. A coupled heat and moisture transport model for arctic soils. Water Resources Research, 10(5): 9951001. DOI:10.1029/WR010i005p00995 
Hansson K, Šimůnek J, Mizoguchi M, et al. 2004. Water flow and heat transport in frozen soil: numerical solution and freezethaw applications. Vadose Zone Journal, 3(2): 693704. DOI:10.2113/3.2.693 
Harlan RL. 1973. Analysis of coupled heatfluid transport in partially frozen soil. Water Resources Research, 9(5): 13141323. DOI:10.1029/WR009i005p01314 
Jame YW, Norum DI. 1980. Heat and mass transfer in a freezing unsaturated porous medium. Water Resources Research, 16(4): 811819. DOI:10.1029/WR016i004p00811 
Konrad JM. 1994. Sixteenth Canadian geotechnical colloquium: frost heave in soils: concepts and engineering. Canadian Geotechnical Journal, 31(2): 223245. DOI:10.1139/t94028 
Konrad JM, Morgenstern NR. 1981. The segregation potential of a freezing soil. Canadian Geotechnical Journal, 18(4): 482491. DOI:10.1139/t81059 
Kruschwitz J, Bluhm J. 2005. Modeling of ice formation in porous solids with regard to the description of frost damage. Computational Materials Science, 32(3–4): 407417. DOI:10.1016/j.commatsci.2004.09.025 
Mikkola M, Hartikainen J. 2001. Mathematical model of soil freezing and its numerical implementation. International Journal for Numerical Methods in Engineering, 52(5–6): 543557. DOI:10.1002/nme.300 
Miller RD, 1978. Frost heaving in noncolloidal soils. In: Proceedings of the Third International Permafrost Conference. Edmonton, Canada, pp. 708–713.

Nishimura S, Gens A, Olivella S, et al. 2009. THMcoupled finite element analysis of frozen soil: formulation and application. Géotechnique, 59(3): 159171. DOI:10.1680/geot.2009.59.3.159 
Nixon JF. 1992. Discrete ice lens theory for frost heave beneath pipelines. Canadian Geotechnical Journal, 29(3): 487497. DOI:10.1139/t92053 
Noborio K, McInnes KJ, Heilman JL. 1996a. Twodimensional model for water, heat, and solute transport in furrowirrigated soil: I. Theory. Soil Science Society of America Journal, 60(4): 10011009. DOI:10.2136/sssaj1996.03615995006000040007x 
Noborio K, McInnes KJ, Heilman JL. 1996b. Twodimensional model for water, heat, and solute transport in furrowirrigated soil: II. Field evaluation. Soil Science Society of America Journal, 60(4): 10101021. DOI:10.2136/sssaj1996.03615995006000040008x 
O'Neill K, Miller RD. 1985. Exploration of a rigid ice model of frost heave. Water Resources Research, 21(3): 281296. DOI:10.1029/WR021i003p00281 
Rempel AW, Wettlaufer JS, Worster MG. 2004. Premelting dynamics in a continuum model of frost heave. Journal of Fluid Mechanics, 498: 227244. DOI:10.1017/S0022112003006761 
Sheng DC, Axelsson K, Knutsson S, 1995. Frost heave due to ice lens formation in freezing soils 1: Theory and verification. Hydrology Research, 26(2): 125–146.

Zhou MM, Meschke G. 2013. A threephase thermohydromechanical finite element model for freezing soils. International Journal for Numerical and Analytical Methods in Geomechanics, 37(18): 31733193. DOI:10.1002/nag.2184 