# Numerical Analysis Method considering Coupled Effects of THMC Multifields on Unsaturated Expansive Soil Subgrade Treated with Lime

Authors: Jie Liu; Pan Chen; Zheng Lu; Hailin Yao
Journal: Geofluids
 (2018)

---

## Abstract

The response model subjected to coupled effect of thermo-hydro-mechano-chemical (THMC) was built in the context of basic theories in the polyporous polyphasic medium mechanics, the mixture theory in the continuum mechanics, and the thermodynamic theories. The finite element discretization of the response model was implemented based on the Galerkin method. The processes of salt leaching and accumulating were analyzed in the numerical results. The results among the numerical results and measured results were compared and discussed. Finally, the solutes migration rule of the soil subjected to the atmosphere eluviations was revealed, and the reasonableness of the coupling model and the finite element method was proved. The agreement between the numerical and the measured results was good, which indicates that the THMC model and finite element program were useful in solving the coupling problems of unsaturated soil. Moreover, the salt dissolution process has a larger effect on the salt movement compared to that of the salt accumulation process. Comparing with the salt leaching effect caused by rainfall, the salt accumulation effect caused by evaporation was smaller.


---

## Body

## 1. Introduction

In China, the stability of subgrade is the basic premise to guarantee highways becoming the leading infrastructure construction. Expansive soil cannot be directly used as the subgrade filling within the subgrade. Engineering properties and mechanical strength indexes of the expansive soil can be increased significantly through modification using the lime. However, the expansive soil subgrade treated by the lime will easily suffer from repetitive drying-wetting cycles and unceasing strength degradation in the rainy and hot weather, which would make the subgrade failure. It is essential to study the effects of solute migration in the expansive soil treated by lime under the atmospheric actions on the subgrade stability. A lot of domestic and overseas experts have developed related research about soil subgrade treatments by lime. Kumar et al. [1] studied, compared, and analyzed variations of compaction and strength properties of the expansive soil before and after treatment by the lime. Nowamooz and Masrouri [2] researched volumetric strains caused by variations in suction or stress of the expansive mixture of limestone soil and flyash. Guney et al. [3] studied the expansibility of expansive soil treated by the lime under drying and wetting cycles. Estabragh et al. [4] researched mechanical behaviors of the expansive soil experienced drying and wetting cycles with different wetting fluids. Tripathy and Subba Rao [5] studied the cyclic swelling and shrinking behaviors of tamped expansive soil through theoretical analysis and indoor experiments. Rajeev et al. [6] researched variation trends of humidity and temperature of the subgrade in the long term under atmospheric actions. Al-Mukhtar et al. [7, 8] studied the micromechanism and mechanical behaviors of the expansive soil treated by the lime. Wang et al. [9] researched variation characteristics of kinetic parameters of the expansive soil treated by the lime. Stoltz et al. [10] performed multiscale analysis of swelling and shrinking of the expansive soil treated by the lime from the microscopic view. Ojuri et al. [11] studied the geotechnical and environmental evaluation of lime-cement stabilized soil-mine tailing mixtures for highway construction. Dash and Hussain [12] researched the influence of lime on shrinkage behavior of soils. Hotineanu et al. [13] studied the effect of freeze-thaw cycling on the mechanical properties of lime-stabilized expansive clays.The above researches show that a lot of experts have performed fruitful researches on physical and mechanical properties of the expansive soil treated by the lime from macroscopic and microscopic views. However, considering atmospheric actions, soil modification, and solute migration rule to study the long-term properties of the expansive soil treated by the lime in the highway subgrade would got closer to the practical results. In this paper, as for the expansive soil treated by the lime under atmospheric actions, the thermo-hydro-mechano-chemical coupled model was built on the basis of the Polyporous Polyphasic Medium Mechanics, the mixture theory in the continuum mechanics, and thermodynamic theories of the irreversible process. The finite element discretization of the response model is implemented based on the Galerkin method. The processes of salt leaching and accumulating are analyzed in the numerical results. The numerical results and measured results are compared and discussed. Finally, the solutes migration rule of the soil subjected to the atmosphere eluviation is revealed, and the reasonableness of the coupling model and the finite element method is proved.


## 2. Coupled Equations for THMC

### 2.1. Basic Assumptions

Some basic assumptions were given before building THMC model to make the problem has universality: (a) the unsaturated soil is polyporous medium, its basic framework is particles, and liquid water and gas exist in the pores of the particle framework. (b) The gas phase is ideal gases with the mixture of dry air and water vapor. The liquid phase consists of water and gases dissolving in water. (c) Mechanical behavior of the polyporous medium (gases, liquids, and solids) satisfies the principle of local stress balance, and the entire framework satisfied the principle of effective stress. (d) The temperature of the soil is the average temperature, which satisfied the local equilibrium hypothesis, heat transfer satisfies the Fourier law, and the framework and water can be considered as isotropic thermoelastic materials. (e) The unsaturated soil satisfies general assumptions and the continuum hypothesis of the mixtures theory. (f) Actions of degradation and stagnant water effect are not taken into account in the solute equilibrium equation.


### 2.2. Balance Equations

