Coupled thermo-hydro-mechanical modeling of frost heave and water migration during artificial freezing of soils for mineshaft sinking

M.Zhelnin, A.Kostin, A.Prokhorov, O.Plekhov, M.Semin, L.Levin

a Institute of Continuous Media Mechanics of the Ural Branch of the Russian Academy of Sciences, Perm, 614013, Russia

b Mining Institute of the Ural Branch of the Russian Academy of Sciences, Perm, 614007, Russia

Keywords:Artificial ground freezing (AGF)Thermo-hydro-mechanical (THM) modeling Frost effects Frozen wall Shaft sinking

ABSTRACT Artificial freezing of water-bearing soil layers composing a sedimentary deposit can induce frost heave and water migration that affect the natural stress-strain state of the soil layers and freezing process.In the present paper, a thermo-hydro-mechanical (THM) model for freezing of water-saturated soil is proposed to study the effects of frost heave and water migration in frozen soils on the formation of a frozen wall and subsequent excavation activity for sinking a vertical shaft.The governing equations of the model are formulated relative to porosity, temperature, and displacement which are considered as primary variables.The relationship between temperature, pore water, and ice pressure in frozen soil is established by the Clausius-Clapeyron equation, whereas the interaction between the stress-strain behavior and changes in porosity and pore pressure is described with the poromechanics theory.Moreover, constitutive relations for additional mechanical deformation are incorporated to describe volumetric expansion of soil during freezing as well as creep strain of soil in the frozen state.The ability of the proposed model to capture the frost heave of frozen soil is demonstrated by a comparison between numerical results and experimental data given by a one-sided freezing test.Also to validate the model in other freezing conditions, a radial freezing experiment is performed.After the validation procedure, the model is applied to numerical simulation of artificial freezing of silt and sand layers for shaft sinking at Petrikov potash mine in Belarus.Comparison of calculated temperature with thermal monitoring data during active freezing stage is presented.Numerical analysis of deformation of unsupported sidewall of a shaft inside the frozen wall is conducted to account for the change in natural stress-strain state of soil layers induced by artificial freezing.

1.Introduction

Construction of vertical mineshafts for mining potash salt deposits is technically complicated due to hydrogeological conditions.The deposits are usually located below water-bearing soft soils which are very unstable.To provide safe shaft sinking through the soil strata,the artificial ground freezing(AGF)method is applied.For vertical shaft sinking,AGF isrealized bya brinenetworkconsistingof a refrigerant plant and freezing pipes(Klein,1981;Andersland and Ladanyi, 2004).The pipes are embedded into the freezing boreholes which are drilled along the contour of projected shaft before excavation works.Circulation of chilled brine through the pipes in the refrigerant plant withdraws heat from the ground.As the ground temperature decreases below the freezing temperature,pore water turns into ice and a wall of frozen soil is formed around the projected shaft.In contrast to their natural state, frozen soils are stable and water impermeable,thus the frozen wall serves as watertight temporary excavation support during the shaft sinking.

AGF is an environmentally friendly technique that allows one to conduct soil stabilization and waterproof sealing in geotechnical and civil engineering construction problems.For decades, the technique has been actively applied for tunnel construction(Jessberger, 1980; Russo et al., 2015), shaft sinking (Vyalov et al.,1979; Levin et al., 2021), and mining (Yun et al., 2017) in complex hydrogeological conditions.Nevertheless, the use of AGF is frequently challenging due to the presence of regional seepage flow in the ground and specific features of freezing of water-saturated soft soils.To extend the applicability of AGF and improve its performance,mathematical models of ground freezing are developed.

One of the first thermo-hydraulic models describing both heat transfer in frozen soil accounting for the latent heat of phase transition and water migration to the freezing front was proposed by Harlan (1973).The mathematical formulation of soil freezing derived by Harlan(1973)has given rise to many subsequent thermohydraulic models (Kurylyk and Watanabe, 2013).Modern thermohydraulic models are successfully utilized to study the influence of heat convection induced by groundwater seepage flow on the AGF process(Vitel et al.,2016;Huang et al.,2018; Li et al.,2020,2021).Theoretical formulation of a large number of thermo-hydraulic models in the context of simulation and optimization of AGF has been recently reviewed by Alzoubi et al.(2020).

AGF of soils can provoke unfavorable uplift of the surface induced by frost heave(Han et al.,2016;Cai et al.,2019;Zheng et al.,2021; Zhou et al., 2021).In particular, during mineshaft sinking,frost heave of soils can lead to a rise in loading acting on excavation sidewall and shaft lining inside the frozen wall (Orlov and Kim,1988; Ji et al., 2019).The rigid ice model developed by Miller(Miller, 1972; O’Neill and Miller, 1985) firstly provided accurate evaluation of frost heave deformation taking into account water migration to the freezing front and ice lens formation.The rigid ice model assumes the existence of a frozen fringe zone between the freezing front and the warmest ice lens.Based on this conception,Gilpin theory (Gilpin, 1980) and segregation potential model(Konrad and Morgenstern, 1981) were developed.The rigid ice model and Gilpin theory are limited only to one-dimensional (1D)problems.For instance, Gilpin theory was used for the elaboration of a separate ice model (Ji et al., 2018; Zhou et al., 2018).On the other hand, segregation potential has been used for developing models that are able to predict frost heave during tunnel construction (Konrad, 2002; Zhou et al., 2021).

Along with the segregation potential, additional volumetric strain tensor was used in a number of works(Michalowski and Zhu,2006;Yang et al.,2006; Cai et al.,2019;Zheng et al.,2021) to take into account the frost heaving phenomenon.Michalowski and Zhu(2006) adapted the porosity rate function and unit growth tensor for estimating the volumetric tensor.Yang et al.(2006) computed the volumetric strain due to the frost heave using the empirical relation for water inflow into the frozen zone.In the thermomechanical model proposed by Cai et al.(2019), the frost heave volumetric strain was provided by the unit growth tensor and frost heave ratio.Also, Zheng et al.(2021) proposed a thermomechanical model and computed the frost heave volumetric strain using the frost heave ratio, which is an integrated approach for predicting the ground surface displacement during AGF.

More complete description of the freezing process of soils can be achieved using thermo-hydro-mechanical (THM) models.The models are able to simulate coupling between heat transfer,water migration to the freezing front due to cryogenic suction,evolution of porosity, and stress-strain state induced by the frost heave during the freezing process of soil.One of the first THM models was developed by Mu and Ladanyi(1987).In the model,the frost heave deformation of frozen soil was described by additional volumetric expansion related to density change due to the phase transition and water accumulation behind the freezing front.The pore pressure was evaluated using the Clausius-Clapeyron equation.Also, the Prandtl-Reuss law was incorporated in the model to determine the creep strain of soil in the frozen state.The effect of frost heave on the stress-strain state of frozen soil was described using an additional volumetric strain in the models proposed by Liu and Yu(2011) and Li et al.(2015, 2018).In the model of Liu and Yu(2011), the volumetric strain caused by water movement was estimated based on the matric potential.Water transfer was described by an extended Richard’s equation with a term related to ice formation that was determined using the Clausius-Clayperon equation.The model was used to study the freezing of pavement in cold region.The model of Li et al.(2015,2018)was developed to explore the frost damage mechanism of a water canal in cold region.The expression for the volumetric strain caused by frost heave was quite similar as that proposed by Mu and Ladanyi (1987).However,the water transport equation was expressed by moisture gradient according to Jame and Norum (1980).Also, constitutive relations were adapted for the computation of viscoplastic strain related to the creep behavior of frozen soil.

Another way to simulate the effect of frost heave on the stressstrain state of frozen soil is given by the principle of effective stress.The early model of Neaupane and Yamabe(2001)takes into account the volumetric expansion of soil due to phase transition induced by changes in pressure head and plastic strain determined using the Mohr-Coulomb criterion.In contrast,a recently proposed model by Bekele et al.(2017) considers frozen soil as a poroelastic material.The model includes Bishop-type effective stress in terms of pore water pressure and ice pressure as well as additional volumetric strain related to the phase transition.Liu et al.(2019)developed a model using the effective stress principle and the porosity rate function similar to that proposed by Michalowski and Zhu(2006) and adapted the model to numerical simulation of AGF in double-barrel culvert sites.Coussy(2005)proposed poromechanics theory which was used by Zhou and Meschke(2013)for simulation of AGF in tunnel construction taking into account the effect of seepage flow on the formation of a frozenwall and cryogenic suction mechanism on the frost heave.Panteleev et al.(2017) adapted the poromechanics theory for simulation of AGF for vertical shaft sinking taking into account thermo-gravitational convection and seepage flow.Another model including the constitutive relations of the poromechanics theory was derived by Tounsi et al.(2019) in a thermodynamically consistent framework according to Vitel et al.(2016) and Rouabhi et al.(2018).The model was applied to a case study of mining operations using AGF.Recently,the poromechanics theory of frozen soils was extended to take into account plastic deformation and viscoplastic damage (Liu et al., 2018; Liu and Lai,2020).However, the models are quite general and have not been applied for computing the stress-strain state of frozen soil in engineering applications.

An alternative way for describing the mechanical behavior of frozen soils was developed by Nishimura et al.(2009).In the model,mechanical constitutive relations with two stress variables were derived based on Barcelona basic model (Alonso et al.,1990).The first variable is similar to the effective stress and the second variable represents the effect of cryogenic suction on the stress-strain state.The constitutive relations were modified by Ghoreishian Amiri et al.(2016a).Modified constitutive relations assume that effective stress depends only on pore water pressure and the effect of cryogenic suction on the stress-strain state is taken into account by elastic and plastic volumetric strains computed using the suction variable.On the basis of the relations,a THM model was developed and applied for numerical simulation of frost heave induced by a chilled pipeline (Haxaire et al., 2017).The extended constitutive relations accounting for viscoplastic strain can be found in Ghoreishian Amiri et al.(2016b).

