Journal of Hydrology 627 (2023) 130401

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Water transfer mechanisms and vapor flow effects in seasonally frozen soils
Ce Zheng a, b, c, Yunfei Chen a, b, c, Wande Gao a, b, c, Xinyang Liang d, Jiří Šimůnek e,
Xiuhua Liu a, b, c, *aSchool of Water and Environment, Chang’an University, Xi’an 710064, China
Key Laboratory of Subsurface Hydrology and Ecological Effects in Arid Region of Ministry of Education, Chang’an University, Xi’an 710064, Chinac
Key Laboratory of Eco-hydrology and Water Security in Arid and Semi-arid Regions of Ministry of Water Resources, Chang’an University, Xi’an 710064, Chinad
Shaanxi Weihe Ecological Zone Protection Center, Xi’an 710043, Chinae
Department of Environmental Sciences, University of California Riverside, Riverside 92521, USAb

A R T I C L E I N F O

A B S T R A C T

This manuscript was handled by Yuefei Huang,
Editor-in-Chief, with the assistance of Jan
Vanderborght, Associate Editor

In seasonally frozen soils, liquid water and water vapor movement are critical in numerous environmental andengineering areas. However, the process of soil water movement becomes complicated due to the coexistence ofliquid water, water vapor, and ice, and the underlying mechanisms are not yet fully understood. Based on thetheory of coupled water, vapor, and heat transport, the modified HYDRUS-1D freezing module was used toinvestigate the mechanisms of soil water transfer in a cold region. Compared with its previous version, the newly
developed model overcomes the numerical instability caused by the phase transformation between liquid waterand ice. The verification results showed that the modified model was numerically stable and highly accurate,
indicating its suitability for the study region. The variations in soil temperature, water content, water storage,
evaporation, and drainage displayed typical seasonal patterns during the analyzed period, reflecting the effects offreeze–thaw cycles and rainfall-infiltration processes. Liquid flow, driven by a water potential gradient, played amajor role under most environmental conditions, except when the soil was frozen when it became negligible. In
comparison, vapor flow, driven by a temperature gradient, dominated soil water flow when ice was present andaccounted for a considerable proportion of the total water flux in dry conditions, proving its significance for
affecting soil water transport. These results clarified the mechanisms of soil liquid water and water vapor
transfers in seasonally frozen soils and provided more insights for deepening the hydrological cycle theory.

Keywords:
Partially frozen unsaturated soil
Coupled water, vapor, and heat transport
Vapor flux
Soil water balance
Modified HYDRUS-1D freezing module

1. Introduction
Regions characterized by permafrost and seasonal frost are exten­
sively distributed in the middle and high latitudes and account for
approximately 55–60 % of the Earth’s surface in the Northern Hemi­
sphere (Hayashi, 2013; Biskaborn et al., 2019). Frozen soils are major
cryosphere components and significantly affect surface processes andthe hydrological cycle in cold regions (Fig. 1), considerably impacting
global climate change (Chadburn et al., 2017; Obu et al., 2019).
Abundant forests, cultivated land, and mineral resources exist inseasonally frozen regions vital to many human activities. Seasonal
freezing-thawing cycles in the abovementioned regions substantially
impact many engineering and agricultural applications (Kreyling et al.,2012; Kurylyk and Watanabe, 2013; Zhao et al., 2013; He et al., 2018).
For instance, the freezing drought stress and infiltration of melting waterin seasonally frozen soils affect and even hinder vegetation growth,which may further change the stability of local ecosystems (Stuurop
et al., 2022). As a result, the freezing-thawing process in cold regions,which leads to numerous environmental and engineering problems, hasbeen drawing increased attention and research interest (Iwata et al.,2013).
As an important component of the Earth’s critical zone, the vadosezone not only serves as an interacting link between the atmosphere,hydrosphere, and lithosphere, but also provides a place for the exchangeof mass and energy in the soil–plant-atmosphere continuums (Banwart
et al., 2011). Soil and water are the vadose zone’s core elements, and soil
water flow processes have long been investigated (Zhao et al., 2020).
Philip and de Vries (1957) pioneered a theory describing the coupled
movement of liquid water, water vapor, and heat in non-isothermal
conditions. Subsequently, many other researchers have been working
during the past several decades on advancing the understanding of these
coupled processes (e.g., Cass et al., 1984; Chung and Horton, 1987; Saito* Corresponding author.
E-mail address: liuxh68@chd.edu.cn (X. Liu).
https://doi.org/10.1016/j.jhydrol.2023.130401
Received 19 May 2023; Received in revised form 28 August 2023; Accepted 13 October 2023
Available online 2 November 2023
0022-1694/© 2023 Elsevier B.V. All rights reserved.

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

designed, user-friendly GUI and various complementary features (e.g.,parameter optimization routines and flexible boundary conditions), the
HYDRUS-1D software package has been broadly adopted to investigate
the coupled transport of water, vapor, and heat. Scanlon et al. (2003)
were the first to consider the effects of vapor flow using the HYDRUS-1D
code, demonstrating the importance of vapor flow on soil water transfer
in arid areas. Saito et al. (2006) further modified HYDRUS-1D byimplementing the surface energy balance conditions, which was subse­
quently included in the standard HYDRUS-1D model (Šimůnek et al.,2008). Since then, extensive numerical research has been conducted on
the coupled movement of liquid water, water vapor, and heat using
HYDRUS-1D. The characteristics of liquid water and water vapor flow,
driven by soil water potential and temperature gradients, as well as thequantitative contribution of vapor flow to soil water flow under dry
conditions, have been well illustrated (e.g., Wang, 2015; Huang et al.,2016; Du et al., 2017; Zhou and Zhao, 2021).
Despite the significant progress in understanding the coupled
movement of liquid water, water vapor, and energy in soils, there is alack of clarity in two research areas on the mechanisms of soil watertransfer in seasonally frozen soils. First, due to complicated interactions
between water in the liquid, gas, and solid phases in unsaturated frozen
soils, studies on soil water movement during an entire hydrological year,
including freezing and non-freezing periods, are limited. The long-termchanges in soil hydrothermal dynamics, water balances, and flux den­sities have not been fully clarified. Second, the specific role of vapor flow
in soil water movement under different conditions (i.e., freezing,
thawing, drought, and rainy conditions) and at different depths should
also be clarified. The implications and effects of vapor flow on local
ecosystems and engineering activities should be explained using thisinformation.
The Mu Us Sandy Land, a typical seasonally frozen area in North­western China, was selected as the research site to address these researchconcerns. The data regarding soil characteristics and temporal changes
in water contents and temperatures were collected during long-term
field observations and experiments. Subsequently, the modified
HYDRUS-1D freezing code simulating the coupled movement of liquidwater, water vapor, and heat was used to explore the related transfer
processes in detail. The specific objectives of this study were (i) to

et al., 2006). A growing number of studies suggested that vapor flow
contributes considerably to mass and energy transfer, particularly inarid and semiarid areas characterized by low soil moisture (Bittelli et al.,2008; Xiang et al., 2012; Yu et al., 2018). Non-isothermal vapor flow is
considered a vital component of total soil water flow. Besides, the phasechange between liquid water and water vapor can release or store
considerable latent energy, and vapor flow can significantly affect thesoil energy balance (Cahill and Parlange, 1998). During the freezing
period, phase changes between ice and liquid water occur when the soil
temperature reaches the freezing point, and ice lenses form in soil pores,
significantly blocking seepage paths and affecting the soil’s permeability(Zhang et al., 2023b). Besides being affected by the matric potential andgravitational gradients, soil water and energy are also transferred fromthe warm end toward the cold end by temperature gradients. Laboratory
tests suggest that thermally driven vapor flow can facilitate ice forma­tion and even dominate moisture transfer when the soil is frozen (Zhang
et al., 2019). It is clear that the transport of liquid water, water vapor,and heat are closely related and strongly influence each other underboth freezing and non-freezing conditions (Yu et al., 2018; Zhang et al.,2023a).
Due to the difficulty and inconvenience of field monitoring and ex­
periments, numerical simulations have been increasingly employed in
relevant research. Algorithms and simulation accuracy have beensignificantly improved and optimized in recent decades with the rapid
development of computer technology (Bao et al., 2016; Hu et al., 2020).
Harlan (1973) first proposed the coupled water and heat transportmodel for unsaturated freezing soil, and the calculation results indicatedthat the freezing process can lead to the redistribution of soil water.
Since then, several numerical models, such as SHAW, CoupModel,
HYDRUS-1D, STEMMUS, and IHTSM, have been developed andemployed to address specific issues (Kurylyk and Watanabe, 2013; Baoet al., 2016). For instance, the improved hydro-thermal-salt-mechanical(IHTSM) model can be used to simulate the coupled water, vapor, heat,and solute migration under unidirectional freezing conditions,providing a reference for the destruction process analysis of the harshgeological environment in cold, arid, and saline areas (Zhang et al.,2021). Various physical phenomena can be investigated with such amodel, better elucidating the hydrological cycle theory. Due to the well-Fig. 1. Conceptual model of hydrological processes in frozen soil regions. Frozen soil affects many soil activities in typical ecosystems, such as snowmelt infiltrationand frozen drought stress.
2

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

