Journal of Hydrology 590 (2020) 125528

Contents lists available at ScienceDirect

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

Research papers

Monitoring and modeling the coupled movement of water, vapor, andenergy in arid areasT

Ce Zhenga,b, Jiří Šimůnekb, Yudong Lua, Xiuhua Liua,c, , Changchun Shid, Huanhuan Lia⁎aSchool of Water and Environment / Key Laboratory of Subsurface Hydrology and Ecological Effects in Arid Region of Ministry of Education, Chang’an University, Xi’an
710064, Chinab
Department of Environmental Sciences, University of California, Riverside, Riverside 92521, USAc
Shaanxi Key Laboratory of Land Consolidation, Chang’an University, Xi’an 710064, Chinad
China State Long-Term Observation and Research Station for Mu Us Desert Ecosystem in Yulin of Shaanxi, Yulin 719000, China

ARTICLE INFOABSTRACTThis manuscript was handled by Emmanouil
Anagnostou, Editor-in-Chief.

In arid areas, vapor flow has been recognized to contribute significantly to the mass and energy transfers and toplay a critical role in maintaining surface vegetation and ecosystems. To better understand the continuous spatial
and temporal variations in liquid water and water vapor contents under different climatic conditions, soil watercontents, temperatures, and micrometeorological variables were observed in-situ in the Mu Us Desert ofnorthwestern China. The collected data were then used to calibrate and validate the Hydrus-1D model simulating either the coupled movement of water, vapor, and energy or only isothermal water flow in soils (the latteras a reference). The results of the coupled model were not only in better agreement with observed data, but also
advanced our understanding of underlying mechanisms of soil water flow. For the coupled model, the isothermal
liquid flux was the most significant component of the total water flux. Three diurnal stages were identified forthermal liquid and vapor fluxes on dry days, while the vapor flux became almost negligible during rainfall
events. The results indicate that isothermal liquid, thermal liquid, and thermal vapor fluxes should be consideredsimultaneously when evaluating soil water flow in arid regions, while the isothermal vapor flux can be neglected. Vapor flow contributed, on average, about 13% of the total water flux in the uppermost soil layer during
the analyzed period at the study area, and this ratio often exceeded 20%.

Keywords:
Vadose zone
Liquid water flow
Vapor movement
Hydrus-1D
Arid region

1. Introduction
Arid and semiarid regions cover more than 40% of the land surface
around the world, and their areas are continuously increasing due tohuman activities and climate change (Gong et al., 2016; Scanlon et al.,2005). In these regions, due to the scarcity of water resources (e.g.,surface water and precipitation), there may exist a series of ecologicalissues, which seriously affects the sustainable development of theseareas and attracts much attention (Kuriqi et al., 2019; Zeng et al.,2009b). During the past several decades, it has been gradually recognized that vapor flow plays a crucial role in maintaining plant
growth and healthy ecosystems, and many researchers have emphasized the role of water vapor in evaluating soil water flow (Banimahd
and Zand-Parsa, 2013; Bittelli et al., 2008). In arid regions, soil water
movement takes place not only in the liquid phase but also in the formof vapor flow, especially under conditions where temperature gradientsare large (Deb et al., 2011; Scanlon et al., 2003). Under very dry soil⁎

conditions, water vapor and liquid water fluxes may be on the sameorder of magnitude near the soil surface, and vapor flow was reported
to account for almost 30% of the total water flux (Cahill and Parlange,
1998; Zhang et al., 2016a). Additionally, vapor flow has a great effecton heat movement, because it carries a substantial amount of the latentenergy of vaporization (Sakai et al., 2009; Zhang et al., 2016b). Since
vapor flow was estimated to account for about 40–60% of the total heat
flux near the soil surface (Cahill and Parlange, 1998), soil temperatures
would be severely underestimated if vapor flow was not considered in
the energy transfer (Kroener et al., 2014). Therefore, it is essential inmany arid-region studies and applications to simultaneously evaluateliquid water, water vapor, and heat transport.
Although the importance of the coupled nature of soil water andheat movement has been long recognized, only isothermal conditions
based on the Richards equation are considered in most studies of soil
water flow. The pioneering theory for the coupled movement of waterand heat in unsaturated soils under non-isothermal conditions was

Corresponding author at: School of Water and Environment, Chang’an University, Xi’an 710064, China.
E-mail address: liuxh68@chd.edu.cn (X. Liu).

https://doi.org/10.1016/j.jhydrol.2020.125528
Received 15 June 2020; Received in revised form 2 September 2020; Accepted 6 September 2020
Available online 15 September 20200022-1694/ © 2020 Elsevier B.V. All rights reserved.

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

developed by Philip and de Vries (1957) (hereafter called the PDV
model). Philip and de Vries (1957) divided the total water flux into four
components: the isothermal liquid, isothermal vapor, thermal liquid,
and thermal vapor fluxes. Subsequently, Milly (1984) modified the PDV
model by adopting the matric head as the driving force for water flow
instead of the total water content. Cass et al. (1984) estimated the enhancement factor of the thermal conductivity for the vapor flux andcompared different approaches. Nassar and Horton (1989) additionally
included in the PDV model the effects of the osmotic pressure on liquidwater and water vapor flow.
The resulting PDV model can simultaneously calculate water flow
driven by soil temperature, pressure head, and osmotic head gradients,and heat flow driven by conduction, convection, and dispersion. Thismodel has been extensively used in many studies (e.g., Deb et al., 2011;
Huang et al., 2016; Wang, 2015; Zeng et al., 2011). As a result of
growing computational capacity and simulation accuracy, numericalmodels based on underlying mechanisms of coupled water, vapor, andheat movement has attracted increasing interest (e.g., Banimahd andZand-Parsa, 2013; Saito et al., 2006; Yu et al., 2018).
Although the studies mentioned above addressed issues related tothe simultaneous movement of liquid water and water vapor, modelingstudies of continuous spatial and temporal variations in driving forces,
and water and vapor fluxes under different climatic conditions, as wellas quantitative evaluations of the contribution of the water vapor flux
to the overall moisture transfer, are still limited. Therefore, the specificobjectives of this study are (i) to monitor the variations of soil watercontents, temperatures, and micro-meteorological variables in the Mu
Us Desert of northwestern China, (ii) to develop and calibrate numerical models for simulating either the coupled movement of water,
vapor, and heat or isothermal water flow (as a reference), (iii) to analyze spatial and temporal variations in soil water contents, temperatures, and water fluxes during the selected typical dry and wet periods,and (iv) to discuss the seasonal movement of soil water and the contribution of vapor flow. The manuscript is organized as follows: Section
2 introduces information about the experimental site and the numericalmodel, Section 3 provides the main modeling results, Section 4 presents
the discussion and interpretation of the main results, and finally,
Section 5 provides conclusions.