In the saturated-unsaturated soil mass, the total stress incrementdσ and effective stress increment dσ′ can be expressed in the incremental forms as follows, respectively:(1)dσ=dσ′-mdp,(2)dσ′=DTdε-dεp-dεT-dεsw,where DT is the tangent elastoplastic modulus matrix; dε is the total strain increment; dεp is the strain increment caused by the pore pressure and dεp=-mdp/3Ks; dεT is the strain increment caused by temperature, and dεT=mβs/3dT; dεsw is the strain increment caused by swelling, and dεsw=-mβswSw; m is the unit array of normal stress; βs is the thermal expansion coefficient of solid particles; βsw is the water swelling expansion coefficient; Ks is the bulk modulus of solid particles; p is the average stress of gases and liquids inside pores, and p=Swpw+Sapa; Sw and Sa are saturability of the pore water and pore gas, respectively, and Sw+Sa=1; pw and pa are pressures of the pore water and pore gas, respectively.The stress equilibrium equation is stated as follows:(3)∂σij∂t-fi=∇·σ-fi=0,where fi is the body force.According to the principle of virtual work, the incremental stress equilibrium equation can be stated as follows:(4)∫ΩδεTdσdΩ-∫ΩδuTdbdΩ-∫ΓδuTdt^dΓ=0,where dσ is the total stress increment; db and dt^ are increments of body force and surface force, respectively; δε and δu are virtual strain and displacement, respectively.Equations (2) and (3) are substituted into (4), and the following is obtained:(5)∫ΩδεTDTdεdtdΩ-∫ΩδεTmSwdpwdt+pwDsndpwdtdΩ+∫ΩδεTDTmβswDsndpwdtdΩ+∫ΩδεTDTm3KsSwdpwdt+pwDsndpwdtdΩ-∫ΩδεTDTmβs3dTdtdΩ-dfdt=0.According to the law of conservation of mass, the fluid continuity equation can be obtained:(6)∂nSφρφπ∂t+∇·nSφρφπu˙π-θπ=0.Similarly, the equation of solid mass conservation can be obtained:(7)∂1-nρs∂t+∇·1-nρsu˙s=0.Equation (7) is the mass conservation equation of the component π in the φ phase. It is the liquid phase as φ is equal to l. It is the gas phase as φ is equal to g. The component is water as π is equal to w, and the component is the gas as π is equal to a. For instance, ρla and ρlw indicate air resolved in the liquid phase and water in the liquid phase, respectively; ρga and ρgw indicate dry air and water steam in the gas phase, respectively; m˙ is the decreasing or increasing rate of unit volume water during evaporation or condensation. u˙s and u˙π are absolute moving velocities of solid particles and fluid, respectively; Sφ is saturability of the water phase or the gas phase.The fluid fluxqφπ, which flows across the unit area per unit time, is stated as follows:(8)qφπ=nSφρφπu˙π=iφπ+nSφρφπu˙rφ+u˙s=qrφπ+nSφρφπu˙s,where iφπ is the diffusion flow of component π in the φ phase; qrφπ is the relative flow when there are movements of the solid phase in the deformed porous medium, and the subscript r indicates that the velocity and flow are relative to the solid phase. According to the definition of relative velocity, it can be obtained that u˙rφ=u˙φ-u˙s. Equation (6) can be stated as follows:(9)∂nSφρφπ∂t+∇·nSφρφπu˙si+∇·qrφπ+θπ=0.According to (6), (7), and the definition of material derivative, D/Dt∗=∂∗/∂t+u˙s·∇∗, it can be obtained that(10)1-nρsDρsDt-DnDt+1-n∇·u˙s=0,(11)nDSφρφπDt+SφρφπDnDt+nSφρφπ∇·u˙s+∇·qrφπ+θπ=0.Substituting (10) into (11), the equilibrium equation of phase component in the deformed porous medium can be obtained as follows:(12)1-nρsDρsDt+nSφρφπDSφρφπDt+∇·u˙s+1Sφρφπ∇·qrφπ+1Sφρφπθπ=0.Under the situation of small deformation, items such as u˙s·∇Sφρφπ and u˙s·∇ρs can be neglected; ∇·u˙s can be presented as the volumetric strain of the framework εv:∇·u˙s=Dεv/Dt≈∂εv/∂t. Equation (12) can be simplified as follows:(13)n∂Sφρφπ∂t+Sφρφπ∂εv∂t+1-nρs∂ρsDt+θπ=-∇·qrφπ.In (13), the right item is the material flow including convective (Darcy) and nonconvective diffusion; the left items indicate the change rates of fluid material caused by the saturability and density (the first item), volumetric strain (the second item), volumetric deformation of soil particles (the third item), and phase change (the fourth item), respectively.Considering the mass conservation of every component,∑kθπ=0, the total equilibrium equation is as follows:(14)n∂∂tSlρla+Sgρga+Slρla+Sgρga∂εv∂t+1-nρs∂ρsDt=-∇·qrla+qrga.As for the water including liquid water and steam in the air (π=w), the total equilibrium equation is as follows:(15)n∂∂tSlρlw+Sgρgw+Slρlw+Sgρgw∂εv∂t+1-nρs∂ρsDt=-∇·qrlw+qrgw.Equation (14) and (15) are continuity equations of the gas and liquid, respectively. For the sake of convenience in writing, in the following derivation process, ρlw=ρl,  ρgw=ρv,  qrlw=ql, and qrgw=qv.Variations of densities of steam, liquid water, and particles with time can be expressed as(16)∂ρl∂t=∂ρw1+γc∂t=1+γc∂ρw∂t+ρwγ∂c∂t=1+γcρwKw∂Pw∂t-ρwβw∂T∂t+ρwγ∂c∂t,∂ρv∂t=∂ρv∂Pw∂Pw∂t+∂ρv∂c∂c∂t+∂ρv∂T∂T∂t,1-nρs∂ρsDt=a-n∂PKs∂t-1-αu˙ii-α-nβsT˙,where ρl=ρw1+γc;  ρw is the density of fresh water; γ is the density ratio of solute to water in the pore water; Kw is the compression modulus of fresh water; βw is the thermal expansion coefficient of fresh water; Ks is the compression modulus of particles; βs is the thermal expansion coefficient of the framework; and a is the Biot coefficient.Considering effects of temperature and solute concentration in the unsaturated soil, the effusion flux of water steam can be described using the revised Fick law as follows:(17)qv=-DvτSa∇ρv,where qv is the effusion flux of water steam, with a unit of kg/m2/s; τ=n-Sw2/3 is the factor considering pore curvature in the soil [14], and Dv=229×10-7T/273.151.75 is steam diffusion coefficient [15].Considering the effects of solute in the unsaturated soil, the density of steamρv can be expressed as(18)ρv=hrρvs=Rhhcρvs.Substituting (16)~(18) into (15), the continuous equilibrium equation of pore water considering combined effects of framework deformation, steam, solute, and temperature can be obtained as follows:(19)Cpe∂ε∂t+Cpp∂pw∂t+CpT∂T∂t+Cpc∂c∂t=∇·Kpp∇pw+ρl2gKkrlμw+∇·KpT∇T+∇·Kpc∇c.The coefficients in the equation are shown in Appendix A.


### 2.3. Energy Conservation Equation

The total internal energy of every fluid phaseMφh and energy exchange flux qφh can be expressed as(20a)Mφh=nSφρφeφ,(20b)qφh=iφh+eφaqrφa+eφwqrφw+nSφρφeφu˙s=iφh+qrφh+nSφρφeφu˙s.Similarly, as for the solid phase, it can be obtained:(21a)Msh=1-nρses,(21b)qsh=ish+esqrs+1-nρsesu˙s=ish+qrsh+1-nρsesu˙s,where eφ is a unit mass of internal energy of phase φ; eφa and eφw are energies of water and gas per unit mass.According to the total energy conservation and the assumption of local thermal equilibrium, it can be obtained:(22)∂∂nSlρlel+nSgρvev+1-nρses-Qh=-∇·imh+qrlh+qrgh+qrsh,where imh is the average heat conduction flux; qrlh,  qrgh, and qrsh are convection fluxes of the liquid phase, gas phase, and solid phase, respectively; Qh is the term of source and convergence. Effects of thermal conduction in deformation of solid phases are neglected.When variations of temperature and pressure are not large enough, the internal energy in the energy balance equation can be replaced by the specific enthalpy. Specific enthalpies of the solid phase, water, and water steam can be expressed as follows:(23a)hs+hl=CsT-To+ClT-To,(23b)hv=hpv+hLo=CpvT-T0+LoTo=ClT-T0+LT,where hs,  hl, and hv are specific enthalpies of the solid phase, water, and water steam, respectively; Cs and Cl are the specific heat of solids and water; Cpv is the specific heat at constant pressure of water steam; Lo is the latent heat of vaporization of water steam.Replacing the internal energy with enthalpy in (23a) and (23b),(24)∂hs∂t=1-nρsCs∂T∂t,(25)∂nSwρlhl∂t=nSwρlCl∂T∂t+nρlClT-T0∂Sw∂t+nSwClT-T0∂ρl∂t,(26)∂nSgρvhpv∂t=nSgρvCpv∂T∂t+nSgCpvT-T0∂ρv∂t+nρvCpvT-T0∂Sg∂t,(27)∂nSgρvhL0∂t=nSgL0∂ρv∂t+nρvL0∂Sg∂t.In the above derivation, effect of the density of solids on energy is neglected, and it is assumed that the specific heat can be considered as the invariant when variations of pressure and temperature are not large.The right part of the energy equation (22) can be abbreviated as follows:(28)qh=-λI∇T+ClT-T0qrl+ClT-T0+LTqrv.By substituting (23a), (23b), and (27) into (22), the energy conservation equation can be obtained as(29)CTp∂pw∂t+CTT∂T∂t+CTc∂c∂t+VTT·∇T=∇·KTp∇pw+∇·KTT∇T+∇·KTc∇c.The coefficients in the equation are shown in Appendix B.


### 2.4. Solute Transport Equation

When the temperature does not vary significantly, it can be considered that migration of solute under the temperature gradient can be neglected; therefore the solute transport equation can be obtained according to the mass conservation law of the solute:(30)∂∂tθc+∂∂tρdcP=∇·Dsh∇c-qc,where cP is the amount of adsorbed solute in soil per mass, which is a function of concentration. Then (30) can be expressed as(31)nSw+ρd∂cp∂c∂c∂t+cDs∂pw∂t=∇·Dsh∇c-qc,where θ is the volumetric water content of pore water; ρd is the density of dry soil; Dsh is the hydrodynamic dispersion coefficient of the solute; q is the water flux; Ds is the moisture capacity.