analyze the temporal variations of soil hydrothermal variables andwater balance components during the entire hydrological year, (ii) toinvestigate the characteristics of water and vapor transfer in isothermal
conditions, and (iii) to clarify the significance and implications of vapor
flow in seasonally frozen soils.

samples collected from the study site. The water content sensors werebased on the Frequency Domain method, given the significantly higher
permittivity of liquid water (≈78) than ice (≈3.4), soil matrix (≈3-5),
and air (≈1) (Cheng et al., 2014). Accordingly, soil water content
measurements in the frozen conditions represent the unfrozen watercontent (Hu et al., 2020). Furthermore, meteorological variables,
including precipitation, air temperature, relative humidity, and windspeed, were monitored at the study site. The temporal changes of thesevariables during the entire hydrological year, starting from November 1,
2017, to October 31, 2018, are shown in Fig. 3. The total rainfall amountduring the analyzed year reached 660.8 mm, far exceeding the average
year amount, and over 95 % of rain fell from April to October when thesoil was not frozen.

2. Materials and methods
2.1. In-situ observations and experiments
2.1.1. Study site
China has been reported as the third largest country worldwide
regarding frozen soil area, of which approximately 72 % is permafrost
and seasonally frozen soil (Wang et al., 2017). As depicted in Fig. 2, mostpermafrost is found in the middle and low latitudes of the Qinghai-Tibet
Plateau region, while seasonally frozen soil is widespread in areas withlatitudes higher than 24◦ . The Mu Us Sandy Land is located between the
Loess Plateau and Ordos Plateau and belongs to a typical seasonallyfrozen region (Gao et al., 2022). The Yu Lin Desert Ecosystem Research
Station of the Mu Us Desert (E: 109◦ 42′29″, N: 38◦ 23′19″N) was selected
as the study site. As indicated by the data collected from a nearby
meteorological station, the climate of this region can be defined as coldand dry, with average annual precipitation and air temperature of
approximately 419 mm and 9.1 ◦ C, respectively (Fig. 3). Generally, soils
are frozen from November till March of the following year, with amaximum frost depth of over 1 m. Given the effects of the seasonal
freeze–thaw process, the soil profile with a depth of 200 cm was selected
to explore the soil liquid water and water vapor transfer mechanisms inthis area.

2.2. Numerical simulations
2.2.1. Theoretical equations
Soil water flow and heat transport were very different under theeffects of the freeze–thaw process than during the non-freezing period.
When the soil is frozen, liquid water, ice, and water vapor coexist in soil
pores, as depicted in Fig. 4 (Mohammed et al., 2018; Yu et al., 2020).
The percentage of liquid water and ice changes according to soil tem­
perature variations. Due to the ice presence, both hydraulic and thermal
conductivities dramatically change, thus affecting soil water and heattransfer. Based on the mass and energy conservation equations, soil
water flow and heat transport in a variably saturated rigid soil are given
as follows (Hansson et al., 2004; Zheng et al., 2021):()∂θl ρv ∂θv ρi ∂θi ∂
∂h
∂T
∂h
∂T
+
+
=
Klh + Klh + KlT + Kvh + KvT(1)∂t ρl ∂t ρl ∂t ∂z
∂z
∂z
∂z
∂z

2.1.2. Soil properties and meteorological variables
At the experimental site, soil physical properties (e.g., particle
composition and bulk density) were determined by collecting soil sam­ples at the 20-cm interval from the ground surface down to a depth of
200 cm. As listed in Table 1, sand is the main textural component at the
study site, and the soil profile consists of three layers. Soil temperatures
and water contents were measured at depths of 10, 20, 50, 100, 130, and200 cm every 10 min using the Stevens Hydro Probe II soil moisturesensors. Before installation, all sensors were calibrated using the soil

(

)

∂Cp T
∂θ i
∂θv ∂
∂T
∂ql T
∂qv T
∂qv
− Lf ρi + L0
=
λ
− Cw
− Cv
− L0
− Cw ST
∂z
∂t
∂t
∂t
∂z
∂z
∂z
∂z
(2)
where θl, θv, and θi are the volumetric liquid water, water vapor, and ice
contents (L3/L-3), respectively, ρl, ρv, and ρi are the densities of liquidwater, water vapor, and ice (M/L-3), respectively, t is time (T), z is thespatial coordinate positive upward (L), Klh (L/T) and KlT (L2 K− 1 T− 1) are

Fig. 2. The site location and the instrumental layout. Adapted from Zheng et al. (2021).
3

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

Fig. 3. Variations of (a) diurnal meteorological variables during the analyzed period and (b) annual precipitation and air temperatures from 1982 to 2021.
Table 1
Measured soil particle composition and soil hydraulic parameters of different layers.
Soil layer (cm)0–80
80–160
160–200

Soil particle composition

Soil hydraulic parameters

Sand (%)

Silt (%)

Clay (%)θr (cm3 cm− 3)

θs (cm3 cm− 3)

α (cm− 1)n (-)

Ks (cm d-1)

95.8
90.1
91.7

3.2
5.4
5.2

1.0
4.5
3.1

0.008
0.015
0.010

0.37
0.42
0.40

0.029
0.039
0.037

1.64
1.51
1.59

1368
811
1440

Fig. 4. Coexistence of ice, water, and vapor in unsaturated frozen soils and the driving forces for liquid and vapor flow.

the isothermal and thermal hydraulic conductivities for liquid water,
respectively, Kvh (L/T) and KvT (L2 K− 1 T− 1) are the isothermal andthermal hydraulic conductivities for water vapor, respectively, h is thepressure head (L), T is the temperature (K), Cp, Cw, and Cv are volumetricheat capacities of moist soil, liquid water, and water vapor (M/L T− 2
K− 1), respectively, Lf is the latent heat of freezing (L2/T− 2), L0 is the4

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

volumetric latent heat of vaporization of liquid water (M/L T− 2), λ is theapparent thermal conductivity (M L T− 3 K− 1), and ql and qv are liquidand vapor fluxes (L/T), respectively.
As shown in Fig. 4, liquid water flow in unsaturated frozen soils is
mainly induced by the pressure head, gravitational, and temperature
gradients (the first three terms on the right-hand side of Eq. (1)), while
the pressure head and temperature gradients drive vapor flow (thefollowing two terms) (Hu et al., 2023). Detailed descriptions of the hy­
draulic and thermal characteristics of the soils are given in Appendix A.

October 31, 2018, was used to evaluate the proposed model. Measured
and simulated temperatures and water contents at various depths areillustrated in Fig. 5. The coefficient of determination (R2), root meansquare error (RMSE), and Nash-Sutcliffe efficiency (NSE) (Nash andSutcliffe, 1970) were the three statistical indicators employed to assess
the model’s accuracy (Table 2).
The simulated soil water contents approximated the observed valuewell, with average values of R2 and RMSE of 0.85 and 0.01 cm3/cm3,
respectively. The fitting accuracy was better when the soil was frozen,
with relatively low soil water contents. The calculated errors were larger
in periods with higher soil water contents, mainly caused by complex
external factors, i.e., intense evaporation and rainfall. The calculated
NSE values were 0.77, 0.89, and 0.71 for the 10, 50, and 100 cm depths,
respectively, indicating that simulated values match well with measured
data in the upper soil layers. Since soil water contents at the 200 cm
depth were less affected by the freeze–thaw cycles and rainfall infiltra­tion processes, only ranging from 0.077 to 0.097 cm3/cm3, the calcu­
lated NSE value was − 2.2. This suggested that NSE is not a useful
criterion when variations in soil water contents are relatively small andshould be more suitable for evaluating dynamic rather than static var­iables, which was also reported by Stumpp et al. (2012). On the other
hand, the correspondence between simulated and measured tempera­ture was significantly better, with an average RMSE of only 0.46 ◦ C.
Overall, the simulation results proved that the proposed model is highly
accurate and applicable to the study area.