properties (Table 1) indicate that the soil profile can be divided into fivelayers. Sand is the dominant soil fraction, accounting for about 95% of
the soil mass in the vadose zone, except in the soil horizon between 80
and 230 cm, where the clay content was slightly higher.
A monitoring well was dug up at the study site down to a depth of
900 cm to install sensors at different depths. Volumetric soil watercontents and temperatures were measured at depths of 10, 20, 50, 100,
130, 200, 400, and 630 cm by the Hydro Probe sensors (Stevens Inc.,
USA). Before monitoring, the sensors were calibrated using gravimetricmeasurements taken from the same soil horizons. Soil heat fluxes weremonitored using the HFP01 probes (Hukseflux Inc., the Netherlands)
installed at depths of 20, 200, and 800 cm. The groundwater level,
which was located at a depth of about 880 cm below the ground in
April, was observed using the Submersible Depth Transmitter probe(Stevens Inc., USA). All sensors in the monitoring well were connectedto the CR1000 data-logger, which recorded data at a 1 min interval andstored them in 10 min averages.
Additionally, micro-meteorological data, including air temperature,
precipitation, relative humidity, and wind speed, were available at the
study site at a 30-minute interval at fixed positions of 20, 20, 200, and240 cm. Fig. 2 shows the measured air temperature and rainfall data at
the study site in 2018. The cumulative rainfall reached a value of
658 mm in 2018, while 97% (about 640 mm) of the total amount occurred between April and October.
2.3. Numerical modeling2.3.1. Water flow
In this study, the Hydrus-1D model was applied to simulate thecoupled transport of water, vapor, and heat in the vadose zone(Šimůnek et al., 2008). For one-dimensional vertical water flow, thegoverning equation for liquid water and water vapor flow in a variablysaturated rigid soil is given by the following mass conservation equation (Saito et al., 2006):

=

t

l

+

t

v

t

=

ql
z

qv
z

S

(1)
3

-3

where θ is the volumetric water content (L L ), which is equal to thesum of the volumetric liquid water θl and water vapor θv contents (L3 L3), t is time (T), ql and qv are the flux densities of liquid water and water
vapor (L T−1), respectively, z is the spatial coordinate positive upward(L), and S is a sink term (T−1), which usually accounts for root water
uptake and can be neglected for bare soil.
Driven by the pressure head, temperature, and gravitational gradients, the total water flux qtotal (L T−1) can be divided into four
components, as shown in Eq. [2]:2. Materials and methods
2.1. Study site
The field observations were conducted at the Yu Lin Desert
Ecosystem Research Station in the Mu Us Desert, northwestern China,
with coordinates of 109°42′29″E and 38°23′19″N (Fig. 1). The soil
surface is flat, and vegetation (mainly Salix Pasmmophila) is quite
sparse, covering less than 5% of the land surface. The climate at the
study site can be characterized as dry and cold, with large diurnal andseasonal temperature differences. The average annual air temperatureis 6.4 °C, with a diurnal range of about 20 °C. The mean annual precipitation and potential evaporation are 360 and 2343 mm, respectively. Rainfall is concentrated mainly in summer, from July to September, when precipitation contributes over 60% of the annual rainfall.
The groundwater level usually fluctuates around 900 cm below the soil
surface during most of the year. Although soil temperatures usually
become negative between November and March, the soil freezingthawing process is not taken into consideration in this study (Zheng
et al., 2020).

qtotal = qlh + qlT + qvh + qvThThT+ 1 + KlT
+ Kvh
+ KvT
z
z
z
z

= Klh

(2)

where qlh, qlT, qvh, and qvT are the isothermal liquid water flux, thermal
liquid water flux, isothermal water vapor flux, and thermal water vapor
flux (L T−1), respectively, Klh (L T−1) and KlT (L2 K−1 T−1) are theisothermal and thermal hydraulic conductivities for liquid water, respectively, Kvh (L T−1) and KvT (L2 K−1 T−1) are the isothermal andthermal hydraulic conductivities for water vapor, respectively, h is thepressure head (L), and T is the temperature (K).
Combining Eqs. [1] and [2], the governing equation for liquid water
and water vapor flow is obtained as follows (Saito et al., 2006):

2.2. In-situ measurementsltSoil samples were collected at a 20-cm interval from the soil surface
to the 300 cm depth and at a 50-cm interval in the rest of the profile(down to a depth of 630 cm) to determine soil physical properties.
Three replicates were taken at each depth. The collected soil physical

+

v

t

=

z

Klh

h
T
h
T
+ Klh + KlT
+ Kvh
+ KvT
z
z
z
z

(3)

In the text below, we will refer to the model that considers all termsof Eq. [3] as the coupled model, and the model that only considers
liquid flow due to isothermal forces as the isothermal model.
2

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

Fig. 1. Location of the study site and the field instrumentation layout.
Table 1
Measured soil particle distribution and optimized soil hydraulic parameters.
Soil layer (cm)Soil particle composition

0–80
80–160
160–230
230–560
560–630

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
94.0
95.4

3.2
5.4
5.2
4.7
3.5

1.0
4.5
3.1
1.3
1.1

0.008
0.015
0.011
0.010
0.011

0.375
0.427
0.415
0.402
0.380

0.027
0.028
0.029
0.021
0.023

1.57
1.51
1.64
1.87
1.76

1368
816
1200
1224
1272

2.3.2. Soil hydraulic properties
The soil water retention curve is the fundamental soil hydrauliccharacteristic required to solve the water flow equation. Using the van
Genuchten (1980) model (Eq. [4]) and the pore-size distribution model(Eq. [5]) of Mualem (1976), the soil water retention function and theisothermal hydraulic conductivity function Klh can be described asfollows, respectively:srr + [1 + | h|n ]m

l =

s h

Klh = Ks Sel [1

Kvh =

KvT =

l

r
r

(5)

Se1/ m )m]2

(8)

d sv
dT

(9)
−1

where D is the vapor diffusivity in soil (L T ), which is defined as
D = Daθaτ, where Da is the vapor diffusivity in air (L2 T−1,
=2.12 × 10−5(T/273.15)2), θa is the air-filled porosity (L3 L−3), τ(=θa7/3/θs2) is the tortuosity factor (-) as defined by Millington and
Quirk (1961), ρw is the density of liquid water (M L-3), ρsv is the saturated vapor density (M L-3), Hr is the relative humidity (-), M is themolecular weight of water (M mol−1, =0.018015 kg mol−1), g is thegravitational acceleration (L T−1, =9.81 m s−2), R is the universal gas
constant of water vapor (J mol−1 K−1, =8.315 J mol−1 K−1), and η isthe enhancement factor (−).
The ρsv can be expressed as a function of temperature:

(4)

-3

(6)

The thermal hydraulic conductivity KlT is defined as follows:

1 d
KlT = Klh hGwT
0 dT

Hr

Mg
RT

2

where θr and θs are the residual and saturated water contents (L L ),respectively, α(L-1), n (-), and m (-) are empirical parameters, Ks is thesaturated hydraulic conductivity (L T−1), l is equal to 0.5 as defined by
Mualem (1976) (-), and Se is the effective liquid saturation (-), which is
described as:
s

D

sv Hr

w

3

Se =

D
w

h<0

