Study on the influence of hydro-thermal-salt-mechanical interaction in saturated frozen sulfate saline soil based on crystallization kinetics

International Journal of Heat and Mass Transfer 146 (2020) 118868

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Study on the influence of hydro-thermal-salt-mechanical interaction in saturated frozen sulfate saline soil based on crystallization kinetics Jing Zhang a,b, Yuanming Lai a,b,⇑, Jifeng Li c, Yanhu Zhao a,b a

State Key Laboratory of Frozen Soil Engineering, Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou, Gansu 730000, China University of Chinese Academy of Sciences, Beijing 100049, China c Qinghai Water Conservancy and Hydroelectric Survey and Design Institute, Xining, Qinghai 810008, China b

a r t i c l e

i n f o

Article history: Received 29 June 2019 Received in revised form 9 October 2019 Accepted 9 October 2019

Keywords: Saturated frozen sulfate saline soil THSM interaction process Crystallization kinetics Mathematical model

a b s t r a c t The freezing of saturated saline soil is a dynamic hydro-thermal-salt-mechanical (THSM) interacting process. One-side freezing experiment of saturated sulfate saline soil in an open system with no-pressure water supplement is carried out. The coupling mechanism of water and salt migration and the soil deformation in the freezing process has been investigated by the one-side freezing experiment. Based on the crystallization kinetics theory, a hydro-thermal-salt-mechanical coupled mathematical model for saturated frozen sulfate saline soil with the effect of phase change is proposed. Moreover, the influence of solute on the physical and mechanical properties of soil during the freezing process is considered. To solve the nonlinear equations, the finite element algorithm is applied to solve the general form of governing equations. Finally, numerical simulation is implemented with the assistance of COMSOL. Validation of the model is illustrated by comparisons between the simulation and experimental results. From this study, it is found that, (1) the phenomenon of macroscopic crystallization can be well illustrated by the microscopic crystallization kinetics theory on the basis of the concepts of the water activity and the solution supersaturation. (2) The pore pressure due to the effect of phase change is main driving force for water and salt migration as well as the deformation of porous medium. It is concluded that the positive pore pressure is the main factor for soil deformation, and the negative pore pressure is the driving force for water migrates to the frozen zone. (3) Salt migrates with water and is rejected into the unfrozen water in the process of ice formation, and the rate of salt rejection gradually increases with the decrease of cooling rate. Therefore, salt crystals with layer distribution are formed in the frozen zone during the freezing process, and the largest salt crystals distribution zone is formed near the freezing front due to the effect of solute diffusion. (4) The calculated results are well agreement with the experimental data, demonstrating that the proposed hydro-thermal-salt-mechanical coupling model can well clarify the mechanism of heat and mass transfer in saturated frozen sulfate saline soil, and predict the deformation due to the effects of frost heave and salt expansion. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Sulfate saline soil is widely distributed in the cold western regions of China, where the temperature changes dramatically. The freezing process of sulfate saline soil involves the interaction of water and salt migration, heat transfer and soil deformation. Especially, the effects of frost heave and salt expansion have caused seriously damage to infrastructures such as roads, railways and other engineering facilities. With the economic development ⇑ Corresponding author at: State Key Laboratory of Frozen Soil Engineering, Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou, Gansu 730000, China. E-mail address: [email protected] (Y. Lai). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118868 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

in western China, roads are increasingly being built in saline soil. Therefore, studying the coupling mechanism of heat and mass transfer is a significant subject to prevent soil deformation due to salt expansion and frost heave of sulfate saline soil. Up to now, many researchers have been using different laboratory, analytical and numerical methods to investigate the frost heave mechanism in a freezing ground. Harlan [1] conducted an analogy between the mechanism of water transport in partially frozen soil and that in unsaturated soil, and proposed a heatfluid transport coupling model. Guymon and Luthin [2] based on an equivalent quasilinear function for the Richards equation, and proposed a one-dimensional heat-moisture transport model. Taylor and Luthin [3] modified the heat-fluid transport model presented by Harlan and presented the frost heave criterion. Jame

2

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

Nomenclature a; b; m aw ; aw;0

testing constants [–] water activity for a given salinity, and at freezing temperature [–] a
; a
;0 average ion activity for a given solution, and saturated solution [–] C p;a heat capacity of a-phase [J/(kg K)] concentration of a-phase [kg/m3] Ca c mass concentration of dissolved salt [kg/m3] Dsh hydrodynamic dispersion coefficient [m2/s] Es compression modulus [MPa] g acceleration of gravity [m/s2] Js ; Jw dissolved salt flux and water flux [kg/(m2 s)] Jq ; JT heat conduction flux and total heat flux [J/(m2 s)] Ki; Kc kinetic parameters [1/s] coefficient of temperature gradient potential [m/K] kT k hydraulic conductivity [m/s] krl relative permeability [–] intrinsic permeability and its initial value [m2] kint ; k0 Lwi ; Lsc ; Lwc latent heats of phase change [J/kg] Ma molar mass of a-phase [kg/mol] ma mass content of a-phase [kg/m3] mass content of crystal water [kg/m3] mwc m molality of the solution [mol/kg] ni ; nc kinetic parameters [–] n; n0 porosity and initial value [–] pressure of a-phase [MPa] Pa Ppor pore pressure [MPa] R ideal gas constant [J/(K mol)] Sa saturation degree of a-phase [W/m3] DSm fusion entropy of ice-water [J/(mol K)] T; T ref temperature and initial temperature [K] T0; Tf freezing temperature of pure water and saline soil [K] U a ; U a;start supersaturation ratio and its trigger value [–] u; ux ; uy ; uz displacement vector, and its component in x, y and z direction [m]

and Norum [4] simulated the phenomena of moisture and heat transport by the Harlan’s model, and considered the influence of soil freezing process on permeability characteristics. It is can be seen that the hydrodynamic models can well illustrate the mechanism of heat and mass transfer. However, the formation of discrete ice lenses is unclear. Gilpin [5] indicated that the rate of frost heave was related to the soil properties and the externally applied boundary conditions, and developed a model to predict ice lensing and frost heave. O’Neill and Miller [6] based on the fundamental thermomechanical considerations, and then proposed a rigid ice model. It has been recognized that the rigid-ice model can describe the formation mechanism of discrete ice lenses, and predict the frost heave, but it contains too many parameters to apply in practical engineering. According to the segregation potential theory, Konrad and Morgenstern [7] established a frost heave model through the concept of segregation-freezing temperature. Considering the limitations of the segregation potential theory, Nixon [8] modified the model presented by Gilpin and gave a discrete ice lenses theory for frost heave in soils. Recognizing the shortcomings of the existing models, Zhou and Zhou [9] built a heat-fluid coupling model through the concept of equivalent water pressure, and illustrated the mechanism of a single ice lens growth. There exists many heat-fluid coupling models, but all of them are not valid to solve the problem induced by stress in the process of freezing. Taking into account the variations of stress and corresponding deformation, some relevant studies have been carried out [10]. Incorporating the effect of phase change, convective heat transfer

Va molar volume of a [m3/mol] ta velocity of a-phase [m/s] trl solution relative velocity [m/s] tM ; tX ; tw stoichiometric number for positive ions M, negative ions X and molecules of water [–] Z Mþ ; Z X  charges of positive ion M and negative ion X b coefficient of thermal stress [K1] c
; c
;0 mean activity coefficient for a given solution, and saturated solution [–] cil surface tension [N/m] et ; ee ; eT total strain, elastic strain and thermal strain [m/m] gl ; g0 dynamic viscosity of solution and its initial value [Pa s] ha volumetric content of a-phase [m3/m3] r a hw ; hw residual unfrozen water and initial moisture content [m3/m3] ka thermal conductivity of a-phase [W/(m K)] l chemical potential [J/mol] qa density of a-phase [kg/m3] r; r0 total stress and effective stress [N/m2] / osmotic coefficient [–] va ratio of a-phase in solution [–] x mass fraction [–] Subscripts w unfrozen water s dissolved salt m material matrix i ice c salt crystal l solution e effective value 0 initial value