2.2.2. Development of the numerical model
Compared with the standard HYDRUS-1D model the research com­munity has widely used, the applications of the existing freezing moduleof HYDRUS-1D developed by Hansson et al. (2004) are minimal, eventhough it was developed many years ago. The major problem is related
to this module’s stability, which often becomes unstable when used forvarious field tests and applications. To overcome these numericalproblems, the available energy approach for the phase change was usedin our modified code to overcome the instability originating from thephase transformation between water and ice (Zheng et al., 2021).
Moreover, a flexible time-stepping method was employed to calculate
the coupled processes. A detailed calculation procedure can be referredto (Zheng et al., 2021).
A numerical model was developed in this study using the modified
HYDRUS-1D freezing code and the hydrogeological conditions andobservation data collected at the study site. The transport domain be­
tween the soil surface and a 200 cm depth was discretized using 2 cm
elements. The analyzed period ranged from November 1, 2017, toOctober 31, 2018, covering an entire hydrological year. The time step of
2 h was used to analyze both diurnal and seasonal variations of soilhydrothermal variables and water fluxes.
The upper boundary condition for water flow was defined using
meteorological variables (i.e., precipitation and evaporation). The lower
boundary was almost unaffected by groundwater due to its distancefrom the groundwater table (around 9 m at the study site), which was
specified as a free drainage condition. The calculated surface tempera­ture (estimated from the temperature measured at the shallow depth andsoil heat flux data, as introduced by Dahiya et al., 2007) and observed
temperature at the 200 cm depth were applied as the upper and lower
boundary conditions for heat transport, respectively. The initial condi­
tions were defined using measured temperatures and water contents atthe beginning of the simulation period, and data between two obser­
vation points were obtained by linear interpolation. In general, theinitial and boundary conditions can be described as follows:
⎧⃒
⃒
⃒
∂h
∂T
∂h
∂T ⃒
⎪
⎪
⎪ ⃒⃒ − Klh − Klh − KlT − Kvh − KvT ⃒⃒ ≤ Em z = 0
⎪
⎪
∂z
∂z
∂z
∂z
⎪
⎪
⎪
⎪
⎨
∂h(z, t)= 0 z = 200 cm
(3)
∂z
⎪
⎪
⎪
⎪
⎪
T(z, t) = T0 (t)z = 0 or z = 200 cm
⎪
⎪
⎪
⎪
⎩
θl (z, t) = θl,i (z), T(z, t) = Ti (z)t = 0

3. Results
3.1. Spatial-temporal variations of soil temperatures and water contents
To display both seasonal and diurnal changes, simulated soil tem­
peratures and unfrozen water contents at 0:00 and 12:00 of each day
during the analyzed period are shown in Fig. 6.
3.1.1. Soil temperature
There were significant temporal and spatial changes in soil temper­
ature (Fig. 6a). Affected by seasonal air temperature variations (Fig. 3),
soil temperature showed a similar trend, first decreasing and thenincreasing. The maximum and minimum surface temperatures occurredon June 26 (39.7 ◦ C) and February 2 (-13.3 ◦ C), respectively. The
location of the maximum soil profile temperature also showed seasonalchanges, at the 200 cm depth at the beginning of the simulation andgradually moving to the soil surface in summer. Since soil temperatures
above and below this point were lower, a divergent zero temperature
gradient plane existed in the soil profile. Driven by temperature gradi­
ents, soil liquid water and water vapor moved away from this point on
both sides.
On the other side, there were significant differences in temperature
variations at different depths. For example, the annual temperaturechange at the 10 cm depth ranged from − 10.1 ◦ C to 33.6 ◦ C, with anamplitude of 43.7 ◦ C. The annual temperature change significantly
decreased as the depth increased to 26.4 ◦ C and 17 ◦ C at depths of 100
and 200 cm, respectively. There was a significant diurnal variation in
the shallow 20 cm soil layer, especially in summer, with the average
daily amplitude reaching approximately 10 ◦ C (Fig. 6a). These consid­
erable temperature changes would result in variations in soil watercontent distributions, driving liquid water and water vapor towards
lower temperature locations.

where Em denotes the maximum flux density at the upper boundary (L/T).
In addition, the soil water retention and isothermal hydraulic con­ductivity are described using the van Genuchten (1980) model and thepore-size distribution model of Mualem (1976) in this numerical model.
The required hydraulic parameters were initially estimated using theneural network method with the measured soil physical properties andthen optimized according to the Levenberg-Marquardt algorithm(Šimůnek et al., 2008). The optimized hydraulic parameters for simu­
lation are listed in Table 1.

3.1.2. Soil water content
Fig. 6b shows that the soil water content in the 80–160 cm layer was
notably higher than at other depths. The formation of a high watercontent layer was closely related to soil physical properties of different
layers (Table 1). Compared with soil temperature variations, diurnalchanges in unfrozen water contents were relatively smaller and mainly
concentrated within the shallow 20 cm depth (e.g., increasing in the

2.2.3. Model validation
The dataset collected at the study site from November 1, 2017, to5

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

Fig. 5. Comparisons of observed and simulated (a) unfrozen water contents and (b) soil temperatures at depths of 10, 50, 100, and 200 cm during theanalyzed period.

To further explore the effects of the water–ice phase change on soil
moisture distributions, Fig. 7 depicts the diurnal variations of liquidwater content, ice content, total water content, and temperature at the 2
cm depth during the freezing-thawing period (from November 1, 2017,
to March 31, 2018). It can be seen that there is a close correlation be­
tween variations in soil water contents and temperatures when the soil isfrozen. Diurnal changes in liquid water and ice contents during theinitial freezing periods were significant (with an average diurnal
amplitude over 0.03 cm3/cm3). Meanwhile, the total water content
increased from 0.065 cm3/cm3 to 0.085 cm3/cm3 in February, whichshould be attributed to the upward soil water movement due to soil
temperature and matric potential gradients. Subsequently, the total soil
water content gradually decreased with increasing temperatures, andthis reduction was mainly caused by evaporation and infiltration ofmelted water.

Table 2
Calculated R2 (coefficient of determination), RMSE (root mean square error),and NSE (Nash-Sutcliffe efficiency) between observed and simulated unfrozen
water contents and soil temperatures.
Depth
(cm)

10
50
100
200

Unfrozenwater
content

Soil
temperature

R2

RMSE
(cm3/
cm3)

NSE

R2

RMSE
(◦ C)

NSE

0.83
0.93
0.75
0.89

0.013
0.009
0.012
0.008

0.77
0.89
0.71
− 2.23

0.98
0.99
0.98
–

0.41
0.44
0.53
–

0.99
0.99
0.98
–

Note: Statistical indicators were not calculated for the 200 cm depth because themeasured temperature at this point was used as the bottom boundary condition.

3.2. Soil water balance

daytime and decreasing at night, as affected by the diurnal cycles of thefreeze–thaw process). Rainfall events and seasonal freeze–thaw pro­
cesses mainly controlled the variation in unfrozen water content.
First, rainfall plays a critical role in soil water content distributions inthe soil profile, primarily determined by the rain intensity. The infil­
tration depth and duration were very different under various rainfall
scenarios. For instance, the maximum rainfall event occurred on August
10, with a total rainfall of 57.1 mm. At this time, the soil water content
reached 0.157 cm3/cm3 at the 10 cm depth, and the infiltration front
eventually reached the bottom of the soil profile (200 cm), increasing
the soil water content there as well. Second, the unfrozen soil watercontent decreased dramatically when the soil was frozen, and thechanges in soil water content were consistent with the temperaturechange (0 ◦ C isoline). The reduction in unfrozen water content closely
followed a decrease in temperature, and the minimum values were only
0.009, 0.011, and 0.041 cm3/cm3 at depths of 10, 50, and 100 cm. When
air temperature increased, the thawing process started from both the soil
surface and the bottom of the frozen zone, rapidly increasing unfrozen
water contents.

The changes in various water balance components provide datacharacterizing the soil water cycle processes (Hou et al., 2016; Brown
et al., 2021). Considering the effect of ice and vapor, the soil (liquid)water storage could be estimated as follows:
S+I +V = P− E− D

(4)where S is the liquid water storage (L), I is the ice storage (L), V is thevapor storage (L), P is precipitation (L), E is evaporation (L), and D isdrainage at the soil profile bottom (L). Fig. 8 depicts temporal changes invarious water balance components during the study period.
3.2.1. Liquid water storage
While the liquid water storage (LWS) fluctuated dramatically during
the simulation period (in the range from 10 to 30 cm), it changed
relatively little over the entire simulated period, decreasing by only 0.9
cm in the whole soil profile, suggesting that annual changes were in a
fairly balanced state (Fig. 8a).
Affected by the freezing process at the beginning of the simulation
6

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

Fig. 6. Spatial-temporal distributions of (a) soil temperature and (b) unfrozen water content from November 1, 2017, to December 31, 2018. The left and rightfigures depict data at 0:00 and 12:00, respectively.

(November through February), the LWS tended to decrease due to thewater–ice phase transformation process (while the ice storage I
increased) and reached its minimum (10.2 cm) on February 16. Whenthe air temperature increased above zero (in March), the LWS rose as theice melted (I decreased to zero). The LWS eventually declined by 2.3 cm
at the end of the melting process compared to the LWS before freezing.
This reduction over these five months was primarily due to melted water
evaporation. Subsequently, the LWS was relatively stable due to limited
precipitation from April to June. From July to September, the study site
was relatively humid, resulting in significant fluctuations in the LWS.
Due to the effects of large rainfall events and intense evaporation, LWS
increases and decreases were relatively rapid at this time. For example,
on August 12, the LWS reached its maximum value of 28.7 cm (almost
three times the minimum value) before declining to 21.5 cm on August28.