Asafoundationfor modernTHMmodels,therigidice modelisalso used.The model enables predicting the frost heave deformation and stress-strain state of frozen soil taking into account the formation and development of discontinuous ice lenses but only 1D freezing processes can be simulated.In a model proposed by Thomas et al.(2009),only water pore pressure was considered,and the impact of the pore pressure on a stress-strain state was given through the Terzaghi type effective stress.To compute the pore pressure, the mass conservation equation including thermal gradient was solved.A new ice lens initiation was predicted by a stress criterion.Anotherapproachfordescribingthe freezing processofsoilwasestablished by Zhou and Li(2012)and Lai et al.(2014).In the model,the deformation of soil was caused by porosity change,which was determined by the massconservationequation.Theequivalentpore pressure includedin the Bishop type effective stress was evaluated by the momentum equilibrium equation.As a criterion for the formation of a new ice lens,the critical porosity was employed.Application of the model to numerical simulation of one-sided freezing experiments has shown a good agreement between computed and measured distributions of porosity and ice lenses, evolution of temperature, and frost heave deformation.Also,the model isableto predict suchessential effects of cryogenic suction as negative pore pressure and consolidation of soil nearthe freezing front.Subsequently,the models have been extended to take into account the influence of pore water salinity(Zhang et al.,2018) and vapor phase (Yin et al., 2018) on the freezing process of soils.

The above literature review has shown that THM models can provide a sufficiently reasonable and adequate description of the freezing process of water-saturated soils.THM models derived within the framework of the poromechanics theory of Coussy have been successfully applied for numerical simulation of AGF in tunnel construction(Zhou and Meschke,2013),mining operations(Tounsi et al.,2019),and vertical shaft sinking(Panteleev et al.,2017)taking into account water movement and frost heave.At the same time,the approach solving the mass conservation equation relative to porosity realized in 1D models by Zhou and Li(2012)and Lai et al.(2014) seems to be promising to achieve adequate prediction of changes in porosity and soil consolidation related to water migration into the frozen zone in frozen soil.

In this paper,a THM model for freezingof water-saturated soil was developed using the poromechanics theory of Coussy to investigate the mechanical behavior of soil.Water transfer in frozen soil is described according to Zhou and Li (2012) and Lai et al.(2014) by solving the mass conservation equation relative to porosity and using theClausius-Clapeyronequation.Changeinequivalent porepressure is expressed by the porosity and the volumetric strain.Furthermore,constitutive relations forplastic volumetric strain induced by the frost heave in the frozen soil and viscoelastic strain related to rheological properties of frozen soil are included.The governing differential equations of the developed model were incorporated in Comsol Multiphysics?and solved numerically using the finite element method(FEM).On the basis of the model,numerical study of AGF for vertical shaft sinking at the Petrikov potash deposit in Belarus was carried out under specific geotechnical conditions.

The paper is organized as follows.In Section 2,the mathematical formulation of the THM model is proposed.In Section 3,the ability of the developed model to describe the freezing process of soil is demonstrated by comparing the results of numerical simulations of laboratory tests with experimental measurements.Further, numerical simulation of artificial freezing of two soil layers at the Petrikov potash deposit in Belarus is performed in Section 4.In Section 5, numerical simulation of deformational process of unsupported circular excavation inside the frozen wall formed by AGF is conducted for the soil layers.Also in Sections 5 and 6,the effect of cryogenic processes on stress distribution and radial displacement is discussed.Finally, the conclusions of the study are presented in Section 6.

2.Theoretical framework

Soil upon freezing is considered as a three-phase porous medium, consisting of solid grains (s), pore water(l), and pore ice (i).Initially,the soil is unfrozen and fully saturated with pore water.A mathematical formulation of the process of soil freezing is derived following mechanics of continuum media on the basis of fundamental conservation laws of mass, energy, and momentum.In order to formulate the governing equations for the freezing process of soil, the following hypotheses are adapted according to Mu and Ladanyi (1987), Zhou and Li (2012), Zhou and Meschke(2013), Lai et al.(2014) and Tounsi et al.(2019):

(1) Local thermal equilibrium is assumed.Therefore, the three phases have the same temperature.

(2) The freezing point of pore water is constant.The freezing point depression due to loading and pore water salinity is negligible.

(3) Density change of each phase is not taken into account during the freezing process.

(4) Movement of pore ice relative to solid grains is not taken into account.

(5) Formation of ice lenses in the frozen soil is ignored.

(6) The soil is an isotropic material with mechanical properties depending on the temperature.The mechanical behavior of the soil is described using the small-deformation theory.

2.1.Water transfer equation

The water transfer equation can be written on the basis of the mass conservation law of pore water and ice as follows:

where ρjSjn is the mass content of pore water(j=l)and ice(j=i),ρjis the density (kg/m3), Sjis the saturation of the phase j, n is the porosity,vlis the velocity of water relative to the solid skeleton(m/s), and t is the time (s).

According to Tice et al.(1976),the ice saturation Siis expressed by an empirical function of temperature T(K)in the following form:

where Tphis the freezing temperature of pore water(K),and α is an experimental parameter.As the soil is supposed to be fully saturated, the pore water saturation Slis estimated as Sl= 1- Si?

The relative velocity of water vlis governed by Darcy’s law as

where k is the hydraulic conductivity(m/s),and ψ is the soil-water potential.The hydraulic conductivity k is given by a power function of the temperature as (Nixon,1991):

where k0is the hydraulic conductivity of unfrozen soil,and β is an experimental parameter.The soil-water potential ψ driving migration of the pore water is defined as

where plis the pore water pressure (Pa), g is the gravitational acceleration (m/s2), and z is the vertical coordinate(m).

According to Zhou and Li (2012) and Lai et al.(2014), the pore water pressure plis expressed by the equivalent pore pressure p and the pore ice pressure pi.The equivalent pore pressure p in thefrozen soil is determined as follows (O’Neill and Miller, 1985;Coussy, 2005):

where χ is a pore pressure parameter,and χ = (1-Si)1?5?The pore ice pressure is estimated using the Clausius-Clapeyron equation as(Thomas et al., 2009; Zhou and Li, 2012):

where p0is the initial pressure(Pa),and L is the latent heat of fusion(J/kg).Then the pore water pressure plcan be expressed from Eqs.(6) and (7) as (Zhou and Li, 2012; Lai et al., 2014):

2.2.Heat transfer equation

The energy conservation law provides the heat transfer equation for frozen soil as

where C is the volumetric heat capacity(J/(K m3)),λ is the thermal conductivity of soil(W/(K m)),clis the specific heat capacity of the water(J/(kg K)),and Qphis the heat source induced by realization of the latent heat during water freezing (W/m3).

The volumetric heat capacity C and the thermal conductivity λ are respectively estimated using the following equations (Mu and Ladanyi,1987):

where cj(j=s,l,i)is the specific heat capacity of the phase j(J/(kg K)); and λjis the thermal conductivity of the phase j.

The heat source Qphcaused by the phase transformation of pore water into ice is written as

2.3.Equilibrium equation

The equilibrium equation is written on the basis of the momentum conservation law and the poromechanics theory of Coussy(Coussy, 2005; Zhou and Meschke, 2013) as follows:

where σ is the total stress tensor of soil,and γ is the unit weight of soil.

The unit weight γ is estimated as

where g is the gravitational acceleration vector.

According to the conception of effective stress,the total stress σ can be written in the following form:

where σ′is the effective stress,b is the effective Biot coefficient,and I is the identity tensor.As the soil is assumed to be isotropic, the relationship between the effective stress σ′and the elastic strain tensor εecan be expressed in the following form:

where K is the effective bulk modulus(Pa),G is the effective shear modulus(Pa),εeis the elastic strain,andis the volumetric part of the tensor εe.

Applying the principle of additive decomposition of the total strain ε, the elastic strain εeis defined as

where εthis the thermal strain, εfhis the plastic volumetric strain related to frost heave of frozen soil,and εveis the viscoelastic strain related to creep of frozen soil.The total strain is expressed by the displacement vector u according to the small strain approximation:

The thermal strain εthis given by the following expression:

where aTis the linear expansion coefficient (1/K), and T0is the initial temperature(K).

Within the poromechanics theory,the equivalent pore pressure can be evaluated from the following expression (Coussy, 2004,2005):

where n0is the initial porosity, and N is the effective Biot tangent modulus (Pa).

The effective mechanical properties are calculated as(Zhou and Meschke,2013):

where X is the effective value;and Xfrand Xunare the values for soil in the frozen and unfrozen states, respectively.

2.4.Mechanical constitutive relations for plastic and viscoelastic strains

During the freezing process, unfrozen pore water in soil transforms into pore ice.Due to the difference in densities of water and ice,the phase transformation can induce frost heave characterized by the volumetric expansion of frozen soil.In fine-grained soils,the frost heave is significantly intensified by cryogenic suction that leads to water migration to the freezing front from the unfrozen zone.In this case, volumetric strain in the frozen zone of finegrained soils can attain considerable values, thus elastic constitutive relations can predict excessive tensile stress.A more reasonable description of the relationship between stress and strain of frozen soil is provided by elastoplastic constitutive models.In the model proposed by Nishimura et al.(2009), part of the volumetric expansion caused by ice pressure is related to plastic volumetric strain, which enables restricting tensile stress.In the model of Ghoreishian Amiri et al.(2016a),volumetric expansion induced byfrost heave is captured by elastic and plastic parts of total volumetric strain, which are expressed in terms of the cryogenic suction.