0
(1

The isothermal and thermal hydraulic conductivities for vapor
fluxes due to pressure head (Kvh) and temperature (KvT) gradients aredefined, respectively, as follows (Saito et al., 2006):

sv =

(7)

(

exp 31.3716

6014.79
TT

7.92495 × 10 3T) × 103(10)

The relative humidity Hr can be calculated using a thermodynamic
relationship between liquid water and water vapor in the soil as follows(Philip and de Vries, 1957):

where Gwt is the gain factor (−), which amends the temperature dependence of the soil water retention curve (defined as 7 in this study), γ
is the surface tension of soil water (M T−2), which can be calculated as
γ = 75.6–0.1425 T-2.38 × 10-4T2, and γ0 is the surface tension at 25 °C(M T−2, =71.89 g s−2).

Hr = exp

3

Mgh
RT

(11)

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

Fig. 2. (a) Measured daily air temperature and precipitation at the study site in 2018; (b), (c), and (d) hourly rainfall during three selected periods.

The enhancement factor, which was derived by Cass et al. (1984) to
describe an increase in the thermal vapor flux due to liquid islands andincreased temperature gradients, can be calculated as follows:= 9.5 + 3

8.5 exp
s

2.6
1+
fc

sensible heat as described by Fourier’s law, sensible heat by convectionof liquid water and water vapor, respectively, and latent heat by vapor
flow:

4

qh =

s

(12)

Q

Cp T
t

=

z

( l)

T
z

Cw

ql T
z

Cv

qv T
z

L0

qv
z

Cw ST

2.3.4. Soil thermal properties
The soil thermal conductivity depends on both the water content
and soil texture, and can be calculated in Hydrus-1D as follows (Chungand Horton, 1987):(14)

where Cn, Cw, Cv, and Cp are volumetric heat capacities (M L T K−1)of solid, liquid water, water vapor, and moist soil, respectively, θn is thevolumetric solid phase fraction (L3 L-3), and L0 is the volumetric latent
heat of vaporization of liquid water (M L-1 T−2), which can be expressed as follows:

L0 = L w wvt

where the last term on the right side represents energy uptake associated with root water uptake that can be neglected for bare soil.

where Sh is the storage of heat in the soil (M L-1 T−2), qh is the total heat
flux density (M T−3), and Q accounts for sources and sinks of energy (ML-1 T−3). The storage of heat is calculated as follows:
-1

+ L0

(17)

(13)

Sh = Cn T n + Cw T l + Cv T v + L 0 v = Cp T + L 0 v

(16)
−1

where λ is the soil thermal conductivity (M L T
K ).
Combining Eqs. [13], [14], and [16], the governing equation forheat movement is described as follows (Nassar and Horton, 1992; Saito
et al., 2006):

2.3.3. Heat transport
The governing equation for the heat movement in a variably-saturated rigid soil is given by the following energy conservation equation(de Vries, 1958):

qhzT+ Cw Tql + Cv Tqv + L 0 qv
z−3

where fc is the mass fraction of clay in the soil (-).

Sh
=
t

( l)

−2

( l ) = b1 + b2 l + b3 l0.5

(18)where b1, b2, and b3 are empirical parameters that depend on soil textural classes (M L T−3 K−1) and are defined in the Hydrus-1D code as
0.228, −2.406, and 4.909 W m−1 K−1 for sands, respectively.

(15)

2.3.5. Initial and boundary conditions
The soil profile was considered to be 630 cm deep, and the spatial
discretization of 1 cm was applied, leading to 631 nodes across the flow
domain. As initial conditions, we applied the measured soil waterwhere Lw is the latent heat of vaporization of liquid water (M L-1 T−2,
=2.501 × 106-2369.2 T J kg−1).
The total heat flux density is defined as the sum of the conduction of4

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

contents and temperatures at 0:00 on April 1, 2018. Due to the limited
number of observation points in the profile, the linear interpolation wasused to obtain the initial values between different measurement depths.
The atmospheric boundary condition was imposed at the soil surface to account for time-dependent data of precipitation and potential
evaporation for liquid water and water vapor flow. Hourly potential
evaporation (Ep) was estimated as follows:

Ep =

vs

va

rv + rs

respectively) were solved using the finite elements and finite differences methods for spatial and temporal discretizations by the Hydrus1D model, respectively. A detailed description of the calculation procedure can be found in Šimůnek et al. (2008).
The measured soil water contents and temperatures data at various
depths (10, 50, 100, 200, 400, and 630 cm) from April 1 to October 31,
2018, were used for model calibration. The optimized soil hydraulicparameters are listed in Table 1. Then the measured data at the same
depths from April 1 to June 30, 2019, were used for model validation.
In addition to the coupled model, the isothermal model was also run asa reference. The performance of both models was evaluated using theroot mean square error (RMSE) and the Nash-Sutcliffe efficiency (NSE,Nash and Sutcliffe, 1970), defined as follows:(19)

where ρvs and ρva are vapor densities (M L-3) at the soil surface and the
20-cm height (where the air temperature is measured), respectively,
and rv and rs are the aerodynamic and soil surface resistances (T L-1) for
vapor movement, respectively. Detailed calculations of the parameters
can be obtained from Bittelli et al. (2008) and Kroener et al. (2014).
Calculated values of Ep were in close agreement with the observed dataat the Yulin meteorological station (a national weather station), which
is approximately 15 km from our study site. Since the groundwaterlevel fluctuated around a depth of 900 cm and there was the lack of
accurate bottom flux data, the free drainage boundary condition wasconsidered at the bottom of the profile at a depth of 630 cm (far above
the groundwater table).
For the soil heat transport, hourly surface temperatures, estimatedfrom measured soil temperatures at shallow depths (Dahiya et al.,2007), and measured soil temperatures at the 630 cm depth were applied as the upper and bottom boundary conditions, respectively.

N

RMSE =

1
(Si
n i=1

NSE = 1

N
(Si
i=1
N
(Oi
i=1

Oi ) 2

(20)

Oi )2
O¯ )2

(21)

where N is the number of observations, Si and Oi are the simulated andobserved values at a given time, respectively, and‾O is the mean ofobserved values. The NSE value close to 1 indicates that simulated
values match very well measured values, while the NSE value close to 0
indicates that the mean of measured data describes the data better thansimulated values (Amin et al., 2014; Wang et al., 2010). Additionally,
the coefficient of determination (R2) was calculated to indicate if there
exists a linear relationship between observed and simulated values(Stumpp et al., 2012).