water content variations at a depth of 200 cm (Fig. 4), confirming anapparent effect of rainfall on drainage. With increased rainfall frequency
and intensity, the maximum drainage flux reached 0.16 cm/d, and thecumulative drainage was 9.1 cm, reflecting intense rainfall events in thisperiod. This rarely occurs in arid and semiarid areas and could be amajor source for local groundwater recharge.
3.3. Seasonal and diurnal soil water transfer
In seasonally frozen soils, the soil temperature and water potential(including matric potential and gravity) gradients are the main driving
forces for the soil liquid water and water vapor movement. The total soil
water flux can be divided into thermal liquid water (qlT), thermal water
vapor (qvT), isothermal liquid water (qlh), and isothermal water vapor(qvh) fluxes, as introduced by Eq. (1). Fig. 9 shows the seasonal anddiurnal variations of soil water fluxes during the analyzed period at
depths of 2, 10, 50, and 100 cm.

3.2.2. Evaporation
The soil surface (the model’s upper boundary) was in direct contact
with the atmosphere, affected by rainfall and evaporation (Fig. 8bcd).
The flux density at the soil surface was dominated by upward evapo­
ration and downward rainfall infiltration. Since evaporation was closely
related to air temperature and radiation, daily evaporation was minimal(cumulation evaporation was only 3.1 cm) during the freezing-thawing
period. Daily evaporation increased rapidly during non-freezing periods,
with significant diurnal fluctuations influenced by environmental fac­
tors. For example, on August 10, the maximum value reached 1.98 cm/
d, far below the average value of 0.16 cm/d. Lastly, the estimated cu­
mulative evaporation reached 57.7 cm, accounting for approximately
87 % of rainfall, reflecting relatively strong evaporation intensity in the
study area.

3.3.1. Thermally driven flow
Driven by a temperature gradient, the qlT and qvT fluxes were always
directed in the same direction, usually upward at night and downward in
the daytime. The liquid water and water vapor fluxes significantly
decreased with increasing depth, with values at 2 cm usually two ordersof magnitude greater than at 100 cm. Meanwhile, while diurnal fluc­tuations of the thermally driven fluxes were apparent in shallow depths,
they became negligible in deeper depths due to decreasing soil tem­
perature gradients and hydraulic conductivities. For example, themaximum upward and downward fluxes in the 2 cm depth reached 7.6e04 and − 0.004 cm/h, while they were only 3.6e-05 and 8.9e-05 cm/h in
the 100 cm depth, respectively. During the analyzed period, their vari­ations were significantly influenced by the freeze–thaw and rainfall
infiltration processes.
When the soil was frozen, liquid water and water vapor flows wererestricted by a reduction in the hydraulic conductivity caused by ice insoil pores. This occurred especially for qlT, which was reduced byapproximately five orders of magnitude compared with the non-freezing
period (i.e., the minimum qlT at the 2 cm depth being 3.3e-10 cm/h).
Driven by a temperature gradient, liquid water and water vapor flowed
mainly toward the surface soil layers during freezing. During thawing,

3.2.3. Drainage
Like evaporation variations, drainage fluxes were relatively small
during the freeze–thaw period (Fig. 8ef). Since the vadose zone was in a
relative water deficit state due to limited precipitation from November
to June, the infiltrated water could hardly migrate down to the 200 cm
depth. As a result, drainage fluxes were constantly smaller than 0.01
cm/d. After higher rainfall events in early July, the bottom flux signif­
icantly increased to 0.07 cm/d. This increase was consistent with the soil
7

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

Fig. 7. Diurnal variations of (a) liquid water content (θl), (b) ice content (θi), (c) total water content (θt), and (d) temperature at the 2 cm depth during the freezingthawing period.

qlT increased rapidly, driving liquid water to flow toward the deep soil
layers.
During the non-freezing period, the variations in qlT and qvT exhibi­
ted completely different response patterns when rainfall occurred. The
qlT values increased instantly due to the increasing hydraulic conduc­tivity for liquid water, and the maximum downward value at the 2 cm
depth reached − 0.045 cm/h on August 11 after continuous rainfall
events, which was significantly higher than the corresponding flux
before rainfall. In comparison, qvT decreased due to reduced tempera­ture gradient and vapor diffusivity. On the other hand, qvT was higher
than qlT during relatively dry periods. For instance, after only 5 mm ofrainfall during the prior 30 days, the maximum value of qvT at the 2 cm
depth reached − 0.004 cm/h on June 22. In the meantime, qlT was only
− 2.8e-04 cm/h, which was far less than qvT. This confirms different
responses of qlT and qvT to rainfall events and indicates that it is crucial
and necessary to consider the vapor flow component of the soil water
movement in dry conditions.

Infiltration water penetrated the soil profile, thus serving as the main
source for groundwater recharge. When rainfall stopped, the upward
flux density generally increased due to evaporation. For instance, theupward qlh reached its maximum (0.026 cm/h) on July 18 after a total
rainfall of 69.8 mm. Due to evaporation, soil water potential gradientstended to increase in the period without rainfall, driving liquid water toflow upward toward the soil surface. The blocking effect severely
affected liquid flow when the soil was frozen, with the flux density
sharply declining (generally 3–7 orders of magnitude) in the frozen
layer.
The qvh flux density in shallow depths was substantially larger than inother depths. The maximum value in the 2 cm depth was 3.4e-05 cm/h,
while merely 2.6e-06 and 8.6e-07 cm/h in the 10 and 50 cm depths,
respectively. Since qvh was significantly smaller than the other flux
components (at least two orders of magnitude) and contributed little tothe total soil water fluxes during the entire period, the effect of qvh could
be neglected in this study.

3.3.2. Isothermally driven flow
Fig. 9 shows that variations in isothermally driven liquid water fluxes
are mainly related to precipitation and evaporation. The downward qlh
flux increased significantly after rainfalls, and the maximum flux density
at the 2 cm depth reached − 0.79 cm/h on July 10 (the same day as qlT).

3.4. Changing characteristics of the vapor flux
Based on the seasonal and diurnal characteristics of vapor flux var­iations, typical freezing (from November 17 to 21, 2017) and drought(from June 14 to 18, 2018) periods were selected to explore the detailed
8

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

Fig. 8. Simulated soil water balance components during the simulation period: (a) liquid water storage, (b) water flux at the soil surface, (c) evaporation, (d)cumulative evaporation, (e) water flux at the bottom boundary, and (f) cumulative water flux at the bottom boundary.

effects of vapor flow. During the selected freezing period, the soil in thesurface layer usually froze at night and thawed in the daytime, and thefreezing depth gradually increased (from 5 to 10 cm below the surface).
In the drought period, there was only 5 mm of rainfall in the prior 30
days, and the soil was relatively dry. Since the influence of theisothermally driven vapor flux could be neglected, only the thermal
vapor flux (qvT) is selected for specific analysis in this section.

soil water flow research when large temperature gradients are present.
3.4.2. The ratio of the water vapor and total water flux
To quantitatively assess the significance of vapor flow, the ratios of
qvT to the total soil water flux during the selected freezing and drought
periods are shown in Fig. 11. Since the flux density magnitudes signif­
icantly decreased with depth, only the variations at 2 and 10 cm depths
are selected.
During the freezing period, the changes in the flux ratios were sig­
nificant and exhibited typical diurnal patterns, which were affected by
the ice formation at night and its melting in the daytime. For example, atthe 2 cm depth, the average flux ratio at night exceeded 90 %, while it
was less than 10 % in the daytime, suggesting a significant influence of
ice on soil water flux variations. Note the flux ratio in the daytime of thefreezing period was smaller than in the non-freezing period, with theformer influenced by the daily melting process. Similar to the wetting
process, the liquid water flux significantly increased when the soil
melted, reducing the vapor flux ratio. As the freezing front moved below
the 10 cm depth, the flux ratio at this depth increased rapidly from 5.1 %on the first day (November 17) to 29.2 % on the fifth day (November
21).
During the drought period, the ratio of the vapor flux to the total
water flux exhibited a significant diurnal trend, which was opposite thatduring the freezing period. The ratios of qvt to the total water flux
increased in the daytime and decreased at night, showing a significant
correlation with diurnal temperature and vapor flux changes. For the 2