### 2.5. Initial Conditions and Boundary Conditions

Initial conditions inside the given analysis areaΩ and on the boundary Γπ(π=u,w,T,c) when t is zero are as follows:(32)u=u0,pw=pw0,T=T0,c=c0,onΩandΓπ.Dirichlet boundary conditions can be expressed as follows:(33)ut=u^t,onΓu;pwt=p^wt,onΓw;Tt=T^t,onΓT;ct=c^t,onΓc.Neumann boundary conditions can be expressed as follows:(34a)ITσ=t-,onΓuq;Kkrlμw∇ψ+DvτSa∇ρv·n=qw+qv+βcρv-ρv∞,onΓwq,(34b)ρwvwL-λM∇T·n=qT+αcT-T∞,onΓTq;Dsh∇c·n=qc,onΓTq,where qw,  qv,  qT,  qc, and t- are assigned liquid water flow, water steam flow, solute flow, and tensile force on the boundary Γπ(π=u,w,T,c). ρv∞ and T∞ are the mass density and temperature of water steam at undisturbed locations far away from boundaries; βc is the mass convection coefficient; αc is the thermal convection coefficient. Mixed boundary conditions may appear during the actual calculation. Therefore boundary conditions should be determined according to the actual situation.The initial concentration distribution is given as(35a)c=c0z,t=0,0≤z≤L,where L is the depth of calculated soil layer.At the earth surface wherez is zero, boundary conditions of solute migration can be one kind of the Dirichlet, Neumann, and Cauchy boundary conditions.(1) As the solute concentration at the earth surface is given, the Dirichlet boundary condition is expressed as follows:(35b)c=cz,z=0,t>0.(2) As the earth surface is in the status of infiltration and the solute concentration cRt or the degree of mineralization of rainfall or irrigation water is given, it is the Cauchy boundary condition expressed as follows since the solute transport flux J=-Dsh∇c+qc:(35c)-Dsh∇c+qc=qcRt,z=0,t>0.As the intensity of water supplyRt is smaller than the soil infiltration capacity, the water movement flux at the earth surface in the above equation is equal to the intensity of water supply Rt; as the intensity of water supply Rt is larger than the soil infiltration capacity, the water movement flux at the earth surface can be calculated through the above equation.(3) As the earth surface is in the status of evaporation, the intensity of evaporationEt adopts the Penman-Wilson formula of the actual evaporation amount at the surface of the unsaturated soil as the evaporation boundary condition:(35d)Et=ΔRn+ηEaΔ+ηA,where Δ is the slope of the vapor pressure-temperature curve of the saturated steam; Rn is the amount of net radiation at soil surface; η is the humidity constant; Ea is the amount of potential evaporation and Ea=fuPgwB-A. Here fu is a function of wind, and fu=0.35×1+0.15u, where u is the wind velocity; Pgw is the vapor pressure of air above the evaporating surface; B is the reciprocal of relative humidity of the air, which is 1/hB; A is the reciprocal of relative humidity of the soil surface, which is 1/Rh.Here it belongs to the Cauchy boundary condition, which can be expressed as follows since the solute flow is zero at the earth surface during evaporation:(35e)Dsh∇c+cEt=0,z=0,t>0.(4) During field or indoor experiments, as it has neither evaporation nor infiltration on the earth surface, it belongs to the Neumann boundary condition, and it can be obtained:(35f)∇c=0.


### 2.6. Finite Element Discretization of the THMC Coupled Governing Equation

Analytical solutions of the THMC coupled mathematical model can be hardly achieved due to the complexity of solution conditions and irregularity of solution region. Therefore, the only way to solve the model is to obtain approximate solutions through numerical methods. The THMC coupled governing equations such as (5), (19), (29), and (31) are spatial discretized as follows:(36)Kdu-dt+Cdp-wdt+IdT-dt=dfdt,Edu-dt+Fp-w+Gdp-wdt+MT-+LdT-dt+Udc-dt=f^,Pdp-wdt+QT-+Jp-w+RdT-dt+Sdc-dt=X,Wc-+Zdc-dt=Y,where(37)K=∫ΩBTDTBdΩ,C=-∫ΩBTm1-DT3KsSw+p-wDsn+DTmβswDsnNpdΩ,I=-∫ΩBTDTmβs3NTdΩ,df=∫ΩNuTdbdΩ+∫ΓNuTdt^dΓ,E=∫ΩNpTSwρw1+γc+1-SwρvmT-mTDT3KsBdΩ,G=∫ΩNpTDsρw1+γc-ρv+n1-Sw∂ρv∂pw+nSw1+γcρwKw+Swρw1+γc+1-Swρv·1-nKw-mTDTm3Ks2·NpdΩ,F=∫Ω∇NpTρlKkrlμw+DvτSaρvMwρlRTNpdΩ,M=∫Ω∇NpT-ρlKkrlμwυϕCRg+DvτSahr∂ρvs∂T+ρgwshc∂Rh∂T+ρvsRh∂hc∂TNTdΩ,L=∫ΩNpTn1-Sw∂ρv∂T-Swρw1+γc+1-Swρv·1-nβs-mTDTm3Ksβs3-nSwρw1+γcβwNTdΩ,U=∫ΩNpTn1-Sw∂ρv∂c+nSwρwγNcdΩ,f^=-∫ΓwNpTqw+qv+βcρv-ρv∞dΓ,P=∫ΩNTTρwCl1+γcT-T0Ds+nSwKw-ρvDsCpvT-T0-L0+n1-Sw∂ρv∂pwCpvT-T0+L0NpdΩ,Q=∫Ω∇NTTλM+LDTVNTdΩ,J=∫Ω∇NTTLDPVNpdΩ,R=∫ΩNTTρcM-nSwClT-T01+γcρwβw+n1-SwCpvT-T0+L0∂ρv∂TNTdΩ,S=∫ΩNTTnSwClT-T0ρwγ+n1-SwCpvT-T0+Lo∂ρv∂cNcdΩ,X=-∫ΓTNTTqT+αcT-T∞dΓ,W=∫Ω∇NcTDsh∇NcdΩ+∫Ω∇NcTqNcdΩ,Z=∫ΩNcTnSw+ρb∂cp∂cNcdΩ,Y=-∫ΓcNcTqcdΓ.Equation (36) can be expressed in the form of matrix as follows:(38)KCI0EGLU0PRS000Zddtu-p-wT-c-+00000FM00JQ0000Wu-p-wT-c-=dfdtf^XY.The Monolithic Augmentation Approach is applied to solve the model. Changes within every time stepΔt are assumed as linear, which means(39)u-p-wT-c-tk+α=N1tN2tu-tkp-wtkT-tkc-tku-tk+Δtkp-wtk+ΔtkT-tk+Δtkc-t+Δtk,where α is a point within the time step Δtk; the item at the left of the equation contains the values of u-,  p-w,  T-, and c- at the time of tk+α, which means α=t-tk/Δtk; N1=1-α and N2=α.Therefore (38) can be transformed into the following iterative solution equation:(40)KCI0EG+αFΔtkL+αFΔtkU0P+αJΔtkR+αQΔtkS000Z+αWΔtktk+αu-p-wT-c-tK+ΔtK=KCI0EG-1-αFΔtkL-1-αMΔtkU0P-1-αJΔtkR-1-αQΔtkS000Z-1-αWΔtktk+αu-p-wT-c-tK+dfdtf^XYtk+αΔtK.The above equation is the spatial discretization of nonlinear equations. Hence the nonlinear governing equation is transformed into linear problem through discretization of space and time. The asymmetric stiffness matrix can be directly solved by using the asymmetric solver.


## 2.1. Basic Assumptions