2.3.6. Water flow and heat transport simulations
The flowchart shown in Fig. 3 summarizes the applied methodology,
which includes four main steps. After collecting experimental data in
the first step, the numerical model was set up in the second step. In
addition to simplifications involved in specifying the boundary conditions, the lack of laboratory experiments to determine the soil hydraulicand thermal properties was another limitation. The parameters representing these properties were instead initially predicted using theneural network method from measured soil physical properties, andthen optimized using the Levenberg-Marquardt algorithm in the Hydrus-1D model (Šimůnek et al., 2012). After the model setup, the soil
water flow and heat transport equations (Eqs. [3] and [17],3. Results
3.1. Model calibration and validation
During both calibration and validation periods, the simulated soil
water contents at selected depths fitted the observed data well (Fig. 4).
Due to complex weather conditions (such as frequent precipitations andintense evaporation), deviations between simulated and measured soil
water contents in the topsoil layers were larger than at deeper depths.

Fig. 3. The conceptual flowchart of the applied methodology. Step 3 describes the calculation procedure during one timestep in the Hydrus-1D model.
5Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

Fig. 4. Comparison of observed and simulated daily (a) soil water contents and (b) soil temperatures at 10, 50, 100, 200, 400, and 630 cm depths during thecalibration (from April 1 to October 31, 2018) and validation (from April 1 to June 30, 2019) periods. There is an interval of 5 months between the calibration andvalidation periods.

6

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

Table 2Model performance statistics for simulated soil water contents (SWC) and temperatures (ST) at depths of 10, 200, and 400 cm.
Criteria

SWC(Coupled model)
SWC(Isothermal model)
ST(Coupled model)

RMSE
NSE
R2
RMSE
NSE
R2
RMSE
NSE
R2

Calibration depths (cm)

Validation depths (cm)

10

200

400

10

200

400

0.017
−0.19
0.69
0.018
−0.35
0.67
0.363
0.98
0.98

0.005
0.05
0.87
0.008
−0.08
0.85
1.394
0.89
0.98

0.002
0.92
0.94
0.003
0.84
0.88
0.817
0.92
0.99

0.014
−0.26
0.68
0.014
−0.26
0.68
1.109
0.96
0.96

0.002
0.31
0.81
0.004
0.02
0.83
0.641
0.96
0.99

0.003
−2.31
0.84
0.004
−6.10
0.82
0.303
0.95
0.97

Note: The units of the RMSE are cm3cm−3 and °C for soil water contents and temperatures, respectively.

Simulated results predicted well both sharp increases due to precipitation and gradual decreases due to evaporation and redistribution.
There exist two main differences in simulation results between theisothermal and coupled models. First, soil water contents simulated by
the isothermal model were relatively lower under continuously dry
conditions compared to those simulated by the coupled model. Second,
the peaks of soil water contents due to rainfall events were overestimated by the isothermal model. The reasons for these discrepancies
are discussed below in detail. In contrast, simulated temperatures fitted
the observed data very well at most depths. The discrepancies in 200
and 400 cm depths are likely caused by selected values of thermal capacities and conductivities, which were not changed from the default
values in the model.
The statistical criteria quantifying the model performance aresummarized in Table 2. For soil water contents, RMSEs were largest at
the shallow depth and decreased with increasing depth. At a certain
depth, RMSEs for the isothermal model were slightly higher than thosefor the coupled model (e.g., 0.008 and 0.005 cm3 cm−3 at the 200 cm
depth for the isothermal and coupled models during the calibration
period, respectively), indicating a better simulation accuracy of thecoupled model. On the other hand, NSE values varied greatly (from−6.1 to 0.92) with depths. Since soil water contents in this study varied
in a very narrow range (especially during the validation period), thesmall denominator of the second term in Eq. [21] resulted in negative
NSE values. This indicates that NSE is not a useful criterion when
variations in observed soil water contents are small (Stumpp et al.,2012). Instead, R2 can be used to evaluate the model prediction precision. Calculated R2 values indicated that simulated water contents
were in good agreement with observed data, particularly in deeper
depths (with R2 > 0.8). For soil temperatures, all three indices indicated a strong correlation between observed and simulated data. For
example, R2 is greater than 0.95 at all depths for both calibration andvalidation periods.

the shallow layer, the matric potential gradually increased with depth
during the three dry days, which caused an upward driving force forwater flow. For isothermal conditions, the matric potential at mostdepths showed a gradually declining trend with time, indicating a soil
drying process. For example, the matric potential decreased from−1134 to −1155 cm at the 10 cm depth on the first day of Period 2.
Small diurnal fluctuations in the matric potential were simulated by the
coupled model as a result of the day/night behavior of liquid water flow
and vapor flow due to temperature gradients (as discussed below indetail). During the rainy period, the matric potential profiles completely changed, and the contours in the top 20 cm depth became instantly dense. This increase in the matric potential quickly propagatedto deeper depths as the infiltration front moved deeper. Once precipitation stopped, the contours of the soil matric potential became
sparse again from top to bottom.
The matric potential below a certain depth showed no obvious
fluctuations (usually less than 5 cm) in each period. Influenced byseasonal climate variations, the matric potential profiles in deep layer
were clearly different in the three periods, with the lowest matric potentials occurring at depths of 330, 600, and 80 cm, respectively. At thesame time, there existed some distinct differences between the simulated results of the two models. For instance, the matric potential simulated by the isothermal model was lower than that simulated by the
coupled model at the same depth in Period 2 (e.g., about −550 and−520 cm for the isothermal and coupled models at the 600 cm depth,
respectively).
In addition to changes in the matric potential, variations of soil
temperature also represent a vital factor affecting soil water movementin the coupled model. Soil temperature variations are shown in Fig. 6,
with contour intervals of 1 and 0.5 °C for shallow (0–80 cm) and deep
(80–600 cm) layers, respectively. In the shallow layer, soil temperatures
exhibited apparent diurnal variations, and these fluctuations decreased
with depth. While the temperature contours were relatively sparse before dawn close to the surface, they became dense during the day time.
Note that each day there exists two contours for the zero-temperature
gradient, indicating the change in the direction of the temperature
gradient. The first zero contour is the divergent plane (which usuallyoccurred during nighttime), with temperature gradients above andbelow this plane being positive (upward) and negative (downward),
respectively. On the other hand, the second zero contour is the convergent plane (which usually occurred during daytime), with temperature gradients above and below this plane being reversed. Soil
temperatures were relatively high during Period 2 due to intense solar
radiation, with an average value of 32.3 °C at the soil surface, while themean surface temperatures in Periods 1 and 3 were 16.9 and 12.8 °C,
respectively. During the rainy period, the contours became more sparse,
indicating that rainfall events decreased diurnal variations in soil
temperature.
Being less affected by solar radiation and precipitation, soil temperatures in deeper layers showed no obvious diurnal changes in both3.2. Variations in soil matrix potential and temperatureThree typical periods in spring (Period 1, from April 21 to 23),summer (Period 2, from August 6 to 8), and autumn (Period 3, from
October 13 to 15) of 2018 were selected to analyze the spatial andtemporal distributions of the soil matric potential and temperature. The
first day of each period represented dry conditions (without rainfall forabout ten days), while rainfall occurred during the next two days (asshown in Fig. 2).
Fig. 5 shows the soil matric potential profiles simulated by both thecoupled and isothermal models for the three selected periods. Counterintervals for shallow and deep soil profiles are 50 and 10 cm, respectively. Affected by evaporation, soil water contents near the soil surface
were very low, which resulted in extremely low soil matric potentials inthis layer (1–2 orders lower than in the rest of the profile). For thisreason, the soil matric potential of the 0–10 cm depth is not depicted. In
7

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

Fig. 5. Spatial and temporal distributions of soil matric potentials simulated by the coupled (a, b, c) and isothermal (d, e, f) models for Period 1 (a, d), Period 2 (b, e),and Period 3 (c, f). Different vertical scales and counter levels are used for the upper (0 to 80 cm) and lower (80–600 cm) parts of the soil profile.

dry and rainy days. However, there existed clear distinctions in soil
temperature profiles during different periods. During the three evaluated periods, the highest temperatures were 13.1, 27.3, and 17.1 °C at
depths of 80, 80, and 306 cm, respectively, while the lowest temperatures were 8.2, 12.7, and 13.5 °C at depths of 320, 600, and 80 cm,
respectively.

reason for the overestimation of the total flux by the isothermal model,
which in turn caused the relatively lower soil water content. However,
the effect of qvT flux became almost negligible during wet days. The
downward and upward isothermal liquid fluxes (qlh) driven by thematric potential gradient dominated flow throughout the soil profile
during (46 h) and after (58 h) the rainfall event, respectively. During
this time, the differences between total water fluxes simulated by thetwo models were smaller than before the rainfall event and were caused
mainly by the thermal liquid flux (qlT). It can be seen that qlT flux in the
top 20 cm layer at 46 h became 1–2 orders of magnitude larger than
before rainfall.
Since the qlh flux dominated soil water flow in deeper layers underboth dry and wet conditions, the total fluxes simulated by the twomodels were similar. The maximum difference occurred at the 250 cm
depth, where the total water flux simulated by the isothermal model(-0.047 cm d-1) was lower than that the one simulated by the coupledmodel (-0.062 cm d-1). This difference could be mainly attributed toneglecting the qlT flux.
Table 3 lists the quantitative contribution of each soil water flux tothe total water flux during Period 2. The qlh flux is the most significant
component throughout the entire profile, followed by the qvT and qlh
fluxes. The contribution of the qvT flux during dry days was considerable in the shallow layer (e.g., 25.9% at 2 cm and 13.2% at 10 cm),while the ratio between the thermal vapor and total water fluxes

3.3. Comparison of water fluxes simulated by the coupled and isothermal
modelsAffected by different driving (isothermal and thermal) forces, soil
water fluxes simulated by the coupled and isothermal models weredifferent. Fig. 7 shows the vertical flux profiles simulated by bothmodels at different times in Period 2, including 0 h (midnight between
dry days), 12 h (noon of a dry day), 46 h (during rain), and 58 h (afterrain).
There are distinct differences between fluxes simulated by the twomodels in the shallow layer during dry days, which are caused mainlyby the thermal vapor flux (qvT). For example, the upward total fluxes atthe 2 cm depth simulated by the isothermal and coupled models were0.094 and 0.039 cm d-1 at 12 h (noon of a dry day), respectively. In
contrast, the downward qvT flux simulated by the coupled model was−0.083 cm d-1 due to a relatively large temperature gradient.
Neglecting the downward qvT flux under dry conditions was the major
8Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

Fig. 6. Spatial and temporal distributions of soil temperatures simulated by the coupled model for Period 1 (a, d), Period 2 (b, e), and Period 3 (c, f). Different verticalscales and counter levels are used for the upper (0 to 80 cm) and lower (80–600 cm) parts of the soil profile. The solid and dotted lines represent convergent anddivergent temperature gradient planes, respectively.

−4.9 cm d-1 in the three periods, respectively, corresponding to applied
rainfall (Fig. 2). The peak of downward flux decreased with time (again
reflecting rainfall) and gradually moved to deeper soil depths. For instance, the depth of the flux peak in Period 1 dropped from 10 to 47 cm
between 50 and 68 h, with the flux ranging from −5.01 to −0.51 cm d1. After rainfall, a divergent plane occurred in the shallow layer andmoved downward, indicating that the qlh flux was again upward toward
the surface due to evaporation. Due to different rainfall amounts, theinfiltration front reached the depth of 98, 71, and 51 cm after 72 h in
the three periods, respectively.
Driven by a temperature gradient, diurnal changes in the qlT and qvTfluxes during three selected dry days are plotted in Fig. 9 and Fig. 10,
respectively. Based on the direction of the qlT fluxes at the soil surface,
three stages can be identified during dry days. For example, in Period 2,
the first stage occurs from midnight to before sunrise (around 7 a.m.)when the qlT flux is directed toward the soil surface. During this stage, adivergent plane continuously dropped from a depth of about 20 cm at
0 h to a depth of about 35 cm at 6 h. In the meantime, the depth withthe maximum upward qlT flux dropped as well, with the flux varying
from 0.005 to 0.007 cm d-1. The qlT flux at the soil surface became
negative (downward) during the second stage (between 8 and 18 h),
and the downward flux increased as the air temperature increased. The
position of the maximum downward flux gradually dropped from adepth of 10 cm to about 20 cm, with the maximum downward flux of−0.021 cm d-1. During the third stage (after 18 h), the surface flux
became again positive (upward), and a divergent plane occurred above
the 20 cm depth. This stage can be regarded as a transitional stage when
the flux profiles change from that of the second stage to that of the firststage. The maximum downward fluxes (at about 35 cm) dropped to−0.017 cm d-1 at 22 h.
Similarly, there exist three stages for the qvT flux as well. Contrary to
the qlT flux, the largest vapor fluxes usually occurred at the soil surface
and decreased sharply with depth. Fluctuations in the qvT flux meantime showed a gradually declining trend. For example, in Period 2, theqvT flux varied from −0.038 to 0.013 cm d-1 and from −0.017 to
0.005 cm d-1 at depths of 10 and 20 cm, respectively, while fluctuationsat depths of 40 and 60 cm were approximately 6 and 12 times smaller.
Below the 60 cm depth, the qvT flux tended to remain stable during a

gradually decreased with increasing depth. In contrast, the contributionof the qlT flux in the deep layer was stable and not negligible. Affectedby the soil texture, the qlh flux decreased in the soil layer between 80
and 160 cm, which resulted in a relatively large contribution of qlT flux
to the total flux (e.g., 31.8% and 27.9% at the 100 cm depth under dry
and wet conditions, respectively). Besides, the isothermal water vapor
flux (qvh), which was always 2–5 orders of magnitude smaller than theother fluxes, can be ignored under both dry and wet conditions.
In general, the results simulated by the coupled model have not only
higher accuracy, but also advance our understanding of the underlying
mechanisms for soil water flow. Therefore, the coupled model should beapplied in evaluating soil water flow in arid regions, especially under
dry conditions with relatively low soil water contents.
3.4. Diurnal soil water dynamics
After comparing the performance of the two models, the results simulated by the coupled model during three selected periods were used
to investigate further the diurnal dynamics of the qlh, qlT, and qvT fluxes
on both dry and wet days. Due to small diurnal variations in water
fluxes in deeper soil layers, only fluxes between the soil surface and a
depth of 80 cm are considered.
Fig. 8 depicts variations of the qlh flux during three selected periods.
During dry days, the qlh flux is upward throughout the entire shallow
soil profile, which indicates that the qlh flux was dominated by thematric potential gradient rather than by the gravitational gradient.
Affected by a relatively high gradient, the qlh flux was large near thesurface, with the maximum values of 0.119, 0.148, and 0.098 cm d-1
during the three periods, respectively. Below the 20 cm depth, the qlh
flux gradually decreased with depth. For instance, the highest qlh values
at depths of 20, 40, and 80 cm in Period 2 were 0.098, 0.082, and0.043 cm d-1, respectively. Additionally, the qlh fluxes were higher insummer compared to those in spring and autumn, which was partly
attributed to the higher matric potential gradient caused by intensive
evaporation in summer.
Due to the negative matric potential gradient, the qlh flux was
downward (i.e., infiltration) during rainy days. The surface flux became
extremely large, with the maximum values of −18.8, −11.1, and9

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

Fig. 7. Simulated vertical flux profiles at different times in Period 2: (a) 0 h, (b) 12 h, (c) 46 h, and (d) 58 h. (qtc and qti are total fluxes simulated by the coupled andisothermal models, respectively; qlh, qlT, qvh, and qvT are isothermal liquid, thermal liquid, isothermal vapor, and thermal vapor fluxes simulated by the coupledmodel).

day. Due to larger temperature gradients, the qvT fluxes were greater insummer. The maximum downward and upward fluxes were −0.119
and 0.038 cm d-1 during Period 2, while they varied from −0.058 to
0.015 cm d-1 and −0.039 to 0.011 cm d-1 during Periods 1 and 3, respectively.
The thermal water flux showed no obvious diurnal variations during
rainy days, and both the qlT and qvT fluxes displayed very different
profiles compared to dry days (Fig. 11). The upward qlT fluxes in theshallow 20 cm layer became very large, with the maximum values of
0.804, 0.512, and 0.202 cm d-1 during the three periods, respectively.
These fluxes were over one order of magnitude larger than the maximum upward qlT fluxes during dry days. On the other hand, the qvT fluxin the uppermost soil layer decreased rapidly and dramatically during
the rainy period. For example, upward surface fluxes at 48 h were 5.6e04, 4.1e-03, and 7.2e-04 cm d-1 during the three periods, respectively,
while they were 9.6e-03, 0.029, and 8.4e-03 cm d-1 at 24 h,
10Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

maximum downward qlT flux (-0.052 cm d-1) was smaller than the qvTflux (−0.144 cm d−1). On the other hand, both the upward qlT and qvTfluxes were relatively small below the convergent plane. In the meantime, a divergent plane of the matric potential gradient, which wasaffected by evaporation and precipitation, fluctuated between depths of
80 and 120 cm. Without considering the influence of rainfall, the instantaneous maximum upward and downward qlh fluxes above andbelow this plane were 0.184 and −0.039 cm d-1, respectively. This
indicates that the soil water flux in the shallow soil profile was dominated by the qvT and qlh fluxes, while the downward qlh flux dominatedthe soil water flux in the deep soil profile (Zeng et al., 2009b).
From June to August, the thermal flux was consistently downward
throughout the soil profile (except for the soil surface layer due todiurnal variations in temperature). Variations in downward qlT and qvTfluxes were from −0.058 to −6.9e-05 cm d-1 and −0.148 to −1.5e05 cm d-1, respectively. Affected by frequent rainfall events in thisperiod, the planes of convergence and divergence of the matric potential gradient appeared alternatively in the shallow layer (Zeng et al.,2009a). For instance, actual downward and upward fluxes reached the
maximum values of −37 and 2 cm d-1, respectively, reflecting the
maximum rainfall and evaporation events, and suggesting that the qlh
flux dominated the total soil water flux during the rainy period.
Meanwhile, various evaporation and precipitation events caused thepreviously divergent water potential gradient plane to fluctuate
quickly. Above this plane, the upward qlh flux was less significant thanthe downward qvT flux in the uppermost soil layer (usually above the
10 cm depth), while it exceeded the thermal flux in the rest of the
profile. Below this plane, the thermal liquid water, thermal water
vapor, and isothermal liquid water fluxes contributed together to thedownward soil water movement, which is a major source for groundwater recharge (Hou et al., 2016).
From September to October, soil temperature started decreasing inthe surface layer due to decreasing air temperature. The daily temperature gradient near the surface became positive in early September,and a divergent plane gradually moved to a depth of about 400 cm at
the end of this period. Without considering the influence of rainfall, theTable 3Contributions of each soil water flux (qlh, qlT, qvh, and qvT are isothermal liquid,
thermal liquid, isothermal vapor, and thermal vapor fluxes) to the total flux
during Period 2.
Depth (cm)

2
10
50
100
200
400
600

Dry conditions

Wet conditions

qlh

qlT

qvh

qvT

qlh

qlT

qvh

qvT

72.8%
80.4%
83.5%
63.3%
84.7%
85.1%
91.2%

1.1%
6.4%
13.2%
31.8%
13.5%
11.3%
7.0%

0.2%
0.0%
0.0%
0.0%
0.0%
0.0%
0.0%

25.9%
13.2%
3.3%
4.8%
1.8%
3.5%
1.8%

90.8%
90.4%
89.9%
67.5%
82.9%
85.3%
91.1%

7.9%
8.7%
8.1%
27.9%
15.1%
11.2%
7.1%

0.0%
0.0%
0.0%
0.0%
0.0%
0.0%
0.0%

1.3%
1.0%
2.0%
4.6%
2.1%
3.5%
1.8%

respectively. This phenomenon could be attributed mainly to two factors. First, the temperature gradient in the shallow layer decreased due
to rainfall events, and second, increased water contents in the soil resulted in lower air contents, and reduced vapor flux.
4. Discussion
4.1. Seasonal soil water migration
In general, the thermal and isothermal liquid water and water vapor
fluxes vary greatly in different months during the year, as shown in
Fig. 12. While some of these variations are in response to short-term
meteorological events (e.g., rainfall), significant seasonal trends can bedetected as well. Since rainfall is mainly concentrated in summer, theinfluence of the rainfall process on soil water flow is only considered in
the period from June to August.
From April to May, soil temperature in the upper part of the vadosezone increased due to increasing air temperature. The convergent planeof the soil temperature gradient occurred in early April at a depth ofabout 200 cm, and then it gradually decreased with time to the bottomof the soil profile. Above the convergent plane, the instantaneousFig. 8. Diurnal patterns of isothermal liquid water flux (qlh) profiles for Period 1 (a, d), Period 2 (b, e), and Period 3 (c, f) on selected dry (a, b, c) and rainy (d, e, f)
days.
11

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

Fig. 9. Diurnal patterns of thermal liquid water flux (qlT) profiles during selected dry days: Period 1 (a, b, c) is in the first row, Period 2 (d, e, f) is in the middle row,and Period 3 (g, h, i) is in the bottom row; the early stage (a, d, g) is in the left column, the middle stage (b, e, h) is in the middle column, and the late stage (c, f, i) is in
the right column.

instantaneous maximum upward qlT and qvT fluxes above this plane
were 0.017 and 0.025 cm d-1, respectively. Meanwhile, the divergent
water potential gradient plane returned to a depth of about 120 cm,
with the maximum upward qlh flux of 0.151 cm d-1. This indicates thatboth the thermal liquid and vapor fluxes and isothermal liquid fluxes
are upward towards the shallow soil layer, supplying soil water forevaporation. Below the divergent plane, the maximum downward qlh
flux was −0.052 cm d-1, which dominated soil water flow in the deep
layer.

sinusoidal variations in the soil water content (Deb et al., 2011; Yu
et al., 2018). The quantitative analysis of vapor flow and its contribution to the total mass transfer is thus of great significance.
In this study, the vapor flux accounted for 13.1%, 6.4%, 6.1%, 4.4%,
and 2.5% of the total soil water flux at depths of 2, 10, 100, 200, and400 cm, respectively, from April to October in 2018. The results of oursimulations are similar to those of Deb et al. (2011). They calculated
that the vapor flux was approximately 10.4% of the total water flux
near the soil surface during 50 days with frequent irrigations. Note that
although ratios of the vapor fluxes to the total fluxes are less significant
in deeper depths, that does not mean that vapor flow can be neglected.
The continuous vapor flux driven by a temperature gradient could affect the overall soil moisture distribution (Fig. 4). Similar behavior was
also reported by Grifoll et al. (2005) and Scanlon et al. (2003). In
contrast, Xiang et al. (2012) indicated that the vapor flux dominated thetotal water fluxes in the deep layer. They observed that liquid flow
could only reach a certain depth (of about 185 cm) and that the liquid
flux became extremely small below this depth. This was likely causedby the fact that liquid flow in arid regions rarely occurs below certain
depths in the sandy and loamy layered profiles. This phenomenon did
not occur in our simulations due to a similar soil texture (sandy) in theentire profile at the study site.
Because vapor transport is impeded during and after rainfall events,
frequent and intensive rainfalls caused lower ratios of vapor fluxes to

4.2. The role of vapor flow
It has been increasingly recognized that in arid regions, vapor flow
plays a crucial role in providing water for vegetation and in ecologicalenvironments. The water vapor fluxes are highest in the surface soil
layer, which can be attributed to (1) the large temperature gradients,(2) the high air-filled porosity due to the low soil water content, and (3)the large vapor density gradients between the atmosphere and the soil
(Bittelli et al., 2008). Although during wet conditions, liquid water
fluxes dominated water flow, water vapor and liquid water fluxes wereon the same order of magnitude in the surface soil layer during dry
conditions (Fig. 7), indicating that the soil water movement occurredboth in the form of water vapor and liquid water (Zeng et al., 2009b).
Besides, the day/night behavior of the vapor flux may induce the12Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

Fig. 10. Diurnal patterns of thermal water vapor flux (qvT) profiles during selected dry days: Period 1 (a, b, c) is in the first row, Period 2 (d, e, f) is in the middle row,and Period 3 (g, h, i) is in the bottom row; the early stage (a, d, g) is on the left, the middle stage (b, e, h) is in the middle column, and the late stage (c, f, i) is on theright.

total fluxes in our simulation. While average annual rainfall at ourstudy site is approximately 360 mm, 640 mm of rain occurred during
the analyzed period, which seldom happens in this arid area. To better
identify general trends of vapor flow at the study site, three additional
rainfall scenarios with rain totals of 512, 320, and 128 mm (accounting
for 80%, 50%, and 20% of rain totals during the simulation period,
respectively) were simulated. The modeling results indicated that thevapor flux at a depth of 2 cm would contribute 15.3%, 20.6%, and31.7% of the total moisture flux for these three rainfall scenarios, respectively. Similarly, Zhang et al. (2016a) also reported larger ratios(i.e., over 20%) of the vapor flux to the total water flux in the uppermost layer during a long term observation in an arid field with little
rainfall. This confirms the importance of vapor flow in arid areas andindicates that it is necessary to consider the effects of vapor flow whenevaluating soil water flow.
The above analysis focuses on the non-freezing period, while vapor
flow may play an important role during the freezing period as well.
Since the liquid water content and the hydraulic conductivity for liquidwater significantly decreases due to the presence of ice, soil water may
be mainly transported by vapor flow in the frozen soil. Even the often
neglected isothermal vapor flux could contribute a considerable part of
the total water flux (Yu et al., 2018). Due to its complexity (e.g., phasechanges between ice and liquid water), the coupled movement of water,
vapor, and heat during freezing periods at the field scale is not fully
understood and needs further investigation.
5. Conclusions
In this study, the Hydrus-1D model simulating the coupled movement of water, vapor, and heat was calibrated and validated using themeasured soil moisture and temperature data. The Hydrus-1D model
simulating only isothermal water flow was also run for comparison.
Since it neglects the effects of vapor and heat on water flow, the isothermal model would underestimate and overestimate soil water contents during dry and rainy conditions, respectively. On the other hand,the simulation results of the coupled model were not only closer to themeasured data, they also advanced our understanding of the underlying
mechanisms of soil water flow. Then, diurnal and seasonal variations ofeach soil water flux were simulated by the coupled model. The results
indicated that the isothermal liquid flux was the most significant
component of the total flux under both dry and wet conditions. Three
diurnal stages were identified for the thermal liquid and vapor fluxes indry days, while the vapor flux became almost negligible during rainfall
events.
13

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

Fig. 11. Diurnal patterns of thermal liquid water (qlT) flux profiles (a, b, c) and thermal water vapor flux (qvT) profiles (d, e, f) for Period 1 (a, d), Period 2 (b, e), andPeriod 3 (c, f) on selected rainy days.

In general, it is necessary to simultaneously consider the isothermal
liquid, thermal liquid, and thermal vapor fluxes in studies evaluating
soil water flow in arid regions. Vapor flow accounted for 13.1% of thetotal moisture transfer near the soil surface during our study. However,this ratio usually exceeds 20% at this study site during more common,drier years. Although this study was carried out in northwestern China,

similar results could be expected in other arid regions. The future study
will focus on the coupled movement of water, vapor, and heat during
the freezing period and the contribution of water vapor to the vegetation growth.

Fig. 12. Schematic illustration of seasonal trends in soil water flows: (a) the thermal liquid water (qlT) and water vapor (qvT) fluxes, (b) the isothermal liquid water
flux (qlh). The solid ellipses show the approximate locations with minimum soil temperatures (a) or water potentials (b), and the empty ellipses show the approximate
locations with maximum soil temperatures (a) or water potentials (b). Seasonal trends in various water fluxes (qlh, qlT, and qvT) are indicated by arrows. The
embedded figure in (b) represents the planes of convergence and divergence occurring in the shallow layer during and after a rainfall event.
14

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.

CRediT authorship contribution statement

influence the work reported in this paper.

Ce Zheng: Conceptualization, Methodology, Writing - originaldraft. Jiří Šimůnek: Writing - review & editing. Yudong Lu: Writing review & editing. Xiuhua Liu: Methodology, Supervision. Changchun
Shi: Investigation, Data curation. Huanhuan Li: Investigation, Data
curation.

Acknowledgments
This research was supported by the National Natural Science
Foundation of China (No. 41877179, 41630634), the Fundamental
Research Funds for the Central Universities, CHD (No. 300102298706,
300102290719), the Shaanxi Water Conservancy Science andTechnology plan project (2019slkj-18) and Fund Project of Shaanxi Key
Laboratory of Land Consolidation. Gratitude is also expressed for thefinancial support from the program of China Scholarship Council.

Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared toAppendix

The list of variables (with units) used in the governing equations for water flow and heat transport.
Cn

Volumetric heat capacity of solid, J m−3 K−1

Cp
Cv
Cw
D
g
h
Klh
KlT
Ks
Kvh
KvT
L0
M
qlh
qlT
qvh
qvT
R
S
T
t
z
γ
λ
θl
θn
θr
θs
θv
ρsv
ρv
ρw

Volumetric heat capacity of moist soil, J m−3 K−1Volumetric heat capacity of water vapor, J m−3 K−1Volumetric heat capacity of liquid water, J m−3 K−1
Vapor diffusivity in air, m2 s−1
Gravitational acceleration, m s−2
Pressure head, mIsothermal hydraulic conductivity for liquid water, m s−1
Thermal hydraulic conductivity for liquid water, m2 K−1 s−1
Saturated hydraulic conductivity, m s−1
Isothermal hydraulic conductivity for water vapor, m s−1
Thermal hydraulic conductivity for water vapor, m2 K−1 s−1Volumetric latent heat of vaporization of liquid water, J m−3
Molecular weight of water, kg mol−1
Isothermal liquid water flux, m s−1
Thermal liquid water flux, m s−1
Isothermal water vapor flux, m s−1
Thermal water vapor flux, m s−1
Universal gas constant, J mol−1 K−1
Sink term, s−1
Temperature, K
Time, sVertical depth from the surface, mSurface tension of soil water, J m−2
Soil thermal conductivity, W m−1 K−1Volumetric liquid water content, m3 m−3Volumetric solid phase fraction, m3 m−3Residual water content, m3 m−3Saturated water content, m3 m−3Volumetric water vapor content, m3 m−3Saturated vapor density, kg m−3
Vapor density, kg m−3
Liquid water density, kg m−3

Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi.org/10.1016/j.jhydrol.2020.125528.

WR023i012p02175.
Dahiya, R., Ingwersen, J., Streck, T., 2007. The effect of mulching and tillage on the waterand 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.
Deb, S.K., Shukla, M.K., Sharma, P., Mexal, J.G., 2011. Coupled liquid water, water vapor,and heat transport simulations in an unsaturated zone of a sandy loam field. Soil Sci.
176 (8), 387–398. https://doi.org/10.1097/ss.0b013e318221f132.
de Vries, D.A., 1958. Simultaneous transfer of heat and moisture in porous media. Trans.
Am. Geophys. Union 39 (5), 909–916. https://doi.org/10.1029/TR039i005p00909.
Gong, J.N., Jia, X., Zha, T.S., Wang, B., Kellomaki, S., Peltola, H., 2016. Modeling theeffects of plant-interspace heterogeneity on water-energy balances in a semiarid
ecosystem. Agr. Forest Meteorol. 221, 189–206. https://doi.org/10.1016/j.
agrformet.2016.01.144.
Grifoll, J., Gasto, J.M., Cohen, Y., 2005. Non-isothermal soil water transport and evaporation. Adv. Water Resour. 28, 1254–1266. https://doi.org/10.1016/j.advwatres.
2005.04.008.
Hou, L.Z., Wang, X.S., 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.

References
Amin, M.G.M., Šimůnek, J., Lægdsmand, M., 2014. Simulation of the redistribution andfate of contaminants from soil-injected animal slurry. Agr. Water Manage. 131,
17–29. https://doi.org/10.1016/j.agwat.2013.09.002.
Banimahd, S.A., Zand-Parsa, S.h., 2013. Simulation of evaporation, coupled liquid water,water vapor and heat transport through the soil medium. Agr. Water Manage. 130,
168–177. https://doi.org/10.1016/j.agwat.2013.08.022.
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.
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.
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/

15

Journal of Hydrology 590 (2020) 125528

C. Zheng, et al.
1016/j.jhydrol.2016.06.036.
Huang, J.T., Hou, R.Z., Yang, H.B., 2016. Diurnal pattern of liquid water and water vapor
movement affected by rainfall in a desert soil with a high water table. Environ. Earth
Sci. 75, 73–88. https://doi.org/10.1007/s12665-015-4938-x.
Kroener, E., Vallati, A., Bittelli, M., 2014. Numerical simulation of coupled heat, liquidwater and water vapor in soils for heat dissipation of underground electrical power
cables. Appl. Therm. Eng. 70 (1), 510–523. https://doi.org/10.1016/j.
applthermaleng.2014.05.033.
Kuriqi, A., Pinheiro, A.N., Sordo-Ward, A., Garrote, L., 2019. Flow regime aspects indetermining environmental flows and maximising energy production at run-of-riverhydropower plants. Appl. Energ. 256, 113980. https://doi.org/10.1016/j.apenergy.
2019.113980.
Milly, P.C.D., 1984. A linear analysis of thermal effects on evaporation from soil. Water
Resour. Res. 20 (8), 1075–1085. https://doi.org/10.1029/wr020i008p01075.
Millington, R.J., Quirk, J.P., 1961. Permeability of porous solids. Trans. Faraday Soc. 57,
1200–1207. https://doi.org/10.1039/tf9615701200.
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, 282–290. https://doi.org/10.1016/00221694(70)90255-6.
Nassar, I.N., Horton, R., 1989. Water transport in unsaturated non-isothermal salty soil:
II. Theoretical development. Soil Sci. Soc. Am. J. 53 (5), 1323–1329. https://doi.org/
10.2136/sssaj1989.03615995005300050005x.
Philip, J.R., de Vries, D.A., 1957. Moisture movement in porous materials under temperature gradient. Trans. Am. Geophys. Union 38 (2), 222–232. https://doi.org/10.
1029/TR038i002p00222.
Saito, H., Š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.
Sakai, M., Toride, N., Šimůnek, J., 2009. Water and vapor movement with condensationand evaporation in a sandy column. Soil Sci. Soc. Am. J. 73 (3), 707. https://doi.org/
10.2136/sssaj2008.0094.
Scanlon, B.R., Keese, K., Reedy, R.C., Šimůnek, J., Andraski, B.J., 2003. Variations in flow
and transport in thick desert vadose zones in response to paleoclimatic forcing (0–90kyr): Field measurements, modeling, and uncertainties. Water Resour. Res. 39 (7),1179–1196. https://doi.org/10.1029/2002wr001604.
Scanlon, B.R., Reedy, R.C., Keese, K.E., Dwyer, S.F., 2005. Evaluation of evapotranspiration covers for waste containment in arid and semiarid regions in thesouthwestern USA. Vadose Zone J. 5 (2), 55–71. https://doi.org/10.2136/vzj2005.
0090.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008. Development and applications of the

HYDRUS and STANMOD software packages and related codes. Vadose Zone J. 7 (2),587–600. https://doi.org/10.2136/vzj2007.0077.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2012. HYDRUS: Model use, calibration, andvalidation. T. ASABE. 55 (4), 1261–1274.
Stumpp, C., Stichler, W., Kandolf, M., Šimůnek, J., 2012. Effects of land cover and fertilization method on water flow and solute transport in five lysimeters: A long-term
study using stable water isotopes. Vadose Zone J. 11 (1), 1–14. https://doi.org/10.
2136/vzj2011.0075.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. https://doi.org/
10.2136/sssaj1980.03615995004400050002x.
Wang, H.Y., Ju, X.T., Wei, Y.P., Li, B.G., Zhao, L.L., Hu, K.L., 2010. Simulation of bromide
and nitrate leaching under heavy rainfall and high-intensity irrigation rates in North
China Plain. Agr. Water Manage. 97, 1646–1654. https://doi.org/10.1016/j.agwat.
2010.05.022.
Wang, X.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.
Xiang, L., Yu, Z.B., Chen, L., Mon, J., Lu, H.S., 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.Y., Zeng, Y.J., Wen, J., Su, Z.B., 2018. Liquid-vapor-air flow in the frozen soil. J.
Geophys. Res. 123 (15), 7393–7415. https://doi.org/10.1029/2018JD028502.
Zeng, Y.J., Su, Z.B., Wan, L., Wen, J., 2011. A simulation analysis of the advective effecton evaporation using a two-phase heat and mass flow model. Water Resour. Res. 47(10), W10529. https://doi.org/10.1029/2011wr010701.
Zeng, Y.J., Su, Z.B., Wan, L., Yang, Z., Zhang, T., Tian, H., Shi, X., Wang, X.H., Cao, W.,
2009a. Diurnal pattern of the drying front in desert and its application for determining the effective infiltration. Hydrol. Earth Syst. Sci. 13, 703–714. https://doi.
org/10.5194/hess-13-703-2009.
Zeng, Y.J., Wan, L., Su, Z.B., Saito, H., Huang, K.L., Wang, X.S., 2009b. Diurnal soil water
dynamics in the shallow vadose zone (field site of China University of Geosciences,
China). Environ. Geol. 58 (1), 11–23. https://doi.org/10.1007/s00254-008-1485-8.
Zhang, M.L., Wen, Z., Xue, K., Chen, L.Z., Li, D.S., 2016a. A coupled model for liquidwater, water vapor and heat transport of saturated–unsaturated soil in cold regions:
Model formulation and verification. Environ. Earth Sci. 75, 701–719. https://doi.
org/10.1007/s12665-016-5499-3.
Zhang, S., Teng, J.D., He, Z.Y., Sheng, D.C., 2016b. Importance of vapor flow in unsaturated freezing soil: a numerical study. Cold Reg. Sci. Technol. 126, 1–9. https://doi.org/10.1016/j.coldregions.2016.02.011.
Zheng, C., Lu, Y.D., Liu, X.H., Šimůnek, J., Zeng, Y.J., Shi, C.C., Li, H.H., 2020. In-situ
monitoring and characteristic analysis of freezing-thawing cycles in a deep vadosezone. Water 12 (5), 1261–1278. https://doi.org/10.3390/w12051261.

16