Following the approach in the present study,plastic volumetric strain εfhis incorporated into the model to simulate the mechanical behavior of frozen soil:

where dλ is the plastic multiplayer,F is the yield surface,and σ′mis the mean effective stress.The mechanical behavior of soil under tension is assumed to be perfectly plastic.The yield surface is determined in the following form:

where the material parameters A and B are expressed by the soil cohesion c (Pa) and angle of internal friction φ as

The plastic volumetric strain εfhis related to volumetric expansion of frozen soil due to the phase transformation of pore water into ice and water migration into the frozen zone.According to the proposed constitutive relations,the strain εfharises when the effective mean stress σ′mexceeds the tensile strength determined by the yield surface F.During the freezing process, the porosity n rises due to frost heave of soil caused by the phase transformation of water into ice and water inflow into the frozen zone.From Eq.(20),porosity growth leads to a rise in the equivalent pore pressure p that contributes to an increase in the effective mean stress as defined by Eq.(15).It is assumed that the plastic volumetric strain εfhis related to plastic dilation of frozen soil mainly induced by the formation of a massive cryogenic structure.If the freezing process of soil is accompanied by intensive ice segregation with the formation of thick ice lenses,then the mechanical behavior of the soil should be described taking into account discrete ice lenses distribution (Thomas et al., 2009; Zhou and Li, 2012; Lai et al., 2014).

Rheological properties of soil in the frozen state are of the utmost importance along with the frost heave, as under external loading,significant creep deformation can evolve.To predict creep deformation of frozen soil, the following constitutive relation for the viscoelastic strain εveis employed (Vyalov et al., 1979;Andersland and Ladanyi, 2004):

where ξ is the total deformation modulus (Pa); m and ω are the material parameters; τ is the dimensionless time; andis the equivalent stress,which is defined as

where devσ′is the deviatoric part of the stress tensor σ′?

The relation given by Eqs.(26)and(27)describes an increase in viscoelastic strain εveof frozen soil with time.The viscoelastic strain εveis included in the model to compute creep deformation in a frozen wall during excavation works(Section 5)and it is not taken into account for numerical simulation of soil freezing.

Capability of the relation(Eq.(26)and(27))to predict the creep deformation of frozen soil is constrained by a decaying creep process.Also, it is assumed that the creep behavior of frozen soil is independent of mean stress and the volumetric creep strain can be neglected.However,due to its simplicity,the relation is widely used for engineering computations (Andersland and Ladanyi, 2004).In particular, the well-known analytical expression for assessing the optimal thickness of the frozen wall by the critical deformation criterion was derived on the basis of the relation (Vyalov et al.,1979).

2.5.Computer implementation of the model

The proposed THM model of frozen soil (Eq.(1)-(27)) was solved in COMSOL Multiphysics?software using the FEM.On the first step,the porosity n,the temperature T and the displacement u were chosen as primary field variables.Weak formulations of the governing equations (Eq.(1), (9), and (13)) relative to these variables are given as

The weak formulations are supplemented by initial and boundary conditions.The initial conditions include initial distributions of temperature T0and porosity n0in the computational domain Ω.The Dirichlet boundary conditions(first type)are given by time-dependent changes in the temperature T, the porosity n,and the displacement u on parts of the boundary ?Ω.The Neumann boundary conditions (second type) are imposed by timedependent changes in the thermal flux -λgradT?n,the mass fluxρlk gradψ?n, and the total stress σ?n on parts of the boundary ?Ω.

Further, the weak formulation of the water transfer equation(Eq.(28)) was incorporated into COMSOL Multiphysics?using the Weak Form PDE (partial differential equation) interface.This interface provides a convenient way to solve user-defined PDEs by representing them in the weak forms and simply putting them into the Weak Expression field of the interface.The software automatically performs the spatial discretization applying user-defined shape functions with the given element order and computes the weak form by element-wise integration over the computational domain Ω using the Gaussian quadrature method(Belytschko et al.,2013).

For the spatial discretization of the weak form(Eq.(28)),linear Lagrange shape functions are employed.Before performing the discretization, the soil-water potential ψ is expressed by the primary field variables combining Eqs.(5), (8) and (20).According to Eq.(5), the soil-water potential ψ depends on the pore water pressure pldetermined from Eq.(20) which includes the temperature T and equivalent pore pressure p defined by Eq.(20).Eqs.(5),(8) and (20) are incorporated into the software as user-defined variables.

The weak formulations of the heat transfer equation (Eq.(29))and the equilibrium equation(Eq.(30))are provided automatically by the Heat Transfer and the Solid Mechanics Interfaces of the software.

In the Heat Transfer Interface, the Heat Source node is used to take into account the heat source Qph(Eq.(12))which is given as a user-defined variable.Linear Lagrange shape functions are employed for the spatial discretization of the weak form of the heat transfer equation.

In the Solid Mechanics Interface, the effective bulk modulus K and shear modulus G of soil are defined according to Eq.(21)which are included in the software as user-defined variables.The External Stress subnode is employed in order to determine the effective stress(Eq.(15)).The External Strain subnodes were used to specify the thermal strain εthand viscoelastic strain εveby entering Eqs.(19) and (26) for components of the strain tensors, respectively.It should be noted that in numerical simulations of the soil freezing process, the viscoelastic strain is not taken into account.Constitutive relations (Eq.(22)-(25)) for the plastic volumetric strain εfhrelated to frost heave are implemented in the Plasticity subnode by entering Eq.(24)as the yield criterion and using an associated flow rule.Quadratic shape functions are used for approximation of the displacement vector u.

Temporal discretization is performed automatically in the software using Backward Differentiation Formula(Brenan et al.,1995).For the solution to nonlinear algebraic equations, a damped Newton method(Deuflhard,1974) with a constant damping factor equal to 0.8 was used.

3.Model validation

To examine the ability of the model to reproduce key features of freezing of saturated soil, numerical simulations of two laboratory tests have been performed, and obtained model predictions have been compared to experimental measurements.The first test is the one-sided freezing of saturated silty clay columns carried out by Lai et al.(2014)in an open system.The second test is the radial freezing of saturated sand in a closed system.

3.1.Numerical simulation of one-sided freezing of silty clay

The capability of the model to predict frost heave deformation induced by the phase transformation of pore water into ice and water migration to the freezing front is studied on the basis of experimental data provided by Lai et al.(2014).Laboratory tests of one-sided freezing of saturated silty clay columns were conducted in an open system for different freezing conditions.The lateral surface of the samples was thermally insulated.The top surface of samples was subjected to constant negative temperature and uniaxial compressive load.The constant positive temperature was maintained on the bottom surface.During the freezing, water supply was provided to the bottom surface.

In the present study, numerical simulation is conducted for samples with a height and a diameter of 10 cm for two freezing conditions.In the first case, cooling temperature Tcon the top surface is -1.6?C and positive temperature Twon the bottom surface is 1.5?C.Overburden pressure Pobon the top surface is 50 kPa.In the second freezing condition, the top surface of soil sample is exposed to the cooling temperature Tcof -4?C and pressure Pobof 100 kPa.On the bottom surface,the temperature Twis kept at 2.5?C.The initial temperature T0of soil is 3?C in the first test and 2.6?C in the second test.The initial porosity n0in both tests is 0.32.

As described, freezing conditions ensure a unidirectional freezing process,and computational domain has a rectangular form representing half of a vertical cross-section of the soil column.The computational mesh consisted of 620 quadrilateral elements.At the top of the computational domain, constant temperature Tc,porosity nb, and pressure Pobare set.According to Zhou and Li(2012), it is supposed that at the start of the tests, pore water contained in the soil near the top surface quickly transforms into pore ice without water migration from the unfrozen zone,thus an increase in porosity at the top end is caused by only volumetric expansion of original water content by 9%.Therefore, the porosity nbis supposed to be 1.09n0.As during the freezing experiment,the water has been supplied to the bottom end of the soil sample, the initial porosity n0and positive temperature Tware imposed to the bottom end of the computational domain.Displacement at the bottom side is fully constrained.On lateral sides of the computational domain, zero-flux thermal and mass conditions are considered, and displacement is admitted only along the vertical direction.

According to experimental data after 14 h of soil freezing in the first condition, temperature in the soil column stabilized and an increase in frost heave was induced only by the growth of the final ice lens.In the second freezing test, stabilization of temperature in soil column was observed after 20 h.As in the proposed model, discrete thick ice lenses are not taken into account,numerical simulation of laboratory experiments is carried out only for a time period of 14 h for the first condition and 20 h for the second one.Thermal properties of water and ice used for the numerical simulation are listed in Table 1.Material properties of silty clay are presented in Table 2.Values of thermal and elastic properties are similar to those in Lai et al.(2014).

Table 1Thermal properties of water and ice.

Table 2Material parameters of silty clay used for the numerical simulation.

Results of the numerical simulation and experimental data for the second freezing test are presented in Figs.1 and 2.Fig.1 shows variations of experimental and calculated temperatures with time at the heights of 20 mm, 40 mm, 60 mm, and 80 mm from the bottom end of the soil column.At each height, the calculated temperature variation with time is quite close to the experimental measurements.

