# Study on Varying Characteristics of Temperature Field and Moisture Field of Shallow Loess in the Freeze-Thaw Period

Authors: Hui Zhang; Bin Zhi; Enlong Liu; Tiehang Wang
Journal: Advances in Materials Science and Engineering
 (2020)

---

## 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.


---

## Body

## 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 Table1.Table 1
Basic physical index.
SitesDry density (g/cm3)Average moisture content (%)Liquid limit (%)Plastic limit (%)Plasticity indexLuochuan1.5321.230.7418.0312.74Tongchuan1.5018.529.7117.7511.96Bin County1.5521.530.918.112.8


### 2.2. Test Method

Figure1 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.Figure 1
Layout of temperature sensors.


## 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 Table1.Table 1
Basic physical index.
SitesDry density (g/cm3)Average moisture content (%)Liquid limit (%)Plastic limit (%)Plasticity indexLuochuan1.5321.230.7418.0312.74Tongchuan1.5018.529.7117.7511.96Bin County1.5521.530.918.112.8


## 2.2. Test Method

Figure1 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.Figure 1
Layout of temperature sensors.


## 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 Figure2–Figure 4. The temperature in Figure 2 is the average value of daily temperature.Figure 2
Temperature variation of shallow loess with time at different sites. (a) Bin County. (b) Luochuan. (c) Tongchuan.
(a)(b)(c)As Figure2 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.Figure3 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.Figure 3
Daily variation of ground temperature with depth in Bin County. (a) December 10th. (b) January 10th. (c) March 10th.
(a)(b)(c)Figure 4
Daily variation of temperature at different depths of the surface in Bin County. (a) December 10th. (b) January 10th. (c) March 10th.
(a)(b)(c)


### 3.2. Moisture Field Variation of Shallow Soil

As Figure5 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.Figure 5
Moisture content variation of shallow loess at different sites with depth. (a) Bin County. (b) Luochuan. (c) Tongchuan.
(a)(b)(c)


## 3.1. Temperature Field Variation of Shallow Soil

According to the test scheme above, the variation of shallow loess with time is presented in Figure2–Figure 4. The temperature in Figure 2 is the average value of daily temperature.Figure 2
Temperature variation of shallow loess with time at different sites. (a) Bin County. (b) Luochuan. (c) Tongchuan.
(a)(b)(c)As Figure2 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.Figure3 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.Figure 3
Daily variation of ground temperature with depth in Bin County. (a) December 10th. (b) January 10th. (c) March 10th.
(a)(b)(c)Figure 4
Daily variation of temperature at different depths of the surface in Bin County. (a) December 10th. (b) January 10th. (c) March 10th.
(a)(b)(c)


## 3.2. Moisture Field Variation of Shallow Soil

As Figure5 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.Figure 5
Moisture content variation of shallow loess at different sites with depth. (a) Bin County. (b) Luochuan. (c) Tongchuan.
(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]:(1)K+NΔtTt=Pt+NΔtTt−Δt,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:(2)kl⋅n=∫−11∫−11λ16|JD2∂Hn∂ξ+D1∂Hn∂ηdξdη+Aαsl,n=i,j,k,m,(3)nl⋅n=∫−11∫−11ρCp|JHlHndξdη,l,n=i,j,k,m,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:(4)λ=4.17w2+1504100.25ρd−3.9,(5)C=ρd1.27+0.021w103.With reference to the relevant literature [30], the array {P} is composed of three terms and calculated according to the following formula:(6)P=P1+P2+P3+P4,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:(7)P1l=∑s=12∑t=12ωsωtρLΔtΔfl|J|Hl|ξs⋅ηt,(8)P2=P2S+P2E+P2A,(9)P2Sl=Qh1−λ1⋅s2,(10)P2El=−ℓ1×4.88×10−8×273+T4⋅s2,(11)P2Al=ℓ2×4.88×10−8×273+Ta4⋅β′⋅s2,(12)P3l=α⋅s⋅Ta2,(13)P4l=U⋅G⋅s2,where l is the node of the boundary of the third kind, Δfl is the variation of solid phase ratio at the node l within the time interval Δt, L is the latent heat of phase change, ωs and ωt are the weighted coefficients, ρ is the density, Qh is the solar radiation volume per square meter of ground level, λ1 is the reflectivity of ground to solar radiation, ℓ1,ℓ2 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. Table2 presents the average solar radiation per ten days in winter in Bin County based on the meteorological data in the past 30 years.Table 2
Solar radiation in Bin County (MJ/m2).
TimeSolar radiationTimeSolar radiationMid-December76.2Late December73.1Early January69.7Mid-January47.9Late January80.8Early February95.3Mid-February84.4Late February99.2Early March128.2The emissivity of the earthℓ1 is taken as 0.68, but it is complex to determine the values of the atmospheric emissivity ℓ2 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 ℓ2and β. Based on the fitting of meteorologically measured natural surface temperature and the consideration of ℓ2 and β in formula (11) as a whole, the expression of ℓ2β′ can be obtained as follows:(14)ℓ2β′=f+0.006Ta+0.004Sd,where Ta is the air temperature (°C), Sd is the relative humidity (%), and f 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:(15)α=5.7+4.1V,where the unit of wind velocity V is m/s.The evaporation of the soil surface can be calculated as follows:(16)U=−0.1+0.118w−0.0043w2uw,where w is the moisture content of surface soil and uw is the evaporation of the water surface.Climatic parameters of Bin County in winter are presented in Table3.Table 3
Climatic parameters of Bin County.
MonthTa/(°C)Sd/(°C)uw/mmV/(m/s)12−3.26541.82.51−5.05945.02.62−2.16257.42.734.563110.02.9