Some basic assumptions were given before building THMC model to make the problem has universality: (a) the unsaturated soil is polyporous medium, its basic framework is particles, and liquid water and gas exist in the pores of the particle framework. (b) The gas phase is ideal gases with the mixture of dry air and water vapor. The liquid phase consists of water and gases dissolving in water. (c) Mechanical behavior of the polyporous medium (gases, liquids, and solids) satisfies the principle of local stress balance, and the entire framework satisfied the principle of effective stress. (d) The temperature of the soil is the average temperature, which satisfied the local equilibrium hypothesis, heat transfer satisfies the Fourier law, and the framework and water can be considered as isotropic thermoelastic materials. (e) The unsaturated soil satisfies general assumptions and the continuum hypothesis of the mixtures theory. (f) Actions of degradation and stagnant water effect are not taken into account in the solute equilibrium equation.


## 2.2. Balance Equations

In the saturated-unsaturated soil mass, the total stress incrementdσ and effective stress increment dσ′ can be expressed in the incremental forms as follows, respectively:(1)dσ=dσ′-mdp,(2)dσ′=DTdε-dεp-dεT-dεsw,where DT is the tangent elastoplastic modulus matrix; dε is the total strain increment; dεp is the strain increment caused by the pore pressure and dεp=-mdp/3Ks; dεT is the strain increment caused by temperature, and dεT=mβs/3dT; dεsw is the strain increment caused by swelling, and dεsw=-mβswSw; m is the unit array of normal stress; βs is the thermal expansion coefficient of solid particles; βsw is the water swelling expansion coefficient; Ks is the bulk modulus of solid particles; p is the average stress of gases and liquids inside pores, and p=Swpw+Sapa; Sw and Sa are saturability of the pore water and pore gas, respectively, and Sw+Sa=1; pw and pa are pressures of the pore water and pore gas, respectively.The stress equilibrium equation is stated as follows:(3)∂σij∂t-fi=∇·σ-fi=0,where fi is the body force.According to the principle of virtual work, the incremental stress equilibrium equation can be stated as follows:(4)∫ΩδεTdσdΩ-∫ΩδuTdbdΩ-∫ΓδuTdt^dΓ=0,where dσ is the total stress increment; db and dt^ are increments of body force and surface force, respectively; δε and δu are virtual strain and displacement, respectively.Equations (2) and (3) are substituted into (4), and the following is obtained:(5)∫ΩδεTDTdεdtdΩ-∫ΩδεTmSwdpwdt+pwDsndpwdtdΩ+∫ΩδεTDTmβswDsndpwdtdΩ+∫ΩδεTDTm3KsSwdpwdt+pwDsndpwdtdΩ-∫ΩδεTDTmβs3dTdtdΩ-dfdt=0.According to the law of conservation of mass, the fluid continuity equation can be obtained:(6)∂nSφρφπ∂t+∇·nSφρφπu˙π-θπ=0.Similarly, the equation of solid mass conservation can be obtained:(7)∂1-nρs∂t+∇·1-nρsu˙s=0.Equation (7) is the mass conservation equation of the component π in the φ phase. It is the liquid phase as φ is equal to l. It is the gas phase as φ is equal to g. The component is water as π is equal to w, and the component is the gas as π is equal to a. For instance, ρla and ρlw indicate air resolved in the liquid phase and water in the liquid phase, respectively; ρga and ρgw indicate dry air and water steam in the gas phase, respectively; m˙ is the decreasing or increasing rate of unit volume water during evaporation or condensation. u˙s and u˙π are absolute moving velocities of solid particles and fluid, respectively; Sφ is saturability of the water phase or the gas phase.The fluid fluxqφπ, which flows across the unit area per unit time, is stated as follows:(8)qφπ=nSφρφπu˙π=iφπ+nSφρφπu˙rφ+u˙s=qrφπ+nSφρφπu˙s,where iφπ is the diffusion flow of component π in the φ phase; qrφπ is the relative flow when there are movements of the solid phase in the deformed porous medium, and the subscript r indicates that the velocity and flow are relative to the solid phase. According to the definition of relative velocity, it can be obtained that u˙rφ=u˙φ-u˙s. Equation (6) can be stated as follows:(9)∂nSφρφπ∂t+∇·nSφρφπu˙si+∇·qrφπ+θπ=0.According to (6), (7), and the definition of material derivative, D/Dt∗=∂∗/∂t+u˙s·∇∗, it can be obtained that(10)1-nρsDρsDt-DnDt+1-n∇·u˙s=0,(11)nDSφρφπDt+SφρφπDnDt+nSφρφπ∇·u˙s+∇·qrφπ+θπ=0.Substituting (10) into (11), the equilibrium equation of phase component in the deformed porous medium can be obtained as follows:(12)1-nρsDρsDt+nSφρφπDSφρφπDt+∇·u˙s+1Sφρφπ∇·qrφπ+1Sφρφπθπ=0.Under the situation of small deformation, items such as u˙s·∇Sφρφπ and u˙s·∇ρs can be neglected; ∇·u˙s can be presented as the volumetric strain of the framework εv:∇·u˙s=Dεv/Dt≈∂εv/∂t. Equation (12) can be simplified as follows:(13)n∂Sφρφπ∂t+Sφρφπ∂εv∂t+1-nρs∂ρsDt+θπ=-∇·qrφπ.In (13), the right item is the material flow including convective (Darcy) and nonconvective diffusion; the left items indicate the change rates of fluid material caused by the saturability and density (the first item), volumetric strain (the second item), volumetric deformation of soil particles (the third item), and phase change (the fourth item), respectively.Considering the mass conservation of every component,∑kθπ=0, the total equilibrium equation is as follows:(14)n∂∂tSlρla+Sgρga+Slρla+Sgρga∂εv∂t+1-nρs∂ρsDt=-∇·qrla+qrga.As for the water including liquid water and steam in the air (π=w), the total equilibrium equation is as follows:(15)n∂∂tSlρlw+Sgρgw+Slρlw+Sgρgw∂εv∂t+1-nρs∂ρsDt=-∇·qrlw+qrgw.Equation (14) and (15) are continuity equations of the gas and liquid, respectively. For the sake of convenience in writing, in the following derivation process, ρlw=ρl,  ρgw=ρv,  qrlw=ql, and qrgw=qv.Variations of densities of steam, liquid water, and particles with time can be expressed as(16)∂ρl∂t=∂ρw1+γc∂t=1+γc∂ρw∂t+ρwγ∂c∂t=1+γcρwKw∂Pw∂t-ρwβw∂T∂t+ρwγ∂c∂t,∂ρv∂t=∂ρv∂Pw∂Pw∂t+∂ρv∂c∂c∂t+∂ρv∂T∂T∂t,1-nρs∂ρsDt=a-n∂PKs∂t-1-αu˙ii-α-nβsT˙,where ρl=ρw1+γc;  ρw is the density of fresh water; γ is the density ratio of solute to water in the pore water; Kw is the compression modulus of fresh water; βw is the thermal expansion coefficient of fresh water; Ks is the compression modulus of particles; βs is the thermal expansion coefficient of the framework; and a is the Biot coefficient.Considering effects of temperature and solute concentration in the unsaturated soil, the effusion flux of water steam can be described using the revised Fick law as follows:(17)qv=-DvτSa∇ρv,where qv is the effusion flux of water steam, with a unit of kg/m2/s; τ=n-Sw2/3 is the factor considering pore curvature in the soil [14], and Dv=229×10-7T/273.151.75 is steam diffusion coefficient [15].Considering the effects of solute in the unsaturated soil, the density of steamρv can be expressed as(18)ρv=hrρvs=Rhhcρvs.Substituting (16)~(18) into (15), the continuous equilibrium equation of pore water considering combined effects of framework deformation, steam, solute, and temperature can be obtained as follows:(19)Cpe∂ε∂t+Cpp∂pw∂t+CpT∂T∂t+Cpc∂c∂t=∇·Kpp∇pw+ρl2gKkrlμw+∇·KpT∇T+∇·Kpc∇c.The coefficients in the equation are shown in Appendix A.


## 2.3. Energy Conservation Equation