3.4.1. Influencing factors
As indicated by Eq. (1), the vapor flux is calculated using the tem­
perature gradient (Tg) and the hydraulic conductivity for vapor flow
(KvT). To explore the correlation between these two variables and qvT,
the qvT variations at 2, 10, 50, and 100 cm depths during the selected
typical freezing and drought periods are depicted in Fig. 10, which also
lists the estimated correlation coefficient (R).
One can see a high correlation between Tg and qvT in both freezing
and drought periods, with an average R of 0.94 for all selected depths.
The Tg range was larger in the surface soil than in a depth of 100 cm. A
temperature gradient reduction resulted in a decrease in the vapor flux.
In contrast, qvT decreased while KvT increased at shallow depths, sug­gesting a negative correlation trend. In deeper layers, there is no sig­
nificant correlation between these two variables. Additionally, theregression analysis results also indicated the Tg’s dominant effect on thevapor flux compared to KvT. The results indicate that temperature gra­dients not only determine the direction of vapor flow, but also dominate
the flux magnitude, suggesting that vapor flow should not be ignored in
9

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

Fig. 9. Variations of (a) thermal liquid, (b) thermal vapor, (c) isothermal liquid, and (d) isothermal vapor fluxes at depths of 2, 10, 50, and 100 cm from November 1,
2017, to October 31, 2018.

and 10 cm depths, the maximum ratio reached 47.1 % and 36.9 % atnoon, respectively, while average values were 16.8 % and 12.8 %. Thisagain clearly shows the importance of vapor flow for soil water flow,
especially for relatively dry conditions (i.e., without rainfall for a longtime).

relatively short time, with the liquid water flux significantly increased(similar to the wetting process). (iii) The relatively stable stage (from
April to June): Less affected by external environmental factors, thechanges in soil hydrothermal dynamics and water fluxes were limited,
and the vadose zone gradually became drier. (iv) The intense fluctuation
stage (from July to the freezing period): Frequent rainfall and intense
evaporation occurred during this stage. Liquid water fluxes significantly
increased after rainfall, whereas vapor flow largely contributed to soil
water flow under continuous drought conditions. As revealed by thesimulation results, the freeze–thaw cycle and the rainfall infiltration
process play a significant role in the soil water cycle.
In addition to affecting soil physical characteristics, soil hydraulicand thermal properties were also significantly affected when the soil wasfrozen (Kurylyk and Watanabe, 2013). The appearance of ice caused thenumber of connected pores in the vadose zone to be reduced. The
blocking effect due to the presence of ice on liquid flow has been
demonstrated through experimental and simulation studies (Zhao et al.,2013; Yu et al., 2018). As a result, vapor flow dominated soil watertransfer, and total soil water fluxes were much lower (including evap­oration and deep percolation) compared with unfrozen conditions

4. Discussion
4.1. Mechanisms of the soil water cycle in seasonally frozen soil
Driven by water potential and temperature gradients, soil liquidwater and water vapor fluxes showed different variations in the entire
hydrological year (Fig. 9). Following the soil water storage variations,the soil water cycle can be divided into the following four stages. (i) The
freezing-induced water accumulation stage (from November to late
February): The soil liquid water and water vapor movement wasimpeded due to the ice formation and mainly directed toward the frozen
layer, increasing the total water content there. (ii) The melting water
infiltration stage (from late February to March): The melting process
occurred in the frozen layer from both the top and bottom sides in a10

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

Fig. 10. Influence of the temperature gradient (Tg) (a, b) and thermal hydraulic conductivity (KvT) (c, d) on the thermal vapor flux (qvT) at different depths (2, 10, 50,
and 100 cm from left to right): (a) qvT(Tg) during the freezing period, (b) qvT(Tg) during the drying period, (c) qvT(KvT) during the freezing period, and (d) qvT(KvT)
during the drying period. The correlation coefficient was not calculated for qvT(KvT) due to its low values.

(Fig. 8). The key issue for clarifying the mechanisms of the soil hydro­
logical cycle during the freeze–thaw period is to explore the specificeffects and influencing patterns of ice presence. Generally, the imped­
ance factor method, instead of the semi-theoretical (based on capillary
and adsorption theory) and empirical methods, is widely employed toestimate the hydraulic conductivity for unsaturated soils since it isconsidered to be relatively reliable (Kurylyk and Watanabe, 2013).
Nevertheless, some researchers pointed out that this approach is not
physically based and provides only a constant impedance value forvaried soil matric potentials (Yu et al., 2018). Proposing an alternative
method for representing and incorporating the blocking effect into themodel would be crucial for further clarifying the mechanisms of soil
water movement during freeze–thaw periods.
On the other hand, rainfall infiltration processes significantly influ­
enced soil liquid water and water vapor flow. Fig. 3 shows thatapproximately 80 % of rainfall occurred from June to September, mostly
in the form of heavy rains or even rainstorms. Given intense radiation inarid and semiarid areas, rainfall water was largely returned to the at­mosphere through surface evaporation (Hou et al., 2016). Moreover,although soil water contents in deeper soil layers showed only small
fluctuations, the cumulative drainage during the analyzed period wassignificant due to the large amount of infiltrated water from the soil
surface.
It should be noted that total precipitation during this hydrologicalyear reached 660.8 mm, far exceeding the average annual precipitation
at the study site. Meteorological data from a nearby weather station
showed that average precipitation was 419.3 mm during the past forty

years, and the maximum and minimum rainfall amount reached 724.9
mm (in 2016) and 248.7 mm (in 2005), respectively (Fig. 3). The
modified HYDRUS-1D freezing module was applied to simulate the soil
water movement during 1982–2021 to evaluate the long-term status of
the soil water balance in the study area. Comparatively, the soil water
cycle has changed only slightly during the freeze–thaw periods in
different years. On the other hand, the soil water storage changed
significantly during different years. Simulated cumulative drainage
fluxes in the years with total rainfall of 724.9, 419.3, and 248.7 mm
were 9.1, 3.2, and 2.2 cm, respectively, indicating that deep drainage
was highly correlated with rainfall intensity.
In this work, the free drainage boundary condition was utilized todescribe the bottom condition, assuming the unit vertical hydraulicgradient condition was satisfied (with a soil matric potential gradientclose to zero) and providing a continuous bottom flux to serve as
groundwater recharge (Hou et al., 2016). The simulation results confirm
the hypothesis that precipitation can recharge local groundwater in the
Mu Us Sandy Land despite a thick vadose zone.
4.2. Implications and effects of vapor flux
Since the liquid–vapor phase transformations occur at all times in thevadose zone, the effect of vapor flow should be considered in most ap­
plications, especially in arid and semiarid regions (Zhou and Zhao,2021). The most direct effect of vapor flow is related to changes in soil
moisture distributions in the vadose zone, such as those discussed in thismanuscript. During the freezing period, vapor flows upward toward the11

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

Fig. 11. Variations in the ratio of the isothermal vapor flux and the total water flux during the selected (a) freezing and (b) dry periods at depths of 2 and 10 cm.

freezing zone, where it condenses, providing an additional water source
for ice formation and increasing the total soil water content in the frozen
layer. Yu et al. (2018) and Zheng et al. (2021) reported this phenomenon
through numerical simulations. Their results suggested that this fraction
could be up to 10 %, which was also confirmed in this study.
On the other hand, with high air-filled soil porosity and large tem­
perature gradients during dry conditions, the thermally driven vapor

flux is on the same order of magnitude as the liquid water flux, thus
affecting the top layer’s soil moisture. Moreover, driven by a tempera­ture gradient, vapor flowed downward toward deep soil horizons during
most of the analyzed period. Although the flux density was limited, thecumulative flux caused by continuous migration was considerable andshould not be ignored since it may become an important source ofgroundwater recharge (Xiang et al., 2012).

Fig. 12. Measured diurnal sap flow of Salix psammophila during freezing (from November to March) and non-freezing (from April to October) periods and its
correlation with different environmental factors. Adapted from Chen et al. (2022).
12

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

