water
ArticleIn-Situ Monitoring and Characteristic Analysis ofFreezing-Thawing Cycles in a Deep Vadose Zone
Ce Zheng 1,2 , Yudong Lu 1 , Xiuhua Liu 1,3, * , Jiří Šimůnek 2 , Yijian Zeng 4 ,Changchun Shi 5 and Huanhuan Li 1
1

2
3
4
5

*

Key Laboratory of Subsurface Hydrology and Ecological Effects in Arid Region of Ministry of Education,
School of Water and Environment, Chang’an University, Xi’an 710054, China; ruudfanschd@139.com (C.Z.);
luyudong18@163.com (Y.L.); 17742499497@163.com (H.L.)
Department of Environmental Sciences, University of California, Riverside, CA 92521, USA;jsimunek@ucr.edu
Shaanxi Key Laboratory of Land Consolidation, Chang’an University, Xi’an 710064, China
Faculty of Geo-Information and Earth Observation (ITC), University of Twente, Hengelosestraat 99,
7514 AE Enschede, The Netherlands; y.zeng@utwente.nl
China State Long-Term Observation and Research Station for Mu Us Desert Ecosystem in Yulin of Shaanxi,
Yulin 719000, China; snsun0509@163.com
Correspondence: liuxh68@chd.edu.cn

Received: 20 February 2020; Accepted: 27 April 2020; Published: 29 April 2020


Abstract: Freeze-thaw cycles play a critical role in affecting ecosystem services in arid regions.
Monitoring studies of soil temperature and moisture during a freeze-thaw process can generate
data for research on the coupled movement of water, vapor, and heat during the freezing-thawing
period which can, in turn, provide theoretical guidance for rational irrigation practices and ecologicalprotection. In this study, the soil temperature and moisture changes in the deep vadose zone wereobserved by in-situ monitoring from November 2017 to March 2018 in the Mu Us Desert. The results
showed that changes in soil temperatures and temperature gradients were largest in soil layers
above the 100-cm depth, and variations decreased with soil depth. The relationship between soil
temperature and unfrozen water content can be depicted well by both theoretical and empiricalmodels. Due to gradients of the matric potential and temperature, soil water flowed from deeper soil
layers towards the frozen soil, increasing the total water content at the freezing front. The vapor flux,
which was affected mainly by temperature, showed diurnal variations in the shallow 20-cm soil layer,
and its rate and variations decreased gradually with increasing soil depths. The freeze-thaw process
can be divided into three stages: the initial freezing stage, the downward freezing stage, and thethawing stage. The upward vapor flux contributed to the formation of the frozen layer during thefreezing process.
Keywords: freeze-thaw process; soil temperature and moisture; unfrozen water content; water vapor;deep vadose zone

1. Introduction
Regions with frozen soil are widely distributed in middle and high latitudes, affecting
approximately 50% of the land around the world [1,2]. Soils in these regions will freeze or thaw inresponse to variations in soil temperatures, resulting in the phase change of soil water among ice,liquid water, and water vapor [3–5]. Abundant forest and mineral resources exist in these areas, and thefrozen soil, representing a vital factor of the biological environment, affects productive activities and thesustainable development of these resources [6]. As the research of the critical zone increases in scope,

Water 2020, 12, 1261; doi:10.3390/w12051261

www.mdpi.com/journal/water

Water 2020, 12, 1261

2 of 18

studies on the dynamics in soil temperature and moisture in the vadose zone during freezing-thawing
periods have attracted a growing interest [7,8].
The freeze-thaw process has a substantial impact on the surface energy balance and soil moisture
distribution, significantly affecting soil properties, such as soil structure, permeability, conductivity,and bulk density, making it difficult and complicated to study water flow, heat transport, and relatedparameters in seasonally frozen soils [9–11]. For instance, the freezing process reduces soil and air
permeabilities, influences the heat exchange, and changes the distribution of soil water. Furthermore,
phase transitions of soil water occur frequently due to freeze-thaw cycles, resulting in variationsin unfrozen water contents at subzero temperatures [12,13]. Physical properties of frozen soil incold regions are strongly dependent on temperature when unfrozen water is present, and both
unfrozen water and soil temperature control the freeze-thaw process and soil water migration [2,14,15].
Also, the presence of salts can alter the freeze-thaw process in soils and their transport during thefreeze-thaw period could lead to increased soil salinization [16,17]. The coexistence of unfrozen waterand ice in frozen soil affects its hydrological and thermal properties, becoming a significant factor
for many engineering and environmental applications [18]. During the freezing process, unfrozen
water flows in the direction of soil matric potential and temperature gradients from unfrozen soil in
deeper soil layers, increasing the water content of the frozen layer [19–21]. The total soil water during
the freezing period can be divided into “freezable” unfrozen water, absorbed water, and unfreezable
water [22,23].
In arid areas with a deep groundwater level, moisture migration in the vadose zone often takes
place in the form of water vapor [24,25], and the transformation between liquid water and water vapor
occurs continuously, which results in water vapor playing a key role in maintaining surface vegetationand ecosystems [26,27]. Water vapor, driven mainly by temperature gradients, flows from warmer
to cooler soil layers and performs a critical function in affecting variations in soil water contents inthe vadose zone [28,29]. The existing literature indicates that vapor flow has a significant influenceon variations in soil water contents during the freeze-thaw process under such conditions [24,30].
Both liquid water and water vapor flow upward towards the freezing front, and it is water vapor thatconnects water transfer below the freezing front and above the evaporation front [30]. Furthermore,since water vapor fluxes can become larger than liquid water fluxes in deep soil layers, only a model
considering the coupled movement of water, vapor, and heat can fully describe critical physicalmechanisms of the hydrological cycle in the vadose zone [31].
Many studies have evaluated coupled interactions between water and heat, and the influenceof freeze-thaw cycles on soil properties. However, most of these studies were carried out either inlaboratories or in areas with shallow groundwater, which have obvious spatial limitations [14,32].
Therefore, in-situ monitoring of the deep vadose zone is highly needed to improve our understandingof dynamic changes in the spatial-temporal distributions of moisture and temperature in frozen soils.
Moreover, monitoring studies of soil water and heat movement under the freeze-thaw process in a
deep vadose zone can provide data required for the development and validation of models simulating
the coupled movement of water, vapor, and heat under such conditions, as well as for establishingpractical guidance for rational irrigation management and ecological protection.
Therefore, the objectives of this study were (i) to monitor the soil temperature and moisturechanges of the 8-m deep vadose zone in the Mu Us Desert during freezing-thawing periods, (ii) toanalyze spatial and temporal distributions of soil temperatures and water contents, and (iii) to evaluate
water vapor fluxes in the soil profile using Fick’s law.
2. Materials and Methods
2.1. Study Site
The Mu Us Desert is located between the Ordos Plateau and the Loess Plateau in northwestern
China (Figure 1), in the north part of the Shaanxi Province, northeast of the Ningxia Hui Autonomous

Water 2020, 12, 1261

3 of 18

Region,of REVIEW
the Inner Mongolia Autonomous Region.
Water
2020,and12, x south
FOR PEER

The climate in the desert can
3 of be18
characterized as dry and cold, with large temperature differences and frequent strong winds. Based on◦ C,36.4
onmeteorologicalannualtemperaturesfrom to−31.4
anthethemeteorological
datadata
sincesince1961,1961,
annualtemperatures
range range
from −31.4
36.4 towith°C,an withaverage◦averageannual temperatureof 6.4°C. annual
Mean annual
precipitationand potential
evaporation
are and360
annual temperatureof 6.4 C.
Mean
precipitationand potential
evaporation
are 360
andrespectively,
over
60%
the
annualrainfall
rainfalloccurs
occursfrom
from July
July to September
23432343
mm,mm,
respectively,
andand
over
60%
of of
the
annualSeptember [33,34].
[33,34].
2
Total annual sunshine time exceeds 3000 h, with total radiation of 608.37 kJ/cm . The
The MuMu Us
Us Desert
Desertisfreezing,
withthethefreeze-thawprocess
fromfrom
Novemberto March.
Theis typically
typicallysubject
subjecttotoseasonalseasonal
freezing,withfreeze-thawprocess
Novemberto March.
frozenlayer
usuallyextends
to
100
cm
below
the
soil
surface.
The frozen layer usually extends to 100 cm below the soil surface.