The total internal energy of every fluid phaseMφh and energy exchange flux qφh can be expressed as(20a)Mφh=nSφρφeφ,(20b)qφh=iφh+eφaqrφa+eφwqrφw+nSφρφeφu˙s=iφh+qrφh+nSφρφeφu˙s.Similarly, as for the solid phase, it can be obtained:(21a)Msh=1-nρses,(21b)qsh=ish+esqrs+1-nρsesu˙s=ish+qrsh+1-nρsesu˙s,where eφ is a unit mass of internal energy of phase φ; eφa and eφw are energies of water and gas per unit mass.According to the total energy conservation and the assumption of local thermal equilibrium, it can be obtained:(22)∂∂nSlρlel+nSgρvev+1-nρses-Qh=-∇·imh+qrlh+qrgh+qrsh,where imh is the average heat conduction flux; qrlh,  qrgh, and qrsh are convection fluxes of the liquid phase, gas phase, and solid phase, respectively; Qh is the term of source and convergence. Effects of thermal conduction in deformation of solid phases are neglected.When variations of temperature and pressure are not large enough, the internal energy in the energy balance equation can be replaced by the specific enthalpy. Specific enthalpies of the solid phase, water, and water steam can be expressed as follows:(23a)hs+hl=CsT-To+ClT-To,(23b)hv=hpv+hLo=CpvT-T0+LoTo=ClT-T0+LT,where hs,  hl, and hv are specific enthalpies of the solid phase, water, and water steam, respectively; Cs and Cl are the specific heat of solids and water; Cpv is the specific heat at constant pressure of water steam; Lo is the latent heat of vaporization of water steam.Replacing the internal energy with enthalpy in (23a) and (23b),(24)∂hs∂t=1-nρsCs∂T∂t,(25)∂nSwρlhl∂t=nSwρlCl∂T∂t+nρlClT-T0∂Sw∂t+nSwClT-T0∂ρl∂t,(26)∂nSgρvhpv∂t=nSgρvCpv∂T∂t+nSgCpvT-T0∂ρv∂t+nρvCpvT-T0∂Sg∂t,(27)∂nSgρvhL0∂t=nSgL0∂ρv∂t+nρvL0∂Sg∂t.In the above derivation, effect of the density of solids on energy is neglected, and it is assumed that the specific heat can be considered as the invariant when variations of pressure and temperature are not large.The right part of the energy equation (22) can be abbreviated as follows:(28)qh=-λI∇T+ClT-T0qrl+ClT-T0+LTqrv.By substituting (23a), (23b), and (27) into (22), the energy conservation equation can be obtained as(29)CTp∂pw∂t+CTT∂T∂t+CTc∂c∂t+VTT·∇T=∇·KTp∇pw+∇·KTT∇T+∇·KTc∇c.The coefficients in the equation are shown in Appendix B.


## 2.4. Solute Transport Equation

When the temperature does not vary significantly, it can be considered that migration of solute under the temperature gradient can be neglected; therefore the solute transport equation can be obtained according to the mass conservation law of the solute:(30)∂∂tθc+∂∂tρdcP=∇·Dsh∇c-qc,where cP is the amount of adsorbed solute in soil per mass, which is a function of concentration. Then (30) can be expressed as(31)nSw+ρd∂cp∂c∂c∂t+cDs∂pw∂t=∇·Dsh∇c-qc,where θ is the volumetric water content of pore water; ρd is the density of dry soil; Dsh is the hydrodynamic dispersion coefficient of the solute; q is the water flux; Ds is the moisture capacity.


## 2.5. Initial Conditions and Boundary Conditions

Initial conditions inside the given analysis areaΩ and on the boundary Γπ(π=u,w,T,c) when t is zero are as follows:(32)u=u0,pw=pw0,T=T0,c=c0,onΩandΓπ.Dirichlet boundary conditions can be expressed as follows:(33)ut=u^t,onΓu;pwt=p^wt,onΓw;Tt=T^t,onΓT;ct=c^t,onΓc.Neumann boundary conditions can be expressed as follows:(34a)ITσ=t-,onΓuq;Kkrlμw∇ψ+DvτSa∇ρv·n=qw+qv+βcρv-ρv∞,onΓwq,(34b)ρwvwL-λM∇T·n=qT+αcT-T∞,onΓTq;Dsh∇c·n=qc,onΓTq,where qw,  qv,  qT,  qc, and t- are assigned liquid water flow, water steam flow, solute flow, and tensile force on the boundary Γπ(π=u,w,T,c). ρv∞ and T∞ are the mass density and temperature of water steam at undisturbed locations far away from boundaries; βc is the mass convection coefficient; αc is the thermal convection coefficient. Mixed boundary conditions may appear during the actual calculation. Therefore boundary conditions should be determined according to the actual situation.The initial concentration distribution is given as(35a)c=c0z,t=0,0≤z≤L,where L is the depth of calculated soil layer.At the earth surface wherez is zero, boundary conditions of solute migration can be one kind of the Dirichlet, Neumann, and Cauchy boundary conditions.(1) As the solute concentration at the earth surface is given, the Dirichlet boundary condition is expressed as follows:(35b)c=cz,z=0,t>0.(2) As the earth surface is in the status of infiltration and the solute concentration cRt or the degree of mineralization of rainfall or irrigation water is given, it is the Cauchy boundary condition expressed as follows since the solute transport flux J=-Dsh∇c+qc:(35c)-Dsh∇c+qc=qcRt,z=0,t>0.As the intensity of water supplyRt is smaller than the soil infiltration capacity, the water movement flux at the earth surface in the above equation is equal to the intensity of water supply Rt; as the intensity of water supply Rt is larger than the soil infiltration capacity, the water movement flux at the earth surface can be calculated through the above equation.(3) As the earth surface is in the status of evaporation, the intensity of evaporationEt adopts the Penman-Wilson formula of the actual evaporation amount at the surface of the unsaturated soil as the evaporation boundary condition:(35d)Et=ΔRn+ηEaΔ+ηA,where Δ is the slope of the vapor pressure-temperature curve of the saturated steam; Rn is the amount of net radiation at soil surface; η is the humidity constant; Ea is the amount of potential evaporation and Ea=fuPgwB-A. Here fu is a function of wind, and fu=0.35×1+0.15u, where u is the wind velocity; Pgw is the vapor pressure of air above the evaporating surface; B is the reciprocal of relative humidity of the air, which is 1/hB; A is the reciprocal of relative humidity of the soil surface, which is 1/Rh.Here it belongs to the Cauchy boundary condition, which can be expressed as follows since the solute flow is zero at the earth surface during evaporation:(35e)Dsh∇c+cEt=0,z=0,t>0.(4) During field or indoor experiments, as it has neither evaporation nor infiltration on the earth surface, it belongs to the Neumann boundary condition, and it can be obtained:(35f)∇c=0.


## 2.6. Finite Element Discretization of the THMC Coupled Governing Equation

