Abstract
The moisture and the temperature variation of shallow loess during the freeze-thaw process have tremendous influences on the freeze-thaw disease and agricultural development. The present paper, based on the field test of the temperature and water content of the shallow loess in the seasonal frozen region, probes into the variation of the temperature and moisture fields of shallow loess in the freeze-thaw period. Field test shows that the daily variation of surface temperature reaches the maximum with the peak temperature appearing around 14:00, and the time when the peak temperature appears gradually postpones as the depth increases. Moreover, the daily variation of soil temperature decreases with the increase of temperature, and the variation can be ignored given the depth of more than 30 cm. In addition, there is an obvious variation of moisture content in shallow soil under the action of freeze-thaw cycles. In this regard, the coupled calculation model is proposed in this paper to simulate the dynamic variation of temperature and moisture in shallow soil in Bin County located in Western China under natural conditions, compare the calculation value with the measurement value, and further verify the rationality of the numerical calculation method and parameter selection. Furthermore, the variation law of the maximum frozen depth with the air temperature is obtained through the simulation of extreme cold climate conditions.
1. Introduction
Most of the loess areas in China, in essence, belong to the seasonal frozen ground region, and the temperature and moisture fields vary with seasons under the influence of natural factors. Especially in winter, the freezing action leads to the migration of moisture to the freezing layer and the increase of moisture content in shallow soil, which accounts for the frequent occurrences of engineering diseases, such as the slope, subgrade collapse, spalling, the frost heave, and thawing settlement of the Loess Plateau [1–5], and affects the agricultural development and the afforestation of mountains and rivers [6–9]. Of a large number of relevant studies, many scholars point out that the diseases are closely related to the variation of the temperature and moisture fields of shallow loess. Based on numerous field investigations, the frost heave and thawing settlement of the subgrade in G214 Expressway are mainly caused by the melting of permafrost [10]. Zhang et al. [11] investigated the law of moisture variation and transfer of the highway subgrade in the seasonal frozen ground region, which proves to be the main factor for the road frost damage. Ye et al. [12], according to the field investigation and experimental study, maintained that the deterioration of loess strength, caused by the moisture transfer under the freeze-thaw effect, plays a vital part in the slope and spalling of loess. Drawing on the numerical simulation analysis, Zhan et al. [13] proposed the close correlation of slope disease in the seasonal frozen ground region with moisture variation. In this regard, it is of great necessity to study the law of temperature and moisture variation of shallow soil and further reveal the occurrence mechanism of engineering diseases.
In recent years, a large number of studies have been conducted in terms of the coupled mechanism of the temperature and moisture in soil based on the laboratory test and numerical simulation. Mao et al. [14], taking aeolian sand as the research object, investigated into the law of moisture transfer with the temperature gradient. Bing et al. [15] investigated the law of moisture transfer in saline soil during the freeze-thaw cycle. Konrad and Morgenstern [16] and Gilpin et al. [17] explored the moisture transfer law and ice formation mechanism in frozen soil with the external load. The temperature and moisture in soil are redistributed during the freezing process (Konrad [18]; Xu [19]), and the distribution law is affected by external environment, initial water content, soil category, soil pore structure, surface vegetation, and so on [20, 21, 22, 23]. The first unsaturated moisture-heat coupling calculation model was established by Harlan [24] and gradually developed into the hydrodynamic model based on the modification by Taylor and Luthin [25]. However, it is difficult to determine the calculation parameters of the model and describe the moisture transfer process [26]. Konrad and Morgenstern [27] concluded that the moisture transfer flux is positively related to the temperature gradient and proposed the segregation potential model according to the experimental test. Zhou et al. [28] put forward the moving-pump model of water transfer in unsaturated soil, which effectively avoids the repartition and recalculation of the moving boundary and simplifies the numerical calculation of the theoretical model. In addition, the method is used to analyze the one-dimensional freezing process and compare the calculated value with the measured value for the verification of rationality. Xu et al. [29] introduced the impedance factor and proposed an improved Harlan model. Moreover, the numerical simulation is carried out for the moisture-heat coupling transmission within the frozen saline loess column, and the comparison between the calculated and measured values verifies the rationality of the method.
At present, most studies concerning the law of temperature and moisture under the freezing action focus on the indoor test with few conducted based on the field test. In this respect, the present paper selects three typical sites in the Loess Plateau in Northern Shaanxi, China, to investigate the variation of temperature field and moisture field of shallow soil in the freeze-thaw cycle and further establish a numerical moisture-heat coupling calculation model for the simulation of natural conditions. Based on the comparison between measured value and numerical value, this paper also verifies the rationality of the calculation model and parameter in an attempt to lay a solid foundation for the exploration of the mechanism of engineering freeze-thaw diseases in the seasonal frozen soil area.
2. Field Test Scheme
2.1. General Situation of Test Sites
In the field test, one test site is selected in Bin County, Tongchuan, and Luochuan, respectively. The test site in Bin County is the flat land on the Loess Plateau with the longitude of E108°4′52″ and the latitude of N35°10′48″ covered with grass on the surface. As for the test site in Luochuan, the platform at the foot of the shady slope on the side of the Loess Plateau, the front zone of gully, is selected with the longitude of E109°24′25″ and the latitude of N35°51′35″. The test site in Tongchuan is a front zone of gully with the longitude of E109°06′57″ and the latitude of N35°7′49″. Of the three sites, a similar vegetation rate is shown, and the shallow soil of the sites is the silty clay, yellowish brown, with a deep level of groundwater. The basic physical indexes are presented in Table 1.
2.2. Test Method
Figure 1 presents the layout of temperature sensors. Temperature sensors are embedded in advance every 10 cm within 1 m of the soil layer and, the temperature for measurement ranges from −30°C to + 85°C with the accuracy of 0.1°C. In an attempt to better reflect the temperature variation of shallow soil in the freeze-thaw period, the test time of ground temperature is from December 15th to March 10th of the next year, and the daily temperature variation is also tested. Soil samples, taken from the excavated pit around the test point at different depths, are used to carry out the experiment for water content in order to investigate the variation law of moisture field of shallow soil.