It can be seen that during 2 h from the start of the test,temperature rapidly decreases at the considered positions of the soil column.After that, a stage of slow cooling can be distinguished.Within the stage, temperature at each height slowly reduces to a minimal value.Due to the small value of the cooling temperature of -1.6?C, temperature decreases below 0?C at heights of 60 mm and 80 mm near the top end.At the positions more distant from the top end, temperature stabilization proceeds more quickly, which can be related to heat influx from the warm bottom end due to thermal conductivity and convective heat transfer induced by the water migration to the freezing front according to Darcy’s law (Eq.(3)) and Eq.(8) for the pore water pressure pl.

Fig.1.Variations of calculated (black lines) and measured (blue lines) temperatures with time at the heights of 20 mm,40 mm,60 mm,and 80 mm from the bottom end of the sample obtained for the first experimental condition(Tc=-1.6 ?C,Tw=1.5 ?C,and Pob = 50 kPa).

Fig.2.Variations of calculated and measured axial displacements on the top end of the sample with time obtained for the first experimental conditions (Tc = -1.6 ?C,Tw = 1.5 ?C, and Pob = 50 kPa).

Fig.2 shows comparison of calculated axial displacement uzat the top end of the sample with experimental measurements.It can be concluded that during the freezing test, uplift of the top end occurs,which indicates frost heave of frozen soil.In the model,the frost heave is described in the following way.The decrease in the temperature T below the freezing temperature Tphleads to the phase transformation of the pore water into ice given by Eq.(2)and the water migration into the frozen zone governed by Darcy’s law(Eq.(3))and Eq.(8)for the pore water pressure pl.According to the water transfer equation(Eq.(1)), the processes lead to an increase in the porosity n in the frozen zone.As the porosity n affects the equivalent pore pressure p according to Eq.(20), the pressure p rises in the frozen zone that contributes by Eq.(15)to the change in the mean effective stress.In turn, the stress-strain constitutive relations (Eq.(16)-(19) and (22)-(25)) eventually give the volumetric expansion of frozen soil and observed increase in the axial displacement uz.

Fig.3.Variations of calculated (black lines) and measured (blue lines) temperatures with time at the heights of 20 mm,40 mm,60 mm,and 80 mm from the bottom end of the sample obtained for the second experimental conditions (Tc = -4 ?C, Tw = 2.5 ?C,and Pob = 100 kPa).

The abovementioned effects illustrate the coupling between thermal, mechanical, and flow processes within the framework of the proposed model and enable us to quite accurately describe not only temperature evolution but also frost heave displacements.It can be seen that the calculated axial displacement presented in Fig.2 is in good agreement with the experimental data during 10 h of the freezing test.After that, the predicted axial displacement is less relative to the measurements, because the growth of the final ice lens cannot be captured by the proposed model(see Fig.2 for ice lenses distribution in Lai et al.(2014)).

Results of the numerical simulation and experimental data for the second freezing test are presented in Figs.3 and 4.Fig.3 shows calculated and experimental curves of temperature variations with time at different positions of the soil column.In contrast to the first freezing condition,in this test,temperature stabilization in the soil column evolves slower, which is related to deeper propagation of the freezing front into the soil column due to lower cooling temperature of -4?C.The calculated curves are in good agreement with the measurements for considered positions except for the height of 60 mm,where predicted temperature smoothly decreases with time,while in the experiment,temperature quickly stabilizes.

Fig.4 shows the calculated and experimentally measured variations of axial displacement uzat the top end of the soil column with time for the second freezing condition.According to the experimental data provided by Lai et al.(2014), soil freezing with lower cooling temperature and larger overburden pressure leads to quicker propagation of the freezing front into the soil column,generally forming a massive cryogenic structure in the absence of thick ice lenses.It can be concluded that in this case,the calculated axial displacements are quite close to the measurements within the considered time period of 20 h.At the beginning stage of the freezing test,the frost heave of the soil column is restricted by the overburden pressure, thus the axial displacement is almost constant.Once the frost heaving pressure overcomes the restricting loading, the axial displacement starts to monotonically rise.After 20 h of the freezing test, the final ice lens arises, thus further freezing is not considered.

Fig.4.Variations of calculated and measured axial displacements at the top end of the sample with time obtained for the second experimental condition (Tc = -4 ?C,Tw = 2.5 ?C, and Pob = 100 kPa).

Within the considered time period, soil freezing in the second test is characterized by dominant massive cryogenic structure without thick ice lenses, which agrees with the assumptions adopted in the developed model.The performed comparison of the calculated transient axial displacement variations to the experimental measurements shows that in this case, the model can accurately describe the frost heave of the soil column with a maximal relative deviation of 3%between calculated and measured values.On the other hand,the predicted axial displacement at the top end of the soil column for the first freezing condition is lower than that measured in the test with a maximum relative deviation of 9% that can be related to growing thick discrete ice lenses.

3.2.Numerical simulation of radial freezing of sand

In Section 3.1,one-sided freezing tests on silty clay columns are considered.Freezing of silty clay is characterized by intensive water migration to the freezing front that significantly promotes frost heave of the soil.However,for underground construction,artificial freezing of saturated sands also has to be performed.Sands are composed of grains of larger size compared to silty clay,therefore,frost heave in frozen sands can be induced only by the volumetric expansion of original water during its transformation into ice.At the same time, water migration to the freezing front due to cryogenic suction in sands is weak or fully absent.

The ability of the proposed model to predict deformation of frozen soil mainly induced by volumetric expansion of original water is analyzed based on the experimental data provided by a laboratory test of radial freezing of saturated quartz sand in a closed system.The freezing test is carried out in the following way.Soil saturated with distilled water is packed into the tough plastic form,as shown in Fig.5a.The plastic form is a cylindrical container with the thickness of the lateral side of 2.8 cm and the thickness of the bottom and top caps of 2 mm.Along the central axis of the cylinder,a copper tube is installed.After soil packing, the resulting soil sample has a hollow cylinder form with a height of 5 mm and the internal and external radii of 1 cm and 5.7 cm, respectively.The radial freezing of the soil is performed by refrigerant circulating through the copper tube.The layout of the freezing system used in the test is shown in Fig.5b.The freezing system consists of a thermally insulated chamber, a freezing chamber, an expansion tank and a pump.During the freezing test,the plastic form is placed inside the thermally insulated chamber to prevent heat influx to the soil sample from the environment.The copper tube is connected to the pump and the expansion tank containing the refrigerant.To control the cooling temperature,the expansion tank is positioned into the freezing chamber where a constant negative temperature of-30?C is maintained.Circulating of the refrigerant along the freezing system is provided by the pump.

In the freezing experiment, measurements of temperature and total radial strain of the frozen soil are carried out using three thermocouples and a fiber-optic strain detector.The location of the thermocouples and the strain detector is shown in Fig.6a.The first thermocouple denoted as TC1 is installed on the surface of the copper pipe to monitor the cooling temperature Tcof soil.The second thermocouple TC2 is located on the lateral side of the plastic container.The third thermocouple TC3 is fixed on the strain detector, which is inserted into the middle of the soil sample at a height of 2.5 mm from the bottom and at a distance of 2.8 cm from the central axis of the plastic form.

The layout of the strain detector is shown in Fig.6b.According to the layout,the strain detector consists of a fiber Bragg grating(FBR)sensor and two small plates fixed at the optical fiber at a distance of 1.2 cm apart from each other.It is supposed that soil deformation in the freezing process causes movement of the plates which can be measured by the FBR sensor.The FBR sensor is placed inside a rigid tube to exclude the mechanical impact of soil particles on it.To perform the temperature compensation of the strain measurements, data of variations of temperature with time registered by the thermocouple TC3 are used.

Fig.5.(a)Water-saturated sand packed in the plastic container;and (b) Layout of the freezing system.1 - Thermal insulation material; 2 - Freezing chamber; 3 - Expansion tank; 4 - Pump; 5 - Plastic container.

Fig.6.(a)Positions of thermocouples TC1,TC2,and TC3 and the strain detector in the plastic container packed by water-saturated sand(all sizes are in mm);and(b)Scheme of the strain detector consisting of optical fiber(1),metallic plates(2),metal protective casing (3), and fiber Bragg grating (FBR) sensor (4).

On the basis of the described experimental conditions,axisymmetric numerical simulation of the freezing process of sand is carried out in a rectangular computational domain with sizes shown in Fig.6a.A computational mesh consisted of 820 quadrilateral elements.Material parameters of the quartz sand are listed in Table 3.Thermal properties of water and ice are the same as those presented in Table 1.The initial temperature T0and porosity n0of the soil are 240?C and 0.35, respectively.The temperature variations measured by the thermocouple TC1 are adopted at the cooled boundary.The measured temperature curve is shown in Fig.7.Also at this boundary, the porosity nbequal to 1.09n0is set and the radial displacement is constrained.To describe confinement of frost heave by sides of the plastic form,the elastic restraint conditions are employed at other boundaries with the stiffness value of 0.12 GPa at the top and bottom sides and of 72 GPa at the lateral side.

Table 3Material parameters of quartz sand.

Fig.8 shows the comparison of variations of the temperature with time at the outer side and the middle of the soil sample obtained in the numerical simulation and measured by the thermocouples TC2 and TC3.It can be seen the calculated temperature curves are in good agreement with the experimental data.A maximum deviation of 2?C between calculated and measured values is reached at an initial stage of the freezing test.Minimal values of the temperature measured by thermocouples TC2 and TC3 are-2?C and-4?C,respectively.At the point corresponding to the position of the thermocouple TC3,decrease in the temperature from 24?C to 0?C is observed for 3 h of freezing.For the remaining 15 h of freezing,the temperature at the point reduces only to-4?C.From the temperature variations measured by the thermocouple TC3, it can be seen that the temperature reaches 0?C for 10.5 h of freezing, and in the remaining 8 h, the temperature reduces to-2?C.The observed retardation of soil freezing can be related to the axisymmetric character of the freezing front propagation.As the freezing front propagates from the cooling tube into the soil, a volume of freezing soil is enlarged, thus the latent heat caused by the phase transition and the heat influx from the unfrozen zone increase.