Analytical solutions of the THMC coupled mathematical model can be hardly achieved due to the complexity of solution conditions and irregularity of solution region. Therefore, the only way to solve the model is to obtain approximate solutions through numerical methods. The THMC coupled governing equations such as (5), (19), (29), and (31) are spatial discretized as follows:(36)Kdu-dt+Cdp-wdt+IdT-dt=dfdt,Edu-dt+Fp-w+Gdp-wdt+MT-+LdT-dt+Udc-dt=f^,Pdp-wdt+QT-+Jp-w+RdT-dt+Sdc-dt=X,Wc-+Zdc-dt=Y,where(37)K=∫ΩBTDTBdΩ,C=-∫ΩBTm1-DT3KsSw+p-wDsn+DTmβswDsnNpdΩ,I=-∫ΩBTDTmβs3NTdΩ,df=∫ΩNuTdbdΩ+∫ΓNuTdt^dΓ,E=∫ΩNpTSwρw1+γc+1-SwρvmT-mTDT3KsBdΩ,G=∫ΩNpTDsρw1+γc-ρv+n1-Sw∂ρv∂pw+nSw1+γcρwKw+Swρw1+γc+1-Swρv·1-nKw-mTDTm3Ks2·NpdΩ,F=∫Ω∇NpTρlKkrlμw+DvτSaρvMwρlRTNpdΩ,M=∫Ω∇NpT-ρlKkrlμwυϕCRg+DvτSahr∂ρvs∂T+ρgwshc∂Rh∂T+ρvsRh∂hc∂TNTdΩ,L=∫ΩNpTn1-Sw∂ρv∂T-Swρw1+γc+1-Swρv·1-nβs-mTDTm3Ksβs3-nSwρw1+γcβwNTdΩ,U=∫ΩNpTn1-Sw∂ρv∂c+nSwρwγNcdΩ,f^=-∫ΓwNpTqw+qv+βcρv-ρv∞dΓ,P=∫ΩNTTρwCl1+γcT-T0Ds+nSwKw-ρvDsCpvT-T0-L0+n1-Sw∂ρv∂pwCpvT-T0+L0NpdΩ,Q=∫Ω∇NTTλM+LDTVNTdΩ,J=∫Ω∇NTTLDPVNpdΩ,R=∫ΩNTTρcM-nSwClT-T01+γcρwβw+n1-SwCpvT-T0+L0∂ρv∂TNTdΩ,S=∫ΩNTTnSwClT-T0ρwγ+n1-SwCpvT-T0+Lo∂ρv∂cNcdΩ,X=-∫ΓTNTTqT+αcT-T∞dΓ,W=∫Ω∇NcTDsh∇NcdΩ+∫Ω∇NcTqNcdΩ,Z=∫ΩNcTnSw+ρb∂cp∂cNcdΩ,Y=-∫ΓcNcTqcdΓ.Equation (36) can be expressed in the form of matrix as follows:(38)KCI0EGLU0PRS000Zddtu-p-wT-c-+00000FM00JQ0000Wu-p-wT-c-=dfdtf^XY.The Monolithic Augmentation Approach is applied to solve the model. Changes within every time stepΔt are assumed as linear, which means(39)u-p-wT-c-tk+α=N1tN2tu-tkp-wtkT-tkc-tku-tk+Δtkp-wtk+ΔtkT-tk+Δtkc-t+Δtk,where α is a point within the time step Δtk; the item at the left of the equation contains the values of u-,  p-w,  T-, and c- at the time of tk+α, which means α=t-tk/Δtk; N1=1-α and N2=α.Therefore (38) can be transformed into the following iterative solution equation:(40)KCI0EG+αFΔtkL+αFΔtkU0P+αJΔtkR+αQΔtkS000Z+αWΔtktk+αu-p-wT-c-tK+ΔtK=KCI0EG-1-αFΔtkL-1-αMΔtkU0P-1-αJΔtkR-1-αQΔtkS000Z-1-αWΔtktk+αu-p-wT-c-tK+dfdtf^XYtk+αΔtK.The above equation is the spatial discretization of nonlinear equations. Hence the nonlinear governing equation is transformed into linear problem through discretization of space and time. The asymmetric stiffness matrix can be directly solved by using the asymmetric solver.


## 3. Assessment and Validation of the Numerical Model

### 3.1. Simulation of Salt Leaching and Accumulating

Combined with the THMC coupled governing equation and initial conditions and boundary conditions (34a)~(34b) and (35a)~(35f), as well as the given initial conditions and boundary conditions of water flow and temperature, the mathematic model of solute migration problem in the soil treated by the lime in the subgrade was then built.We choose a calculation example to study the salt leaching and accumulating processes in field to investigate the migration rule of salinity in the unsaturated soil, and the applicability and accuracy of the model are verified: dynamic test was developed to study the changes of the water and salt of the soil in the condition of fresh water and salt water flushing [16]. The test was conducted on a plot (6.1 m × 6.1 m). The average water content of the soil was 0.2. The observation and sample point was set in the depth of 30, 60, 90, 120, 150, and 180 cm, respectively. Based on the test, the calculated results of the THMC model we raised will be contrasted with the test results, and the model will be verified.According to the observation data, Bresler [17] has calculated the experimental process above, and the soil water characteristic curve is generalized as follows:(41)hθ=-1300e-10.5θcm,θ≤0.35,hθ=-1.59×107e-36.6θcm,θ≤0.35.The relationship between hydraulic conductivity and water content is(42)Kθ=3.24×10-8e35.8θcm/min.During calculation, the hydrodynamic dispersion coefficientDshν,θ adopts the following equation:(43)Dshν,θ=λν+Daaebθ,where the value of Da is 0.04 cm2/min; a ranges from 0.001 to 0.005; b equals 10; λ ranges from 0.28 to 0.55 and equals 0.28, 0.39, and 0.55 in the calculation, respectively.Since there are no calculating parameters such as temperature and volumetric deformation, the model was simplified as the water-solute coupled model, and the dynamic of soil water and salt was numerically simulated using the Bresler hydraulic characteristic parameters. Figure1 showed the comparison between numerically calculated results and measured results, which indicated that distributions of water content and concentration were consistent with measured data. After two hours of irrigation using the saline solution, salinity migrated downward with water at a depth up to 30 cm, and the concentration decreased when the depth increased from zero (at the earth surface) to 30 cm. Then the salt was leached using the fresh water. After nine hours and eleven hours of salt leaching, wetting fronts occurred at depths of 35 cm and 110 cm, respectively, while solute fronts (the maximal concentration) appeared at depths of 40 cm and 50 cm, with peak values decreased from 210 to 170 and 150, respectively. After thirteen hours of salt leaching, solute front appeared with a depth of 60 cm with a peak value of 140. After seventeen hours of salt leaching, wetting front appeared at a depth of 170 cm.Figure 1
The comparison between numerically simulated and measured results of distributions of water content and concentration in the soil.
(a)
The comparison between simulated and measured results of water content (b)
Simulated results of solute distribution distribution under different values ofλ (c)
Comparison between simulated and measured results of solute distribution after 9 hours of salt leaching (d)
Comparison between simulated and measured results of solute distribution after 11 hours of salt leachingIt is assumed that the underground water level is deep, and the initial salt concentration of the surface soil is up to 130 g/L. The initial salt concentration was shown in Figure1. Then the surface soil was leached using fresh water for 2000 minutes. After stopping the irrigation, the surface soil was evaporated according to the following rules:(44)Et=Eo=0.65cm/d,θ>0.27,Et=Eoaθ+b=2θ+0.11cm/d,θ≤0.27,⁡where E(t) is the evaporation intensity, which is not a constant and relates to the volumetric water content of the surface soil θ.In order to study the process of salt accumulating during evaporation, in this research, initial conditions and boundary conditions provided by Zhang [18] were adopted to calculate the above problem. Distributions of solute and volumetric water content at 250 min, 500 min, 700 min, 3000 min, and 21 days, separately, and the calculation results were showed in Figure 2.Figure 2
Numerically simulated results of salt leaching and accumulating.
(a)
Distributions of concentration at different time (b)
Distributions of water content at different timeFigure2 shows that, after being leached for 500 minutes, the concentration peak moved from the earth surface to the location with a depth of about 50 cm, and the maximum concentration decreased from 130 g/L to 25 g/L. After stopping the leaching, salt content of the surface soil gradually increased. However, due to the deep level of underground water, salinity continued to migrate downward with the downward infiltration of water. With the decreasing of water content, the evaporation intensity of surface soil rapidly reduced; therefore the salt accumulating became slow. When it was 21 days after leaching, salt content concentration in surface soil obviously rose to 20 g/L, which was because of the fact that when the water content in surface soil was very low, salt concentration had an obvious variation when the salt content was not large.Simulation results show that processes of salt leaching and accumulating were different from each other. During the 2000 minutes of salt leaching, the solute front kept moving downward, and its peak value kept decreasing. After stopping the salt leaching and following evaporating for 21 days, the amount of evaporation was determined through (44), which gradually reduced with the decreasing of water content in the surface soil, while the solute front did not move upward obviously due to the evaporation. Simulation results indicate that processes of salt leaching and accumulating were different. Relative to the process of salt accumulating, the process of salt leaching affected the downward migration of salinity more, which was at some extent related to the fact that evaporation amount reduced with the decreasing of water content in surface soil.