3. Result and Analysis of the Field Test
3.1. Temperature Field Variation of Shallow Soil
According to the test scheme above, the variation of shallow loess with time is presented in Figure 2–Figure 4. The temperature in Figure 2 is the average value of daily temperature.

(a)

(b)

(c)
As Figure 2 shows, the temperature of surface soil in winter decreases with time, and freezing depth gradually increases. The freezing depth in Bin County on January 10th is 26 cm, Luochuan on January 25th is 61 cm, and Tongchuan on January 26th is 28 cm. By March 10th, the temperature returns to the positive level, and the surface soil temperature also rises without the existence of frozen soil. However, the temperature is relatively low at the underground depth of 50 cm below in Bin County, 90 cm below in Luochuan, and 60 cm below in Tongchuan since the heat conduction of energy in soil requires time with the increase of surface soil temperature and the heat exchange between soil and atmosphere. Therefore, the temperature in the deeper soil continues to decrease under the original conditions, which indicates the hysteresis from the macroperspective.
Figure 3 and Figure 4 indicate that the variation of daily temperature gradually decreases with deeper soil deposit. The daily variation of temperature in surface soil reaches the maximum with the peak temperature appearing around 14:00, and the time when the peak temperature appears gradually postpones as the depth increases. With the depth of more than 20 cm, the daily temperature varies slightly, and the temperature variation almost reaches 0 given the depth of more than 30 cm. Under such circumstances, the daily temperature variation can be neglected with the only consideration of the temperature variation with seasons. In Figure 3(b) and Figure 4(b), crossover phenomenon appears in the temperature curve within 10 cm to 20 cm below the surface. In other words, when the ground temperature increases to the maximum, the soil temperature decreases in the range of 10–20 cm, which is caused by the lag of the daily change of temperature in the process of heat conduction.

(a)

(b)

(c)

(a)

(b)

(c)
3.2. Moisture Field Variation of Shallow Soil
As Figure 5 indicates, there is a significant increase in the moisture content of the freezing layer at three sites, while an opposite trend is shown in the moisture content of unfrozen soil under the freezing layer. It can also be seen from Figure 3(a) and Figure 4(a) that, from December 15th to January 10th, the surface temperature decreases, the freezing front moves downward, and the thickness of the freezing layer increases. In addition, an increasing trend is shown in the moisture content of both whole freezing layer and freezing front since the suction force is formed on the freezing front, and the moisture in the unfroze area migrates to the front continuously, which leads to the increase of moisture content in the freezing layer when the soil temperature drops below the freezing point in winter. With regard to the unfrozen soil in the deeper parts, the moisture content moves upward under the influence of temperature gradient and moisture difference. However, the low moisture transfer speed fails to supply the upward part, which results in the decrease of moisture content in the unfrozen soil of the adjacent freezing layer. When the shallow soil melts with the increase of temperature in March, part of moisture in the original frozen soil will infiltrate downward, while part of it will evaporate upward. Owing to different evaporation capacities, the moisture content is also distributed differently in the shallow soil, which further accounts for the smaller evaporation capacities in Bin County and Tongchuan than those in Luochuan. Therefore, as shown in the curve of March 10th in Figures 5(a) and 5(c), the moisture content of the original frozen soil in Bin County and Tongchuan presents a decreasing trend, and the content increases at about 30 cm below the original frozen layer. Due to the large evaporation in Luochuan, most of melted and infiltrated moisture evaporates, so there is no obvious increase in the curve under the freezing layer.