After clarifying the vapor transfer mechanisms and their influence on
soil hydrothermal dynamics, it is possible to explore the implications
and effects of vapor flow on actual activities. With abundant desertifi­
cation and sparse vegetation, the Mu Us Sandy Land has long been
considered a fragile ecological area. Windbreak and sand fixation
vegetation, such as Salix psammophila and Artemisia aragana, is widelydistributed at the Mu Us Sandy Land (Huang et al., 2016; Hayat et al.,2020). It is necessary to investigate the sources of vegetation water
consumption since the balance between vegetation growth and water
consumption needs to be clarified for successful local vegetation resto­
ration and environmental protection. Through long-term monitoring ofvegetation physiological characteristics, Chen et al. (2022) indicatedthat the sap flow of Salix psammophila during the freezing period still
represents a considerable proportion of the total annual sap flow
(although smaller than during the non-freezing period, as shown in
Fig. 12). Since the liquid flow is severely impeded when the soil isfrozen, condensation of water vapor could be a dominant water sourcefor vegetation growth during the dormant (freezing) season. Meanwhile,
Salix psammophila exhibits significant differences in daily water con­sumption patterns during different periods, which could also be related
to diurnal variations in vapor flow. By coupling the process-based root
water uptake model (e.g., Couvreur et al., 2012) with the HYDRUS-1D
freezing module, the hydrological cycle of the “soil-plant-atmosphere”
system of the seasonally frozen regions could be clarified in future
research, and the specific contribution of vapor flow to root water up­
take and vegetation growth during both the growing and dormant sea­sons could be evaluated.
In addition to influencing agricultural and ecological applications,vapor flow also plays a significant role in numerous engineering appli­
cations. For instance, frost heaves are commonly observed on theembankment of highways and railways in cold and arid areas. This is
primarily caused by the impermeable cover layer (Wu et al., 2016; Heet al., 2018), with vapor condensed and stored in the topsoil layer as ice.
However, previous research has mainly focused on the frost heaveprocess and waterproofing projects related to liquid water migration.
The impact of vapor migration has been frequently overlooked, leading
to engineering problems (Zhang et al., 2021). Thus, by conducting field
monitoring and experiments, as well as numerical simulations, designsof engineering constructions in cold areas must clarify the mechanismsof vapor transfer below different impeding layers and propose specificprevention measures.

high accuracy, and can be applied to the study region.
As indicated by the simulation results, the variations in soil tem­
perature, water content, water storage, evaporation, and drainage dis­played typical seasonal patterns. The soil water cycle in seasonallyfrozen regions can be divided into four stages: freezing-induced water
accumulation stage, melting water infiltration stage, relatively stablestage, and intense fluctuation stage. Driven by soil water potential
gradients, liquid water flow played a major role in affecting total soil
water flow, and it generally became negligible when the soil was frozen.
In contrast, water vapor fluxes induced by temperature gradients
showed significant diurnal changes and were relatively larger thanliquid water fluxes under freezing and drought conditions, representing
up to over 90 % and, on average, approximately 20 % of the overallwater flux.
In general, these results clarified the detailed mechanisms of soilhydrological processes in cold regions. Since the implications of liquidwater and water vapor transfers according to the current monitoring of
the “soil-plant-atmosphere” system are significant, the effects and spe­cific contributions of soil water flow to vegetation and ecosystems need
to be quantitatively determined, which could be crucial for local vege­
tation restoration and environmental protection of fragile ecologicalareas.

5. Conclusions

Acknowledgments

Following the theory of the coupled transport of water, vapor, andenergy, the modified HYDRUS-1D freezing code was used to develop a
numerical model for the Mu Us Sandy Land using the in-situ monitoringdata to investigate the mechanisms of soil water transfer. The simulated
results indicated that the proposed model is numerically stable, provides

This research was supported by the National Natural Science Foun­dation of China (No. 42202279, 42372288), the Fundamental ResearchFunds for the Central Universities, CHD (No. 300102292105,
300102292904), and China Postdoctoral Science Foundation (No.
2022M720536).

CRediT authorship contribution statement
Ce Zheng: Conceptualization, Methodology, Writing – original draft.
Yunfei Chen: Investigation, Data curation, Writing – review & editing.
Wande Gao: Investigation, Data curation, Writing – review & editing.
Xinyang Liang: Investigation, Methodology. Jiří Šimůnek: Writing –review & editing, Supervision. Xiuhua Liu: Methodology, Supervision,
Writing – review & editing, Funding acquisition.
Declaration of competing interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Data availability
Data will be made available on request.

Appendix A:. Evaluation of soil properties
A.1 Soil hydraulic properties
Following the pore-size distribution model of Mualem (1976), the hydraulic conductivity under isothermal conditions could be described using thevan Genuchten model (1980):
θ s − θrm
[1 + |αh|n ]

(A1)

[
(
)m ]2
Klh = Ks Sel 1 − 1 − Se1/m(A2)

θ l = θr +where θr and θs are the residual and saturated water contents (L3/L-3), respectively, α (L-1), n (-), and m (-) are empirical coefficients, Ks is the saturated
hydraulic conductivity (L/T), and Se is the effective liquid saturation (-), that can be defined as follows:13

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401(θl − θr )(θs − θr )

Se =

(A3)

With the ice formation, water flow in the soil pores is subject to significantly enhanced resistance (a blocking effect from the ice), which graduallystrengthens with an increase in the ice fraction. The hydraulic conductivity under isothermal conditions for water flow in the freezing soil (Kflh, L/T) is
described below (Lundin, 1990):(A4)Kflh = 10− ψ Q ⋅Klh

where ψ is the water flow impedance factor owing to ice formation (-) and Q is the solid ice phase’s weight ratio to all three water (-) phases, which canbe calculated as follows:
Q=

ρi θi
ρi θ i + ρl θ l(A5)

The thermal hydraulic conductivity can be determined using the following equation:()
1 dg
KlT = Klh hGwT
g0 dT

(A6)

where Gwt denotes the gain factor (-), which is used for correcting the temperature dependence in a soil water retention curve, γ0 is the surface tensionof water at 25 ◦ C (i.e., 71.89 g s− 2), and γ is the surface tension of water in soil (M/T− 2), which can be calculated as follows:(A7)

γ = 75.6 − 0.1425T − 2.38 × 10− 4 T 2

The isothermal hydraulic conductivity and thermal hydraulic conductivity for vapor flow are defined as follows, respectively (Saito et al., 2006):
D

Kvh =

KvT =

ρ Hr
ρw sv RT

Mg

(A8)

dρsv
dT

(A9)

D

ρw

ηH r

where M is the water molecular weight (0.018 kg mol− 1), g is the gravitational acceleration (9.8 m s− 2), R is the water vapor universal gas constant(8.315 J mol− 1 K− 1). Besides, η is the enhancement factor (-), and D is the water vapor diffusivity in the soil (L2/T), defined as:(A10)

D = Da θ a τ-32

where θa is the air-filled porosity (L3/L ), and Da is the water vapor diffusivity in the air (L /T) expressed as:()2T
Da = 2.12×10 - 5
273.15(A11)

and τ is the tortuosity factor (–) as follows:

τ=
as:

θ7/3
a
θ2s

(A12)

The vapor density can be expressed as aproduct of relative humidity and saturated vapor density. The saturated vapor density (ρsv, M/L-3) is written(
exp 31.3716 −)
6014.79
−
T

ρsv =

T

7.92495 × 10− 3 T(A13)

× 10− 3and the relative humidity (Hr) could be expressed with a thermodynamic correlation between water and vapor in soils (Philip and de Vries, 1957):()
Mgh
Hr = exp(A14)
RT
Under the effect of the liquid islands as well as the increased gradients of temperature, an enhancement factor for the increasing thermal vapor flux
is described below (Cass et al., 1984):{ [() ]4 }
θ
2.6 θ
(A15)
η = 9.5 + 3 − 8.5exp −1 + √̅̅̅̅
θs
fc θs
where fc is the clay mass fraction in soil (–).
A.2 Soil freezing curve
Based on the theory of thermodynamic equilibrium and the hypothesis of ice pressure being zero for frozen soil, the matric potential of soil can be
estimated for subfreezing temperature by using the generalized Clapeyron equation described below (Hansson et al., 2004; Yu et al., 2018):

14

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401

Lf T
ln
g Tf

(A16)

h=