The variations of measured and predicted total radial strains with time are presented in Fig.9.Comparison between the results of numerical simulation and measurements demonstrates that the model can quite correctly predict the soil deformation in the freezing test.According to the presented numerical and experimental data, some conclusions about the mechanical behavior of the frozen soil can be drawn.At the initial stage of the experiment,negative values of the total radial strain indicate that unfrozen soil contracts at the location of the strain detector due to a negative thermal strain given by Eq.(19) and a mechanical impact induced by the volumetric expansion of the frozen zone caused by a rise in the equivalent pore pressure p according to Eq.(15).As the freezing front attains the strain detector, volumetric expansion of soil caused by the phase transformation of water into ice results in an increase in the total radial strain.At the final stage of the experiment, volumetric expansion of soil near the lateral side of the plastic form leads to contraction of frozen soil at the position of the strain detector.However, owing to the high stiffness of the frozen soil described by Eq.(21), decrease in the total radial strain is less than that at the initial stage of the experiment.Similar mechanical behavior has been detected by strain gages during the freezing experiment conducted by Tounsi et al.(2020) for limestone.

The proposed scheme of the laboratory experiment approximately reproduces soil freezing conditions around a vertical freezing borehole drilled for AGF.Restraint of frost heave of soil in the test by the caps of the plastic container can be related to a situation when artificial freezing soil layer is bounded from the top and the bottom by two rigid soil strata.Restraint of frost heave from the lateral side of the plastic container simulates an effect of unfrozen soil surrounding the area of the freezing activity.

4.Numerical simulation of AGF for shaft sinking

Potash salt at the Petrikov deposit is bedded under watersaturated soil strata.To perform vertical shaft sinking, AGF was carried out from the ground surface to the depth of 265 m.For this purpose,41 freezing boreholes were drilled along the circle with a radius of 8.25 m.According to the project documentation, the freezing borehole radius was equal to 7.3 cm.Distance between the axes of the two neighboring boreholes was equal to 1.1 m.Mineshaft radius in design was equal to 5.25 m.Fig.10 shows the location of the freezing boreholes around the intended mineshaft.

To control the formation of a frozen wall during AGF, thermal monitoring was carried out up to the depth of 255 m.Temperature measurements for two thermal monitoring boreholes are considered.The first borehole denoted as TM1 is located between two freezing boreholes at a distance of 1.25 m from the freezing contour.The second borehole denoted as TM2 is drilled at a distance of 2.25 m from the nearest freezing borehole.Temperature measurements were carried out along the depth of the monitoring boreholes by fiber-optic sensors with an accuracy of 0.3?C.A detailed description of the implementation of AGF and thermal monitoring at the Petrikov potash mine can be found in Levin et al.(2021).

The following hypotheses were supposed to carry out numerical simulation of AGF:

Fig.7.Variation of temperature with time measured by the thermocouple TC1.

Fig.8.Comparison of calculated temperatures with measurements provided by thermocouples TC2 and TC3.

Fig.9.Variations of calculated and measured total radial strain with time.

(1) Heat transfer between neighboring soil layers has a minor effect on frozen wall formation.

(2) Mechanical interaction between neighboring soil layers is neglected.

(3) Water migration between neighboring soil layers is neglected.

(4) Deviation of borehole axis from the vertical direction is not taken into account.

According to the first three hypotheses, each layer can be considered separately.The fourth assumption enables us to use symmetry conditions and reduce the size of the computational domain.These hypotheses considerably reduce numerical costs and enable us to carry out calculations without supercomputers.

Two water-saturated layers are considered.The first one is silt located at the depths of 50-58 m and the second one is sand located at the depths of 65-85 m.

4.1.Geometry and boundary conditions

The geometry of the computational domain and boundary conditions are shown in Fig.11.The area is half of a sector of the circle located between two freezing boreholes.Distance from the center of the mineshaft to the outer boundary of the area is equal to 16.5 m.Other sizes correspond to the project documentation.The upper boundaries of the silt and sand layers correspond to the depths of 50 m and 65 m, respectively.

Fig.11.The geometry of the computational domain and the layout of the boundary conditions.

A computational mesh consists of 38,219 triangular prism elements.The base of each prism is parallel to the horizontal crosssection of the domain.Around the freezing borehole, mesh refinement is performed such that elements adjacent to the borehole boundary have a minimum size.With an increase in distance from the borehole, the element size rises to a maximum value on the outer side of the computational domain.

According to the scheme presented in Fig.11, coolant temperature Tbhis given on the borehole boundary as well as a constant value of porosity nbequal to 1.09n0.Only vertical displacements are allowed on the borehole boundary and the outer boundary of the domain.Constant temperature T and porosity n equal to initial values T0and n0,respectively,are applied on the outer boundary of the domain.Zero-flux thermal and mass conditions are given on the other boundaries.The bottom boundary of the domain is fixed in the vertical direction.The top surface is loaded by the overburden pressure Pob.The symmetry boundary condition is set onthe lateral sides of the area.Also, horizontal displacements at the inner edge of the domain are restricted due to symmetry.

Fig.10.The layout of freezing and thermal monitoring boreholes (TM1 and TM2) in cross-section of soil stratum at the Petrikov potash deposit.Red lines denote the area of interest.

The initial temperature T0,porosity n0,hydrostatic pressure p0and overburden pressure Pobfor silt and sand are given in Table 4.Dependence of the coolant temperature Tbhon the time is presented in Fig.12.The coolant temperature was maintained at the inlet of the freezing pipes during AGF at the Petrikov mine.

4.2.Identification of material parameters

Material parameters used for the numerical simulation of artificial freezing of the silt and sand layers were defined on the basis of experimental results provided by the standard experimental program which was performed before AGF at the Petrikov mine.Elastic, thermo-physical, and filtration characteristics of the considered soils were experimentally obtained for positive and negative temperature values.Identification of other parameters was carried out on the basis of laboratory tests for measuring frost heave ratio which corresponded to Russian union standard GOST 28622-90 (2005).According to this experiment, one-sided freezing of water-saturated cylindrical sample was carried out with a constant temperature in an open system.The diameter and height of the sample were equal to 10 cm.The temperature at the top Tcwas negative and equal to-6?C,and the temperature at the bottom was positive and equal to 1?C.The bottom surface was in permanent contact with water.The lateral sides of the sample were thermally insulated.Freezing of the sample was continuing up to the depth of 9 cm.After that, vertical displacement uz(absolute strain) of the top surface was measured.Axial displacements uzof the top surface in silt and sand samples at the end of the experiment were equal to 6.5 mm and 4.1 mm, respectively.

Numerical simulation of laboratory tests was carried out using the same experimental condition.Due to the symmetry, the simulated domain was a rectangle with dimensions corresponding to half of the cylinder cross-section.The computational mesh consisted of 620 quadrilateral elements.The top surface had a constant cooling temperature Tcand porosity nbequal to 1.09n0.The bottom boundary had a constant positive temperature Twand initial porosity n0.The lateral sides of the sample were imposed to zero thermal and mass fluxes.Horizontal displacements at the lateral surfaces were constrained and the bottom boundary was fixed while the top surface was free of any loadings.

Material parameters of silt and sand used in the simulation are listed in Table 5.Thermo-physical properties of water and ice are given in Table 1.Elastic,thermal and filtration characteristics of the soils were obtained experimentally.Other parameters were identified by a series of calculations to achieve a match between numerical and experimental values of the axial displacement.

Table 4Initial values for the soil layers.

Table 5Material parameters for silt and sand.

Fig.13 shows variations of axial displacement uzat the top end of sample with time and typical distribution of porosity obtained by the numerical simulation for silt and sand samples.The calculated values of uzfor both soils correspond to the measured values at the end of the test.However,the graphs are qualitatively different.Silt freezing is slower than sand freezing due to its thermo-physical properties, higher initial water content, and convective heat transfer induced by water migration.As a result,freezing of the silt sample up to the depth of 9 cm takes 20 h while the freezing of the sand sample takes 11 h.Another factor affecting the shape of the frost heave curve is the water migration to the phase transition front.In the silt sample,frost heave is induced by freezing of initial pore water as well as water supplied from the bottom surface due to the cryogenic suction.As a consequence,the uplift of the silt top surface is uniform and slows down gradually with the increase in the freezing depth.

Frost heave of the sand sample occurs only due to the initial water content.This is the reason why the vertical displacement uzrises sharply at the beginning of the experiment when freezing of most pore water takes place.After that, surface uplift slows down rapidly and uzbecomes almost constant.

Porosity distribution in the silt sample has three different areas.The first area is characterized by a porosity increase under the temperature T lower than Tph.The temperature in the second area is also lower than Tph, but the porosity in this area reaches its minimum value.The third area corresponds to the unfrozen zone.The increase in the porosity in the first area is 12% with regard to the initial value.Therefore, there is an additional water inflow to the frozen zone induced by cryogenic suction.On the other hand, a decrease in porosity in the second and third areas shows consolidation of the soil near the freezing front.

Unlike silt, an increase in porosity for the frozen zone of sand corresponds to 9%, which is induced by volumetric expansion of water initially contained in pores.Porosity distribution shows that it grows monotonously behind the freezing front and is equal to the initial value before the freezing front.