(a)

(b)

(c)
4. Numerical Simulation of Moisture-Heat Coupling of Shallow Soil in Winter
4.1. Finite Element Equation and Parameter Determination of Temperature Field
Climate factors, such as radiation, evaporation, humidity, and wind speed, vary with time, under whose influence the temperature field of shallow soil belongs to the unsteady temperature field with phase change. The finite element equation is expressed as follows [30]:where [K] is the temperature stiffness matrix, [N] is the unsteady temperature changing matrix, {T} is the column vector of temperature value, Δt is the time step, {P} is the synthetic array, and subscript t is the time.
Based on the isoparametric quadrilateral element, the parameters of the stiffness matrix and the temperature changing matrix are determined as follows:where I, j, k, and m are the element nodes, λ is the thermal conductivity, ρCp is the specific heat capacity of volume, │J│, D1, D2, and H are the functions of the coordinate of integral base points (ζ,η) related to the element nodes, α is the convective heat transfer coefficient, and s is the length of the heat exchanging boundary. A is taken as 0 for the boundary of the nonthird kind. For the third kind of boundary, A is taken as 1/3 when l equals to n; otherwise, A is 1/6. The values of thermal conductivity and specific volume heat capacity are determined based on the research findings of Wang et al. [31] in the following formula:
With reference to the relevant literature [30], the array {P} is composed of three terms and calculated according to the following formula:where {P1} is the array with phase change, {P2}is the radiation heat exchange array, consisting of the solar shortwave radiation array {P2S}, the ground longwave radiation array {P2E}, and the atmospheric longwave radiation array {P2A}, {P3} is the convective heat transfer array, and {P4} is the evaporation heat consumption array:where l is the node of the boundary of the third kind, is the variation of solid phase ratio at the node l within the time interval , L is the latent heat of phase change, ωs and ωt are the weighted coefficients, ρ is the density, is the solar radiation volume per square meter of ground level, 1 is the reflectivity of ground to solar radiation, are the emissivity of earth and atmosphere, respectively, is the absorptivity of the earth to atmospheric radiation, T and Ta are the temperatures of the earth and atmosphere, respectively, U is the evaporation of surface soil, and G is the latent heat of moisture vaporization, which is taken as 2,459.1 kJ/kg.
According to the meteorological data, the solar radiation per square meter of the surface level can be obtained. Table 2 presents the average solar radiation per ten days in winter in Bin County based on the meteorological data in the past 30 years.
The emissivity of the earth is taken as 0.68, but it is complex to determine the values of the atmospheric emissivity and the radiation absorptivity of the earth to atmosphere , which are closely related to temperature, cloud amount, humidity, dust content, and so on. Temperature and humidity can reflect not only the amount of water vapor in the air but also the amount of cloud and air quality. In this regard, the present paper selects temperature and humidity as the characteristic indexes of climate for the determination of and . Based on the fitting of meteorologically measured natural surface temperature and the consideration of and in formula (11) as a whole, the expression of can be obtained as follows:where is the air temperature (°C), is the relative humidity (%), and is the regional coefficient considering the influence of other factors, which is taken as 0.25 in this paper.
The convective heat transfer coefficient consists of natural convection and forced convection, which can be calculated with the following formula:where the unit of wind velocity V is m/s.
The evaporation of the soil surface can be calculated as follows:where is the moisture content of surface soil and is the evaporation of the water surface.
Climatic parameters of Bin County in winter are presented in Table 3.
4.2. Finite Element Equation and Parameter Determination of Moisture Field
In winter, the moisture migration in shallow loess is unsteady under the freezing action with the finite element equation expressed as follows [32]:where [D] is the stiffness matrix, {F} is the flow vector array reflecting boundary conditions, [E] is the capacity matrix, {h} is the water head array, and Δt is the time step. With regard to the frozen areas, the permeability coefficient is taken as 0 owing to the small amount of moisture content, which can be neglected, and the calculation of water head has no significance for moisture transfer. For the unfrozen area, the value of water head is calculated according to the following formula:where h is the gravitational potential head, hu is the matric suction head, and hc is the phase change head.
hu is calculated according to formula (19) based on the existing literature [33]:where , , , , m = 1–1/n, and Tt is the temperature of soil; ρd is the dry density (g/cm3).
The phase change head is related to the density and moisture content of soil, which is expressed based on the research conclusion of the author according to the water transfer test of indoor frozen loess:where hc is the phase change head (m), ρd is the dry density (g/cm3), and Sr is the saturation.
4.3. Calculation and Result of Moisture-Heat Coupling
The calculation process of moisture-heat coupling is to couple the two equations, equations (1) and (17), and for each time step, no fixed values are adopted for hydrothermal parameters, such as thermal conductivity, specific heat capacity, latent heat with phase change, and water head. First of all, according to the initial moisture content, the temperature field is calculated within this time step, and then the calculation of moisture content is performed through equation (17) based on the temperature field above. After convergence, the process moves to the next time step, and the calculation results of the previous time step are used for the hydrothermal parameters, and so forth. In order to improve the convergence performance, the iteration method is employed for the calculation of moisture-heat coupling, and the matrix elimination method is used between the two iterations for the direct solution. Then, the next iteration value is determined by the low relaxation factor, modified by the following formula:where and are the ith iteration values, and are the values for the next iteration, and ω is the relaxation factor. The low relaxation iteration is adopted with ω < 1, and the ω value is adjusted according to the iteration times until the calculation results meet the required accuracy.
Based on the idea above, Fortran language is adopted to develop finite element calculation software. With the soil temperature and moisture in Bin County on December 15th as the initial conditions, this paper employs the temperature-moisture coupling model to calculate the variation of temperature and moisture in shallow soil, as shown in Figure 6 and Figure 7. As Figure 6(a) shows, the calculated value of the temperature field on January 10th and March 10th is basically consistent with the measured value, which is also true of the moisture field in Figure 7. In view of the calculation results, the rationality of the temperature-moisture coupling numerical calculation method is verified, which further validates the feasibility of the application of equation (20) concerning the water head at the phase change interface, obtained from the indoor test, in the sites. It can be found from Figure 6(b) that the maximum freezing depth is 42 cm, and an increasingly decrease in the surface soil temperature is shown from December. The temperature reaches its lowest level in late January and then begins to rise. By March 10th, the temperature in soil has returned to the normal level. Similar to the measured variation, the temperature in the deeper part decreases with the increase of air temperature owing to the lag of heat conduction.