and the development of ice lenses, a coupled thermo-hydromechanical model was submitted [11]. According to the principle of effective stress, Nishimura et al. [12] put forward a coupled thermo-hydro-mechanical model to characterize the interaction in saturated soil during the process of freeze-thaw cycle. Combined with the Clapeyron equation, Lai et al. [13] introduced a hydrothermal-mechanical coupling model for the frost heave, and suggested that the negative pore water pressure was the main driving force for water migrates from the unfrozen zone to the frost front. With the consideration of the interaction of heat and mass transfer as well as stress development in the freezing soil, Zhou et al. [14] advanced a hydro-thermal-mechanical coupling model by extending the coupled heat-fluid model presented by Zhou and Zhou [9]. With the development of volume averaging theory [15], more and more coupling models have been propounded using the innovative approach. Tong et al. [16] combined with the volume averaging theory, and gave a numerical method for modeling coupled thermo-hydro-mechanical process of geomaterials with multiphase fluid flow. Koniorczyk [17] established a mathematical model for heat and water transport of deformable materials considering the kinetics of water phase change. Na and Sun [18] applied the volume averaging theory, built a stabilized thermohydro-mechanical finite element model to investigate the freezethaw action of frozen porous media. All of these abovementioned thermo-hydro-mechanical coupled models, especially the models established through the volume averaging theory, can well simulate the interaction of water, heat and stress in the

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

freezing process. However, neither of them takes into account the effect of solute, and cannot describe the phenomenon of salt migration and crystallization. Generally salt exists in two phases in porous medium: dissolved salt and salt crystals. It has been recognized that salt in porous medium will reduce the service life of buildings [19]. Most of the pores are blocked by the precipitated salt, which changes the internal structure and reduces the hydraulic conductivity [20,21]. In addition, the influence of salt on freezing temperature and unfrozen water content is obvious [22]. It is found that the crystallization pressure has a great destructive effect on the porous materials, which is related to the supersaturation of the solution [23]. Salt migrates with water, and the crystallization and accumulation of salt have significant effects on heat and mass transfer. At the same time, the deformation of soil is caused by salt expansion. According to experimental and theoretical studies, many of the mathematical and numerical models have been established to illustrate the mechanism of water and salt transfer, and to predict the soil deformation. It has been discovered that dissolved salt would be rejected into the unfrozen water during the freezing process. Moreover, a larger salt crystals distribution zone will be formed near the freezing fringe due to the diffusion and migration of salt [24–26]. Panday and Corapcioglu [26] propounded a hydrothermal-salt mathematical model to elucidate the phenomenon of salt rejection in frozen saline soil. Pavlik et al. [27] advanced a coupled water and salt transport model including advection and diffusion effects. However, the phase change of salt is neglected. Nguyen et al. [28] did not take into account the heat flux while describing salt transport. As is well-known, substantial amount of heat will be released and consumed during the processes of crystallization and dissolution, respectively. Therefore neglecting heat transfer might be erroneous. In view of the heat flux during salt crystallization, Koniorczyk [19] proposed a coupled heat, moisture and salt transport model by employing the kinetics of salt phase change, but ignored the impact of crystallization pressure on the stress field. Castellazzi et al. [29] extended the coupled heat, moisture and salt model presented by Castellazzi et al. [30] including the crystallization/dissolution and the hydration/dehydration of salts with N crystallized forms, like sulphates. However, the extended model is only valid in the range of positive temperature. Zhang et al. [31] introduced solute conservation equation into the traditional three-field coupling model, and gave a hydro-thermalsalt-mechanical mathematical model. However, the effect of salt on the moisture field and the soil deformation is neglected. According to the above-mentioned analysis results, scholars have preliminarily explored the coupling mechanism of saline soil, and achieved certain achievements. The disadvantage is that, in many cases, the moisture and temperature fields in the hydro-thermalsalt-mechanical mathematical coupling model are still dominated by that of the traditional three-field coupling model, ignoring the influence of salt on the heat and mass transfer. Moreover, the effects of salt crystallization and accumulation are rarely considered in the solute transport equation. Generally speaking, the study on the interaction of water and salt migration, heat transfer, salt crystallization and accumulation and soil deformation during the freezing process of sulfate saline soil is insufficient. Therefore, providing a more complete hydro-thermal-salt-mechanical coupled mathematical model for saturated frozen sulfate saline soil is extremely significant. In this paper, one-side freezing experiment of saturated sulfate saline soil in an open system with no-pressure water supplement is implemented. Based on experimental results and theoretical analyses, a coupled hydro-thermal-salt-mechanical mathematical model for saturated frozen sulfate saline soil considering the effect of phase change is presented. The governing equations consist of

3

mass balance equation, energy balance equation and stress balance equation. In this model, the crystallization kinetics models of ice and salt crystals growth were presented respectively by the activity of water and the supersaturation of solution. The generalized Clapeyron equation is adopted to solve ice pressure. Additionally, the influence of salt on the physical and mechanical properties of soil during the freezing process is considered. To solve the mathematical model, finite algorithm is introduced to the general form of equations. Finally, numerical simulation of the experimental soil column is performed to validate this model. Comparisons between the simulated and experimental results illustrate that the model simulating results match well with the test results. 2. Crystallization kinetics 2.1. Ice-water phase change kinetics For ice-water phase change, Bronfenbrener and Korin [32] presented a kinetic model of solidification in a fine-grained porous medium. The rate of ice-water phase change is clarified by the crystallization kinetics model. Koop et al. [33] concluded that the water activity was the driving force for homogeneous ice nucleation, and established the water-activity-based kinetics. In this paper, the water activity is the dominant factor for ice crystallizai tion. Therefore, the phase change rate @m of water freezing and @t thawing is given by [34]:

@mi ¼ @t