Fig.14 shows distributions of vertical displacement uzalong the depth of the silt and sand samples after 2 h and 1 h,respectively.It can be seen that the presented plots are in qualitative agreementwith porosity distributions.The curve corresponding to silt shows negative displacement uznear the freezing front, which means compression of the solid skeleton in this area.According to Eq.(15)for the effective stress,the equivalent pore pressure p in this area is negative.This is the reason for water migration to the freezing front in the silt sample.Displacement uzin the frozen zone increases monotonously, which corresponds to soil tension under the influence of positive pore pressure.It should be noted that the presented distribution of uzis in agreement with the results predicted by modern rigid ice models (Thomas et al., 2009; Zhou and Li, 2012).Distribution of uzalong the sand sample shows only a monotonous rise behind the freezing front, which displays an absence of cryogenic suction in sand and emergence of the positive pore pressure in the frozen zone.

Fig.12.Variation of coolant temperature Tbh at the inlet of the freezing pipe with time.

4.3.Numerical results of AGF of water-saturated soils

Figs.15 and 16 demonstrate temperature(T)distributions in the sand stratum after 16 d, 36 d, 52 d and 126 d of artificial freezing and in the silt stratum after 16 d, 38 d, 70 d and 210 d of artificial freezing.For clearer representation,the distributions are shown for the whole soil layer.These results were obtained by mirroring the numerical solution for the computational domain (Fig.11) with regard to the symmetry planes.

Two stages of the formation of the frozen wall can be distinguished during artificial freezing.In the initial stage, cylinders of frozen soil initiate and grow around freezing boreholes.Further freezing provides their advance up to the middle plane.Then, the coalescence of the frozen cylinders forms a closed continuous frozen wall.During the second stage, there is an increase in frozen wall thickness caused by freezing of the soil located inside the contour of the freezing boreholes as well as soil adjacent to the outside of the frozen wall.Propagation of the inner boundary of the frozen wall up to the mineshaft contour is achieved after 126 d of freezing for the sand and after 210 d for the silt.Slower freezing of silt is related to higher water content,lower thermal diffusivity, and convective heat transfer to the freezing front.

Fig.17 shows the calculated and measured variations of the temperature with time at the thermal monitoring boreholes TM1 and TM2 at the depths of sand and silt strata.The data are presented for the active freezing stage.It can be seen that the calculated temperature curves are close to the experimental measurements.Due to smaller distance between the thermal monitoring borehole TM1 and the freezing pipe,the temperature at the borehole TM1 decreases more quickly than that at the borehole TM2.The artificial freezing of the sand stratum evolves more rapidly than that of the silt, which is in agreement with the presented temperature distributions.

Fig.13.Variations of the vertical displacement uz at the top surface of silt(a)and sand(c)samples with time(blue dashed line corresponds to the experimental values measured at the end of the tests).Porosity distribution n after 2 h of freezing for the silt sample(b)and after 1 h of freezing for the sand sample(d)(the white line corresponds to the freezing front where T = Tph).

Fig.14.Distribution of the vertical displacement uz along the depth of the silt sample after 2 h of freezing and of the sand sample after 1 h of freezing.

Figs.18 and 19 present distributions of porosity n after 16 d,36 d,52 d and 126 d of artificial freezing for the sand stratum and after 16 d, 38 d, 70 d and 210 d of freezing for the silt stratum.

In the sand, the porosity is homogenously distributed throughout the frozen wall thickness.In the frozen zone,there is an increase of porosity by 9%with regard to the initial value.It can be explained by relating the frost heave of frozen sand only to water expansion initially contained in pores without additional water inflow induced by the cryogenic suction.

The qualitatively different character of porosity distribution is observed in the silt layer.Artificial freezing of this layer is accompanied by intensive frost heave in the frozen zone and consolidation of the soil close to the freezing front.Porosity near the freezing borehole rises by 21% with regard to the initial value, which demonstrates intensive water inflow into frozen wall induced by the cryogenic suction.An increase in the porosity around the freezing pipes due to water migration was also predicted in the numerical simulation of AGF for tunnel construction performed by Li et al.(2021).On the other hand, the outflow of water from unfrozen soil leads to its consolidation and reduction in porosity.Consolidation of the soil rises with the propagation of the freezing front to the center of the mineshaft because there is no external water supply in the unfrozen ground inside the frozen wall.For instance,reductions in the porosity near the internal and external boundaries of the frozen wall are 10%and 8%,respectively,on the 70th day of freezing.After 210 d of freezing,the porosity reduction near the inner boundary of the frozen wall becomes more significant and attains 14%.A similar situation was observed in a one-sided freezing test for fine-grained soil cylinders in a closed system(Hansson et al.,2004).In the test,the porosity in the unfrozen zone near the bottom end of soil samples reduces with time of freezing due to water migration from the unfrozen part of the samples to the freezing front.Also consolidation of unfrozen soil due to water migration was predicted during tunnel construction using AGF by a frost heave prediction method developed in Zhou et al.(2021).

Fig.20a shows porosity distribution and the water velocity vlalong the frozen wall thickness for the silt after 70 d of artificial freezing.It can be seen that the porosity distribution is inhomogeneous and the velocity vector is directed to the frozen wall boundaries due to the cryogenic suction.Water flow stops only behind the freezing front at the temperature below Tph, which agrees with the experimental observations (Miller, 1972; Thomas et al., 2009; Zhou and Li, 2012; Lai et al., 2014).

Water inflow to the external boundary of the frozen wall is higher than that to the internal boundary which leads to more intensive frost heave.This is because the external side of the frozen wall is surrounded by unfrozen soil, which is supplied by thegroundwater.At the same time,a quantity of pore water inside the frozen wall is restricted.As a result,the porosity near the external boundary of the frozen wall is 3%higher than that near the internal boundary.

Fig.15.Temperature T (?C) in sand layer after (a) 16 d, (b) 36 d, (c) 52 d and (d) 126 d of freezing (white line corresponds to the freezing front where T = Tph).

Fig.16.Temperature T (?C) in silt layer after (a) 16 d, (b) 38 d, (c) 70 d and (d) 210 d of freezing (white line corresponds to the freezing front where T = Tph).

Moreover,Fig.20a illustrates an area of lower porosity near the middle plane.During the first stage of the frozen wall formation when the freezing front propagates from the freezing borehole to the middle plane, there is a water outflow from the area to the freezing front due to the cryogenic suction.Additional soil compression in the area is induced by the frost heave evolving in the increasing frozen cylinder.As a result, there is a significant reduction in porosity near the middle plane and soil consolidation.Fig.20b illustrates a profile of volumetric ice content along the line connecting borehole axis with the middle plane.It can be seen that in this area,volumetric ice content is 17%lower than that near the freezing borehole.

During AGF,the natural stress-strain states of the soil layers are changing along with the decrease in temperature and porosity evolution.

Fig.21 shows distributions of the total volumetric strain εvoland the effective mean stress σ′min the sand stratum after 52 d and 126 d of artificial freezing.As water migration to the freezing front in sand is absent,it can be concluded that volumetric expansion of water upon freezing leads to increases in the total volumetric strain and the effective mean stress.It can be explained by the porosity growth of 9%, which induces the increases in the equivalent pore pressure p and the effective mean stress according to Eqs.(15)and(20).In spite of the volumetric expansion of frozen water,the total volumetric strain and the effective stress remain negative in the frozen zone, which means that the overburden pressure Pobsuppresses the frost heaving pressure and compresses the frozen soil.

Besides,Fig.21 illustrates that in unfrozen soil inside the frozen wall,the total volumetric strain and the effective stress are less than those in unfrozen soil surrounding the frozen wall.It can be concluded that due to the frost heaving pressure in the frozen zone,rise in the thickness of the frozen wall leads to contraction of the unfrozen soil inside the frozen wall, which is characterized by decreases in the total volumetric and the mean effective stress.At the same time, the unfrozen soil inside the frozen wall restraints the frost heaving pressure, thus the total volumetric strain and the effective mean stress reduce along the thickness of the inner part of the frozen wall relative to the contour of freezing boreholes.

According to Eq.(20), contraction of unfrozen soil inside the frozen wall leads to an increase in the pore water pressure.Fig.22 shows the variation of pore water pressure plwith time relative to the initial hydrostatic pressure at a point located at a distance of 1.3 m from the center of the intended mineshaft.It can be observed that within the first stage of the formation of frozen wall, the change in the pore pressure is pretty low.After 36 d of freezing,when the closed continuous frozen wall is formed, the pore pressure curve changes qualitatively and the pore pressure quickly rises with time.

Fig.23 presents the total volumetric strain εvoland the effective mean stress σ′min the silt stratum after 70 d and 210 d of the artificial freezing.Similar to the sand, increases in the total volumetric strain and the effective mean stress are observed in the frozen zone of the silt stratum.However, due to water inflow into the frozen zone, artificial freezing of the silt stratum causes a significant frost heaving pressure, which, in spite of the effect of the overburden pressure, causes a change in the initial compressive stress into tensile stress in most part of the frozen wall.Plastic volumetric strain arises where the effective mean stress satisfiesthe yield criterion (Eq.(24)).From Fig.23a and c, it can be concluded that the most considerable plastic volumetric strain evolves near the freezing borehole.At the same time,in the area of lower porosity,the frozen silt remains in a compressive state.In the area,in spite of the phase transformation of water into ice,the total volumetric strain and the effective mean stress take negative values.Another effect of water flow from the unfrozen zone into the frozen zone is soil contraction near the freezing front.From Fig.23b and d,it can be seen that the effective mean stress reduces near the freezing front that contributes to the soil contraction.Also inside the frozen wall,the effective mean stress is less than that outside.It means that the mechanical impact on the unfrozen soil inside the wall induced by the extension of the frozen wall is higher.