(a)

(b)

To further probe into the effect of extreme temperature on the maximum freezing depth in Bin County, the temperatures of −8°C, −11°C, and −15°C in January are selected, respectively, for the calculation of the maximum freezing depth, as shown in Figure 8.

As Figure 8 indicates, the maximum freezing depth increases with the decrease of air temperature, but the increment gradually drops. The fitting of the relationship between the maximum freezing depth and air temperature is presented in equation (23) with the correlation coefficient of 0.9561. Therefore, equation (23) is feasible to fit the relationship between the temperature and the maximum freezing depth:where h is the maximum freezing depth (cm) and Ta is the absolute value of air temperature (°C), which is the air temperature at 1.5 m above the ground.
5. Conclusions
(1)With the time going by in winter, the temperature of surface soil decreases, and the freezing depth increases. Moreover, as the air temperature rises after winter, there is a lag of temperature variation in deeper parts. With the increase of depth, the daily temperature variation gradually narrows.(2)The daily temperature variation reaches the maximum with the peak value appearing around 14:00, and as the depth increases, the time of peak temperature moves backward.(3)As the temperature of surface soil decreases, the freezing front moves downward, the thickness of the freezing layer increases, and an increasing trend is also presented in the moisture content in the freezing layer. However, there is a significant decrease of moisture content in the unfrozen soil below the freezing layer.(4)The present paper establishes the numerical calculation model of moisture-heat coupling for the simulation of sites and calculates the temperature field and moisture field of shallow loess in Bin County. The calculation results are basically consistent with the measured data, which verifies the rationality of the calculation model and parameter selection, as well as the feasibility of the model in the simulation of the hydrothermal dynamic law of shallow loess in winter in the Loess Plateau.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by Shaanxi Technology Innovation Guidance Project (Fund) (Project no. 2019QYPY-148) and the National Natural Science Foundation of China (Project no. 41702346).