### 3.2. Numerical Simulation of Migration of Ca2+ Ions in the Soil Subgrade Treated by the Lime under Conditions of Rainfall and Evaporation

The adsorption process in soil is usually described by the adsorption isotherm, which is an empirical equation describing relationship between the adsorption amount of solute and the solute concentration in the solution under the constant temperature. Comparison study shows that, among the adsorption isotherms (linear equation, Freundlich equation and Langmuir equation, etc.), the Freundlich adsorption isotherm has a better effect in fitting the experimental data. Results of fitting experiments are shown in Figure3. The expression of adsorption isotherm is as follows:(45)S=KfCN,where C is the solute concentration; Kf and N are fitting parameters, with values of 0.1522 and 0.5848, respectively.Figure 3
Fitting of the Freundlich adsorption isotherm of Ca2+ ions.According to quantitative analysis results of mineral compositions using the X-ray diffraction, it can be obtained that the effective montmorillonite content was about 6%~18%, and its average value of 12% was taken as the reference value. Fitting parametersKf and N of the adsorption isotherm of the expansive soil in Nanning were 0.0228 and 0.5848, respectively.Situations in temperature variation and meteorological parameters under conditions of rainfall and evaporation were listed in Tables1 and 2. One assumes that, under the eluviations of atmosphere, shear strength and elasticity modulus of the soil treated by the lime share the same attenuation law, which is(46)EcE0=τcτ0,where E is the elasticity modulus of soil treated by the lime, with a unit of kPa; c is the solute concentration, with a unit of kg/kg. The diffusion coefficient of Ca2+ ion in the expansive soil was selected as 3.8 × 10−6 cm2/s according to references. The soil water characteristic curve is showed in Figure 4. According to results of the compaction test of soil treated by the lime, control indexes of the subgrade filling were 95% of the maximum dry density and 2% larger than optimal moisture, which means the initial volumetric water content of the subgrade was 0.31, and the saturability was 86%. Indoor experiment tested that the content of Ca2+ ions was 4.803%. According to the fitting curve of the adsorption isotherm shown in Figure 3, the equilibrium concentration was 104.5 kg/m3 when the solution was in the equilibrium.Table 1
Calculation parameter.
Volume specific heat coefficient(J·(m3·°C)−1) Thermal conductivity(J·(s·m·°C)−1) Permeability coefficient/m·s−1 Elasticity modulusMPa Poisson ratio Bulk modulus(kPa) Thermal expansion coefficient(K−1) C s C w λ s λ w k s E μ K w K s β w β s 2.24 × 106 4.15 × 106 6.0 0.57 1.25 × 109 60 0.3 4.3 × 109 1.4 × 109 6.3 × 105 0.9 × 106Table 2
Weather data of some region in China.
Temperature/°C Relative humidity% Wind speedm/s Net radiationMJ/(m2·d) Rainfall MJ/mm Highest Lowest Average Highest Lowest Average 35.04 19.86 24.68 99.99 46.74 73.37 0.482 14.2 0Figure 4
The soil water characteristic curve of soil treated by the lime.According to simulation and calculation of rainfall, infiltration and evaporation of soil subgrade treated by the lime, trends of subgrade deformation, and distributions of the volumetric water content in soil and Ca2+ ion content along the depth were obtained, as shown in Figures 5 and 6. As can be seen from Figure 5, under the action of rainfall, water continuous seeping into the soil layer of the subgrade, and rain for 20 days later, the impact depth is about 1.5 m. With the evaporation of the atmosphere, the surface of subgrade soil is constantly losing water, and the influence depth of evaporation for 20 days is about 1.0 m. From the calculation result of Figure 6, in the action of the eluviation, Ca2+ ions migrated downward. After continuous rainfall, in the region within 10 cm from the surface soil, Ca2+ ion content reduced significantly and accumulated downward due to the absorption of Ca2+ ions by the expansive soil and the low diffusivity of Ca2+ ions in the expansive soil. When evaporation began, the soil layer contained a lot of Ca2+ ions, and therefore salt accumulation occurred at the soil surface under the condition of continuous evaporation. Relative to the effect of salt leaching caused by rainfall, the effect of salt accumulating caused by evaporation was smaller.Figure 5
Distributions of water content under conditions of continuous rainfall and evaporation.
(a)
Profile evolution in distributions of water content under the condition of rainfall content (b)
Profile evolution in distributions of water content under the condition of evaporationFigure 6
Distributions of Ca2+ ion content under conditions of continuous rainfall and evaporation.
(a)
Profile evolution in distributions of Ca2+ ion content under the condition of rainfall (b)
Profile evolution in distributions of Ca2+ ion content under the condition of evaporation


## 3.1. Simulation of Salt Leaching and Accumulating

Combined with the THMC coupled governing equation and initial conditions and boundary conditions (34a)~(34b) and (35a)~(35f), as well as the given initial conditions and boundary conditions of water flow and temperature, the mathematic model of solute migration problem in the soil treated by the lime in the subgrade was then built.We choose a calculation example to study the salt leaching and accumulating processes in field to investigate the migration rule of salinity in the unsaturated soil, and the applicability and accuracy of the model are verified: dynamic test was developed to study the changes of the water and salt of the soil in the condition of fresh water and salt water flushing [16]. The test was conducted on a plot (6.1 m × 6.1 m). The average water content of the soil was 0.2. The observation and sample point was set in the depth of 30, 60, 90, 120, 150, and 180 cm, respectively. Based on the test, the calculated results of the THMC model we raised will be contrasted with the test results, and the model will be verified.According to the observation data, Bresler [17] has calculated the experimental process above, and the soil water characteristic curve is generalized as follows:(41)hθ=-1300e-10.5θcm,θ≤0.35,hθ=-1.59×107e-36.6θcm,θ≤0.35.The relationship between hydraulic conductivity and water content is(42)Kθ=3.24×10-8e35.8θcm/min.During calculation, the hydrodynamic dispersion coefficientDshν,θ adopts the following equation:(43)Dshν,θ=λν+Daaebθ,where the value of Da is 0.04 cm2/min; a ranges from 0.001 to 0.005; b equals 10; λ ranges from 0.28 to 0.55 and equals 0.28, 0.39, and 0.55 in the calculation, respectively.Since there are no calculating parameters such as temperature and volumetric deformation, the model was simplified as the water-solute coupled model, and the dynamic of soil water and salt was numerically simulated using the Bresler hydraulic characteristic parameters. Figure1 showed the comparison between numerically calculated results and measured results, which indicated that distributions of water content and concentration were consistent with measured data. After two hours of irrigation using the saline solution, salinity migrated downward with water at a depth up to 30 cm, and the concentration decreased when the depth increased from zero (at the earth surface) to 30 cm. Then the salt was leached using the fresh water. After nine hours and eleven hours of salt leaching, wetting fronts occurred at depths of 35 cm and 110 cm, respectively, while solute fronts (the maximal concentration) appeared at depths of 40 cm and 50 cm, with peak values decreased from 210 to 170 and 150, respectively. After thirteen hours of salt leaching, solute front appeared with a depth of 60 cm with a peak value of 140. After seventeen hours of salt leaching, wetting front appeared at a depth of 170 cm.Figure 1
The comparison between numerically simulated and measured results of distributions of water content and concentration in the soil.
(a)
The comparison between simulated and measured results of water content (b)
Simulated results of solute distribution distribution under different values ofλ (c)
Comparison between simulated and measured results of solute distribution after 9 hours of salt leaching (d)
Comparison between simulated and measured results of solute distribution after 11 hours of salt leachingIt is assumed that the underground water level is deep, and the initial salt concentration of the surface soil is up to 130 g/L. The initial salt concentration was shown in Figure1. Then the surface soil was leached using fresh water for 2000 minutes. After stopping the irrigation, the surface soil was evaporated according to the following rules:(44)Et=Eo=0.65cm/d,θ>0.27,Et=Eoaθ+b=2θ+0.11cm/d,θ≤0.27,⁡where E(t) is the evaporation intensity, which is not a constant and relates to the volumetric water content of the surface soil θ.In order to study the process of salt accumulating during evaporation, in this research, initial conditions and boundary conditions provided by Zhang [18] were adopted to calculate the above problem. Distributions of solute and volumetric water content at 250 min, 500 min, 700 min, 3000 min, and 21 days, separately, and the calculation results were showed in Figure 2.Figure 2
Numerically simulated results of salt leaching and accumulating.
(a)
Distributions of concentration at different time (b)
Distributions of water content at different timeFigure2 shows that, after being leached for 500 minutes, the concentration peak moved from the earth surface to the location with a depth of about 50 cm, and the maximum concentration decreased from 130 g/L to 25 g/L. After stopping the leaching, salt content of the surface soil gradually increased. However, due to the deep level of underground water, salinity continued to migrate downward with the downward infiltration of water. With the decreasing of water content, the evaporation intensity of surface soil rapidly reduced; therefore the salt accumulating became slow. When it was 21 days after leaching, salt content concentration in surface soil obviously rose to 20 g/L, which was because of the fact that when the water content in surface soil was very low, salt concentration had an obvious variation when the salt content was not large.Simulation results show that processes of salt leaching and accumulating were different from each other. During the 2000 minutes of salt leaching, the solute front kept moving downward, and its peak value kept decreasing. After stopping the salt leaching and following evaporating for 21 days, the amount of evaporation was determined through (44), which gradually reduced with the decreasing of water content in the surface soil, while the solute front did not move upward obviously due to the evaporation. Simulation results indicate that processes of salt leaching and accumulating were different. Relative to the process of salt accumulating, the process of salt leaching affected the downward migration of salinity more, which was at some extent related to the fact that evaporation amount reduced with the decreasing of water content in surface soil.