Figure 1.
1. Location
Location ofof thethe study
study site
site andand thethe field
field instrumentation
instrumentation layout.
layout.
FigureThe MuMu Us
Us Desert
Desert isis oneone ofof thethe most
most vulnerable
vulnerable ecologicalecological areas
areas inin China.
China. With
With rapid
rapid economic
economicThedevelopment
andlarge-scaleconstructionofenergybases,
groundwater
resources
have
been
extensively
development and large-scale construction of energy bases, groundwater resources have beendeveloped, developed,
resulting inresulting
a significant
in thegroundwaterlevel and anlevel
increasein increase
the vadoseextensively
in a decline
significant
declinein the groundwaterand aninzone
thickness[35–37].
Duringthefreezing
period,diurnaltemperaturechangeshaveaconsiderable
the vadose zone thickness [35–37]. During the freezing period, diurnal temperature changes have aeffect on soil effectwateronmovement,
dynamics
soil waterincontents
andtemperatures
are closely
considerable
soil water andmovement,
andindynamics
soil watercontents
and temperatures
relatedto vegetation
growth, landdesertification,and geologic
hazards.
areclosely
related to vegetation
growth,land desertification,andgeologic hazards.
2.2. In-Situ Observations
2.2. In-Situ Observations
The field work was conducted at the Yu Lin Desert Ecosystem Research Station (longitude
The
field work was conducted at the Yu Lin Desert Ecosystem Research Station (longitude◦
109 420 2900 E, latitude 38◦ 230 1900 N, at an altitude of 1157 m above the sea level). At the study
109°42′29′′ E, latitude 38°23′19′′ N, at an altitude of 1157 m above the sea level). At the study site, sandsite, sand is the dominant soil fraction in the active soil layer and the land cover is dominated byis the dominant soil fraction in the active soil layer and the land cover is dominated by Salix
Salix Pasmmophila and Artemisia. Vegetation is sparse during the freezing period, covering less than 5%
Pasmmophila and Artemisia. Vegetation is sparse during the freezing period, covering less than 5% ofof the land surface. As shown in Figure 1, a monitoring well with a diameter of 150 cm was dug to a
the land surface. As shown in Figure 1, a monitoring well with a diameter of 150 cm was dug to adepth of 900 cm, which is below the groundwater table (the groundwater level was located at a soil
depth of 900 cm, which is below the groundwater table (the groundwater level was located at a soil
depth of 852 cm on 1 November 2017). Soil samples were collected to conduct the particle size analysis
depth of 852 cm on 1 November 2017). Soil samples were collected to conduct the particle size analysis(Table 1). Sand and clay accounted for approximately 95% and 1% of the soil mass in the vadose zone,(Table 1). Sand and clay accounted for approximately 95% and 1% of the soil mass in the vadose zone,
respectively, except for the soil horizon between 80 and 230 cm, where the clay content was slightly
respectively, except for the soil horizon between 80 and 230 cm, where the clay content was slightly
higher. For the measurements not to be influenced by external meteorological conditions, the well was
higher. For the measurements not to be influenced by external meteorological conditions, the wellcovered by a 15 cm thick piece of polyethylene foam with the thermal conductivity of 0.03 W/m/K and
was covered by a 15 cm thick piece of polyethylene foam with the thermal conductivity of 0.03 W/m/K
a 10 cm thick well cover.
and a 10 cm thick well cover.
Monitoring sensors were installed at different depths of the soil profile through the wall of
the well (Table 2). To minimize the influence of the well itself on measured values, the average
distance of sensors from the wellbore was 50 cm. Monitoring started 6 months after the sensors
were installed to allow soil to settle and establish better contact with the sensors. The Hydra Probe
II sensors, which are based on the frequency domain (FD) measurement principles, were calibratedusing gravimetric measurements (samples were oven-dried at 105 ◦ C for 24 h; gravimetric water

Water 2020, 12, 1261

4 of 18

contents were converted to volumetric water contents using the soil bulk density) on soil samples
collected from the study site. The concept of unfrozen water content measurements is based on theconfirmed similarity between soil drying-wetting and freezing-thawing processes [3,38]. When testingat 100 MHz, the permittivity of liquid water (≈78) is significantly higher than those of the soil matrix(≈3–5), ice (≈3.4) and air (≈1) [15,38,39]. Moreover, the soil permittivity of the liquid phase is much
more sensitive to soil temperature changes than those of solids and air [40]. Therefore, observed
soil water contents were assumed to represent unfrozen water contents when the soil was frozen.
Soil samples collected from different depths during the freezing period were taken to the laboratory
to determine the total water content by an oven drying method. All sensors in the monitoring well
were linked to a solar-powered automatic data-logger (CR1000), which recorded data at a 10 min
interval. An automatic micrometeorological station was installed to record related meteorologicalvariables at a 30 min interval, such as wind speed, air temperature, relative humidity, and precipitation.
The observed data at the study site, from 1 November 2017 to 31 March 2018, was collected to analyze
the variations in soil temperature and moisture during the freeze-thaw process.
Table 1. Measured soil particle composition and bulk density at different depths.
Soil Layer (cm)
0–10
10–50
50–80
80–160
160–230
230–560
560–800

Soil Particle Composition (%)
Sand
Silt
Clay
96.2
95.7
95.5
90.1
91.7
94.0
95.4

3.1
3.2
3.3
5.4
5.2
4.7
3.5

0.7
1.1
1.2
4.5
3.1
1.3
1.1

Bulk Density(g/cm3 )
1.57
1.59
1.57
1.51
1.52
1.57
1.57

Table 2. Monitoring instruments at the study site. Positive heights mean that sensors were installed
aboveground, while negative heights indicate that sensors were installed in the monitoring well below
the surface. Soil temperatures and water contents measured at a depth of 2 cm were considered torepresent the changes at the soil surface.
Variables

Sensors

Manufacturers

Height/Depth (cm)Wind speed (m/s) andwind direction

Davis Cup
Anemometer

Decagon240

Precipitation (mm)ECRN-100

Decagon200

VP-3

Decagon20

Soil temperature (◦ C) andwater content (cm3 /cm3 )Hydra Probe II

Stevens−2, −10, −20, −50, −100,−130, −200, −400, −630,−660 −800

Soil matric potential (hPa)TensioMark

Stevens−2, −10, −20, −50, −100,−130, −200, −400, −630,−660 −800

Soil heat flux (W/m2 )HFP01

Hukseflux−20, −200, −800

Groundwater level (cm)

Submersible Depth
Transmitter

Stevens−880

Air temperature (◦ C) andrelative humidity

Figure 2 shows the mean, maximum, and minimum air temperatures during the observation
period. Air temperatures varied greatly during each day, with a mean diurnal temperature range of
22.8 ◦ C, which dropped sharply in November and then increased in March. The mean temperature of
the coldest month (January) and warmest month (March) during the study period reached −11.5 and6.1 ◦ C, respectively. Rainfall and snow during the study period were only 20.8 mm. It is clear that

Water 2020, 12, x FOR PEER REVIEW

5 of 18