Fig.17.Comparison between the calculated and measured temperatures at the positions of the thermal monitoring boreholes TM1 and TM2 in(a)sand and(b)silt during the active stage of freezing.

Fig.24 shows the variation of pore water pressure plwith time in the silt stratum.Within the first stage of a frozen wall formation,the pore water pressure slightly increases.After 38 d,when a closed frozen wall is formed,the pore water pressure quickly rises due to compression of the unfrozen soil inside the wall.At 110 d,the pore pressure attains maximum value and starts to decrease due to the effect of cryogenic suction.Thus,after the closure time,a significant frost heaving pressure evolving in the frozen zone induces contraction of unfrozen soil inside the frozen wall and rise in the pore pressure.However, the cryogenic suction leads to water flow from the unfrozen soil to the inner side of the frozen wall.Inside the closed frozen wall,there is no external water supply,therefore,the water migration leads to decrease in the pore water pressure.

Fig.18.Porosity n in the sand stratum after (a) 16 d, (b) 36 d, (c) 52 d and (d) 126 d of freezing (white line corresponds to the freezing front where T = Tph).

Fig.19.Porosity n in the silt stratum after (a) 16 d, (b) 38 d, (c) 70 d and (d) 210 d of freezing (white line corresponds to the freezing front where T = Tph).

Fig.20.(a) Porosity n and water velocity vl (green arrows) near the freezing front and inside the frozen zone after 70 d of freezing (white line corresponds to the freezing front where T = Tph); and (b) Distribution of volumetric ice content along the line connecting borehole axis with the middle plane.

5.Displacement of unsupported shaft sidewall inside a frozen wall in silt and sand layers

Vertical shaft sinking under the protection of a frozen wall is performed by low depth stopes with subsequent reinforcement of internal excavation boundary by a tubbing lining (Vyalov et al.,1979).Therefore, the internal unsupported boundary of the excavation deforms only on a small part of the area located between the bottom surface of the lining and the bottom of the stope.This unsupported sidewall of the excavation is deformed under the lateral loading from the surrounding unfrozen soil up to the installation of the shaft lining.In this period of time,significant creep deformation of the frozen soil can evolve that may complicate the installation of the shaft lining and induce a collapse of the excavation sidewall.To account for the creep deformation of frozen soil, constitutive relations(Eq.(26)and(27))for the viscoelastic strain are included in the proposed model.To evaluate the impact of frost heave and cryogenic suction evolving the displacement of an unsupported excavation sidewall inside the frozen wall,numerical simulation of stress-strain state in silt and sand layers subjected to excavation works has been carried out.

The following studies for the partially frozen silt layer have been made.Firstly, displacement of the excavation sidewall has been defined taking into account the impact of frost heave and cryogenic suction on the stress-strain state of the silt layer.Secondly, a similar problem has been solved considering the effect of the frost heave only and neglecting the effect of cryogenic suction.Finally, calculation without any effect of frost heave and cryogenic suction has been carried out.For the partially frozen sand layer, only two cases have been considered which correspond to calculations with and without frost heave.The mechanical part of the developed THM model together with the constitutive equation for viscoelastic strain has been used for numerical simulation.

Fig.21.Total volumetric strain εvol(a,c)and effective mean stress σ′m (MPa)(b,d)in sand after 52 d(a,b)and 126 d(c,d)of freezing(white line corresponds to the freezing front where T = Tph).

5.1.Geometry and boundary conditions

The geometry of the simulated domain and boundary conditions are shown in Fig.25.The domain corresponds to the one presented in Fig.11 which takes into account the excavated soil layer.The excavation sizes correspond to the project of mineshaft construction at the Petrikov potash mine (Kostina et al., 2020).The radius and height of the excavation are equal to 5.25 m and 5 m,respectively.

Calculations for each layer are carried out for the thickness of the frozen wall corresponding to the case that the soil has been frozen up to the shaft sidewall.Boundary conditions correspond to those applied for numerical simulation of AGF.Vertical displacement is allowed at the freezing borehole boundary, outer boundary,and inner boundary of the domain.Symmetry condition is set at the lateral sides of the area.The top surface is loaded by the overburden pressure Pob.The bottom of the stope is free from any loadings.The lateral side of the stope is also free and fixed in the horizontal direction at the upper edge which simulates a restraint effect of the lining.

Hydrostatic pressure p0and overburden pressure Pobfor the soil layers are given in Table 4.Distribution of mechanical properties(Table 5)corresponds to temperature(T)distribution obtained from numerical simulation of AGF.Initial stress and pore pressure distributions defined in Section 4.4 are set for the problem which considers effects of both frost heave and cryogenic suction.These distributions are applied for the whole domain excluding unfrozen soil below the bottom of the stope because of the unloading during shaft sinking.Numerical simulation has been carried out for 24 h.According to the project documentation, this time is required to install a tubing lining.

5.2.Identification of rheological properties for silt and sand in frozen state

Creep strain of frozen silt and sand is defined by constitutive equations (Eq.(26) and (27)).Rheological parameters in Eq.(26)were obtained by numerical simulation of mechanical tests for long-term uniaxial loading of cylindrical samples performed as a part of the standard experimental program before AGF at the Petrikov mine.The diameter of the samples was equal to 5 cm and the height was 8 cm.The mechanical tests were performed according to the Russian union standard GOST 12248-96(2005).The uniaxial loading was applied to the top surface for 24 h.Temperature of the samples was maintained at -8?C.

Material parameters used for the modeling of creep curves are presented in Table 6.These parameters were obtained according to the experimental data of axial strain εzfor two different loadings(Pob= 1.13 MPa and Pob= 2.26 MPa).Comparison between calculated and measured axial strain εzwith time is shown in Fig.26.It can be seen that calculated curves are in good agreement with experimental results.For both soils, the maximum deviation does not exceed 4%.

Table 6Creep properties of frozen silt and sand.

Data presented in Fig.26 illustrate the decaying creep behavior of considered frozen soils.An initial increase in creep rate corresponds to the primary stage.After that,a nearly constant creep rate is observed which is referred to as the secondary stage.Although silt and sand samples were subjected to the same loadings, the value of vertical strain for silt is higher than that for sand.Therefore,frozen silt is characterized by a more pronounced rheological flow under long-term loading.

5.3.Numerical analysis of displacement of unsupported shaft sidewall

Fig.27 presents the distribution of radial displacement urin the partially frozen silt layer subjected to excavation activity after 24 h of the deformation process of unsupported excavation wall.The distribution was obtained by numerical simulation of the mechanical behavior of the soil stratum in three different cases.In the first case,numerical simulation was performed taking into account a change in the stress state of the soil layer induced by the frost heave and cryogenics suction evolving during the frozen wall formation.In the second case,only the impact of the frost heave on the stress state was taken into account.In the third case, only nonhomogeneous mechanical properties were given in the soil layer according to the temperature distribution.In all cases,the frozen wall is adjacent to the excavation contour that is reached after 210 d of artificial freezing.From the presented distributions, it can be concluded that the lateral pressure induces bending of the excavation sidewall and uplifting of the excavation bottom due to compression of unfrozen soil inside the frozen wall.Maximum radial displacement is reached on the excavation sidewall,and it is equal to 2.8 cm for the first case,1.9 cm for the second case, and 1.5 cm for the third case.Therefore,the cryogenic processes lead to an increase in the radial displacement of the excavation sidewall.Unfrozen soil surrounding the frozen wall restrains the frost heaving pressure evolving in artificial freezing of soil, thus the lateral pressure acting on the outer side of the wall rises.Water inflow into the frozen zone contributes to the volumetric expansion of frozen soil, therefore, in the first case, the maximum radial displacement is higher than that in the second one.If the cryogenic processes are not taken into account, the radial displacement is underestimated by 47%.

Fig.28 demonstrates the distribution of the radial displacement urin the partially frozen sand after 24 h of deformation process of unsupported excavation wall.The distributions were obtained for two cases.In the first case,the stress state of the soil layer is affected by the frost heave evolving during artificial freezing.In the second case, only nonhomogeneous mechanical properties are considered which are given according to temperature distribution.It can be seen that the presented distributions are in qualitative agreement with the numerical results for the silt stratum.The maximum displacement is attained at the excavation sidewall and it is equal to 1.1 cm for the first case and 0.8 cm for the second computation.Thus, if the effect of frost heave is not taken into account, the radial displacement is underestimated by 27%.

Fig.29 shows profiles of the radial displacement uralong the line segment which crosses the soil layers in the middle plane and is located at the height of 2.5 m from the excavation bottom.The mostsignificant radial displacement on the excavation sidewall in the silt stratum is observed for the computations with consideration of frost heave and cryogenic suction and in the sand stratum for the computation considering frost heave.The absolute value of the radial displacement in the unfrozen soil surrounding the frozen wall almost linearly increases to the outer side of the frozen wall.The frozen wall restraints the lateral pressure applied by the surrounding soil mass and prevents soil inrush into the excavation.The outer part of the frozen wall relative to the contour of freezing boreholes is compressed by the lateral pressure coming from the surrounding soil mass.In the outer part of the frozen wall, the absolute value of the radial displacement reduces to its minimum at the contour of the freezing boreholes.Due to redistribution of the loading throughout the thickness of the frozen wall,the inner part of the frozen wall moves into the excavation.In the inner part of the frozen wall,absolute value of the radial displacement increases and attains its maximum on the excavation sidewall.

Fig.23.Total volumetric strain εvol(a, c) and effective mean stress (MPa)(b, d) in silt after 70 d (a, b)and 210 d (c, d) of freezing(white line corresponds to the freezing front where T = Tph).