## 3.2. Numerical Simulation of Migration of Ca2+ Ions in the Soil Subgrade Treated by the Lime under Conditions of Rainfall and Evaporation

The adsorption process in soil is usually described by the adsorption isotherm, which is an empirical equation describing relationship between the adsorption amount of solute and the solute concentration in the solution under the constant temperature. Comparison study shows that, among the adsorption isotherms (linear equation, Freundlich equation and Langmuir equation, etc.), the Freundlich adsorption isotherm has a better effect in fitting the experimental data. Results of fitting experiments are shown in Figure3. The expression of adsorption isotherm is as follows:(45)S=KfCN,where C is the solute concentration; Kf and N are fitting parameters, with values of 0.1522 and 0.5848, respectively.Figure 3
Fitting of the Freundlich adsorption isotherm of Ca2+ ions.According to quantitative analysis results of mineral compositions using the X-ray diffraction, it can be obtained that the effective montmorillonite content was about 6%~18%, and its average value of 12% was taken as the reference value. Fitting parametersKf and N of the adsorption isotherm of the expansive soil in Nanning were 0.0228 and 0.5848, respectively.Situations in temperature variation and meteorological parameters under conditions of rainfall and evaporation were listed in Tables1 and 2. One assumes that, under the eluviations of atmosphere, shear strength and elasticity modulus of the soil treated by the lime share the same attenuation law, which is(46)EcE0=τcτ0,where E is the elasticity modulus of soil treated by the lime, with a unit of kPa; c is the solute concentration, with a unit of kg/kg. The diffusion coefficient of Ca2+ ion in the expansive soil was selected as 3.8 × 10−6 cm2/s according to references. The soil water characteristic curve is showed in Figure 4. According to results of the compaction test of soil treated by the lime, control indexes of the subgrade filling were 95% of the maximum dry density and 2% larger than optimal moisture, which means the initial volumetric water content of the subgrade was 0.31, and the saturability was 86%. Indoor experiment tested that the content of Ca2+ ions was 4.803%. According to the fitting curve of the adsorption isotherm shown in Figure 3, the equilibrium concentration was 104.5 kg/m3 when the solution was in the equilibrium.Table 1
Calculation parameter.
Volume specific heat coefficient(J·(m3·°C)−1) Thermal conductivity(J·(s·m·°C)−1) Permeability coefficient/m·s−1 Elasticity modulusMPa Poisson ratio Bulk modulus(kPa) Thermal expansion coefficient(K−1) C s C w λ s λ w k s E μ K w K s β w β s 2.24 × 106 4.15 × 106 6.0 0.57 1.25 × 109 60 0.3 4.3 × 109 1.4 × 109 6.3 × 105 0.9 × 106Table 2
Weather data of some region in China.
Temperature/°C Relative humidity% Wind speedm/s Net radiationMJ/(m2·d) Rainfall MJ/mm Highest Lowest Average Highest Lowest Average 35.04 19.86 24.68 99.99 46.74 73.37 0.482 14.2 0Figure 4
The soil water characteristic curve of soil treated by the lime.According to simulation and calculation of rainfall, infiltration and evaporation of soil subgrade treated by the lime, trends of subgrade deformation, and distributions of the volumetric water content in soil and Ca2+ ion content along the depth were obtained, as shown in Figures 5 and 6. As can be seen from Figure 5, under the action of rainfall, water continuous seeping into the soil layer of the subgrade, and rain for 20 days later, the impact depth is about 1.5 m. With the evaporation of the atmosphere, the surface of subgrade soil is constantly losing water, and the influence depth of evaporation for 20 days is about 1.0 m. From the calculation result of Figure 6, in the action of the eluviation, Ca2+ ions migrated downward. After continuous rainfall, in the region within 10 cm from the surface soil, Ca2+ ion content reduced significantly and accumulated downward due to the absorption of Ca2+ ions by the expansive soil and the low diffusivity of Ca2+ ions in the expansive soil. When evaporation began, the soil layer contained a lot of Ca2+ ions, and therefore salt accumulation occurred at the soil surface under the condition of continuous evaporation. Relative to the effect of salt leaching caused by rainfall, the effect of salt accumulating caused by evaporation was smaller.Figure 5
Distributions of water content under conditions of continuous rainfall and evaporation.
(a)
Profile evolution in distributions of water content under the condition of rainfall content (b)
Profile evolution in distributions of water content under the condition of evaporationFigure 6
Distributions of Ca2+ ion content under conditions of continuous rainfall and evaporation.
(a)
Profile evolution in distributions of Ca2+ ion content under the condition of rainfall (b)
Profile evolution in distributions of Ca2+ ion content under the condition of evaporation


## 4. Conclusions

(1) The processes of salt leaching and accumulating in the soil treated by the unsaturated lime were different. Comparing with the process of salt accumulating, the influence of the process of salt leaching to the downward migration of salinity was larger, which was closely related to the fact that evaporation amount reduced to the decreasing of water content in soil surface.(2) Under the action of eluviation, the Ca2+ ion content in deeper soil increased due to the Ca2+ ion content accumulating and moving down, and the influence depth reached 1.5 m. Under the condition of continuous evaporation, salt accumulation occurred on the surface of the expensive soil subgrade. Comparing with the effect of salt leaching caused by rainfall, the effect of salt accumulating caused by evaporation was smaller.(3) Comparing the experimental results with the theoretical results, the following can be concluded: the thermo-hydro-mechano-chemical (THMC) model which was built based on the basic theories in the polyporous polyphasic medium mechanics, the mixture theory in the continuum mechanics, and thermodynamic theories of the irreversible process can accurately calculate and solve the coupling problem of solute migration during the treatment of expansive soil using the unsaturated lime.



---

*Source: 4089612-2018-03-21.xml*