22.8 °C, which dropped sharply in November and then increased in March. The mean temperature
5 of 18of the coldest month (January) and warmest month (March) during the study period reached −11.5
and 6.1 °C, respectively. Rainfall and snow during the study period were only 20.8 mm. It is clearthat snow
had effectlittle effecton variationsin the unfrozen
water content
at shallow
depthstheduringthesnow
had little
on variationsin the unfrozen
water content
at shallow
depths duringfreezing
freezing
period2b,c).
(Figure 2b,c).
period(FigureWater 2020, 12, 1261

Figure(a) Measured
Measuredair
airtemperature,temperature,rainfall,
rainfall,andandsnow
snowduring
duringtheperiodfromNovember
2017
Figure 2.
2. (a)
theperiodfromNovember
2017
to
to
March
2018
at thestudy
unfrozenwater
content
cm depth
during
two
March
2018
at thestudy
site;site;
andand(b,c)(b,c)
unfrozenwater
content
at theat10thecm10depth
duringtwo snowsnow
periods.
periods.

2.3. Mathematical Calculations
2.3. Mathematical Calculations
2.3.1. Soil Heat Transport
2.3.1. Soil Heat Transport
Many analytical solutions of the soil energy balance equation coupled with Fourier’s law exist.
analyticalsolutionsof the soil
energybalanceequation
coupledwithboundary
Fourier’s condition
law exist.
One Manyof themholds forhomogeneous
soils,
no flow
conditions,andan upperOne ofbythem
holds for homogeneous
soils,
no flow
conditions,
upper boundary
condition
given
a temperaturesine wave withaveragetemperature
T0 and(◦ C),ana temperatureamplitude
A0
given
by
a
temperature
sine
wave
with
averagetemperatureT
0 (°C), a temperature amplitude A0 (°C),◦( C), the time period τ (e.g., 24 h or 365 day), and a constant (both temporally and spatially) thermal
the time period
τ (e.g., 24 h or 365 day), and a constant (both temporally and spatially) thermal
diffusivity
KT [41]:diffusivity KT [41]:zz
T (z, t) = T0 + A0 exp(− ) sin ωt + φ −
(1)
d
d
z
z


z , t )t =is Ttime) sin (cm),where T is the temperatureT(◦(C),z is− depth
0 + A(h),0 exp(
 ω t +dφis− the damping depth (cm), ω(1)is
d
d

−1
the angular frequency (h ), and ϕ is the phase constant (-). The damping depth d and the diurnal◦ C) at depth z can be calculated as follows:
amplitude
where T is A
the(°C), t is time (h), z is depth (cm), d is the damping depth (cm), ω is thez ( temperature
angular frequency (h−1), and φ is the phase constant(-). The damping depth d and the diurnalrKτ
amplitude Az (°C) at depth z can be calculated as follows:
T
d=
(2)
π
z
Az = A0 exp(− )d(3)Water 2020, 12, 1261

6 of 18

2.3.2. Unfrozen Water Content
Based on the analogy between freezing and drying processes and the assumption that the ice
pressure is equal to zero at freezing conditions, the relationship between soil subfreezing temperatureand the liquid water matric potential can be expressed using the generalized Clapeyron equation [13]:
h=

Lf
g

ln

T
T0

(4)

where h is the liquid water matric potential (m), Lf is the latent heat of vaporization (approximately
3.34 × 105 J/kg), g is the gravitational acceleration (=9.81 m/s2 ), T is the temperature (K), and T0 is thefreezing point of liquid water (=273.15 K).
According to the van Genuchten [42] model, the soil–water characteristic curve, which relates thewater content with the pressure head, is as follows:θl (h, T ) = θr + h

θs − θr
1 + −αh(T )

(5)

n im

where θl , θr , and θs are the liquid, residual, and saturated water contents (m3 /m3 ), respectively, and α(m−1 ), n (-), and m (-) are empirical parameters. The soil hydraulic parameters can be estimatedfrom soil physical properties using neural network-based pedotransfer functions implemented in
the numerical model HYDRUS-1D [43]. Equations (4) and (5) can be used to characterize the soil
freezing characteristics curve to predict the unfrozen soil water content with measured subfreezing
temperatures instead of measured pressure heads.
Additionally, an empirical method proposed by Xu et al. [44] can estimate the variations of
unfrozen water contents with subfreezing temperatures:θl = a(−T )b(6)

where coefficients a and b are empirical constants that depend on soil texture and the initial watercontent. These two coefficients can be estimated using two (subscripts 1 and 2) freezing temperatures
Tf (◦ C) and corresponding unfrozen water content θ (m3 /m3 ) data pairs as follows:ln θ1 − ln θ2
ln T f 1 − ln T f 2(7)

a = θ1 T f 1 b

(8)

b=−

2.3.3. Thermal Vapor Flux
In general, the vapor density ρv (kg/m3 ), which is used to calculate water vapor fluxes, can becalculated as a function of the saturated vapor density ρsv (kg/m3 ) and relative humidity Hr as
follows [45]:ρv = ρsv × Hr
(9)The saturated vapor density can be expressed as a function of temperature:
ρsv =



− 7.92495 × 10−3 T
exp 31.3716 − 6014.79
TT× 10−3(10)

Water 2020, 12, 1261

7 of 18

and the relative humidity Hr can be calculated using a thermodynamic relationship between liquidwater and vapor in the soil as follows:Mgh
Hr = exp
RT

!
(11)

where M is the molecular weight of water (=0.018015 kg/mol) and R is the universal gas constant ofwater vapor (=8.315 J/mol/K).
The thermal vapor flux qv (m/s) can be calculated using Fick’s Law:
qv = −KvT

∂T
∂z

(12)

where KvT is the thermal vapor hydraulic conductivity (m2 /K/s), which can be expressed as follows:
KvT =

dρsv
D
ηHr
ρw
dT

(13)

where D is the vapor diffusivity in soil (m2 /s), which is defined as D = Da θa τ, where Da is the vapor
diffusivity in air (=2.12 × 10−5 (T/273.15)2 m2 /s), θa is the air-filled porosity (m3 /m3 ), τ is the tortuosity
factor (=θa 7/3 /θs 2 ) [46], ρw is the density of liquid water (kg/m3 ), and η is the enhancement factor (-)
that can be calculated as follows [47]:  4 
θ
2.6  θ   

η = 9.5 + 3 − 8.5 exp−1 + p   θs
fc θs (14)

where fc is the mass fraction of clay in the soil (-).
3. Results
3.1. Measured Spatial-Temporal Distribution of Soil Temperature
The spatial-temporal distribution of soil temperature measured during the observation period
is shown in Figure 3. At the beginning of monitoring, soil temperatures in the upper part of thevadose zone were relatively low, while the maximum temperature (15.4 ◦ C) was observed at a depth of
400 cm, resulting in an upward temperature gradient. As air temperatures kept falling after November,soil temperatures decreased accordingly and were invariably higher than air temperatures, resulting
in the soil being exothermic. Due to a persistent upward heat flux, the depth of the maximum soil
temperature in the vadose zone continued to increase, eventually reaching the 800 cm depth on
1 January. During this time period, the heat flux recorded by HFP01 at a depth of 800 cm decreased froma positive (downward) value of 0.76 W/m2 to a negative (upward) value of −0.02 W/m2 , which wasconsistent with observed soil temperatures. Measured data indicated that soil temperatures above
a depth of 200 cm reached their lowest values gradually in February and subsequently rose withan increase in air temperatures, especially in the top 50 cm layer, where they exceeded 10 ◦ C. As aresult, soil temperatures between depths of 100 and 200 cm became the lowest throughout the vadosezone. It should be noted that soil temperatures in the lower part of the vadose zone retained a slowly
decreasing trend, even at the end of the observation period.

Water 2020, 12, 1261
Water 2020, 12, x FOR PEER REVIEW

8 of 18
8 of 18

Figureof ofsoilsoiltemperaturesduringthe monthsof (a)ofNovember
2017,
Figure 3.
3. Spatial
Spatialandandtemporal
temporalvariations
variations
temperaturesduringthe months(a) November(b)December
2017,(c)January2018,(d)February
2018,
and
(e)
March
2018.
2017, (b) December 2017, (c) January 2018, (d) February 2018, and (e) March 2018.

It is evident that the spatial distribution of temperature gradients (Figure 4) is consistent with soil
It is evident that the spatial distribution of temperature gradients (Figure 4) is consistent withtemperatures, with temperature gradients decreasing with increasing soil depths. Note that variations
soil temperatures, with temperature gradients decreasing with increasing soil depths. Note thatin temperature gradient at soil depths above the 20 cm depth differed greatly (about 1–2 orders ofvariations in temperature gradient at soil depths above the 20 cm depth differed greatly (about 1–2
magnitude) from those at deeper soil layers. During the freezing period, the temperature gradientorders of magnitude) from those at deeper soil layers. During the freezing period, the temperature
varied from 0.248 to 0.055 ◦ C/cm between the soil surface and a depth of 100 cm, reflecting obvious
gradient varied from 0.248 to 0.055 °C/cm between the soil surface and a depth of 100 cm, reflecting
variations in soil temperature in this layer. On the contrary, temperature gradients were relatively
obvious variations in soil temperature in this layer. On the contrary, temperature gradients were
low in the lower parts of the vadose zone. During the thawing period, the soil layer above a depth ofrelatively low in the lower parts of the vadose zone. During the thawing period, the soil layer above
100 cm displayed negative (downward) temperature gradients, especially in the top 20 cm depths,
a depth of 100 cm displayed negative (downward)◦ temperature gradients, especially in the top 20 cm
with the gradients ranging from −0.102 to −0.051 C/cm.
depths, with the gradients ranging from −0.102 to −0.051 °C/cm.

Water 2020, 12, 1261
Water 2020, 12, x FOR PEER REVIEW

9 of 18
9 of 18

FigureSpatial andandtemporal
temporalvariations
variationsofofsoil
soiltemperaturetemperaturegradients
gradients
duringmonthsFigure 4.
4. Spatial
duringthethemonths
of of
(a)
(a)
November
2017,
December
2017,(c) January
2018,(d) February(e) March
November2017,(b) (b)December
2017,(c) January
2018,(d) February
2018,2018,
and and(e) March2018.2018.
The
The
temperature
gradient
was calculated
as ∆T/∆z= (Ti+1where
where
Ti are soil
− Ti )/∆z,temperature
gradient
was calculated
as ΔT/Δz=(Ti+1 − Ti)/Δz,
Ti+1 andTi areTi+1
soil andtemperaturesattemperaturesatdepthsiandi+1,respectively,
and
∆z
is
the
verticaldistance
betweendepths i and i + 1, respectively, and Δz is the vertical distance between depths i and i + 1. depths i
and i + 1.

Table 3 provides a list of quantitative variables characterizing variations of soil temperatures atTable 3 provides a list of quantitative variables characterizing variations of soil temperatures atdifferent depths. During the observation period, soil temperatures at depths of 2, 10, and 20 cm
different depths. During the observation period, soil temperatures at depths of 2, 10, and 20 cm ranged
ranged from −12.4 to 15.9 °C, −10.1 to 14.0 °C, and −8.6 to 13.1 °C, respectively, and were significantly
from −12.4 to 15.9 ◦ C, −10.1 to 14.0 ◦ C, and −8.6 to 13.1 ◦ C, respectively, and were significantly affectedaffected by changes in air temperatures (with correlation coefficients of 0.952, 0.938, and 0.899,
by changes in air temperatures (with correlation coefficients of 0.952, 0.938, and 0.899, respectively).
respectively). The range of temperature variations decreased with depth, and the temperature
The range of temperature variations decreased with depth, and the temperature changes becamechanges became relatively small below the 200 cm depth (less than 10 °C). Furthermore, therelatively small below the 200 cm depth (less than 10 ◦ C). Furthermore, the coefficient of variation (Cv ),coefficient of variation (Cv), as expected, decreased as the soil depth increased, indicating a weakening
as expected, decreased as the soil depth increased, indicating a weakening trend in soil temperature
trend in soil temperature variations with depth.
variations with depth.
As the soil in this study is relatively dry (as discussed below in detail), the effect of the specific
As the soil in this study is relatively dry (as discussed below in detail), the effect of the specificheat of freezing/thawing can be neglected, and the Equation (3) can be used to approximate soil
heat of freezing/thawing can be neglected, and the Equation (3) can be used to approximate soil
temperatures with depth at subzero temperatures as well. Figure 5 illustrates measured andtemperatures with depth at subzero temperatures as well. Figure
5 illustrates measured and calculated
calculated (with T0 = −2.1 °C, A0 =5.4 °C, and KT =20.7 cm2/h) diurnal soil temperature amplitudes as(with T0 = −2.1 ◦ C, A0 = 5.4 ◦ C, and KT = 20.7 cm2 /h) diurnal soil temperature amplitudes as a function
a function of depth. The results show that calculated and measured temperature amplitudes agreedof depth. The results show that calculated and measured temperature amplitudes agreed well, with awell, with a correlation coefficient of 0.98. The observed diurnal temperature amplitude at the soil
correlation coefficient of 0.98. The observed diurnal temperature amplitude at the soil surface, and at
surface, and at depths of 10 and 20 cm, were 5.4, 2.7, and 1.4 °C, respectively, while the amplitudedepths of 10 and 20 cm, were 5.4, 2.7, and 1.4 ◦ C, respectively, while the amplitude became much
became much smaller (less than 0.1 °C) below the 100 cm depth. Calculated values were always
slightly lower than measured values, which was mainly because the measured amplitude at a depthWater 2020, 12, x FOR PEER REVIEW

10 of 18

Water
2020,
12, 1261
of 18
of 2 cmwasadopted to represent temperature variations at the soil surface, resulting in a lower10
valueof the diurnal temperature amplitude at the surface (A0).
smaller (less than 0.1 ◦ C) below the 100 cm depth. Calculated values were always slightly lower
Table 3. Daily temperature characteristics calculated for different depths. The coefficient of variation
than measured values, which was mainly because the measured amplitude at a depth of 2 cm wasis calculated based on the Kelvin temperature scale.
adopted to represent temperature variations at the soil surface, resulting in a lower value of the diurnal
Average
Minimum
MaximumStandard
Coefficient of
temperature amplitude
at the surface(A0 ).
Depth
Temperature
Temperature
Temperature
Deviation
Variation(cm)
Table 3. Daily temperature
characteristics
calculatedfordifferent
depths.
The
coefficientofvariation
is
(°C)
(°C)
(°C)
(°C)
calculated based on the Kelvin temperature scale.

2
10Depth (cm)
20
2
50 10
100 2050130 100
200 130400 200400
630 630
660 660
800
800

−2.1
Average−1.5
Temperature (◦ C)
−0.8
−2.1
0.7−1.5
2.7−0.8
0.7
4.32.7
6.94.3
6.9
10.710.7
12 12
12.112.1
12.3
12.3

−12.4
Minimum−10.1
Temperature (◦ C)
−8.6
−12.4
−5.3
−10.1
−1.7−8.6−5.3
0.6−1.7
3.5 0.6
7.9 3.5
7.9
10.410.4
10.710.7
11.6
11.6

15.9
Maximum
14.0
Temperature (◦ C)
13.1
15.9
10.614.0
11.913.1
10.6
13.011.9
14.513.0
15.414.5
15.4
14.114.1
13.813.8
13.3
13.3

6.67
Standard
5.90
Deviation (◦ C)
5.19
6.67
4.03
5.90
5.19
3.66
4.03
3.48
3.66
3.48
3.35
3.35
2.43
2.43
1.26
1.26
1.06
1.06
0.61
0.61

0.0246
Coefficient of
0.0217
Variation

0.01910.02460.01470.02170.01910.01330.01470.01250.01330.01250.01200.01200.00860.00860.00440.00440.00370.00370.00210.0021

Changes in the diurnal soil temperature amplitude as a function of depth.
Figure 5. Changes

3.2.
Measured Spatial-Temporal
Spatial-Temporal Distribution
3.2. Measured
Distribution ofof Soil
Soil Water
Water Content
Content
Figurethethespatial-temporalvariationsof the unfrozen
water content
the observation
Figure 66showsshows
spatial-temporalvariationsof the unfrozen
waterduringcontent
during theperiod.
Thechangesintheunfrozenwater
content
weredistinctive
betweenthesoil
surfaceandaobservation period. The changes in the unfrozen water content were distinctive between
the
soil
depth
100acm.
Dueofto100
intensive
radiationsolar
and radiation
sparse rainfall,
the average
watercontent waterin thesurfaceofanddepth
cm. Duesolar
to intensive
and sparse
rainfall,the average
3 /cm3 . During the freezing period, the unfrozen watersoil
profile
was
very
low:
approximately
0.06
cm
content in the soil profile was very low: approximately 0.06 cm3/cm3. During the freezing period, thecontent
decreasedwith time
elapsed.with
The time
measureddataTheshoweda significant
decreasea
unfrozengraduallywater content
gradually
decreased
elapsed.
measureddata showedinthe unfrozen
water
content
at depthsof 10, 20,
50, andof100
cm50,
onand
24 November,
December,significant
decreasein theunfrozenwater content
at depths
10, 20,100 cm on 241November,
10
December,and26
January,
respectively,
demonstrating
the
phase
changeofliquidwater.
the1 December, 10 December, and 26 January, respectively, demonstrating the phase change ofOnliquidothera decreasein theunfrozen
content water
belowcontent
the 100 below
cm depthwascm
only
about,
oronlylesswater.hand,On theother hand,a decrease
inwaterthe unfrozen
the 100depth
was
3 /cm3 . As temperatures increased, the frozen layer gradually thawed after late February.
than,
0.01
cm
3
3
about, or less than, 0.01 cm /cm . As temperatures increased, the frozen layer gradually thawed after
Comparedwiththe period
beforefreezing,
the soilwaterthecontentchangedafterchanged
the frozen
layerlate February.
Comparedwiththe period
beforefreezing,
soil watercontent
after
thecompletely
melted.
For
example,adecrease
in
the
soil
water
content
at
a
depth
of
100
cm
was
about
frozen layer completely melted. For example, a decrease in the soil water content at a depth of 100
3 /cm3 . Since water contents were not affected by external factors, such as
16%,
from
0.11416%,
to 0.096
cm wasabout
fromcm
0.114
to 0.096 cm3/cm3. Since water contents were not affected by external
rainfall,thisobservation
provesthat the freeze-thawprocess
contributed process
to the redistribution
factors, such as rainfall, this observation
proves thatthe freeze-thaw
contributed of
to soilthemoisture
in theofsoil
Meanwhile,the correlation
coefficients
betweensoil temperatures
redistribution
soilprofile.
moisture
in the soil profile.
Meanwhile,the correlation
coefficients
between andsoil
unfrozenwaterand
contents
werewater
generally
highwerewith generally
an averagehighvalue
of about
0.9 (asvalue
shownofinabout
Table0.9
4),temperatures
unfrozen
contents
withan average
indicating
trendsdecreasing
and trends
then increasing)in soil temperatures
and water contents(as shownsimilar
in Table4), (firstindicating
similar(first decreasing
and then increasing)in soilatcorresponding
depths.
The
minimum
value(0.754)occurred
at
the
100
cm
depth,
which occurred
was mainly
temperatures and water contents at corresponding depths. The minimum value (0.754)
atWater 2020, 12, 1261
Water 2020, 12, x FOR PEER REVIEW

11 of 18
11 of 18

the 100
depth,decrease
which wasmainlydue toa suddendecrease
thewhile
unfrozenwater contentdropped(Figuredue
to acmsuddenin theunfrozenwater
content(Figurein6c)
soil temperature6c) while soil
temperature
dropped gradually (Figure 3c).
gradually(Figure3c).

Figure 6. Spatial and temporal variations of the soil unfrozen water contents during the months of
Figure 6. Spatial and temporal variations of the soil unfrozen water contents during the months of (a)(a) November 2017, (b) December 2017, (c) January 2018, (d) February 2018, and (e) March 2018.
November 2017, (b) December 2017, (c) January 2018, (d) February 2018, and (e) March 2018.
Table 4. Correlation coefficients between soil temperatures and unfrozen water contents at various
Table 4.All
Correlation
coefficientssoil temperatures
depths.
values pass
the test atbetween
a 0.01 significancelevel. and unfrozen water contents at various
depths. All values pass the test at a 0.01 significance level.
0 cm
10 cm
20 cm
50 cm
100 cm 130 cm 200 cm 410 cm 630 cm 660 cm 800 cm

0 cm
10 0.918
cm 20 cm
50 cm
100 0.754
cm 1300.912
cm
0.907
0.906
0.887

200
cm
0.957

410
cm
0.919

630
cm
0.914

660 cm 0.906
800 cm

0.915

0.907 0.918
0.906
0.887
0.754
0.912
0.957
0.919
0.914
0.915
0.906
A wetter soil layer with the soil water content from 0.08 to 0.12 cm3 /cm3 existed between depthsof 80 and 230 cm, which may be attributed to two factors. First, due 3to thetemperature gradient,
A wetter soil layer with the soil water content from 0.08 to 0.12 cm /cm3 existed between depths
both liquid water and water vapor would flow upward and accumulate in this layer. Second,of 80 and 230 cm, which may be attributed to two factors. First, due to the temperature gradient, bothalthough the entire profile was mainly composed of sandy soil, the clay and silt fraction represented aliquid water and water vapor would flow upward and accumulate in this layer. Second, although theconsiderable proportion (approximately 10%) of this soil layer, resulting in higher water retention andentire profile was mainly composed of sandy soil, the clay and silt fraction represented a considerable
more restricted soil water movement in this layer.
proportion (approximately 10%) of this soil layer, resulting in higher water retention and more
Figure 7 shows the relationship between the unfrozen water content and subfreezing temperature,
restricted soil water movement in this layer.
and Table 5 lists the fitted parameters for both theoretical (Equation (5)) and empirical (Equation (6))
Figure 7 shows the relationship between the unfrozen water content and subfreezing
models. It is apparent that both models fitted observed data well at two different depths. The unfrozen
temperature, and Table 5 lists the fitted parameters for both theoretical (Equation (5)) and empiricalwater content began decreasing when soil temperature dropped below 0 ◦ C, and the temperature(Equation (6)) models. It is apparent that both models fitted observed data well at two different
range between 0 and −2 ◦ C can be regarded as an apparent phase transformation temperature interval
depths. The unfrozen water content began decreasing when soil temperature dropped below 0 °C,
and the temperature range between 0 and −2 °C can be regarded as an apparent phase transformation

Water 2020, 12, 1261

12 of 18

Water 2020, 12, x FOR PEER REVIEW

12 of 18

forthe Mu Us
sandyforsoil.
temperaturea decreasein thea decrease
unfrozeninwatercontent
temperatureinterval
the InMuthis
Us sandy
soil. In interval,
this temperatureinterval,
the unfrozen
accounted
for
over
75%
of
the
total
water
content.
Most
free
water
in
the
soil
matrixfroze
in thiswater content accounted for over 75% of the total water content. Most free water in the soil matrix
temperatureinterval, andinterval,
unfrozenwater
remainedin veryonly
small
poressmall
where
ice cannotbefroze in this temperatureandunfrozenwateronly
remainedin very
pores
where ice
easilyformed.
The
downward
trend
in
the
unfrozenwater
content
slowedwhentemperature
rangedcannot be easily formed. The downward trend in the unfrozen water content slowed when
◦ C when unfrozen water consisted mainly of film water and absorbed water.
between−2 and−4 between
temperature
ranged−2 and −4 °C when unfrozen water consisted mainly of film water and
Owing
to
this
restriction,
remaining
unfrozenremaining
water (close
to the residual
water tocontent)
cannot water
easilyabsorbed water. Owing to this restriction,
unfrozenwater (close
the residual
freeze,andthechangeintheunfrozenwater
content
becamerelatively
low
when
soil
temperature
content) cannot easily freeze, and the change in the unfrozen water content became relatively low
was
below
−4 ◦ C. Also, wasthe slopeof−4the°C.
decliningtrend
at aofdepthof 10 cm trend
depthatwas
steeperatwhensoil temperature
below
Also, theslope
the declining
a depthofthan
10 cm
a
depth
of
20
cm.
The
minimum
unfrozenwater
content
was
lower
at
a
depth
of
10
cm,
which
was
depth was steeper than at a depth of 20 cm. The minimum unfrozen water content was lower at amainlyto slightly
soildue
textures.
depth ofdue
10 cm,
which different
was mainly
to slightly different soil textures.

Figure 7. The relationship between the unfrozen soil water content and soil temperature at depths of
10 (a) and 20 (b) cm. Figures
Figures show
show experimental data, as well as the van Genuchten (Equation (5)) andempirical (Equation (6)) functions. R
R is
is the correlation coefficient, and RMSE is the root mean squareerror (cm33/cm/cm33).).
Table 5. The fitted parameters for the van Genuchten and empirical models.
Table 5. The fitted parameters for the van Genuchten and empirical models.
Depth (cm)Van Genuchten Model

Empirical Modelvan Genuchten Model
Empirical Model
θs
α
n
m
a
b
θr
θs
α
n
m
a
b
0.01
0.35
0.034
1.55
0.35
0.023
−0.2940.01 0.35
0.35 0.034
0.034 1.55
0.35
0.023
−0.294
0.014
1.61
0.38
0.024
−0.1410.014 0.35 0.034 1.61 0.38 0.024−0.141

Depth (cm) θr
10
10
20

20

Table
the freezing
freezing
Table 66 lists
lists thethe measured
measured total
total water
water content
content andand liquid
liquid water
water content
content data
data during
during theperiod.
The
results
indicated
that
the
total
soil
water
content
above
the
depth
of
50
cm
displayedanperiod. The results indicated that the total soil water content above the depth of 50 cm displayed anincreasing
trend, especially
especiallyin
inthe
theshallow
shallowdepth
depthofof2020
cm,
and
ratio
of liquid
water
content
toincreasing trend,cm,
andthetheratio
of liquid
water
content
to the
the
total
water
content
gradually
decreasedwithtimeelapsed.
This
phenomenon
occurredmainlytotal water content gradually decreased with time elapsed. This phenomenon occurred mainly
becausethefreezing
front
resultedin ainsharp
declinein thebecause aa decrease
decreasein
inthe
theunfrozen
unfrozenwater
watercontent
contentatat
thefreezing
front
resulteda sharpdecline
in
soil
water
potential,
causingtheunfrozenwater
fromdeepersoil
layerstoflow
upward,towardsthethe soil water potential, causing the unfrozen water from deeper soil layers to flow upward, towards
freezing
front,
andandthento freeze
there.
Similarly,
water
the freezing
front,
thento freeze
there.
Similarly,
watervapor
vaporwaswasflowing
flowingfrom
fromdeeper
deeper soil
soil layers
layersupwarddue
tothetemperature
gradient(asdiscussed
below
indetail).
upward due to the temperature gradient (as discussed below in detail).

Water 2020, 12, 1261

13 of 18

Table 6. The total soil water content θt (cm3 /cm3 ) measured using the drying method and the liquidsoil water content θl (cm3 /cm3 ) measured using the Hydra Probe II sensors at different soil depths andat different times.
Depth (cm)
10
20
30
40
50

19 December 2017θtθlθl /θt

θt

0.055
0.052
0.050
0.051
0.051

0.061
0.056
0.053
0.055
0.055

0.012
0.019
0.034

0.22
0.37
0.67

9 January 2018θl
θl /θt
0.013
0.018
0.033

0.21
0.32
0.60

θt

5 February 2018θl
θl /θt

0.063
0.059
0.055
0.056
0.054

0.011
0.018
0.027

0.17
0.31
0.50

3.3. Water Vapor Flux
The calculation by Equation (9) showed that the vapor density remained in near-saturated
conditions in the vadose zone during the non-freezing period, with the relative humidity reaching or
acceding 99%. When the soil was frozen, the relative humidity declined with a decrease in the soilpressure head and temperature (e.g., it was about 90% at −10 ◦ C), causing a reduction in the vapor
density. Moreover, the results indicate that variations in soil temperature had a great influence on thevapor density. For example, a temperature increase of 1 ◦ C would produce an increase of over 6% in
the vapor density.
Figure 8 shows the distribution of the vapor density in the soil profile at three typical dates:
before freezing (on 16 November), during the stable freezing stage (on 1 February), and after melting(on 25 March). Before the freezing period started, the vapor density was much smaller in the shallow
soil layer than below it, while the maximum value of 12.5 × 10−6 g/cm3 occurred at the 400 cm depth.
During the freezing period, the pressure head and soil temperatures in the frozen layer decreased
sharply, resulting in a low vapor density and an upward vapor density gradient in the vadose zone.
By comparison, the vapor density in the shallow layer increased significantly after melting, and the soil
depth between 100 and 200 cm had the lowest vapor density (around 6.7 × 10−6 g/cm3 ). Additionally,
the vapor density in the top 20 cm soil layer showed diurnal variations in all three cases, with an
upward gradient of the vapor density during nighttime and a downward gradient of the vapor density
duringWater the2020,daytime.
12, x FOR PEER REVIEW
14 of 18

Figure8. Thedepth
distribution
of ofthethevapor
densityonon
differentdays:(a)(a)
beforefreezing(16
November),Figure8. Thedepth
distribution
vapor
densitydifferent
days:
beforefreezing(November(b)16),
duringthestablefreezingstage(1
February),
and
(c)
after
melting(25
March).
(b) during the stable freezing stage (February 1), and (c) after melting (March 25).

Water 2020, 12, 1261

14 of 18

Figure 9 shows thermal vapor fluxes in the 10–20 and 50–100 cm soil layers on three typical days
calculated using Equation (12). It is apparent that the vapor flux is affected by temperature changes
in the shallow 20 cm depth. In this soil layer, on 16 November, water vapor flowed upward (about
0.003 cm/day at a depth of 10 cm) during most of the day, which was similar to variations on 1 February.
However, with an increase in air temperature during the daytime, soil temperatures at the 10 cm
depth were higher than at the 20 cm depth from 2:00 p.m. to 7:00 p.m. on 16 November and from1:00 p.m. to 8:00 p.m. on 1 February, respectively, resulting in downward vapor flow during thistime interval. Due to a rapid increase in air temperature after melting, the downward vapor flux
sustained from 11:00 a.m. to 10:00 p.m. on 25 March. The maximum downward vapor flux was−0.017 cm/day, which was almost five times higher than before. For depths between 50 and 100 cm,
soil temperatures fluctuated only slightly during the day, resulting in small vapor fluxes. While vapor
fluxesFigure
were8.upwardondistribution
16 November
0.002 oncm/day)and1 February(about (November
0.001 cm/day),
The depthof the (aboutvapor density
differentdays:(a) before freezing
they were
downward
on
25
March
(about−0.001
cm/day).
16), (b) during the stable freezing stage (February 1), and (c) after melting (March 25).

Figure 9. Diurnal variations of the thermal vapor flux at two layers (10–20 and 50–100 cm) on different
days: (a)(a) November
16 November,(b)February
1 February,and(c)(c)March
25 March.
qv10–20
qv50–100are calculated
and andqv50–100are calculated
vapor vapor
fluxes16, (b)1, and25. qv10–20
fluxessoil layers
10–20
and 50–100
cm, respectively.
Positive
indicate
upwardflow,
in soil inlayers
10–20 and
50–100cm, respectively.
Positivevaluesvalues
indicate
upwardvapor vapor
flow, while
while
negative
indicate
downward
negativevaluesvalues
indicate
downward
vaporvapor
flow. flow.

Assoil layer,
layer, theflux showed
showed aa decreasing
decreasing trend
trend with
with increasing
increasing soil
soil depths.
depths.
As forfor thethe deep
deep soil
the vaporvapor flux
For
example,
owing
to
the
relatively
lowtemperature
gradients,
thethermalvapor
fluxesat 800the
For example, owing to the relatively low temperature gradients, the thermal vapor fluxes at the−4 , 1.6 × 10−4 , and 2.9 × 10−4 cm/day on 16 November, 1 February,
800
cm
depth
were
−1.4
×
10
−4
−4
−4
cm depth were −1.4 × 10 , 1.6 × 10 , and 2.9 × 10 cm/day on November 16, February 1, and March
and25 March, respectively.
25, respectively.
4. Discussion
4.1. Characteristics of the Freeze-Thaw Process
The freezing depth, defined as the depth with zero temperature, is often used to characterize thefreeze-thaw process [9,26], as shown in Figure 10. According to the measured soil temperature andmoisture data, the total freeze-thaw process, with the maximum freezing depth to be at a depth of
118 cm, lasted 112 days (from 19 November 2017, to 10 March 2018), and could be divided into threemain stages: the initial freezing stage, the downward freezing stage, and the thawing stage.

The freezing depth, defined as the depth with zero temperature, is often used to characterize thefreeze-thaw process [9,26], as shown in Figure 10. According to the measured soil temperature andmoisture data, the total freeze-thaw process, with the maximum freezing depth to be at a depth of
118 cm,Water
2020,lasted
12, 1261112 days (from 19 November 2017, to 10 March 2018), and could be divided into15threeof 18
main stages: the initial freezing stage, the downward freezing stage, and the thawing stage.

FigureFigure 10.
10. The freezing depth
depth curve
curve during
during thethe freeze-thaw
freeze-thaw process.
process.

DuringDuring thethe initial
initial freezing
freezing stage,stage, thethe soil
soil above
above thethe 2020 cm
cm depth
depth experienced
experienced several
several freeze-thaw
freeze-thaw
cycles,
which
could
be
characterizedbynightly
freezing
anddaily
thawing.
The
diurnalvariations
cycles, which could be characterized by nightly freezing and daily thawing. The diurnal variations ofof◦ C) were larger than at a depth of 20 cm (0.6 ◦ C). The diurnal
temperaturetemperature at
at aa depth
depth of
of 10
10 cm
cm (2.5(2.5 °C)
were larger than at a depth of 20 cm (0.6 °C). The diurnal
amplitude
amplitude ofof thethe unfrozen
unfrozen water
water content
content at
at aa depth
depth of
of 10
10 cm
cm waswas considerably
considerably larger
larger than
than at
at aa depth
depth
of
20
cm,
indicatingthattheprocess
of
the
water
phase
transition
was
more
intenseinthe
topof 20 cm, indicating that the process of the water phase transition was more intense in the top 10
10 cm
cm
depths.
Duringthe
thedownward
downward
freezing
stage,
freezing
process
accelerated
andfreezing
the freezing
depths. Duringfreezing
stage,
the thefreezing
process
accelerated
and thedepth
depth
gradually
deepened.
Unliketheshallower
soil in shallower
depths,the soilthe 20 which
cm depth,
gradually
deepened.
Unlike thesoil indepths, thesoil below
the below
20 cm depth,
was
which
was
less
influencedbysolar
radiation,did
notundergomultiplefreeze-thaw
cycles,butinsteadless influenced by solar radiation, did not undergo multiple freeze-thaw cycles, but instead directly
directly
frozea witha freezingrate
1.3 cm/day.
daily
range
shallowsoil
soil temperaturestemperatures further
froze withfreezing
rate of
1.3ofcm/day.
TheThedaily
range
ofofshallow
further◦ C at a depth of 10 cm), and the unfrozen water content continuously declined
increased(e.g.,
5.1
increased (e.g., 5.1 °C at a depth of 10 cm), and the unfrozen water content continuously declined asastemperatures
dropped.
At thelowestpoint,
unfrozenwater
content
decreasedby 79%
temperatures
dropped.
At thelowestpoint,
the theunfrozenwater
content
decreasedby 79%
and and65%
3
3
3 /cm3 ),
65%
in
soil
depths
of
10
cm
(from
0.053
to
0.011
cm
/cm
)
and
20
cm
(from
0.052
to
0.018
cm
3
3
3
in soil depths of 10 cm (from 0.053 to 0.011 cm /cm ) and 20 cm (from 0.052 to 0.018 cm /cm3),
respectively,
100 cmcm (51%).
(51%).
respectively,while
while thethe decrease
decrease was
was relatively
relatively small
small inin soilsoil depths
depths of
of 50
50 cm
cm (53%)(53%) andand 100
The
both
the
top
and
thethebottomof thefrozenlayer,
andandthe shallow
soil
Thethawing
thawingprocessprocessappeared
appearedfrom
from
both
the
top
and
bottomof thefrozenlayer,
the shallow
experienced
freeze-thaw
cyclesonce
again.
The
entirethawing
processlasted
only
25
days,
which
was
soil experienced freeze-thaw cycles once again. The entire thawing process lasted only 25 days, which
strikingly
shortershorter
than thefreezing
stage. Dueto a significant
increase in springtime
temperature,
was strikingly
than
the freezing
stage.
Due to a significant
increase inair
springtime
air
soil
temperatures
rapidly,
resultingin a melting
rate ofto 4.7 cm/day,temperature,soilincreased
temperatures
increased
rapidly,
resultingin up
a meltingrate ofwhich
up tosubstantially
4.7 cm/day,exceeded
the freezingexceeded
rate [3]. the freezing rate [3].
which substantially4.2. Vapor Migration during the Freezing-Thawing Period
4.2. Vapor Migration during the Freezing-Thawing Period
Due to generally low liquid water contents in arid areas, vapor flow performs a critical function
Due to generally low liquid water contents in arid areas, vapor flow performs a critical functionof the hydrological cycle in the vadose zone during the seasonal freezing-thawing process. At theof the hydrological cycle in the vadose zone during the seasonal freezing-thawing process. At the
beginning of the observation period, the layer between the soil surface and the 200 cm soil depth
beginning of the observation period, the layer between the soil surface and the 200 cm soil depth
became a low soil temperature zone (Figure 3a). As the freezing process started, temperatures in thisbecame a low soil temperature zone (Figure 3a). As the freezing process started, temperatures in thislayer kept a declining trend. The decreasing rate (from November 1 to 1 February) of soil temperatures
layer kept a declining trend. The decreasing rate (from November 1 to February 1) of soil
at depths of 50, 100, and 200 cm were 0.16, 0.14, and 0.11 ◦ C/day, respectively. Driven by the upward
temperatures at depths of 50, 100, and 200 cm were 0.16, 0.14, and 0.11 °C/day, respectively. Driven
temperature gradient, water vapor flowed upward towards the soil surface and evaporated into theby the upward temperature gradient, water vapor flowed upward towards the soil surface andatmosphere when the frozen layer was thin. Subsequently, it condensed gradually in the seasonallyevaporated into the atmosphere when the frozen layer was thin. Subsequently, it condensed
frozen layer, contributing to the ice formation [5,23]. The cumulative upward vapor flux at the 10 cm
gradually in the seasonally frozen layer, contributing to the ice formation [5,23]. The cumulative
depth was about 0.56 cm during this period, which affected the soil moisture distribution at such a low
upward vapor flux at the 10 cm depth was about 0.56 cm during this period, which affected the soil
liquid water content.
moisture distribution at such a low liquid water content.
As air temperature increased in late February, soil temperature in the shallow layer increased
rapidly, causing the vapor flux to move mainly downward. The cumulative downward vapor flux atthe 10 cm depth was approximately 0.17 cm during the thawing period, and the flux maintained arapidly increasing trend with a gradual increase in temperature gradients. At this time, temperatures
at soil depths between 100 and 200 cm became relatively low (Figure 3e), and water vapor migrated
from both above and below to the relatively low-temperature layer, which contributed to the formationof the wet layer.

Water 2020, 12, 1261

16 of 18

5. Conclusions
Based on in-situ observations in the Mu Us Desert, the changes in soil temperatures and watercontents during the freeze-thaw period were studied in this manuscript. The results showed that soil
temperature displayed a decreasing trend during the freezing period and the depth with the maximum
soil temperature in the vadose zone kept increasing from 400 cm to 800 cm. On the contrary, the soil
layer above the 100 cm depths displayed negative (downward) temperature gradients during themelting period. The soil water content in the profile was generally low (only 0.06 cm3 /cm3 ) before
freezing, except for the soil horizon between 80 and 230 cm. Both theoretical and empirical models
captured the relationship between the unfrozen water content and subfreezing temperature well.
The total water content in the frozen layer increased due to the upward soil water flux from deeper soil
layers. The entire freeze-thaw process can be divided based on the measured data into three stages,
including the initial freezing stage, the downward freezing stage, and the thawing stage. According to
Fick’s Law, the thermal vapor flux in the shallow 20 cm depth showed markedly diurnal variations,with vapor flowing upward during the nighttime and downward during the daytime, while themagnitude of the vapor flux gradually decreased with increasing soil depths. Calculation results
indicated that vapor moved upward towards the frozen layer and contributed to the ice formation
during the freezing process, while it flowed downward during the thawing process and contributed to
the formation of the wet layer.
Although this study was carried out in northwestern China, similar results could be expected
for other regions with similar soil and climate conditions. Further studies will focus on quantitative
calculations of the fully coupled movement of water, vapor, and heat during the freeze-thaw process
in the deep vadose zone, and will evaluate the influence of gradients of the matric potential andtemperature on transport processes of liquid water and water vapor.
Author Contributions: C.Z. and X.L. were responsible for the original concept, and C.Z. wrote the original
manuscript. C.Z., C.S., and H.L. conducted experiments and data analysis. Y.L., J.Š., and Y.Z. provided valuablesuggestions, and J.Š. edited the manuscript. All authors have read and agreed to the published version of
the manuscript.
Funding: This research was supported by the National Natural Science Foundation of China (Nos. 41877179
and 41630634), the Fundamental Research Foundation of the Central Universities, CHD (Nos. 300102298706 and300102290719), the Shaanxi Water Conservancy Science and Technology plan project (2019slkj-18), the Fund projectof Shaanxi Key Laboratory of Land Consolidation and the program of China Scholarship Council.
Acknowledgments: The authors are grateful to the editor and reviewers for providing positive and constructive
comments for this paper.
Conflicts of Interest: The authors declare no conflict of interest.

References
1.
2.

3.

4.
5.

Chen, J.; Gao, X.; Zheng, X.; Miao, C.; Liu, P.; Du, Q.; Xu, Y. Transformation between phreatic water and soil
Water during freeze–thaw periods. Water 2018, 10, 376. [CrossRef]
Wen, Z.; Ma, W.; Feng, W.; Deng, Y.; Wang, D.; Fan, Z.; Zhou, C. Experimental study on unfrozen watercontent and soil matric potential of Qinghai-Tibetan silty clay. Environ. Earth Sci. 2012, 66, 1467–1476.
[CrossRef]
Cheng, Q.; Sun, Y.; Jones, S.B.; Vasilyev, V.I.; Popov, V.V.; Wang, G.; Zheng, L. In situ measured and simulatedseasonal freeze-thaw cycle: A 2-year comparative study between layered and homogeneous field soil profiles.
J. Hydrol. 2014, 519, 1466–1473. [CrossRef]
Wang, W.; Li, H.; Wang, J.; Hao, X. Water vapor from western Eurasia promotes precipitation during thesnow season in northern Xinjiang, a typical arid region in central Asia. Water 2020, 12, 141. [CrossRef]
Zhang, Y.; Cheng, G.; Li, X.; Han, X.; Wang, L.; Li, H.; Chang, X.; Flerchinger, G.N. Coupling of a simultaneousheat and water model with a distributed hydrological model and evaluation of the combined model in acold region watershed. Hydrol. Process. 2012, 27, 3762–3776. [CrossRef]

Water 2020, 12, 1261

6.
7.
8.
9.
10.
11.

12.
13.

14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.

25.
26.
27.
28.
29.
30.

17 of 18

Hayashi, M. The cold vadose zone: Hydrological and ecological significance of frozen-soil processes. Vadose
Zone J. 2013, 12, 37–44. [CrossRef]
Yin, X.; Liu, E.; Song, B.; Zhang, D. Numerical analysis of coupled liquid water, vapor, stress and heattransport in unsaturated freezing soil. Cold Reg. Sci. Technol. 2018, 155, 20–28. [CrossRef]
Zheng, D.; Rogier, V.D.V.; Su, Z.; Wen, J.; Wang, X.; Yang, K. Impact of soil freeze-thaw mechanism on therunoff dynamics of two Tibetan rivers. J. Hydrol. 2018, 563, 382–394. [CrossRef]
Monrabal-Martinez, C.; Scibilia, E.; Maus, S.; Muthanna, T.M. Infiltration response of adsorbent amended
filters for stormwater management under freezing/thawing conditions. Water 2019, 11, 2619. [CrossRef]
Steelman, C.M.; Endres, A.L.; Kruk, J.V.D. Field observations of shallow freeze and thaw processes using
high-frequency ground-penetrating radar. Hydrol. Process. 2010, 24, 2022–2033. [CrossRef]
Zhang, Y.; Ma, W.; Wang, T.; Cheng, B.; Wen, A. Characteristics of the liquid and vapor migration ofcoarse-grained soil in an open-system under constant-temperature freezing. Cold Reg. Sci. Technol. 2019, 165,
102793. [CrossRef]
Guo, D.; Yang, M.; Wang, H. Characteristics of land surface heat and water exchange under different soil
freeze/thaw conditions over the central Tibetan Plateau. Hydrol. Process. 2011, 25, 2531–2541. [CrossRef]
Yi, J.; Zhao, Y.; Shao, M.; Zhang, J.; Cui, L.; Si, B. Soil freezing and thawing processes affected by the different
landscape in the middle reaches of the Heihe River Basin, Gansu, China. J. Hydrol. 2014, 519, 1328–1338.
[CrossRef]
Romanovsky, V.E.; Osterkamp, T.E. Effects of unfrozen water on heat and mass transport processes in theactive layer and permafrost. Permafr. Periglac. 2000, 11, 219–239. [CrossRef]
Watanabe, K.; Wake, T. Measurement of unfrozen water content and relative permittivity of frozen unsaturated
soil using NMR and TDR. Cold Reg. Sci. Technol. 2009, 59, 34–41. [CrossRef]
Lu, C.; Qin, W.; Zhao, G.; Zhang, Y.; Wang, W. Better-fitted probability of hydraulic conductivity for a silty
clay site and its effects on solute transport. Water 2017, 9, 466. [CrossRef]
Zhao, X.; Xu, S.; Liu, T.; Qiu, P.; Qin, G. Moisture distribution in sloping black soil farmland during thefreeze–thaw period in northeastern China. Water 2019, 11, 536. [CrossRef]
Hansson, K.; Šimůnek, J.; Mizoguchi, M.; Lundin, L.C.; Th. van Genuchten, M. Water flow and heat transport
in frozen soil: Numerical solution and freeze-thaw applications. Vadose Zone J. 2004, 3, 527–533. [CrossRef]
Iwata, Y.; Hirota, T. Monitoring over-winter soil water dynamics in a freezing and snow-covered environment
using a thermally insulated tensiometer. Hydrol. Process. 2005, 19, 3013–3019. [CrossRef]
Iwata, Y.; Hirota, T.; Hayashi, M.; Suzuki, S.; Hasegawa, S. Effects of frozen soil and snow cover on cold-season
soil water dynamics in Tokachi, Japan. Hydrol. Process. 2010, 24, 1755–1765. [CrossRef]
Kozlowski, T.; Nartowska, E. The Unfrozen water content in representative bentonites of different origin
subjected to cyclic freezing and thawing. Vadose Zone J. 2013, 12, 1196–1206. [CrossRef]
Henry, H.A.L. Soil freeze-thaw cycle experiments: Trends methodological weaknesses and suggested
improvements. Soil Biol. Biochem. 2007, 39, 977–986. [CrossRef]
Kozlowski, T. Some factors affecting supercooling and the equilibrium freezing point in soil-water systems.
Cold Reg. Sci. Technol. 2009, 59, 25–33. [CrossRef]
Zeng, Y.; Su, Z.; Wan, L.; Yang, Z.; Zhang, T.; Tian, H.; Shi, X.; Wang, X.; Cao, W. Diurnal pattern of the drying
front in desert and its application for determining the effective infiltration. Hydrol. Earth Syst. Sci. 2009, 13,
703–714. [CrossRef]
Zhou, Y.; Wang, X.; Han, P. Depth-dependent seasonal variation of soil water in a thick vadose zone in the
Badain Jaran Desert, China. Water 2018, 10, 1719. [CrossRef]
Zeng, Y.; Su, Z.; Wan, L.; Wen, J. Numerical analysis of air-water-heat flow in the unsaturated soil: Is it
necessary to consider air flow in land surface models? J. Geophys. Res. 2011, 116, D20107. [CrossRef]
Zhang, S.; Teng, J.; He, Z.; Sheng, D. Importance of vapor flow in unsaturated freezing soil: A numerical
study. Cold Reg. Sci. Technol. 2016, 126, 1–9. [CrossRef]
Saito, H.; Šimůnek, J.; Mohanty, B. Numerical analyses of coupled water, vapor, and heat transport in thevadose zone. Vadose Zone J. 2006, 5, 784–800. [CrossRef]
Zeng, Y.; Su, Z.; Wan, L.; Wen, J. A simulation analysis of the advective effect on evaporation using a
two-phase heat and mass flow model. Water Resour. Res. 2011, 47, W10529. [CrossRef]
Yu, L.; Zeng, Y.; Wen, J.; Su, Z. Liquid-vapor-air flow in the frozen soil. J. Geophys. Res. 2018, 123, 7393–7415.
[CrossRef]

Water 2020, 12, 1261

31.
32.
33.

34.

35.

36.

37.
38.

39.
40.

41.
42.
43.
44.
45.
46.
47.

18 of 18

Xiang, L.; Yu, Z.; Chen, L.; Mon, J.; Lu, H. Evaluating coupled water, vapor, and heat flows and their influenceon moisture dynamics in arid regions. J. Hydrol. Eng. 2012, 17, 565–577. [CrossRef]
Zhao, Y.; Huang, M.; Horton, R.; Liu, F.; Peth, S.; Horn, R. Influence of winter grazing on water and heat flow
in seasonally frozen soil of Inner Mongolia. Vadose Zone J. 2013, 12, 105–115. [CrossRef]
Bai, C.; He, X.; Tang, H.; Shan, B.; Zhao, L. Spatial distribution of arbuscular mycorrhizal fungi, glomalinand soil enzymes under the canopy of Astragalus adsurgens Pall. In the Mu Us sandland, China. Soil Biol.
Biochem. 2009, 41, 941–947. [CrossRef]
Cheng, D.; Li, Y.; Chen, X.; Wang, W.; Hou, G.; Wang, C. Estimation of groundwater evaportranspiration
using diurnal water table fluctuations in the Mu Us Desert, northern China. J. Hydrol. 2013, 490, 106–113.
[CrossRef]
Li, H.; Lu, Y.; Zheng, C.; Yang, M.; Li, S. Groundwater level prediction for the arid oasis of northwest China
based on the artificial bee colony algorithm and a back-propagation neural network with double hidden
layers. Water 2019, 11, 860. [CrossRef]
Liu, X.; Guo, C.; He, S.; Zhu, H.; Li, J.; Yu, Z.; Qi, Y.; He, J.; Zhang, J.; Muller, C. Divergent gross nitrogentransformation paths in the topsoil and subsoil between abandoned and agricultural cultivation land in
irrigated areas. Sci. Total Environ. 2020, 716, 137148. [CrossRef]
Zheng, C.; Lu, Y.; Guo, X.; Li, H.; Sai, J.; Liu, X. Application of HYDRUS-1D model for research on irrigation
infiltration characteristics in arid oasis of northwest China. Environ. Earth Sci. 2017, 76, 785–794. [CrossRef]
Zhang, M.; Wen, Z.; Xue, K.; Chen, L.; Li, D. A coupled model for liquid water, water vapor and heattransport of saturated–unsaturated soil in cold regions: Model formulation and verification. Environ. Earth
Sci. 2016, 75, 701–719. [CrossRef]
Cheng, Q.; Sun, Y.; Xue, X.; Guo, J. Characteristics for estimation of soil moisture characteristics using adielectric tube sensor. Soil Sci. Soc. Am. J. 2014, 78, 133–138. [CrossRef]
Sun, Y.; Cheng, Q.; Xue, X.; Fu, L.; Chai, J.; Meng, F.; Lammers, P.S.; Jones, S.B. Determining in-situ soil
freeze-thaw cycle dynamics using an access tube-based dielectric sensor. Geoderma 2012, 189–190, 321–327.
[CrossRef]
Radcliffe, D.E.; Šimůnek, J. Soil Physics with HYDRUS: Modeling and Applications; CRC Press: Boca Raton, FL,
USA, 2010; p. 373.
Van Genuchten, M.T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.
Soil Sci. Soc. Am. J. 2010, 44, 892–898. [CrossRef]Šimůnek, J.; van Genuchten, M.T.; Šejna, M. Development and applications of the HYDRUS and STANMOD
software packages and related codes. Vadose Zone J. 2008, 7, 587–600. [CrossRef]
Xu, X.; Wang, J.; Zhang, L. Frozen Soil Physics; Science Press: Beijing, China, 2010; p. 43.
Camillo, P.J.; Gurney, R.; Schmugge, T. A soil and atmospheric boundary layer model for evapotranspirationand soil moisture studies. Water Resour. Res. 1983, 19, 371–380. [CrossRef]
Millington, R.J.; Quirk, J.P. Permeability of porous solids. Trans. Faraday Soc. 1961, 57, 1200–1207. [CrossRef]
Cass, A.; Campbell, G.S.; Jones, T.L. Enhancement of thermal water vapor diffusion in soil. Soil Sci. Soc. Am.
J. 1984, 48, 25–32. [CrossRef]© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open accessarticle distributed under the terms and conditions of the Creative Commons Attribution(CC BY) license (http://creativecommons.org/licenses/by/4.0/).