(

C w K i ð1  aw Þni

if aw < 1 and C w > 0;

C i K i ðaw  1Þni

if aw  1 and C i > 0;

ð1Þ

where K i and ni are kinetic parameters. C w and C i are the concentrations of water and ice, respectively. aw is the water activity at arbitrary temperature, and aw < 1 is the threshold value of ice crystal growth. It can be seen the first term in Eq. (1) indicates the process of ice crystal growth (aw < 1 and C w > 0), while the second term in Eq. (1) refers to the melting process of ice crystals (aw  1 and C i > 0). 2.2. Solution-crystal phase change kinetics Given a salt of composition M tM X tX  tw H2 O consisting of tM positive ions M of charge Z Mþ , tX negative ions X of charge Z X  and tw molecules of water, and the dissolution reaction is defined as:

MtM X tX  tw H2 O ¡

tM MZMþ þ tX X ZX þ tw H2 O

ð2Þ

where tM , tX and tw are stoichiometric number for positive ions M, negative ions X and molecules of water, respectively. Espinosa et al. [35] assumed that the driving force for salt crystallization was the difference between the chemical potentials of solution and salt crystal, Dl

Dl ¼ ln RT



a
a
;0

v  v  aw w ¼ lnU a  aw;0

ð3Þ

where a
and a
;0 are the average ion activities of given solution and saturated solution, respectively. t ¼ tM þ tX is the total ions mole number per mole solution. aw and aw;0 are the water activities at any temperature and equilibrium state, respectively. U a is the supersaturation ratio of solution, R is gas constant and T is the temperature. The supersaturation is the driving factor for salt crystallization. Based on the crystallization kinetics equation established by Espic nosa et al. [35], the phase change rate @m for crystallization and @t dissolution can be expressed as:

4

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

@mc ¼ @t

(

C s K c ðU a  U a;start Þnc

if U a  U a;start and C s > 0;

C c K c ð1  U a Þnc

if U a < 1 and C c > 0;

ð4Þ

where K c and nc are kinetic parameters. C s and C c are the concentrations of solute and salt crystal, respectively. U a;start  1 is the threshold value of salt crystal growth. Usually U a;start  1 for primary crystallization but U a;start ¼ 1 once the crystallization process has begun. It is noted that the first term in Eq. (4) denotes the process of salt crystallization (U a  U a;start and C s > 0), whereas the second term in Eq. (4) represents the dissolution process of salt crystals (U a < 1 and C c > 0). 3. Mathematical model Frozen saline soil is composed of solid, liquid and gas phases, which is a multi-phase continuous porous medium. When the soil is saturated, only the interaction between solid and liquid phases is considered, and the influence of gas phase is neglected. As is known that sodium sulfate is more destructive than other salts. Therefore, only one type of salt and salt crystal is considered, i.e., mirabilite (Na2SO410H2O). It is assumed that the liquid phase is an ideal solution consisted of unfrozen water and dissolved salt. Additionally, the solid phase consists of soil particles, ice and salt crystals. Moreover, the following assumptions are made to establish the model: (1) The soil is an isotropic and homogeneous porous medium. (2) The change of mass flux caused by soil deformation is considered. (3) Soil particles are incompressible and porosity changes can lead to soil deformation. (4) Moisture and salt migration are dominated by liquid form, without considering the movement of ice and salt crystals. In the following, unfrozen water and dissolved salt are described by subscripts w and s, respectively. Soil particles, ice and salt crystals are denoted by subscripts m, i and c, respectively. The total porosity, defined as the volume of voids per unit volume of porous medium, is defined by n. The content of each component can be described by the concentrationC a , defined as the mass of a in per unit volume of porous medium, or by the corresponding saturation degreeSa , defined as the pore volume occupied bya:

C a ¼ ðnSa Þqa ¼ ha qa

ð5Þ

where qa is the mass density of a, and ha ¼ nSa is the volume of a in per unit volume of porous medium. The relationship of the degrees of saturation all the components can be represented as:

Sw þ Ss þ Si þ Sc ¼ 1

ð6Þ

Additionally, the volumes of all the components satisfy the relation:

hm þ hw þ hs þ hi þ hc ¼ 1

ð7Þ

For saturated soil, the liquid in the porous medium is the mixture of water and dissolved salt. Thus, the saturation of solution is formulated as follows:

S l ¼ Sw þ S s

ð8Þ

The ratios of dissolved salt and moisture in the solution can be written as:

S Sl

vs ¼ s ; vw ¼

Sw Sl

ð9Þ

The concentration of salt in the solution is defined as:

c ¼ vs q s

ð10Þ

Based on the volume averaging theory, the coupled model, including mass balance equation, energy balance equation and stress balance equation, was presented. At the same time, the effects of ice and salt crystallization on porous media are considered by employing the crystallization kinetics approach. 3.1. Mass conservation equation 3.1.1. Moisture conservation equation Based on the concept of equivalent water content, it is suggested that moisture in pores consist of unfrozen water, ice and crystal water during the freezing process of saturated sulfate saline soil. Therefore, neglecting the inertial and viscous effects [36], the moisture conservation equation can be established from the equation of linear momentum balance for the liquid and the mass balance equation of water, which can be expressed as [18,19]:

@ ðmw þ mi þ mwc Þ þ rðnSw qw tl Þ ¼ 0 @t

ð11Þ

where mw ¼ nSw qw dV, mi ¼ nSi qi dV and mwc ¼ tw MMwc nSc qc dV are mass contents of unfrozen water, ice crystals and crystal water in per unit volume of porous medium, respectively. tl is the absolute velocity of liquid solution. Based on Darcy’s law, the equation of linear momentum balance for the liquid is given by:

trl ¼ kru

ð12Þ

where trl is the relative velocity of liquid solution, k denotes the hydraulic conductivity, and u is the total soil-water potential. Moreover, the relative velocity of solution with respect to a moving solid is given by [36]:

trl ¼ nSl ðtl  tm Þ

ð13Þ

where tm ¼ @u is the solid skeleton velocity. @t Substituting Eqs. (12) and (13) into Eq. (11), the moisture balance equation can be written as:

    tr @ Mw nSw qw þ nSi qi þ tw nSc qc þ r nSw qw l @t Mc nSl   @u ¼0 þ r nSw qw @t

ð14Þ

where M w and Mc are the molar mass of water and mirabilite, respectively. u denotes displacement vector. Based on the crystallization kinetics theory, Eq. (14) is further simplified as:

    tr @Sw @n @u þ r nSw qw l þ r nSw qw þ Sw qw @t @t @t nSl @mi @mwc  ¼ @t @t

nqw

where

@mi @t

¼ @t@ ðnSi qi Þ and

@mwc @t

¼ @t@



ð15Þ



tw MMwc nSc qc are the crystalliza-

tion rate of ice and the formation rate of crystal water, respectively. Lai et al. [13] concluded that the total soil-water potential in freezing soil consists of the matric potential up ¼ qPwg, the gravity w

ug ¼ z and the temperature gradient potential uT ¼ kT DT. It has been found that the osmotic potential us ¼ MRTc s ql g potential

induced by the presence of dissolved salt cannot be neglected [20]. Therefore, the total soil-water potential of the moisture migration in frozen saline soil is defined as:

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

u ¼ up þ uT þ ug þ us ¼

Pw

qw g

þ z þ kT DT þ

RTc Ms ql g

ð16Þ

where Pw is pore water pressure, kT is the coefficient of the temperature gradient potential, Ms is solute molar mass, ql ¼ vs qs þ vw qw denotes the equivalent density of liquid solution. The traditional Clapeyron equation is not valid in frozen saline soil due to the presence of solute. Loch [37] extended the general Clapeyron equation by introducing the osmotic pressure P s . In addition, it is assumed that the liquid solution is an ideal dilute solution whose density is approximately equal to that of pure water. The extension of Clapeyron equation is given as [26]:

RTc Pi T  T0   ¼ Lwi qw Ms qw qi T0 Pl

ð17Þ

where Pl , Pi are the pressures of liquid solution and ice, respectively. Lwi is the latent heat of ice-water phase change, T 0 is the freezing temperature of water. Eq. (17) describes the thermodynamic equilibrium relationship among ice pressure (Pi ), osmotic pressure (P s ), temperature (T) and liquid solution pressure (P l ), and reflects that temperature and liquid solution pressure is dependent. Thus, it is suggested that kT ¼ 0 [13]. The relative seepage velocity of liquid solution can be given by:

trl ¼ 

krl kint

gl

ðrpl þ ql grzÞ

ð18Þ

where kint denotes the intrinsic permeability of soil. krl and gl represent the relative permeability and the dynamic viscosity of liquid solution, respectively. 3.1.2. Solute conservation equation Salt exists in two phases in porous medium: dissolved salt and salt crystals. The hygrothermal properties of porous materials are seriously affected by salt migration and crystallization [29]. Based on the crystallization kinetics approach, the solute conservation equation, including the effect of salt migration and crystallization, can be given as:

@ ðms þ mc Þ þ rðnSl J S Þ ¼ 0 @t

ð19Þ

where ms ¼ nSs qs dV and mc ¼ nSc qc dV are mass contents of dissolved salt and salt crystals in per unit volume of porous medium, respectively. J S denotes the flux of solute. Salt migrates through advection and hydrodynamic dispersion. The flux of salt is defined as:

J S ¼ ctl  Dsh rc

ð20Þ

where Dsh is the coefficient of hydrodynamic dispersion. Substituting Eqs. (12), (13) and (20) into Eq. (19), the solute balance equation can be expressed as:

    tr @Ss @n @u þ r nSs qs l þ r nSs qs þ Ss qs @t @t @t nSl @mc  r½Dsh rðnSs qs Þ ¼  @t

nqs

where

@mc @t

5

3.1.4. Soil particles mass balance During the freezing process, ice and salt crystals in porous media will lead to the increase of porosity and micro cracks. The crystallization of ice and salt will change the inner structure of pore and cause the deformation of soil skeleton [26]. Assuming that the density of soil particle remains unchanged, the mass conservation equation of soil particles is given as:

  @ @u ¼0 ½qm ð1  nÞ þ r qm ð1  nÞ @t @t

ð24Þ

3.2. Heat transfer equation The energy conservation equation, considering the latent heat released by ice and salt crystallization and the variation of energy caused by the mass transfer of each component, has the following form [29]:



qC p

@T @mi @mc @mwc þ rJ T ¼ Lwi þ Lsc þ Lwc e @t @t @t @t

ð25Þ

 where qC p e is the effective specific heat capacity; Lsc andLwc are the latent heats of solution-salt crystal and free water-crystal water phase change, respectively. J T is total heat flux, which can be defined as:

J T ¼ J q þ J w C p;w T þ J s C p;s T

ð26Þ

where C p;a denotes the specific heat capacity of a. Jq and Jw represent the heat flux and the water flux, respectively.

J q ¼ ke rT J w ¼ nSw qw

ð27Þ 

trl nSl

þ

 @u @t

ð28Þ

where ke is the effective thermal conductivity. Taking into account the influence of the content variation of each component on the effective specific heat capacity and the effective thermal conductivity, the parameters mentioned above can be described as [13]:





qC p e ¼ nSw qw C p;w þ nSs qs C p;s þ nSi qi C p;i þ nSc qc C p;c þ ð1  nÞqm C p;m

ke ¼ ðkw Þhw ðks Þhs ðki Þhi ðkc Þhc ðkm Þð1nÞ

ð29Þ ð30Þ

where ka denotes the thermal conductivity of a component: Substituting Eqs. (20), (27) and (28) into Eq. (25), the energy conservation equation can be given as:



@T   rðke rT Þ þ Jw C p;w þ J s C p;s rT e @t @mi @mc @mwc þ Lsc þ Lwc ¼ Lwi @t @t @t

qC p

ð31Þ

ð21Þ 3.3. Stress balance equation

¼ @t@ ðnSc qc Þ denotes the growth rate of salt crystals.

3.1.3. Ice and salt crystals conservation equations Based on the hypotheses mentioned above, the movement of ice and salt crystals is neglected. The mass balance equations of ice and salt crystals can be described as [34]:

  @Si @n @u @mi nqi þ r nSi qi ¼ þ Si qi @t @t @t @t

ð22Þ

  @Sc @n @u @mc ¼ þ r nSc qc þ Sc q c nqc @t @t @t @t

ð23Þ

The stress equilibrium equation based on the effective stress principle can be expressed as [17]:

rr þ c ¼ 0

ð32Þ

where r is the total stress, andc is unit weight. The unit weight varies with the change of each component in soil, which can be obtained as:

c ¼ qe g ¼ ½nSw qw þ nSs qs þ nSi qi þ nSc qc þ ð1  nÞqm g

ð33Þ

where qe is the equivalent density of soil, which is related to the density of each component.

6

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

Based on the effective stress theory, the total stress is defined as:

r ¼ r0 þ Ppor

ð34Þ

where r0 represents the effective stress, and P por denotes the pore pressure. The pore pressure consists of pore water pressure (Pw ), pore ice pressure (Pi ), osmosis pressure (P s ) and salt crystal pressure (P c ). Thus, the pore pressure can be written as:

Ppor ¼ Sw Pw þ Ss Ps þ Si P i þ Sc P c

ð35Þ

The total strain et for freezing soil, including elastic strain e and thermal strain eT , can be written as: e

et ¼ e þ e e

ð36Þ

T

The relationship between strain and displacement is defined as:

et ¼

@ux @uy @uz þ þ @x @y @z

ð37Þ

where ux , uy and uz represent displacements in x, y and z directions, respectively. The elastic stress-strain relationship of freezing soil can be given as:

dr0 ¼ Es dee

ð38Þ

where Es denotes the compressive modulus. Lai et al. [38] studied the relationship between the compressive modulus and temperature, and proposed the following empirical formula:

Es ¼ a þ bjTjm

ð39Þ

where a and b are test parameters. For the silty clay, a ¼ 2:8 MPa and b ¼ 2:6. Especially, when temperature of soil is above 0 °C, then b ¼ 0, and m is a nonlinear index which is smaller than 1, and in this study m ¼ 0:6. The thermal strain induced by the temperature variation is commonly given by



eT ¼ b T  T ref



ð40Þ

where b is the linear thermal dilatation coefficient, T ref represents the reference temperature, which is belong to the initial temperature of soil. Substituting Eqs. (33) and (34) into Eq. (32), the stress balance equation can be described as:  r r0 þ Ppor þ ½nSw qw þ nSs qs þ nSi qi þ nSc qc þ ð1  nÞqm g ¼ 0 ð41Þ

tionship between the intrinsic permeability and the effective porosity can be represented as [40]:

kint ¼ k0

ð1  n0 Þ2 n0 3 ð1  ne Þ ne 3

ð43Þ

2

where k0 and n0 are the initial values of intrinsic permeability and porosity, respectively. Ophori [41] studied the influence of temperature and salt concentration on the dynamic viscosity of solution. Koniorczyk [42] summarized the research results, and proposed the dynamic viscosity equation as follows:



gl ¼ g0 1 þ A1 x þ A2 x2 þ A3 x3



ð44Þ

where g0 denotes the dynamic viscosity of pure water at a given temperature and A1 ¼ 1:85, A2 ¼ 4:1, A3 ¼ 44:5 are material coefficients. x ¼ C s =ðC w þ C s Þ is the mass fraction. 4.2. The solubility of sodium sulfate The solubility of sodium sulfate depends on temperature. Koniorczyk [42] combined with the previous experimental data, and built a fitting formula to describe the relationship between the solubility of sodium sulfate and temperature. However, the expression can only illuminate salt crystallization over the range of positive temperature. According to the crystallization kinetics, the solubility of sodium sulfate is related to the supersaturation. Marion and Farren [43] modified the Pitzer-equation parameters in the SMW model presented by Spencer et al. [44], and proposed the FREZCHEM model to illuminate the solubility characteristics of sulfate minerals over the temperature range from 60 °C to 25 °C, i.e., sodium sulfate (Na2SO4). Basing on the FREZCHEM model, we chose the threshold value of U a as the criterion to judge the saturation state of the solution at a specific temperature. In addition, it is assumed that when U a equals 1, the solution is considered to be saturated. The relationship among the supersaturation of solution, the mean activity coefficient (c
), the molality of solution (m) and the water activity (aw ) can be written as [45]:



a
Ua ¼ a
;0

v 

aw aw;0

v w

c
¼ c
;0

!v 

m m0

v 

aw aw;0

v w ð45Þ

4. Constitutive relations

where c
and c
;0 are the mean activity coefficients at any temperature and equilibrium state, respectively. m and m0 are the molalities of solution at the given temperature and equilibrium state, respectively. The water activity related to the osmotic coefficient / can be given by [46]

4.1. Moisture migration

lnaw ¼ /tm

It is well-known that ice and salt crystals have adverse influences on the hydraulic properties of porous medium. The growth of ice and salt crystals in porous medium may partially or completely block the pore, leading to change in the internal structure of porous materials, and ultimately to a decrease in the effective porosity and permeability [21]. Koniorczyk and Gawin [39] gave the relationship between salt crystal content and the effective porosity. Wu et al. [34] found that the effective porosity depends on ice crystal content in the porous material. Based on the above researches, the effective porosity can be defined as:

According to the FREZCHEM model, the variation of molality of Na2SO4 solution over the temperature range 12 °C to 25 °C is calculated. From Fig. 1, the calculated results are in agreement with the experimental data [47], illustrating that the above method can be applied to calculate the change of molality of Na2SO4 solution at negative temperature.

ne ¼ nð1  Si  Sc Þ

ð42Þ

where ne denotes the effective porosity. The clogging effect induced by ice and salt crystals seriously reduces the seepage properties of porous medium [20]. The rela-

Mw 1000

ð46Þ

4.3. Crystallization pressure Salt crystallization in pores is one of the main factors leading to the destruction of porous materials, which is closely related to the durability of underground structures, roads and other buildings. Additionally, a better understanding of the crystallization pressure can clarify the salt expansion phenomenon on the surface of finegrained materials. The stress generation is related to the supersat-

7

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

vature radius, and h is the contact angle. It can be seen that the second term in Eq. (50) indicates the freezing temperature decrease with the increase of the molality of solution in porous medium, and the third term in Eq. (50) represents the freezing temperature decrease with the decrease of pore size. The variation of freezing temperature of saline soil with the molality of Na2SO4 solution was calculated by Eq. (50). From Fig. 2, the agreement between the calculated results and the experimental data [49] is seen to be good. 4.5. Unfrozen water content As a result of the influence of specific surface area and salinity, there is still a certain amount of unfrozen water in the pores below freezing temperature. Michalowski and Zhu [50] developed a three-parameter model to illuminate the variation of unfrozen water content during the freezing process. The relationship between unfrozen water content and temperature is as follows:

 hw ðT Þ ¼ hrw þ haw  hrw eC 1 ðTT f Þ

Fig. 1. The molality of sodium sulfate solution versus temperature.

uration at the interface between the solution and salt crystal when the latter is about to be confined between the pore walls. The crystallization pressure equation is given by [23]:

Pc ¼

tRT Vc

m c ln þ ln
m0 c
;0

!

ð47Þ

where V c denotes the molar volume of salt crystal. Neglecting the non-ideal behavior of the solution, it is assumed that the mean activity coefficient c
is equal to 1 [45]. Thus, Eq. (47) can further be simplified as:

tRT

m Pc ¼ ln Vc m0

ð48Þ

Based on the microscopic distribution of crystallization pressure, Castellazzi et al. [30] presented a simplified equation to calculate the macroscopic tensile stress through the pore model. It is assumed that the stress is related to the content of salt crystal. Thus, the macroscopic tensile stress is mainly resulted from ice pressure and salt crystal pressure. Additionally, the influence of pore water pressure and osmotic pressure on the stress is neglected, because the content of unfrozen water in frozen zone is very small. The macro tensile stress can be calculated by:

r 0 ¼ Si P i þ S c P c where

ð49Þ

r0 signifies the macro tensile stress.

hrw

ð51Þ haw

where is the residual unfrozen water, is the initial moisture content in unfrozen soil and C 1 is experimental parameter reflecting the rate of decay. Considering the influence of freezing degree, Wu et al. [51] extended the above model to elucidate the variation of unfrozen water content in non-equilibrium state. However, the extension model cannot clarify the change of unfrozen water content caused by salt crystallization before ice-water phase change. Consequently, the extension model of unfrozen water content in this paper can be presented as:

8 Tf < T  Tc > < n  hc ; hw ¼ n  hc  hi ; ðn  hc  hi Þ > hw ðT Þ > : hw ðT Þ; ðn  hc  hi Þ  hw ðT Þ

ð52Þ

where hc and hi are the volumes of salt and ice crystals, respectively. T c is the temperature of solution-salt crystal phase change. In this paper, hrw ¼ 0:001, haw ¼ 0:209 and C 1 ¼ 0:16oC1 are calculation parameters. 5. Model implement 5.1. Finite element algorithm The chosen primary variables of the developed model are the solution pressure pl , saturation degrees of dissolved salt Ss , ice crystal Si and salt crystal Sc , porosity n, temperature T and soil dis-

4.4. The freezing temperature of saline soil It is well-known that the freezing temperature of saline soil is lower than that of bulk pure water due to the physicochemical effect. Sufficient experimental data have shown that the freezing temperature depression in saline soil is related to pore size, unfrozen water content and solute content [22,48]. Based on the phase transition theory, Xiao et al. [48] established an equation to calculate the freezing point of saline soil, including the main influencing factors of pore size and solute content. The equation can be written as:

Tf ¼ T0 þ

RT 0 lnaw 2cil cosh þ r DSm V i DSm

ð50Þ

where T 0 and T f are the freezing temperatures of pure water and saline soil, respectively. V i is molar volume of ice crystal, DSm denotes the entropy of per molar volume during phase change, cil is the surface energy of ice-liquid interface, while r is its mean cur-

Fig. 2. Relationships between freezing temperature and sodium sulfate content.

8

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

placement u. The governing Eqs. (15), (21)–(24), (31) and (41) are highly nonlinear because the coefficients in these equations vary with time and they will affect each other. These equations are solved with the assistance of COMSOL. The above-mentioned equations can be written as:

da

@n þrC¼f @t

ð53Þ

where n indicates variable in governing equations, n ¼ fpl ; Ss ; Si ; Sc ; n; T; ug. da denotes damping coefficient. C is flux vector, and f is source term. Generalized solution to Eq. (53) could be given in the form

      @n  da ; n þ r  C; n ¼ f ; n X X @t X

ð54Þ



where n is virtual displacement. Eq. (54) could be written as:

        @n  da ; n  C; r n ¼ f ; n þ n  C; n X X @X @t X

ð55Þ

Using backward difference, the time terms are discretized by

@n ni  ni1 ¼ @t Dt

ð56Þ

where the superscript i refers time i, i  1 for time i  1. X is the calculated domain. @ X is the boundary of the calculating domain, and n is the outward normal unit vector of boundary. Substituting Eq. (56) into Eq. (55) gives

ni  ni1  da ;n Dt

!

       C; r n ¼ f ; n þ n  C; n X

X

X

ð57Þ

@X

Then, Eq. (57) could be rewritten as

          da ni ; n  Dt C; r n ¼ Dt f ; n þ da ni1 ; n þ Dt n  C;n X

X

X

X

@X

ð58Þ To stabilize the variable n, the following scheme is used:

n¼

ni þ ni1 2

ð59Þ

Finally, Eq. (59) could be rewritten as:

        2 da ni ; n  Dt C; r n ¼ 2Dt f ; n þ 2 da ni1 ; n X X X   X    Dt da ni1 ; n þ 2Dt n  C; n ð60Þ X

@X

The following form is defined in the above equation

Z

ðf 1 ; f 2 ÞX ¼

X

ðf 1 ; f 2 ÞdX

l is the matrix of Lagrange multiplier. Especially, P ¼ 0 denotes the Neumann boundary condition. 6. Experimental and numerical validation Based on the experimental data of the freezing process of saturated sulfate saline soil, the numerical simulation of the freezing process is carried out by using the hydro-thermal-saltmechanical coupled model developed in this paper. Validation of the model is illuminated by the comparisons between the calculated results and experimental data. 6.1. Test and calculation conditions One-side freezing experiment for saturated sulfate saline soil is carried out in an open system, and the validity of the proposed model is illuminated by comparing with the experimental data. Experiment is carried out in the freezing and thawing test chamber of XT5405B series in State Key Laboratory of Frozen Soil Engineering. Operating temperature of the test chamber ranges from 35 °C to +35 °C. The desalinized silty clay is taken from Qinghai-Tibet Plateau. The accumulative curve of particle size grading of soil sample is shown in Fig. 3. In the process of sample preparation, 1.0% content of sodium sulfate was added to the desalinized silty clay to make the soil column with 20.00 cm in diameter and 21.00 cm in height. After the vacuum saturation process, the water content closes to 20.9%. The one-side freezing experimental device for testing the soil column is shown in Fig. 4a. It consists of upper and lower cooling plates, whose temperatures can be controlled by a separate cooling circuit. The sample container is positioned inside the freezing and thawing test chamber to control environmental temperature. Eleven temperature sensors (thermistor, precision: ±0.05 °C) were assigned with interval 2.0 cm and 0.5 cm above the bottom of the soil sample to measure the temperature variations. Moreover, a displacement sensor (precision: ±0.01 mm) was installed at the top of upper plate to measure the displacement of the soil column. Insulation material with the thickness of 2.0 cm was used around the sample container to prevent the impact of outside temperature. The water supplement system consists of a Mariotee flash connected with the bottom plate through plastic tube. Additionally no-pressure water supplement process is performed. From Fig. 4a, the top of the sample is cold end subjected to negative temperature and the bottom is warm end subjected to positive temperature. Data Taker 80 was employed to collect the test data such as temperature, displacement with 5-min time interval. At the end of

ð61Þ

where f 1 and f 2 are two factors of each term in equation, such as 

f 1 ¼ da ni and f 2 ¼ n in the first term at the left of Eq. (60). 5.2. Initial and boundary conditions The initial condition can be given as [51]:

nðx; y; z; 0Þ ¼ n0

ð62Þ

where n0 is a known function of time or a prescribed value of variable n at the initial time. Boundary conditions in above model, including the Dirichlet and Neumann boundary conditions, can be expressed as [13]:

n  C ¼ G þ ð@ P=@nÞTl where n is the outward normal unit vector of boundary, G is source term at boundary, ð@ P=@nÞT is a matrix of the elastic constraint, and

Fig. 3. Accumulative curve of particle size grading of soil sample.

9

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

Fig. 4. Experimental device and finite element model. The symbols in (a) respectively represent a Displacement sensor; b and i Coolant inlet and outlet; c and h Upper and lower plates, respectively; d Sample container; e Insulation material; f Soil sample; g Temperature sensor; j Freezing and thawing test chamber; k Solution supply system.

the experiment, soil column was divided into 21 parts with a thickness of 1 cm, and the moisture and salt contents in the soil column were measured. An axisymmetric calculation model is shown in Fig. 4b, in which AB is symmetry axis. In the process of numerical calculation, AD and CD are impervious boundary conditions, and sodium sulfate solution is supplied through the permeable boundary of BC. The radial displacements of CD and BC are constrained, and the vertical displacement of BC is set to 0. At the beginning of the test, the temperature of the sample is uniformly distributed, Tjt¼0 ¼ 25 °C. CD is thermal isolation boundary, and the temperatures of AD and BC remain unchanged after decreasing from initial values to predetermined temperatures. TjAD ¼ 11:2 °C, TjBC ¼ 4 °C are used for numerical calculation. The main calculation parameters are shown in Table 1. 6.2. Validation of numerical result Three cases of grid system, named as Grid 1, 2 and 3, have been considered for the mesh independence test, respectively. The central temperature of points A and B as well as the displacement of point A are selected as references to evaluate the precision of the  and D , respectively. The values of adopted grids, named as T AB A

Table 1 Parameters of soil column for numerical simulation.

qw qs qi qc qm C p;w C p;s C p;i C p;c C p;m kw ks ki kc km Lwi Lwc Lsc b Ki Kc n0

m

k0

Value

Units

Source

1000 2700 917 1460 2700 4180 891 2090 1743 850 0.58 0.14 2.22 1.35 1.5 6.01 73.04 2.34 2  106 5.82  103 1.80  103 0.357 0.33 5.06  1019

kg/m3 kg/m3 kg/m3 kg/m3 kg/m3 J kg1 K1 J kg1 K1 J kg1 K1 J kg1 K1 J kg1 K1 W m1 K1 W m1 K1 W m1 K1 W m1 K1 W m1 K1 kJ/mol kJ/mol kJ/mol K1 s1 s1 – – m2

Ref. [13] Ref. [29] Ref. [13] Ref. [29] Ref. [51] Ref. [34] Ref. [29] Ref. [34] Ref. [29] Ref. [34] Ref. [13] Literature Ref. [13] Ref. [51] Ref. [13] Ref. [13] Ref. [51] Ref. [51] Ref. [13] Ref. [34] Ref. [51] Calculated Ref. [13] Ref. [34]

 and D T AB A are given in Table 2. The percentage changes of the

quantities are written inside the bracket. It is observed that the maximum change of 3.61% occurs in the value of DA as the mesh changes from G1 to G2, whereas no significant variations in the calculated results is observed when the mesh is changed from G2 into G3. To reduce the computational time, according to the test results it is found that the Grid 2 can well satisfy the calculated precision. A whole finite element mesh with 840 elements and 903 nodes, in the Grid 2, is shown in Fig. 4b. 6.3. Results and analyses 6.3.1. Temperature The transient temperature variations at the different positions of the soil sample are shown in Fig. 5. It can be seen from Fig. 5 that although there are obvious differences during the initial stage, the distribution of calculated temperatures is consistent with that of the experimental data. Moreover, temperature change rate varies with position, indicating that the closer to the cold end, the faster

the cooling rate. As mentioned above, the heat transfer in freezing soil depends on many factors, such as conduction and convection heat transfer, latent heat released by phase change, and heat sources and sinks. It is well known that the whole temperature variation process can be divided into three stages: fast cooling, slow cooling and steady stage. In addition, the temperature near the freezing front reaches the steady stage at the last. The temperature distributions along the height of soil column at different time are shown in Fig. 6. As shown in Fig. 6, the calculated temperatures along the column have obvious multi-linear distribution characteristics at the early stage. It is well acknowledged that the thermodynamic parameters of saline soil are related to the content of each component in the freezing process. Therefore, the slopes of these lines are different near the phase change temperature since latent heat changes the temperature distribution of the phase change zone. At the end of the experiment, the temperature

10

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

Table 2  , D T AB A values of the soil column for the grid independence test. Grid system

No. of elements

No. of nodes

T

G1 G2 G3

576 840 1200

627 903 1275

4.321 4.458 (3.23%) 4.464 (0.13%)



AB

/°C

DA /mm 15.079 15.623 (3.61%) 15.631 (0.05%)

6.3.2. Moisture migration and ice crystallization The comparison of the calculated result for the total water content (equivalent to the contents of unfrozen water, ice crystals and

crystal water) with the experimental data is shown in Fig. 7. As shown in Fig. 7, the total water content increased significantly above the freezing front while slightly decreased below it. It is seen that the redistribution of total water content is the result of the migration of unfrozen water from the unfrozen zone to the freezing fringe. In the early stage, the freezing front advances downwards rapidly, owing to a high thermal gradient near the cold end, and there is not sufficient time for water migration. Therefore, the total water content near the cold end is smaller than that near the freezing front. Furthermore, as the freezing front propagation rate decreases, the unfrozen water has enough time to migrate toward the freezing fringe so that the total water content near the freezing front is the maximum. It can be seen from Fig. 8 that the unfrozen water content decreases gradually with the temperature depression in the soil sample. The slope of unfrozen water content curve in the kinetic zone (defined as the range in which supercooling occurs) is very steep while it is relatively gentle far from the phase change zone. Additionally, the ice content increases with the decrease of unfrozen water in the soil column, indicating that the unfrozen water in the pores phase into ice crystals with the propagation of frost. Moreover, the temperature of salt primary crystallization is larger than the freezing temperature, which leads to the consumption of unfrozen water and the decrease of that before the development of freezing. It is found that the ice-water phase change kinetics can well illustrate the variations of unfrozen water and ice contents with temperature in saline soil during the freezing process. The variations of ice crystallization rate are shown in Figs. 9 and 10. It can be seen that the crystallization rate increases first and then decreases with the decreasing of unfrozen water content and cooling rate. It is well known that the supplement of unfrozen water near the cold end is insufficient, which leads to the crystallization time near the cold end is relatively short. Meanwhile, the peak value of ice crystallization rate near the freezing front decreases gradually, indicating that it is related to the cooling rate. From Eq. (1), the rate of ice crystallization is taken as a function of the unfrozen water content and water activity [32].

Fig. 6. Temperature distribution versus height.

Fig. 7. Distribution of total water content of soil column.

Fig. 5. Comparisons between experimental and calculated temperatures at different heights.

distribution is approximately linear owing to the cooling rate gradually reaches a steady state. The phase change temperature was taken as a critical index for the variations of moisture and salt contents during the freezing process [49]. In addition, crystallization temperature is closely related to salt content. Based on the FREZCHEM model [44], it is concluded that the primary crystallization temperature of solution is significantly larger than that of ice-water phase change under the experimental condition of 1.0% salt content. Therefore, it is important to study the effect of salt crystallization during the freezing process.

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

Fig. 8. The contents of unfrozen water and ice versus temperature at different time.

Fig. 9. Ice crystallization rates versus time.

6.3.3. Salt migration and crystallization Salt migration and crystallization have great influence on the internal structure and mechanical properties of soil. Based on the above theoretical analysis, the mechanism of salt migration and crystallization was studied by numerical calculation and experimental investigation. The variations of the molality of the solution along the column with freezing time are shown in Fig. 11. The initial salt content is 1.0%, which is 0.4 mol/kg converting to molality. As shown in Fig. 6, the primary crystallization temperature of solution is significantly larger than that of ice-water phase change. Meanwhile, the temperature drops rapidly in the early stage. Therefore, the sodium sulfate solution in the pores near the cold end first reached a supersaturated state, resulting in the solution concentration decreasing with the precipitation of salt crystals. Subsequently, the concentration of pore solution successively satisfies the supersaturated state and decreases gradually along the direction of frost propagation. In addition, the supersaturation of solution increases gradually with the freezing of unfrozen water, resulting in secondary crystallization. Consequently, it can be seen from Fig. 11 that the solution concentration starts to decrease rapidly after maintaining a slow rate of reduction for a period of time. Generally speaking, the variations of temperature and the unfrozen water content can lead to the change of solution concentration. As shown in Fig. 12, the salt content near the cold end does not change significantly compared with the initial value since the supplement of solution is insufficient. It is well acknowledged that salt migrates through advection and hydrodynamic dispersion [19]. Convection leads to the increase of salt content in the frozen zone, and diffusion effect aggravates salt accumulation at the edge of discrete ice lenses. Therefore, the accumulation of salt in the frozen zone increases gradually with the increasing of water migration. As is well known that the ice grows only by association with water molecules, and molecules of solute are rejected into the unfrozen water [26]. Consequently, the salt distribution during the freezing process has obvious zonal distribution characteristics. Meanwhile, the amount of salt rejection decreases with the increase of cooling rate, making the phenomenon of salt redistribution unobvious near the cold end. Moreover, the effect of salt rejection near the freezing fringe is significantly enhanced by the influences of the decrease of cooling rate and the formation of discrete ice lenses. Finally, the largest salt distribution zone is formed closes to the freezing front. The variations of soil temperature, salt content and unfrozen water content will cause the change of solution concentration. When the pore solution is supersaturated, it begins to crystallize

Fig. 10. Ice crystallization rates along the column at different time.

Therefore, the main factors affecting the ice crystallization rate are the unfrozen water content and cooling rate in the frozen zone. Through the comprehensive analysis of the calculated results, it is recognized that the peak value of ice crystallization rate is determined by cooling rate, and its duration is related to the unfrozen water content in the frozen zone with the propagation of frost.

11

Fig. 11. The molality variations of solution along the column versus time.

12

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

Fig. 12. Distribution of salt content of soil column.

and aggregate. The distribution properties of volumetric salt crystal content at different time are shown in Fig. 13. As shown in Fig. 13a that there is a small amount of salt crystals in pores at the initial stage, and the salt content near the cold end is relatively large. Subsequently, the content of salt crystals increases with the propagation of frost. Salt crystals are layered and distributed discontinuously along the height of soil column toward the direction of frost propagation. Meanwhile, when the temperature reaches the steady stage, the salt crystal content near the freezing front is obviously higher than that in other zones. From Figs. 12 and 13, it is seen that unfrozen water and dissolved salt crystallize alternately during the freezing process, forming discrete ice lenses and discontinuously-distributed-layered salt crystals zones. Moreover, the content of salt crystals near the freezing front is obviously larger than that in other regions, owing to the migration of salt from the unfrozen zone to the frozen fringe and the diffusion effect during the forming process of discrete ice lenses. 6.3.4. Pore water pressure and pore pressure The distributions of pore water pressure in the frozen zone, the freezing fringe and the unfrozen zone at different time are shown in Fig. 14. It is well-known that the variation of the pore water pressure depends on the unfrozen water content and temperature. The pore water pressure first drops rapidly and then increases with

the decrease of temperature due to the effects of phase change and mechanical equilibrium. It is shown that the pore water pressure at the frozen fringe is less than that in other regions due to the effect of ice-water phase change. Moreover, it is seen that the variation of temperature in positive temperature range has no significant influence on pore water pressure. The variation of pore water pressure is controlled by the effects of stress and hydraulic boundary at the freezing front. From Fig. 14, it is deduced that the negative pore water pressure is the main driving force for water and salt migration from the frozen zone to the frozen fringe. From Eq. (35), the pore pressure depends on a variety of substances and their corresponding contents in pores, especially ice and salt crystals. As shown in Fig. 15, the pore pressure in the frozen zone increases significantly due to the ice and salt crystallizations at the early stage. It is well acknowledged that the crystallization pressure is related to the supersaturation at the interface between the solution and salt crystal when the latter is confined between the pore walls. Therefore, the positive pore pressure increases with the decrease of temperature and the propagation of frost. Meanwhile, it can be seen that the maximum negative pore pressure exists at the freezing front. The negative pore pressure along the column is approximately linear in the unfrozen zone [11]. It has been recognized that positive pore pressure is the dominant factor for the deformation of soil sample when it exceeds the tensile strength of the porous medium, while the negative pore pressure is the driving force for water and salt migration. 6.3.5. Soil displacement The variations of soil displacement and freezing front with freezing time are shown in Fig. 16. From Fig. 16, it can be seen that the calculated results are agreement with the experimental data. The calculated freezing front is around 60.0 mm from the bottom of the soil sample at 96 h, and that of the test is 62.9 mm at the same time. During the freezing process, water and salt migrate continuously from the unfrozen zone to the frozen zone, driven by the negative pore pressure, which provides material source for the formation of ice and salt crystals. As mentioned above, the positive pore pressure is the dominant factor for the deformation of soil sample when it exceeds the tensile strength of the porous medium. The alternating formation of ice and salt crystals in pores leads to the deformation of frost heave and salt expansion. As shown in Fig. 16, the calculated soil displacement is 15.6 mm at 96 h, and that of test is 15.3 mm at the same time. By comparing the calculated and experimental results, it is found

Fig. 13. The distribution of volumetric salt crystal content at different time. (a) 1 h, (b) 24 h and (c) 90 h.

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

13

7. Conclusions

Fig. 14. Pore water pressure distribution versus height.

Fig. 15. Pore pressure distribution versus height.

The process of salt migration and precipitation is very complex, which is closely related to heat and mass transfer and soil deformation. A coupled hydro-thermal-salt-mechanical model with phase change is proposed to simulate the one-side freezing process of saturated sulfate saline soil. The crystallization kinetics theory is applied to illuminate the mechanism of ice and salt crystallization. Meanwhile, the water activity and the solution supersaturation are the driving forces for ice-water and solution-salt crystal phase change, respectively. The generalized Clapeyron equation is employed to calculate the ice and water pressure for phase equilibrium condition of water and ice in saline soil. Additionally, the influence of solute on the physical and mechanical properties of soil during the freezing process is investigated. The general form of the governing equations and boundaries are applied to describe the mechanism of saturated sulfate saline soil during the one-side freezing process. Finally, the comparisons of calculated results with experimental data are made to validate the model. The following conclusions are made: (1) Based on the crystallization kinetics, the growth mechanism of ice and salt crystals in the freezing process is clarified through the concepts of the water activity and the solution supersaturation, and it is demonstrated that the microscopic crystallization kinetics theory can well elucidate the phenomenon of macroscopic crystallization. (2) The pore pressure induced by the effect of phase change is the main driving force for water and salt migration as well as the deformation of porous medium. The positive pore pressure is the dominant factor for the deformation of the soil sample when it exceeds the tensile strength of porous material, while the negative pore pressure is the driving force for water migrates to the freezing front. (3) Solute migrates with water and is rejected into the unfrozen water in the process of ice formation. Additionally, the rate of salt rejection gradually increases with the decrease of cooling rate. Therefore, salt crystals with layer distribution are formed in the frozen zone during the freezing process. As a result of the effect of solute diffusion, the largest salt crystals distribution zone is formed near the freezing front. (4) The crystallization of ice and salt in the porous medium will cause the change of macroscopic crystallization pressure on the pore walls, which is the dominant factor for soil deformation. Therefore, preventing the migration and accumulation of water and salt, and reducing the phase change of water and salt crystals can effectively reduce the soil deformation caused by frost heave and salt expansion. (5) The calculated results are agreement with the experimental data, demonstrating that the proposed hydro-thermal-saltmechanical coupling model can well illustrate the mechanism of heat and mass migration of saturated frozen sulfate saline soil, and predict soil deformation due to the effects of frost heave and salt expansion.

Declaration of Competing Interest Fig. 16. Variations of soil displacement and freezing front of soil column with time.

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.

that the proposed hydro-thermal-salt-mechanical coupling model can well elucidate the mechanism of heat and mass migration of saturated frozen sulfate saline soil, and predict the deformation induced by frost heave and salt expansion.

Acknowledgements This research was supported by National Key Research and Development Program of China (Grant No. 2018YFC0809605),

14

J. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118868

Key Research Program of Frontier Sciences of Chinese Academy of Sciences (Grant No. QYZDY-SSW-DQC015), National key Basic Research Program of China (973 Program No. 2012CB026102), the National Natural Science Foundation of China (Grant No. 41230630), the Program of the State Key Laboratory of Frozen Soil Engineering (Grant No. SKLFSE-ZT-23), the National Natural Science Foundation of China (41601074).

References [1] R. Harlan, Analysis of coupled heat-fluid transport in partially frozen soil, Water Res. Res. 9 (5) (1973) 1314–1323. [2] G.L. Guymon, J.N. Luthin, A coupled heat and moisture transport model for arctic soils, Water Res. Res. 10 (5) (1974) 995–1001. [3] G.S. Taylor, J.N. Luthin, A model for coupled heat and moisture transfer during soil freezing, Can. Geotech. J. 15 (4) (1978) 548–555. [4] Y.W. Jame, D.I. Norum, Heat and mass transfer in freezing unsaturated porous media, Water Res. Res. 16 (4) (1980) 811–819. [5] R.R. Gilpin, A model for the prediction of ice lensing and frost heave in soils, Water Res. Res. 16 (5) (1980) 918–930. [6] K. O’Neill, R.D. Miller, Exploration of a rigid ice model of frost heave, Water Res. Res. 21 (3) (1985) 281–296. [7] J.M. Konrad, N.R. Morgenstern, The segregation potential of a freezing soil, Can. Geotech. J. 18 (4) (1981) 482–491. [8] J.F. Nixon, Discrete ice lens theory for frost heave in soils, Can. Geotech. J. 28 (6) (1991) 843–859. [9] Y. Zhou, G.Q. Zhou, Intermittent freezing mode to reduce frost heave in freezing soils-experiments and mechanism analysis, Can. Geotech. J. 49 (6) (2012) 686–693. [10] S.Y. Li, M.Y. Zhang, W.S. Pei, et al., Experimental and numerical simulations on heat-water-mechanics interaction mechanism in a freezing soil, Appl. Therm. Eng. 132 (2018) 209–220. [11] H.R. Thomas, P. Cleall, Y.C. Li, et al., Modeling of cryogenic processes in permafrost and seasonally frozen soils, Geotechnique 59 (3) (2009) 173–184. [12] S. Nishimura, A. Gens, S. Olivella, et al., THM-coupled finite element analysis of frozen soil: formulation and application, Geotechnique 59 (3) (2009) 159–171. [13] Y.M. Lai, W.S. Pei, M.Y. Zhang, et al., Study on theory model of hydro-thermalmechanical interaction process in saturated silty soil, Int. J. Heat Mass Transf. 78 (2014) 805–819. [14] G.Q. Zhou, Y. Zhou, K. Hu, et al., Separate-ice frost heave model for onedimensional soil freezing process, Acta Geotech. 10 (2018) 1–11. [15] M. Hassanizadeh, W.G. Gray, General conservation equations for multi-phase systems: 1. Averaging procedure, Adv. Water Resour. 2 (1979) 131–144. [16] F.G. Tong, L.R. Jing, R.W. Zimmerman, A fully coupled thermo-hydromechanical model for simulating multiphase flow, deformation and heat transfer in buffer material and rock masses, Int. J. Rock Mech. Min. Sci. 47 (2) (2010) 205–217. [17] M. Koniorczyk, Coupled heat and water transport in deformable porous materials considering phase change kinetics, Int. J. Heat Mass Transf. 81 (2015) 260–271. [18] S.H. Na, W.C. Sun, Computational thermo-hydro-mechanics for multiphase freezing and thawing porous media in the finite deformation range, Comput. Methods Appl. Mech. Eng. 318 (2017) 667–700. [19] M. Koniorczyk, Salt transport and crystallization in non-isothermal, partially saturated porous materials considering ions interaction model, Int. J. Heat Mass Transf. 55 (2012) 665–679. [20] N. Mokni, S. Olivella, E.E. Alonso, Swelling in clayey soils induced by the presence of salt crystals, Appl. Clay Sci. 47 (2010) 105–112. [21] A.H. Alizadeh, M. Akbarabadi, E. Barsotti, et al., Salt precipitation in ultratight porous media and its impact on pore connectivity and hydraulic conductivity, Water Res. Res. 54 (2018) 2768–2780. [22] J.Z. Zhou, C.F. Wei, Y.M. Lai, et al., Application of the generalized clapeyron equation to freezing point depression and unfrozen water content, Water Res. Res. 54 (11) (2018) 1–20. [23] A. Naillon, P. Joseph, M. Prat, Ion transport and precipitation kinetics as key aspects of stress generation on pore walls induced by salt crystallization, Phys. Rev. Lett. 120 (3) (2018) 034502.

[24] S.M. Hassanizadeh, T. Leijnse, On the modeling of brine transport in porous media, Water Res. Res. 24 (3) (1988) 321–330. [25] G.C. Baker, T.E. Osterkamp, Salt redistribution during freezing of saline sand columns at constant rates, Water Res. Res. 25 (8) (1989) 1825–1831. [26] S. Panday, M.Y. Corapcioglu, Solute rejection in freezing soils, Water Res. Res. 27 (1) (1991) 99–108. [27] Z. Pavlik, L. Fiala, J. Madera, et al., Computational modelling of coupled water and salt transport in porous materials using diffusion-advection model, J. Franklin Instit. 348 (7) (2011) 1574–1587. [28] T.Q. Nguyen, J. Petkovic, P. Dangla, Modelling of coupled ion and moisture transport in porous building materials, Constr. Build. Mater. 22 (11) (2008) 2185–2195. [29] G. Castellazzi, C. Colla, S. de Miranda, et al., A coupled multiphase model for hygrothermal analysis of masonry structures and prediction of stress induced by salt crystallization, Constr. Build. Mater. 41 (2013) 717–731. [30] G. Castellazzi, S. de Miranda, L. Grementieri, Multiphase model for hygrothermal analysis of porous media with salt crystallization and hydration, Mater. Struct. 49 (3) (2016) 1039–1063. [31] X.D. Zhang, Q. Wang, T.W. Yu, et al., Numerical study on the multifield mathematical coupled model of hydraulic-thermal-salt-mechanical in saturated freezing saline soil, Int. J. Geomech. 18 (7) (2018) 04018064. [32] L. Bronfenbrener, E. Korin, Kinetic model for crystallization in porous media, Int. J. Heat Mass Transf. 40 (5) (1997) 665–679. [33] T. Koop, B. Luo, T. Peter, Water activity as the determinant for homogeneous ice nucleation in aqueous solutions, Nature 406 (2000) 611–614. [34] D.Y. Wu, Y.M. Lai, M.Y. Zhang, Heat and mass transfer effects of ice growth mechanisms in a fully saturated soil, Int. J. Heat Mass Transf. 86 (2015) 699– 709. [35] R.M. Espinosa, L. Franke, G. Deckelmann, Phase changes of salts in porous materials: Crystallization, hydration and deliquescence, Constr. Build. Mater. 22 (2008) 1758–1773. [36] N. Khalili, M.A. Habte, S. Zargarbashi, A fully coupled flow deformation model for cyclic analysis of unsaturated soils including hydraulic and mechanical hysteresis, Comput. Geotech. 35 (6) (2008) 872–889. [37] J.P.G. Loch, Thermodynamic equilibrium between ice and water in porous media, Soil Sci. 126 (2) (1978) 77–80. [38] Y.M. Lai, M.Y. Zhang, S.Y. Li, Theory and Application of Cold Regions Engineering, Science Press, Beijing, 2009. [39] M. Koniorczyk, D. Gawin, Heat and moisture transport in porous building materials containing salt, J. Build. Phys. 31 (2008) 279–300. [40] J. Choo, W.C. Sun, Cracking and damage from crystallization in pores: coupled chemo-hydro-mechanics and phase-field modeling, Comput. Methods Appl. Mech. Eng. 335 (2018) 347–379. [41] D.U. Ophori, The significance of viscosity in density-dependent flow of groundwater, J. Hydrol. 204 (1998) 261–270. [42] M. Koniorczyk, Modelling the phase change of salt dissolved in pore waterequilibrium and non-equilibrium approach, Constr. Build. Mater. 24 (7) (2010) 1119–1128. [43] G.M. Marion, R.E. Farren, Mineral solubilities in the Na-K-Mg-Ca-Cl-SO4-H2O system: a re-evaluation of the sulfate chemistry in the Spencer-Møller-Weare model, Geochim. Cosmochim. Acta. 63 (9) (1999) 1305–1318. [44] R.J. Spencer, N. Møller, J.H. Weare, The prediction of mineral solubilities in natural waters: A chemical equilibrium model for the Na-K-Ca-Mg-Cl-SO4H2O system at temperatures below 25°C, Geochim. Cosmochim. Acta. 54 (3) (1990) 575–590. [45] M. Steiger, Crystal growth in porous materials-I: the crystallization pressure of large crystals, J. Crystal Growth 282 (2005) 455–469. [46] K.S. Pitzer, Ion interaction approach: theory and data correlation, in: K.S. Pitzer (Ed.), Activity Coefficients In Electrolyte Solutions, second ed., CRC Press, Berlin, 1991, pp. 75–153. [47] X.Z. Xu, J.C. Wang, L.X. Zhang, Frozen Soil Physics, Science Press, Beijing, 2001. [48] Z. Xiao, Y.M. Lai, M.Y. Zhang, Study on the freezing temperature of saline soil, Acta Geotech. (2017). [49] X.S. Wan, Y.M. Lai, C. Wang, Experimental study on the freezing temperatures of saline silty soils, Permafrost Periglac. 26 (2) (2015) 175–187. [50] R.L. Michalowski, M. Zhu, Frost heave modeling using porosity rate function, Int. J. Numer. Anal. Methods Geomech. 30 (8) (2006) 703–722. [51] D. Wu, Y. Lai, M. Zhang, Thermo-hydro-salt-mechanical coupled model for saturated porous media based on crystallization kinetics, Cold Reg. Sci. Technol. 133 (2017) 94–107.