### 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]:(17)D+EΔtht+Δt=EΔtht+F,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:(18)h=hg+hu+hc,where hg 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]:(19)θh=θr+θs−θr1+αhnm,where θr=−0.38+0.36ρd, θs=1−0.38ρd, α=e8.89−0.08Tt+−9.36+0.07Ttρd, n=10.51−0.78Tt+0.015Tt2+−14+1.14Tt−0.022Tt2ρd+5.56−0.41Tt+0.008Tt2ρd2, 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:(20)hc=0.00025e6.4ρd1−0.5Sr,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:(21)T¯i=ωTi+1−ωT¯i−1,(22)h¯i=ωhi+1−ωh¯i−1,where T¯i−1 and h¯i−1 are the ith iteration values, T¯i and h¯i 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.Figure 6
Numerical calculation results of temperature field.
(a)(b)Figure 7
Numerical calculation results of moisture field.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 Figure8.Figure 8
The relationship between the temperature and the maximum frozen depth.As Figure8 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:(23)h=18.3×Ta+3.7,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.


## 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]:(1)K+NΔtTt=Pt+NΔtTt−Δt,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:(2)kl⋅n=∫−11∫−11λ16|JD2∂Hn∂ξ+D1∂Hn∂ηdξdη+Aαsl,n=i,j,k,m,(3)nl⋅n=∫−11∫−11ρCp|JHlHndξdη,l,n=i,j,k,m,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:(4)λ=4.17w2+1504100.25ρd−3.9,(5)C=ρd1.27+0.021w103.With reference to the relevant literature [30], the array {P} is composed of three terms and calculated according to the following formula:(6)P=P1+P2+P3+P4,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:(7)P1l=∑s=12∑t=12ωsωtρLΔtΔfl|J|Hl|ξs⋅ηt,(8)P2=P2S+P2E+P2A,(9)P2Sl=Qh1−λ1⋅s2,(10)P2El=−ℓ1×4.88×10−8×273+T4⋅s2,(11)P2Al=ℓ2×4.88×10−8×273+Ta4⋅β′⋅s2,(12)P3l=α⋅s⋅Ta2,(13)P4l=U⋅G⋅s2,where l is the node of the boundary of the third kind, Δfl is the variation of solid phase ratio at the node l within the time interval Δt, L is the latent heat of phase change, ωs and ωt are the weighted coefficients, ρ is the density, Qh is the solar radiation volume per square meter of ground level, λ1 is the reflectivity of ground to solar radiation, ℓ1,ℓ2 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. Table2 presents the average solar radiation per ten days in winter in Bin County based on the meteorological data in the past 30 years.Table 2
Solar radiation in Bin County (MJ/m2).
TimeSolar radiationTimeSolar radiationMid-December76.2Late December73.1Early January69.7Mid-January47.9Late January80.8Early February95.3Mid-February84.4Late February99.2Early March128.2The emissivity of the earthℓ1 is taken as 0.68, but it is complex to determine the values of the atmospheric emissivity ℓ2 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 ℓ2and β. Based on the fitting of meteorologically measured natural surface temperature and the consideration of ℓ2 and β in formula (11) as a whole, the expression of ℓ2β′ can be obtained as follows:(14)ℓ2β′=f+0.006Ta+0.004Sd,where Ta is the air temperature (°C), Sd is the relative humidity (%), and f 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:(15)α=5.7+4.1V,where the unit of wind velocity V is m/s.The evaporation of the soil surface can be calculated as follows:(16)U=−0.1+0.118w−0.0043w2uw,where w is the moisture content of surface soil and uw is the evaporation of the water surface.Climatic parameters of Bin County in winter are presented in Table3.Table 3
Climatic parameters of Bin County.
MonthTa/(°C)Sd/(°C)uw/mmV/(m/s)12−3.26541.82.51−5.05945.02.62−2.16257.42.734.563110.02.9


## 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]:(17)D+EΔtht+Δt=EΔtht+F,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:(18)h=hg+hu+hc,where hg 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]:(19)θh=θr+θs−θr1+αhnm,where θr=−0.38+0.36ρd, θs=1−0.38ρd, α=e8.89−0.08Tt+−9.36+0.07Ttρd, n=10.51−0.78Tt+0.015Tt2+−14+1.14Tt−0.022Tt2ρd+5.56−0.41Tt+0.008Tt2ρd2, 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:(20)hc=0.00025e6.4ρd1−0.5Sr,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:(21)T¯i=ωTi+1−ωT¯i−1,(22)h¯i=ωhi+1−ωh¯i−1,where T¯i−1 and h¯i−1 are the ith iteration values, T¯i and h¯i 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.Figure 6
Numerical calculation results of temperature field.
(a)(b)Figure 7
Numerical calculation results of moisture field.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 Figure8.Figure 8
The relationship between the temperature and the maximum frozen depth.As Figure8 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:(23)h=18.3×Ta+3.7,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.



---

*Source: 1690406-2020-05-19.xml*