Fig.24.Variation of pore water pressure pl relative to the initial hydrostatic pressure p0 with time in silt at a point located at a distance of 1.3 m from the center of the intended mineshaft.

Fig.25.Geometry of the domain and layout of the boundary conditions.

In spite of higher overburden pressure in the sand, the radial deformation is less than that in the silt stratum.It can be explained by the higher stiffness and better creep properties of the sand in the unfrozen and frozen states compared to the silt.Nevertheless, the computed radial displacement of the excavation sidewall does not exceed the critical value of 10 cm in both layers.

In addition to the estimation of the radial displacement of the excavation sidewall,the stress-strain state of the partially frozen silt and sand layers subjected to the excavation works can be analyzed.Fig.30 shows distributions of the equivalent stress σ′eqand the radial strain εrin the silt and sand layers obtained taking into account frostheave and cryogenic suction in the silt and frost heave in the sand evolving during artificial freezing.It can be seen that the equivalent stress in the frozen zone is larger than that in the unfrozen zone due to the higher stiffness of the frozen soil.As the horizontal displacement on the boundary of the freezing borehole is constrained, the equivalent stress reaches maximum values on the boundary.The radial strain is positive in the unfrozen soil surrounding the frozen wall and in the inner part of the frozen wall bounded by the excavation sidewall and the freezing borehole.At the same time, in the outer part of the frozen wall and unfrozen soil inside the frozen wall,the radial strain is negative.The change of sign of the radial strain throughout the thickness of the frozen wall means that the boundary of the freezing borehole affects the mechanical behavior of the frozen wall.The outer part of the frozen wall supported by the rigid boundary of the freezing borehole endures the loading from the surrounding unfrozen soil and constraints its inrush into the excavation.The outer part of the frozen wall is only constrained at the upper edge by the shaft lining (Fig.25), therefore, the frozen soil moves into the excavation under the loading that is characterized by an increase in the radial displacement and contraction of the unfrozen soil inside the frozen wall.

Fig.26.Variations of absolute value of axial strain εz with time for frozen silt (a) and sand (b).Markers correspond to experimental data, and solid lines correspond to calculated results.The curves were obtained for two different values of vertical compressive loads, i.e.Pob = 1.13 NPa and Pob = 2.26 NPa.

Fig.27.The radial displacement ur(cm)in the silt stratum obtained after 24 h of the deformation process for the computations(a)with frost heave and cryogenic suction,(b)with only frost heave, and (c) without frost heave and cryogenic suction (the blue lines correspond to the frozen wall boundaries determined on isotherm where T = Tph, and the displacement scaling factor is χ = 11).

Fig.28.The radial displacement ur(cm)in the sand stratum obtained after 24 h of the deformation process for the computations(a)with and(b)without frost heave(the blue lines correspond to the frozen wall boundaries determined on isotherm where T = Tph, and the displacement scaling factor is χ = 23).

Fig.29.Profiles of the radial displacement ur along the line segment lying in the middle plane of (a) silt and (b) sand strata located at the height of 2.5 m from the excavation bottom.The presented curves correspond to the computations with frost heave and cryogenic suction (1), only frost heave(2), and without frost heave and cryogenic suction (3).

6.Conclusions

In this study, a coupled THM model of freezing of watersaturated soil is proposed.Artificial freezing of water-bearing soil layers for vertical shaft sinking,and deformation of an unsupported sidewall of mineshaft excavated in the partially frozen soil layers are investigated.The proposed model includes the water transfer equation,the heat transfer equation,and the equilibrium equation which are incorporated in COMSOL Multiphysics software to be solved numerically relative to primary variables of porosity, temperature, and displacement using the FEM.The phase transformation of water upon freezing is governed by an empirical temperature-dependent ice saturation function.Interaction between the equations is established by constitutive relations of poromechanics of freezing porous media proposed by Coussy, the effective stress of the Bishop type, and the Clausius-Clapeyron equation.Moreover, simple constitutive relations are used to describe the volumetric expansion of frozen soil due to frost heave and creep of frozen soil under long-term loading.The key feature of the proposed model is that material parameters of soil required for numerical simulation can be identified using experimental data provided by a series of standard tests usually performed beforeapplication of the AGF technique.Results of numerical simulation of laboratory freezing tests demonstrate that the proposed model is able to predict the mechanical deformation and temperature of frozen soil both in one-sided and radial freezing conditions.The applicability of the developed model to studies of AGF in geotechnical problems was illustrated by numerical simulation of freezing and shaft sinking activities in silt and sand strata at the Petrikov potash mine.To perform the simulation, material parameters of the soils were calibrated based on experimental data of laboratory freezing and creep tests conducted on the soil samples.According to the obtained results,the following conclusions can be drawn:

Fig.30.The equivalent stress (MPa) (a, c) and the radial strain εr (b, d) in the excavated silt (a, b) and sand (c, d) layers.

(1) Comparison of calculated results with existing experimental data from one-sided freezing tests for silty-clay soil samples in an open system shows that the proposed model is able to predict variations of frost heave and temperature with time up to the moment of temperature stabilization in the soil samples.If the freezing process of soil evolves with the main growth of massive cryogenic structure without thick ice lenses,then the predicted values of the frost heave are highly consistent with measurements.In the case of the formation of thick ice lenses, the calculated values of frost heave are underestimated compared to measurements.Results of numerical simulation of radial freezing test on sand samples are in good agreement with the experimental data.The model allows one to describe both contraction and expansion of the soil due to effect of frost heaving pressure evolving during the freezing test.

(2) Numerical simulation of artificial freezing of water-bearing silt and sand layers at the Petrikov potash mine allows one to conclude that the calculated variations of temperature with time at the position of thermal monitoring boreholes are close to the data of thermal monitoring carried out during the active stage of artificial freezing.Artificial freezing of the sand proceeds without water migration, thus the porosity in the frozen zone increases by 9%.The volumetric expansion in the frozen zone leads to contraction of unfrozen soil inside a closed frozen wall that entails a monotonous increase in the pore water pressure with the extension of the frozen wall.In the silt stratum during artificial freezing,water flows from the unfrozen to frozen zone due to negative pore pressure arising near the freezing front.The water migration to the freezing front significantly contributes to volumetric expansion in the frozen zone that is related to a 21% growth of porosity near the freezing borehole.At the same time, in the area adjacent to the middle plane, the porosity is reduced due to water outflow to the freezing front during the first stage of freezing when the freezing front propagates from the freezing borehole to the middle plane.The volumetric ice content in the area is 17% less than that near the freezing borehole.Free water supplement from unfrozen soil mass surrounding the frozen wall causes a rise in the porosity on the outside of the frozen wall.Near the freezing front, the soil contracts due to negative pore pressure.Significant frost heaving pressure near the freezing borehole leads to the change in stress state from initial compression to tension.Under the mechanical impact from the frozen zone, unfrozen soil inside the frozen wall is compressed.The pore water pressure inside the frozen wall quickly rises after the closure time of the frozen wall, but water migration to the inner boundary of the wall eventually causes a decrease in the pore water pressure.

(3) During shaft sinking in the silt and sand layers under the protection of the formed frozen wall, the unsupported excavation sidewall is bending and the excavation bottom is uplifted under the pressure from the surrounding unfrozen mass.If the effects of frost heave and water migration on the stress state of the silt stratum are not taken into account,then the radial displacement of the excavation sidewall is underestimated by 47%.In the sand stratum, the frost heaving pressure induced by artificial freezing during shaft sinking leads to rise of the radial displacement of the excavation sidewall by 27% during shaft sinking.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research was supported by 17-11-01204 project (Russian Science Foundation).

List of symbols

φ Angle of internal friction

uzAxial displacement

εzAxial strain

bunBiot coefficient of soil in the unfrozen state

bfrBiot coefficient of soil in the frozen state

NunBiot coefficient of soil in the unfrozen state

NfrBiot coefficient of soil in the frozen state

KunBulk modulus of soil in the unfrozen state

KfrBulk modulus of soil in the frozen state

x, y, z Coordinate

c Cohesion of soil

ρjDensity of the phase j

u Displacement vector

τ Dimensionless time

b Effective Biot coefficient of freezing soil

N Effective Biot tangent modulus of freezing soil

K Effective bulk modulus of freezing soil

G Effective shear modulus of freezing soil

p Equivalent pore pressure

X Effective mechanical parameter

εeElastic strain

σ′ Effective stress

TphFreezing temperature of pore water

g Gravity acceleration

QphHeat source due to the phase change

k Hydraulic conductivity of freezing soil

k0Hydraulic conductivity of soil in the unfrozen state

I Identity tensor

n0Initial porosity

p0Initial pressure

T0Initial temperature

L Latent heat

PobOverburden pressure

n Porosity

χ Pore water pressure coefficient

pjPressure of phase j

XunProperty of soil in the unfrozen state

XfrProperty of soil in the frozen state

urRadial displacement

δvarRelative error

SjSaturation of the phase j

χ Scale factor

ψ Soil water potential

GunShear modulus of soil in the unfrozen state

GfrShear modulus of soil in the frozen state

cjSpecific heat capacity of the phase j

T Temperature

α,β Testing constants for saturation and hydraulic conductivity

ω, m Testing constants for creep

λjThermal conductivity of the phase j

εthThermal strain

t Time

ξ Total deformation modulus

ε Total strain

σ Total stress

vlVelocity of water

εveViscoelastic strain

C Volumetric heat capacity of soil

ClVolumetric heat capacity of pore water

aTVolumetric thermal dilation coefficient

F Yield surface