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.
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 increment and effective stress increment can be expressed in the incremental forms as follows, respectively:where is the tangent elastoplastic modulus matrix; is the total strain increment; is the strain increment caused by the pore pressure and ; is the strain increment caused by temperature, and ; is the strain increment caused by swelling, and ; is the unit array of normal stress; is the thermal expansion coefficient of solid particles; is the water swelling expansion coefficient; is the bulk modulus of solid particles; is the average stress of gases and liquids inside pores, and ; and are saturability of the pore water and pore gas, respectively, and ; and are pressures of the pore water and pore gas, respectively.
The stress equilibrium equation is stated as follows:where is the body force.
According to the principle of virtual work, the incremental stress equilibrium equation can be stated as follows:where is the total stress increment; and are increments of body force and surface force, respectively; and are virtual strain and displacement, respectively.
Equations (2) and (3) are substituted into (4), and the following is obtained:According to the law of conservation of mass, the fluid continuity equation can be obtained:Similarly, the equation of solid mass conservation can be obtained:Equation (7) is the mass conservation equation of the component in the phase. It is the liquid phase as is equal to . It is the gas phase as is equal to . The component is water as is equal to , and the component is the gas as is equal to . For instance, and indicate air resolved in the liquid phase and water in the liquid phase, respectively; and indicate dry air and water steam in the gas phase, respectively; is the decreasing or increasing rate of unit volume water during evaporation or condensation. and are absolute moving velocities of solid particles and fluid, respectively; is saturability of the water phase or the gas phase.
The fluid flux , which flows across the unit area per unit time, is stated as follows:where is the diffusion flow of component in the phase; is the relative flow when there are movements of the solid phase in the deformed porous medium, and the subscript indicates that the velocity and flow are relative to the solid phase. According to the definition of relative velocity, it can be obtained that . Equation (6) can be stated as follows:According to (6), (7), and the definition of material derivative, , it can be obtained thatSubstituting (10) into (11), the equilibrium equation of phase component in the deformed porous medium can be obtained as follows:Under the situation of small deformation, items such as and can be neglected; can be presented as the volumetric strain of the framework . Equation (12) can be simplified as follows: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, , the total equilibrium equation is as follows:As for the water including liquid water and steam in the air (), the total equilibrium equation is as follows: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, , , , and .
Variations of densities of steam, liquid water, and particles with time can be expressed aswhere ; is the density of fresh water; is the density ratio of solute to water in the pore water; is the compression modulus of fresh water; is the thermal expansion coefficient of fresh water; is the compression modulus of particles; is the thermal expansion coefficient of the framework; and 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:where is the effusion flux of water steam, with a unit of kg/m2/s; is the factor considering pore curvature in the soil [14], and is steam diffusion coefficient [15].
Considering the effects of solute in the unsaturated soil, the density of steam can be expressed as
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:The coefficients in the equation are shown in Appendix A.
2.3. Energy Conservation Equation
The total internal energy of every fluid phase and energy exchange flux can be expressed asSimilarly, as for the solid phase, it can be obtained:where is a unit mass of internal energy of phase ; and 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:where is the average heat conduction flux; , , and are convection fluxes of the liquid phase, gas phase, and solid phase, respectively; 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:where , , and are specific enthalpies of the solid phase, water, and water steam, respectively; and are the specific heat of solids and water; is the specific heat at constant pressure of water steam; is the latent heat of vaporization of water steam.
Replacing the internal energy with enthalpy in (23a) and (23b),
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:
By substituting (23a), (23b), and (27) into (22), the energy conservation equation can be obtained asThe 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:where is the amount of adsorbed solute in soil per mass, which is a function of concentration. Then (30) can be expressed aswhere is the volumetric water content of pore water; is the density of dry soil; is the hydrodynamic dispersion coefficient of the solute; is the water flux; is the moisture capacity.
2.5. Initial Conditions and Boundary Conditions
Initial conditions inside the given analysis area and on the boundary when is zero are as follows:
Dirichlet boundary conditions can be expressed as follows:
Neumann boundary conditions can be expressed as follows:where , , , , and are assigned liquid water flow, water steam flow, solute flow, and tensile force on the boundary . and are the mass density and temperature of water steam at undisturbed locations far away from boundaries; is the mass convection coefficient; 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 aswhere is the depth of calculated soil layer.
At the earth surface where 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:(2) As the earth surface is in the status of infiltration and the solute concentration 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 :
As the intensity of water supply 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 ; as the intensity of water supply 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 evaporation adopts the Penman-Wilson formula of the actual evaporation amount at the surface of the unsaturated soil as the evaporation boundary condition:where is the slope of the vapor pressure-temperature curve of the saturated steam; is the amount of net radiation at soil surface; is the humidity constant; is the amount of potential evaporation and . Here is a function of wind, and , where is the wind velocity; is the vapor pressure of air above the evaporating surface; is the reciprocal of relative humidity of the air, which is ; is the reciprocal of relative humidity of the soil surface, which is .
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:(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:
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:whereEquation (36) can be expressed in the form of matrix as follows:
The Monolithic Augmentation Approach is applied to solve the model. Changes within every time step are assumed as linear, which meanswhere is a point within the time step ; the item at the left of the equation contains the values of , , , and at the time of , which means ; and .
Therefore (38) can be transformed into the following iterative solution equation:
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:
The relationship between hydraulic conductivity and water content is
During calculation, the hydrodynamic dispersion coefficient adopts the following equation:where the value of is 0.04 cm2/min; ranges from 0.001 to 0.005; 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. Figure 1 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.
(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 leaching
It 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 Figure 1. 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:where 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.
(a) Distributions of concentration at different time
(b) Distributions of water content at different time
Figure 2 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 Figure 3. The expression of adsorption isotherm is as follows:where is the solute concentration; and are fitting parameters, with values of 0.1522 and 0.5848, respectively.
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 parameters and 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 Tables 1 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 iswhere is the elasticity modulus of soil treated by the lime, with a unit of kPa; 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.
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.
(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 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.
Appendix
A. The Coefficients of the Continuous Equilibrium Equation of Pore Water (See (21a) and (21b))
B. The Coefficients of the THMC Coupled Governing Equation (See (30))
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the Youth Innovation Promotion Association CAS, the outstanding youth fund of Hubei Province (2017CFA056), the Technology Service Network Initiative (no. KFJ-STS-ZDTP-037), and the National Natural Science Foundation of China (nos. 41472286, 41472290, and 41672312).