where Tf is the water freezing point (K), which varies depending on the soil texture and the initial water content. Integration of Eq. (A16) with Eq. (A1)
leads to the correlation between the content of liquid water and the temperature (i.e., SFCC - soil freezing characteristic curve) shown below.
θ − θr⃒ (s)⃒ ]θ l = θr + [⃒ Lf T ⃒n m
1 + ⃒⃒α g ln Tf ⃒⃒(A17)

A.3 Soil thermal properties
For wet soil, its volumetric heat capacity comprises the individual corresponding values of all phases, including the solid (Cn), liquid (Cw), vapor(Cv), and ice (Ci) phases:(A18)

Cp = Cn θ n + Cw θ l + Cv θ v + Ci θ iwhere θn denotes the solid phase volumetric fraction (L3/L-3) that can be estimated as follows:(A19)θ n = 1 − θs

By combining the apparent thermal conductivity of soil and the porous medium thermal conductivity with the relative thermal dispersivity, thefollowing equation for the unfrozen condition could be obtained (Šimůnek et al., 2008):(A20)

λ = λ0 (θl ) + βl Cw |qw |where βl is the longitudinal thermal dispersivity (L) and λ0 is the porous medium thermal conductivity (M L T− 3 K− 1), which can be described as(Chung and Horton, 1987):(A21)

λ0 (θl ) = b1 + b2 θl + b3 θ0.5lwhere b1, b2, and b3 represent the empirical coefficients that depend on the soil texture (M L T− 3 K− 1). Because the ice could affect the heat conduction
as well, the method proposed by Hansson et al. (2004) is adopted for determination of the frozen soil thermal conductivity, which is written as:(A22)

λ0 (θl ) = b1 + b2 (θl + Fθi ) + b3 (θl + Fθi )0.5
where F denotes the discrepancy in the thermal conductivity between liquid water and the soils (-).
Additionally, in Eq. (2), the volumetric latent heat of liquid water vaporization can be calculated as:(A23)

L0 = ρw Lw
where Lw denotes the latent heat of liquid water vaporization (L2/T− 2), which can be estimated as follows:(A24)

Lw = 2.501 × 106 − 2369.2T

References

Bittelli, M., Ventura, F., Campbell, G.S., Snyder, R.L., Gallegati, F., Pisa, P.R., 2008.
Coupling of heat, water vapor, and liquid water fluxes to compute evaporation inbare soils. J. Hydrol. 362, 191–205. https://doi.org/10.1016/j.jhydrol.2008.08.014.
Brown, S., Wagner-Riddle, C., Debruyn, Z., Jordan, S., Berg, A., Ambadan, J.T.,
Congreves, K.A., Machado, P.V.F., 2021. Assessing variability of soil water balancecomponents measured at a new lysimeter facility dedicated to the study of soil
ecosystem services. J. Hydrol. 603, 127037 https://doi.org/10.1016/j.
jhydrol.2021.127037.
Cahill, A.T., Parlange, M.B., 1998. On water vapor transport in field soils. Water Resour.
Res. 34 (4), 731–739. https://doi.org/10.1029/97WR03756.
Cass, A., Campbell, S., Jones, T.L., 1984. Enhancement of thermal water vapor diffusion
in soil. Soil Sci. Soc. Am. J. 48 (1), 25–32. https://doi.org/10.2136/
sssaj1984.03615995004800010005x.
Chadburn, S., Burke, E., Cox, P., Friedlingstein, P., Hugelius, G., Westermann, S., 2017.
An observation-based constraint on permafrost loss as a function of global warming.
Nature Clim. Change 7, 340–344. https://doi.org/10.1038/nclimate3262.
Chen, Y., He, J., He, Y., Gao, W., Zheng, C., Liu, X., 2022. Seasonal hydrological traits in
Salix psammophila and its responses to soil moisture and meteorological factors indesert areas. Ecol. Indic. 136, 108626 https://doi.org/10.1016/j.
ecolind.2022.108626.
Cheng, Q., Sun, Y., Jones, S.B., Vasilyev, V.I., Popov, V.V., Wang, G., Zheng, L., 2014. In
situ measured and simulated seasonal freeze-thaw cycle: A 2-year comparative study
between layered and homogeneous field soil profiles. J. Hydrol. 519, 1466–1473.
https://doi.org/10.1016/j.jhydrol.2014.09.023.

Banwart, S.A., Bernasconi, S., Bloem, J., Blum, W., Brandao, M., Brantley, S.L.,
Chabaux, F., Duffy, C., Kram, P., Lria, G., Lundin, L., Nikolaidis, N.P., Novak, M.,
Panagos, P., Ragnarsdottir, K.V., Reynolds, B., Rousseva, S., de Ruiter, P.C., vanGaans, P., van Riemsdijk, W., White, T., Zhang, B., 2011. Soil processes andfunctions in critical zone observatories: hypotheses and experimental design. Vadose
Zone J. 10 (3), 974–987. https://doi.org/10.2136/vzj2010.0136.
Bao, H., Koike, T., Yang, K., Wang, L., Shrestha, M., Lawford, P., 2016. Development ofan enthalpy-based frozen soil model and its validation in a cold region in China.
J. Geophys. Res. Atmos. 121 (10), 5259–5280. https://doi.org/10.1002/
2015JD024451.
Biskaborn, B.K., Smith, S.L., Noetzli, J., Matthes, H., Vieira, G., Streletskiy, D.A.,
Schoeneich, P., Romanovsky, V.E., Lewkowicz, A.G., Abramov, A., Allard, M.,
Boike, J., Cable, W.L., Christiansen, H.H., Delaloye, R., Diekmann, B., Drozdov, D.,
Etzelmüller, B., Grosse, G., Guglielmin, M., Ingeman-Nielsen, T., Isaksen, K.,
Ishikawa, M., Johansson, M., Johannsson, H., Joo, A., Kaverin, D., Kholodov, A.,
Konstantinov, P., Kröger, T., Lambiel, C., Lanckman, J.P., Luo, D., Malkova, G.,
Meiklejohn, I., Moskalenko, N., Oliva, M., Phillips, M., Ramos, M., Sannel, A.B.K.,
Sergeev, D., Seybold, C., Skryabin, P., Vasiliev, A., Wu, Q., Yoshikawa, K.,
Zheleznyak, M., Lantuit, H., 2019. Permafrost is warming at a global scale. Nat.
Commun. 10, 264. https://doi.org/10.1038/s41467-018-08240-4.

15

C. Zheng et al.

Journal of Hydrology 627 (2023) 130401
Philip, J.R., de Vries, D.A., 1957. Moisture movement in porous materials undertemperature gradient. Trans. Am. Geophys. Union 38 (2), 222–232. https://doi.org/
10.1029/TR038i002p00222.
Saito, H., Šimůnek, J., Mohanty, B., 2006. Numerical analyses of coupled water, vapor,and heat transport in the vadose zone. Vadose Zone J. 5 (2), 784–800. https://doi.
org/10.2136/vzj2006.0007.
Scanlon, B.R., Keese, K., Reedy, R.C., Šimůnek, J., Andraski, B.J., 2003. Variations inflow and transport in thick desert vadose zones in response to paleoclimatic forcing(0–90 kyr): Field measurements, modeling, and uncertainties. Water Resour. Res. 39(7), 1179–1196. https://doi.org/10.1029/2002wr001604.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.T., 2008. The HYDRUS-1D
Software Package for Simulating the One-Dimensional Movement of Water, Heat,and Multiple Solutes in Variably Saturated Media, Version 4.0, HYDRUS Software
Series 3. Department of Environmental Sciences. University of California Riverside,
Riverside, California, USA.
Stumpp, C., Stichler, W., Kandolf, M., Šimůnek, J., 2012. Effects of land cover andfertilization method on water flow and solute transport in five lysimeters: A longterm study using stable water isotopes. Vadose Zone J. 11 (1), 1–14. https://doi.org/
10.2136/vzj2011.0075.
Stuurop, J.C., van der Zee, S.E.A.T.M., French, H.K., 2022. The influence of soil texture
and environmental conditions on frozen soil infiltration: A numerical investigation.
Cold Reg. Sci. Technol. 194, 103456 https://doi.org/10.1016/j.
coldregions.2021.103456.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. https://doi.
org/10.2136/sssaj1980.03615995004400050002x.
Wang, X., 2015. Vapor flow resistance of dry soil layer to soil water evaporation in arid
environment: An overview. Water 7 (12), 4552–4574. https://doi.org/10.3390/
w7084552.
Wang, L., Zhou, J., Qi, J., Sun, L., Yang, K., Tian, L., Lin, Y., Liu, W., Shrestha, M.,
Xue, Y., Koike, T., Ma, Y., Li, X., Chen, Y., Chen, D., Piao, S., Lu, H., 2017.
Development of a land surface model with coupled snow and frozen soil physics.
Water Resour. Res. 53 (6), 5085–5103. https://doi.org/10.1002/2017WR020451.
Wu, Q., Zhang, Z., Gao, S., Ma, W., 2016. Thermal impacts of engineering activities andvegetation layer on permafrost in different alpine ecosystems of the Qinghai-Tibet
Plateau, China. Cryosphere 10 (4), 1695–1706. https://doi.org/10.5194/tc-101695-2016.
Xiang, L., Yu, Z., Chen, L., Mon, J., Lu, H., 2012. Evaluating coupled water, vapor, andheat flows and their influence on moisture dynamics in arid regions. J. Hydrol. Eng.
17 (4), 565–577. https://doi.org/10.1061/(asce)he.1943-5584.0000474.
Yu, L., Zeng, Y., Wen, J., Su, Z., 2018. Liquid-vapor-air flow in the frozen soil.
J. Geophys. Res. 123 (15), 7393–7415. https://doi.org/10.1029/2018JD028502.
Yu, L., Zeng, Y., Su, Z., 2020. Understanding the mass, momentum, and energy transfer
in the frozen soil with three levels of model complexities. Hydrol. Earth Syst. Sci. 24(10), 4813–4830. https://doi.org/10.5194/hess-24-4813-2020.
Zhang, Y., Ma, W., Wang, T., Cheng, B., Wen, A., 2019. Characteristics of the liquid andvapor migration of coarse-grained soil in an open-system under constanttemperature freezing. Cold Reg. Sci. Technol. 165, 102793 https://doi.org/10.1016/
j.coldregions.2019.102793.
Zhang, X., Zhai, E., Wu, Y., Sun, D., Lu, Y., 2021. Theoretical and numerical analyses onhydro-thermal-salt-mechanical interaction of unsaturated salinized soil subjected totypical unidirectional freezing process. Int. J. Geomech 21 (7), 04021104. https://doi.org/10.1061/(ASCE)GM.1943-5622.00020.
Zhang, X., Shu, C., Fujii, M., Wu, Y., Sun, D., Ye, P., Bao, Y., 2023a. Numerical andexperimental study on water-heat-salt transport patterns in shallow bare soil withvarying salt contents under evaporative conditions: A comparative investigation.
J. Hydrol. 621, 129564 https://doi.org/10.1016/j.jhydrol.2023.129564.
Zhang, X., Shu, C., Wu, Y., Ye, P., Du, D., 2023b. Advances of coupled water-heat-salttheory and test techniques for soils in cold and arid regions: A review. Geoderma
432, 116378. https://doi.org/10.1016/j.geoderma.2023.116378.
Zhao, Y., Nishimura, T., Hill, R., Miyazaki, T., 2013. Determining hydraulic conductivity
for air-filled porosity in an unsaturated frozen soil by the multistep outflow method.
Vadose Zone J. 12 (1), 1–10. https://doi.org/10.2136/vzj2012.0061.
Zhao, M., Wang, W., Wang, Z., Chen, L., Ma, Z., Wang, Q., 2020. Water use of Salix in thevariably unsaturated zone of a semiarid desert region based on in-situ observation.
J. Hydrol. 591, 125579 https://doi.org/10.1016/j.jhydrol.2020.125579.
Zheng, C., Šimůnek, J., Zhao, Y., Lu, Y., Liu, X., Shi, C., Li, H., Yu, L., Zeng, Y., Su, Z.,
2021. Development of the Hydrus-1D Freezing module and its application insimulating the coupled movement of water, vapor, and heat. J. Hydrol. 598, 126250
https://doi.org/10.1016/j.jhydrol.2021.126250.
Zhou, H., Zhao, W., 2021. Evolution of soil-water states in the vadose zone of a desert soil
after an extreme rainfall event and its impact on the ecosystem. Hydrogeol. J. 29,
2127–2147. https://doi.org/10.1007/s10040-021-02372-6.

Chung, S.O., Horton, R., 1987. Soil heat and water flow with a partial surface mulch.
Water Resour. Res. 23 (12), 2175–2186. https://doi.org/10.1029/
WR023i012p02175.
Couvreur, V., Vanderborght, J., Javaux, M., 2012. A simple three-dimensional
macroscopic root water uptake model based on the hydraulic architecture approach.
Hydrol. Earth Syst. Sci. 16, 2957–2971. https://doi.org/10.5194/hessd-9-49432012.
Dahiya, R., Ingwersen, J., Streck, T., 2007. The effect of mulching and tillage on thewater and temperature regimes of a loess soil: Experimental findings and modeling.
Soil till. Res. 96, 52–63. https://doi.org/10.1016/j.still.2007.02.004.
Du, C., Yu, J., Wang, P., Zhang, Y., 2017. Analysing the mechanisms of soil water andvapour transport in the desert vadose zone of the extremely arid region of northern
China. J. Hydrol. 558, 592–606. https://doi.org/10.1016/j.jhydrol.2017.09.054.
Gao, W., Zheng, C., Liu, X., Lu, Y., Chen, Y., Wei, Y., Ma, Y., 2022. NDVI-based
vegetation dynamics and their responses to climate change and human activities
from 1982 to 2020: A case study in the Mu Us Sandy Land, China. Ecol. Indic. 137,
108745 https://doi.org/10.1016/j.ecolind.2022.108745.
Hansson, K., Šimůnek, J., Mizoguchi, M., Lundin, L.C., van Genuchten, M.T., 2004. Water
flow and heat transport in frozen soil: Numerical solution and freeze-thaw
applications. Vadose Zone J. 3 (2), 693–704. https://doi.org/10.2113/3.2.693.
Harlan, R.L., 1973. Analysis of coupled heat-fluid transport in partially frozen soil. Water
Resour. Res. 9 (5), 1314–1323. https://doi.org/10.1029/WR009i005p01314.
Hayashi, M., 2013. The cold vadose zone: Hydrological and ecological significance offrozen-soil processes. Vadose Zone J. 12, 37–44. https://doi.org/10.2136/
vzj2013.03.0064.
Hayat, M., Zha, T., Jia, X., Iqbal, S., Qian, D., Bourque, C.P.A., Khan, A., Tian, Y., Bai, Y.,
Liu, P., Yang, R., 2020. A multiple-temporal scale analysis of biophysical control ofsap flow in Salix psammophila growing in a semiarid shrubland ecosystem ofnorthwest China. Agr. Forest Meteorol. 288–289, 107985 https://doi.org/10.1016/j.
agrformet.2020.107985.
He, Z., Zhang, S., Teng, J., Yao, Y., Sheng, D.C., 2018. A coupled model for liquid watervapor-heat migration in freezing soils. Cold Reg. Sci. Technol. 148, 22–28. https://doi.org/10.1016/j.coldregions.2018.01.003.
Hou, L., Wang, X., Hu, B.X., Shang, J., Wan, L., 2016. Experimental and numericalinvestigations of soil water balance at the hinterland of the Badain Jaran Desert forgroundwater recharge estimation. J. Hydrol. 540, 386–396. https://doi.org/
10.1016/j.jhydrol.2016.06.036.
Hu, G., Zhao, L., Zhu, X., Wu, X., Wu, T., Li, R., Xie, C., Hao, J., 2020. Review ofalgorithms and parameterizations to determine unfrozen water content in frozen
soil. Geoderma 368, 114277. https://doi.org/10.1016/j.geoderma.2020.114277.
Hu, G., Zhao, L., Li, R., Park, H., Wu, X., Su, Y., Guggenberger, G., Wu, T., Zou, D.,
Zhu, X., Zhang, W., Wu, Y., Hao, J., 2023. Water and heat coupling processes and its
simulation in frozen soils: Current status and future research directions. Catena 222,
106844. https://doi.org/10.1016/j.catena.2022.106844.
Huang, J., Zhou, Y., Wenninger, J., Ma, H., Zhang, J., Zhang, D., 2016. How water use of
Salix psammophila bush depends on groundwater depth in a semi-desert area.
Environ. Earth Sci. 75 (7), 556–568. https://doi.org/10.1007/s12665-016-5376-0.
Iwata, Y., Hayashi, M., Suzuki, S., Hirota, H., Hasegawa, S., 2013. Effects of snow coveron soil freezing, water movement, and snowmelt infiltration: A paired plot
experiment. Water Resour. Res. 46 (9), W09504. https://doi.org/10.1029/
2009WR008070.
Kreyling, J., Peršoh, D., Werner, S., Benzenberg, M., Wöllecke, J., 2012. Short-term
impacts of soil freeze-thaw cycles on roots and root-associated fungi of Holcus
lanatus and Calluna vulgaris. Plant Soil 353, 19–31. https://doi.org/10.1007/
s11104-011-0970-0.
Kurylyk, B., Watanabe, K., 2013. The mathematical representation of freezing andthawing processes in variably-saturated, non-deformable soils. Adv. Water Resour.
60, 160–177. https://doi.org/10.1016/j.advwatres.2013.07.016.
Lundin, L.C., 1990. Hydraulic properties in an operational model of frozen soil.
J. Hydrol. 118 (1-4), 289–310. https://doi.org/10.1016/0022-1694(90)90264-X.
Mohammed, A.A., Kurylyk, B.L., Cey, E.E., Hayashi, M., 2018. Snowmelt infiltration andmacropore flow in frozen soils: Overview, knowledge gaps, and a conceptual
framework. Vadose Zone J. 17, 180084 https://doi.org/10.2136/vzj2018.04.0084.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res. 12 (3), 513–522. https://doi.org/10.1029/
wr012i003p00513.
Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models. Part IA discussion of principles. J. Hydrol. 10 (3), 282–290. https://doi.org/10.1016/
0022-1694(70)90255-6.
Obu, J., Westermann, S., Bartsch, A., Berdnikov, N., Christiansen, H.H., Dashtseren, A.,
Delaloye, R., Elberling, B., Etzelmueller, B., Kholodov, A., Khomutov, A., Kaab, A.,
Leibman, M.O., Lewkowicz, A.G., Panda, S.K., Romanovsky, V., Way, R.G.,
Westergaard-Nielsen, A., Wu, T., Yamkhin, J., Zou, D., 2019. Northern Hemispherepermafrost map based on TTOP modelling for 2000–2016 at 1 km2 scale. EarthScience Rev. 193, 299–316. https://doi.org/10.1016/j.earscirev.2019.04.023.

16