Journal of Hydrology 598 (2021) 126250

Contents lists available at ScienceDirect

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

Research papers

Development of the Hydrus-1D freezing module and its application insimulating the coupled movement of water, vapor, and heat
Ce Zheng a, b, Jiří Šimůnek b, Ying Zhao c, Yudong Lu a, Xiuhua Liu a, d, *, Changchun Shi e,Huanhuan Li a, Lianyu Yu f, Yijian Zeng f, Zhongbo Su fa

School 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, China
Department of Environmental Sciences, University of California Riverside, Riverside 92521, USAc
College of Resources and Environmental Engineering, Ludong University, Yantai 264025, Chinad
Shaanxi Key Laboratory of Land Consolidation, Chang’an University, Xi’an 710064, ChinaeChina State Long-Term Observation and Research Station for Mu Us Desert Ecosystem in Yulin of Shaanxi, Yulin 719000, Chinaf
Faculty of Geo-information Science and Earth Observation (ITC), University of Twente, Enschede 7514 AE, Netherlandsb

A R T I C L E I N F O

A B S T R A C T

This manuscript was handled by Corrado Cor­
radini, Editor-in-Chief, with the assistance ofMohsen Sherif, Associate Editor

In cold regions, freeze-thaw cycles play a critical role in many engineering and agricultural applications andcause soil water flow and heat transport studies to be much more complicated due to phase changes involved. A
fully coupled numerical module for simulating the simultaneous movement of water, vapor, and heat duringfreezing-thawing periods was developed and incorporated in the Hydrus-1D software in this study. To avoid
numerical instabilities caused by a sudden increase in the apparent heat capacity during a phase change, a new
approach based on the available energy concept was adopted to adjust soil temperature when the freezing
temperature is reached. The proposed freezing module’s performance was then validated using experimental
data collected at three field sites with typical seasonal freezing/thawing processes. Results showed that themodel could efficiently obtain a convergent solution and that simulated soil moisture and temperature variations
captured the observed data well. Driven by soil matric potential and temperature gradients, both liquid water
and water vapor flowed towards the freezing front. The isothermal liquid flux was the most significant
component of overall flow in most soil depths except in the frozen layer, where it decreased by 1–5 orders ofmagnitude from values before freezing. Instead, the thermal vapor flux was the dominant moisture transfer
mechanism in the frozen layer and contributed about 10% to the ice formation. These results indicate that themodel, which considers the coupled movement of water, vapor, and heat, can better describe the physicalmechanisms of the hydrological cycle in the vadose zone during the freezing-thawing periods.

Keywords:
Freeze-thaw process
Unfrozen water content
Soil temperature
Vapor flow
Hydrus-1D

1. Introduction
The areas with frozen soil are widely distributed in middle and highlatitudes, occupying 55–60% of the land surface in the Northern
Hemisphere and approximately 72% of the Chinese territory (Wanget al., 2017; Zhang et al., 2012). Freezing-thawing cycles in these regions
substantially impact many engineering and agricultural applications,such as frost heave damage in high-speed railway embankments andmelting water infiltration (Kelleners et al., 2016; Kurylyk and Watanabe,
2013; Zhao et al., 2013). Meanwhile, phase changes during freezingthawing cycles significantly affect soil hydraulic and thermal proper­
ties, as well as soil water flow and heat transport (Lundin, 1990;

Watanabe and Wake, 2009; Yi et al., 2014; Zhang et al., 2012). For
example, due to the blocking effect of the ice present at subfreezing
temperatures, the isothermal hydraulic conductivity for liquid water isestimated to decrease by approximately 1 and 2 orders of magnitude
when about 15% and 20% of soil water is frozen, respectively (Zhaoet al., 2013). Therefore, the freezing-thawing process causes studies of
the coupled movement of soil water and heat to be more complex andcrucial for many environmental and engineering problems in cold
regions.
Vapor flow has been recognized to contribute to transfers of massand energy, especially in arid areas, such as northwestern China (Debet al., 2011; Zeng et al., 2009a). In these regions, soil moisture* 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.2021.126250
Received 10 December 2020; Received in revised form 13 March 2021; Accepted 19 March 2021
Available online 26 March 2021
0022-1694/© 2021 Elsevier B.V. All rights reserved.

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

movement takes place in the form of both liquid water and water vapor
flows, and the transformation between liquid water and water vapor
occurs continuously (Li et al., 2010; Sakai et al., 2009; Scanlon et al.,
2003; Zeng et al., 2009b; Zhao et al., 2013). Since phase changes be­
tween liquid water and water vapor could release or store a substantial
amount of latent energy, vapor flow has a considerable impact on thesoil energy balance (Cahill and Parlange, 1998). When soil is frozen,
vapor flow may play a dominant role in moisture transfer (Teng et al.,2018). Also, laboratory freezing experiments have shown that vapor
flow promotes ice formation, and its contribution increases withdecreasing water content before freezing (Yu et al., 2018).
The theory of the coupled water, vapor, and heat transport in porous
media under non-isothermal conditions was pioneered by Philip and deVries (1957). Since then, many researchers have been working during
the past several decades on advancing the understanding of these
coupled processes (Cass et al., 1984; Milly, 1984; Nassar and Horton,
1989; Saito et al., 2006). With increasing computational capacity andimproving simulation accuracy, underlying mechanisms of coupledwater and heat transfers under unfrozen soil conditions have attracted
growing interest and are widely reported (Deb et al., 2011; Saito et al.,2006; Sakai et al., 2009; Zeng et al., 2009b).
Multiple numerical models have been developed for simulating thesoil water and heat movement during freeze-thaw periods, and signifi­
cant progress has been achieved over recent decades in modeling theseprocesses (Flerchinger and Saxton, 1989; Guymon and Luthin, 1974;
Hansson et al., 2004; Jansson, 2012; Šimůnek et al., 2016; Yu et al.,2018, 2020). Due to different original concerns, these models differ inthe formulations of water flow and heat transport equations and areapplied to address specific issues (Bao et al., 2016; Kurylyk and Wata­nabe, 2013; Li et al., 2010). However, relative contributions of vapor
flow and its phase transition to the total water and heat balances areusually considered less significant in frozen soils and neglected (Li et al.,2010; Slater et al., 1998). This simplification may be practical for moist
soils but not for arid areas (Bao et al., 2016; Guymon and Luthin, 1974).
By neglecting the role of vapor flow, the physical mechanism behind soil
water flow and heat transport during the freezing period cannot be fully
understood (Flerchinger and Saxton, 1989; Hansson et al., 2004; Kur­
ylyk and Watanabe, 2013).
Meanwhile, existing studies of the coupled water, vapor, and heattransport during the freeze-thaw process contain some shortcomings
that need to be overcome. Most models were calibrated and validated
using data from laboratory freezing experiments or literature instead offield experiments (He et al., 2018; Teng et al., 2018, 2020). Usually,
actual field conditions are more complicated (e.g., time-variable atmo­
spheric boundary and multi-layered soil profiles), resulting in thesetypes of models being unsuitable for practical applications. Additionally,
the fixed freezing point approach, which converts the entire liquid watercontent into the ice content at the freezing point (e.g., 0℃), is adopted insome of these models to estimate the ice content (Li et al., 2010; Niu andYang, 2006). However, since some liquid water always exists in thefrozen soil no matter how low the temperature is, this approach repre­
sents a severe simplification of the soil freezing process (Sun et al., 2012;
Teng et al., 2020). While this approach may avoid numerical instabilities
caused by phase changes, it results in low precision and unrealistic
simulation results.
This study’s main objective is to develop a fully coupled numericalmodel for simulating the simultaneous movement of water, vapor, andheat in complicated field situations during freeze-thaw cycles that isboth physically realistic and computational stable and efficient. Since
the Hydrus-1D model (Šimůnek et al., 2008, 2016) has been widely
shown to perform well in evaluating the coupled movement of water,
vapor, and heat under nonfreezing conditions (e.g., Deb et al., 2011;
Saito et al., 2006; Sakai et al., 2009; Zeng et al., 2009a, 2009b), theproposed numerical model is incorporated into the Hydrus-1D code. The
previous version of the Hydrus-1D Freezing module was developed andsuccessfully used to evaluate laboratory column freezing experiments by

Hansson et al. (2004). The model was then also used by Zhao et al.
(2016) to simulate the variations in soil moisture and temperature in afield experiment. However, the current model becomes significantly lessstable and may encounter numerical problems in more complicatedsituations and neglects the effects of vapor flow. Thus, the specific ob­
jectives of this study are (i) to develop a rigorous and stable numericalfreezing module incorporated into the Hydrus-1D software for simu­
lating the coupled movement of water, vapor, and heat under freezing
conditions, (ii) to evaluate the model performance using in-situ obser­
vation data of soil water contents and temperatures, and (iii) to analyze
spatial and temporal variations in soil water contents, temperatures, andwater fluxes during the freezing-thawing process.
2. Model development
2.1. Mathematical description
Soil water flow in seasonally frozen soils with both above and below
freezing temperatures is described by the following mass conservation
equation (Hansson et al., 2004):∂θl ρv ∂θv ρi ∂θi∂++=−(q + qlT + qvh + qvT ) - S∂z lh∂t ρl ∂t ρl ∂t
(

=

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

)

(1)

-S

where θl, θv, and θi are the volumetric liquid 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), qlh and qlT are the isothermal andthermal liquid water fluxes (L T− 1), respectively, qvh and qvT are theisothermal and thermal vapor fluxes (L T− 1), respectively, Klh (L T− 1)
and KlT (L2 K− 1 T− 1) are the isothermal and thermal hydraulic conduc­tivities for liquid water, respectively, Kvh (L T− 1) and KvT (L2 K− 1 T− 1) are
the isothermal and thermal hydraulic conductivities for water vapor,respectively, h is the pressure head (L), T is the temperature (K), and S isa sink term (T− 1), which usually accounts for root water uptake. The first
three terms on the right-hand side of Eq. (1) represent liquid water flow
driven by the pressure head gradient, gravity, and temperature gradient,and the next two terms are vapor flow driven by the pressure head andtemperature gradients.
According to the pore-size distribution model of Mualem (1976), theisothermal hydraulic conductivity for liquid water can be describedusing the van Genuchten (1980) model as:
θs − θ rm[1 + |αh|n ]

(2)

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

θl = θr +where θr and θs are the residual and saturated water contents (L3 L− 3),
respectively, α (L− 1), n (− ), and m (− ) are empirical parameters, l (− ) isa pore connectivity parameter, which normally is set to 0.5 (Mualem,1976), Ks is the saturated hydraulic conductivity (L T− 1), and Se is theeffective liquid saturation (− ). Meanwhile, the resistance to water flow
in soil pores increases significantly due to the presence of ice. The
blocking effect gradually increases with an increase in the ice content
portion of the total water content. The isothermal hydraulic conduc­tivity for liquid water flow in freezing soil (Kflh, L T− 1) can then be
expressed as:
Kflh = 10− ψ Q ⋅Klh

(4)

where ψ is the impedance factor for liquid water flow due to the ice
formation (− ), which usually ranges from 0.5 to 6.5 depending on soil
properties (Zhao et al., 2016), and Q is the mass ratio of the ice content
2

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

T− 1), respectively. The components of soil heat fluxes, on the right-handside of Eq. (9), include heat flow by conduction, convection of sensible
heat with liquid water flow, convection of sensible heat with vapor flow,
and convection of latent heat with vapor flow.
The apparent soil thermal conductivity, which combines the thermal
conductivity of the porous medium in the absence of flow with thethermal dispersivity, can be expressed for unfrozen conditions as follows(Hansson et al., 2004):

to the total water content (− ).
The thermal hydraulic conductivity for liquid water flow andisothermal and thermal hydraulic conductivities for vapor flow are
defined as follows (Saito et al., 2006):()
1 dγ
KlT = Klh hGwT
(5)
γ0 dT
Kvh =

KvT =

D

ρw

Mg
RT

(6)

dρsv
dT

(7)

Hr ρsv

D

ρw

ηH r

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

where GwT is the gain factor (–), which corrects the temperaturedependence of the soil water retention curve and is defined as 7 for sands(Sakai et al., 2009), γ0 is the surface tension at 25 ◦ C (M T− 2, =71.89 gs− 2), γ is the surface tension of soil water (M T− 2), Hr is the relative
humidity (–), ρsv is the saturated vapor density (M L− 3), M is the mo­lecular weight of water (M mol− 1, =0.018 kg mol− 1), g is the gravita­
tional acceleration (L T− 2), R is the universal gas constant (J mol− 1 K− 1,=8.315 J mol− 1 K− 1), η is the enhancement factor (–), and D is the vapor
diffusivity in soil (L2 T− 1), which can be calculated as follows (Campbell,1985):λ0 (θl ) = b1 + b2 θl + b3 θ0.5l

λ0 (θl ) = b1 + b2 (θl + Fθi ) + b3 (θl + Fθi )0.5

where Cp is the volumetric heat capacity of moist soil (M L− 1 T− 2 K− 1),which is defined as the sum of the volumetric heat capacities of the solid(Cn, =1.9 MJ m− 3 K− 1), liquid water (Cw, =4.2 MJ m− 3 K− 1), water
vapor (Cv, =6.3 KJ m− 3 K− 1), and ice (Ci, =2.1 MJ m− 3 K− 1) phases
multiplied by their respective volumetric fractions (Campbell, 1985), Lf
is the latent heat of freezing (L2 T− 2, ≈3.34 × 105 J kg− 1), Lw is the latent
heat of vaporization of liquid water (L2 T− 2), λ is the apparent thermal
conductivity (M L T− 3 K− 1), and ql and qv are liquid and vapor fluxes (Lh=θl = θr + [Equation

References

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

Hillel (1971)

Da = 2.12 × 10− 5 (T/273.15)2

Campbell(1985)2τ = θ7/3
a /θs

Saturated vapour
density ρsv, (ML− 3 )
Relative humidity
Hr, (–)
Enhancement
factor η, (–)ρsv = exp(31.37 − 6014.79/T − 7.92 ×Millington
and Quirk(1961)
Hillel (1971)

Latent heat of water
vaporization Lw,(L2 T− 2)

10− 3 T)/T × 10− 3

Hr = exp(Mgh/RT)

η = 9.5 +3θl /θs − 8.5exp(− ((1 + 2.6/

Lf T
ln
g Tf

√̅̅̅̅
fc )θl /θs )4 )Lw = 2.501 × 106 − 2369.2T(13)

where Tf is the freezing point of liquid water (K), which varies with soil
texture and the initial water content. This indicates that the liquid watercontent is a function of temperature and can be obtained by the soil
water freezing curve. When soil is frozen, substituting Eq. (13) into Eq.
(2) leads to the relationship between the liquid water content andtemperature (the soil freezing characteristics curve) as follows:

Table 1
Parameters used in the mathematical model.
Surface tension ofsoil water γ, (MT− 2)Vapor diffusivity inthe air Da, (L2
T− 1)
Tortuosity factor τ,
(–)

(12)

where F accounts for the difference between the thermal conductivities
of liquid water and in soils (–) and can be defined as suggested by
Hansson et al. (2004).
For unfrozen soils, the liquid water content is a function of thepressure head. It can be estimated from the soil water characteristics
curve as a function of only two unknown parameters (h and T). Due toanother unknown parameter (θi) under freezing conditions, one addi­
tional equation is required to solve this new parameter (Kurylyk and
Watanabe, 2013; Sun et al., 2012). Based on the thermal dynamic
equilibrium theory and the assumption that the ice pressure is equal tozero when soil is frozen, the soil matric potential can be estimated forsubfreezing temperatures using the generalized Clapeyron equation as
follows (Hansson et al., 2004):

where θa is the air-filled porosity (L3 L− 3), Da is the vapor diffusivity inthe air (L2 T− 1), and τ is the tortuosity factor (–). Detailed descriptions ofrelated parameters are presented in Table 1.
The governing equation for heat transport in a variably-saturated
rigid soil is given by the following energy conservation equation (Nas­sar and Horton, 1989):
(
)
∂Cp T
∂θi
∂θv ∂
∂T
∂ql T
∂qv T
∂qv
− Lf ρi + Lw ρw
=
λ
− Cw
− Cv
− L0
− Cw ST
∂z
∂t
∂t
∂t
∂z
∂z
∂z
∂z
(9)

Parameter(11)

where b1, b2, and b3 are empirical parameters that depend on soil texture(M L T− 3 K− 1) and equal to 0.228, − 2.406, and 4.909 W m− 1 K− 1 forsands, respectively. Since ice can also influence heat conduction, thethermal conductivity of freezing soils can be calculated as follows:(8)

D = Da θ a τ(10)θ − θr⃒ (s)⃒ ]⃒ Lf T ⃒n m
1 + ⃒⃒α g lnTf ⃒⃒(14)

The mathematical model of the coupled movement of water, vapor,and energy in soils under freezing-thawing conditions described above
has surprisingly few additional parameters in addition to those requiredby the standard water flow and heat transport models (i.e., soil hydraulic(retention and conductivity) and thermal (capacity and conductivity)
parameters. The values of most additional parameters (e.g., latent heatof freezing and vaporization, vapor density, vapor diffusivity) and theirtemperature dependences are known and can be found in the literature(see Table 2). The two main additional new parameters are the imped­
ance factor ψ in Eq. (4) and the F parameter in Eq. (12), which affect thehydraulic and thermal conductivities, respectively, in frozen soils. Thus,
the coupled model is still sensitive mainly to the soil hydraulic andthermal parameters as in the standard water flow and heat transportmodels.

Philip and deVries (1957)
Cass et al.
(1984)
Sakai et al.
(2009)

fc: mass fraction of clay in the soil (–).
3

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

Table 2
Measured soil particle composition and soil hydraulic parameters at the Yu Lin Desert (YD), Inner Mongolia Grassland (IG), and Tibetan Plateau Maqu Grassland (MG)sites.
Study site

Soil layer (cm)Soil particle composition
Sand (%)
Silt (%)

Clay (%)

Soil hydraulic parametersθs (cm3 cm− 3)
θr (cm3 cm− 3)

α (cm− 1)n (–)

Ks (cm d− 1)

YD

0–80
80–160
160–200
0–15
15–25
25–35
35–100
0–10
10–40
40–160

95.8
90.1
91.7
62.0
73.0
78.1
78.7
44.1
44.3
65.6

1.0
4.5
3.1
15.8
12.4
10.9
9.7
9.0
10.1
5.6

0.008
0.015
0.010
0.042
0.045
0.042
0.069
0.035
0.039
0.045

0.029
0.039
0.037
0.026
0.025
0.021
0.019
0.041
0.042
0.073

1.64
1.51
1.59
1.73
2.05
1.75
1.93
1.33
1.36
1.58

1368
811
1440
165
133
113
72
12.5
8.1
5.9

IG

MG

3.2
5.4
5.2
22.2
14.6
11.0
11.6
46.9
45.6
28.8

2.2. Calculation procedure

0.375
0.427
0.415
0.572
0.472
0.453
0.452
0.5
0.45
0.41

a linear equation system, i.e., the matrix equation. After incorporating
the boundary conditions, the Gaussian elimination method is used tosolve the matrix equation and to obtain the new pressure head (and thecorresponding liquid water content) at the new time level. The Picard
iteration procedure is used to solve this highly nonlinear equation. The
iterative process continues until absolute changes in pressure heads andsoil water contents in all nodes are less than the prescribed tolerance(Kelleners, 2013; Šimůnek et al., 2008).
Next, the Galerkin finite element and finite difference methods are

A flexible time-stepping method is used in the Hydrus-1D model for
calculating the coupled movement of liquid water, water vapor, andheat. The main computational procedure at each timestep is shown in
Fig. 1 and is described below in detail.
First, the soil water flow equation is solved using the finite difference
method for both spatial and temporal discretizations (Šimůnek et al.,2008). The modified Richards equation is discretized and assembled into

Fig. 1. Computational flowchart of the Hydrus-1D Freezing module during one timestep. Variables are listed in the upper right corner of the figure. The superscript k+ 1 indicates the current iteration in this timestep. The procedure loops twice between water flow and heat transport routines within one timestep (iLoop).
4

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

applied for calculating the spatial and temporal derivatives for the heattransport equation, respectively (Šimůnek et al., 2008). Nodal values ofsoil water contents and fluxes calculated by solving the water flow
equation are used as known values in the heat transport equation. A
similar solution strategy is applied for calculating the new soil temper­
atures until achieving convergency. The steps mentioned above arelooped twice at each timestep (i.e., iLoop = 1, 2). The main reason forthe second loop is to update the solution of the water flow equation
using updated values of soil temperatures and ice contents (if ice exists).
These are obtained by the heat transport equation solution, which re­flects the effects of the freezing process and temperature changes on soil
water flow within the current timestep.
Note that there are two unknowns when solving each governing
equation: the liquid water content θl and the pressure head h in the waterflow equation and the ice content θi and the temperature T in the heattransport equation. A mass-conservative method proposed by Celia et al.
(1990) is used in the time difference scheme to eliminate θl and θi fromthe two equations, as shown below (Hansson et al., 2004):∂θl θj+1,k+1− θjl θj+1,k+1− θj+1,kθj+1,k − θjl
l
≈ l
= l
+ l
∂t
Δt
Δt
Δt
= Cj+1,kLf ρi

hj+1,k+1 − hj+1,k θj+1,k
− θjl
+ l
Δt
Δt

∂θi
θj+1,k+1 − θjiθj+1,k+1 − θj+1,kθj+1,k − θjii≈ L f ρi i
= Lf ρi i+ Lf ρi i
∂t
Δt
Δt
Δt
=−

(
)j+1,kLf2 ρiT j+1,k+1 − T j+1,kθj+1,k − θji⋅C+ Lf ρi igT
Δt
Δt

2.4. Model stability
The apparent heat capacity (Ca, M L− 1 T− 2 K− 1), which includes thelatent heat as an additional term in the heat capacity, can be calculatedusing the generalized Clapeyron equation and the soil water capacity as
follows:
Ca = Cp − Lf ρi

L2f ρidθi= Cp +
C
dT
gT

(18)

In fully saturated soils, the apparent heat capacity increases by about
2 and 5 orders of magnitude for silty clay and sand, respectively, when
the freezing point (0 ◦ C) is reached. The freezing point for unsaturated
soils is lower than 0 ◦ C and can be estimated as follows:( )gh
Tf = T0 ⋅exp
(19)
− T0
Lf
where T0 (273.15 ◦ C) is the absolute temperature. When the freezing
point is reached, the significant difference in the apparent heat capacity
for temperatures above and below the freezing point can lead to aconsiderable underestimation of the soil temperature. This may result in
oscillations in calculated temperatures between negative and positive
values within a timestep (Hansson et al., 2004; Wang et al., 2017). A
method based on available energy for a phase change is adopted in our
simulations to overcome this instability caused by a phase change be­
tween liquid and solid states of water. The new temperature at the first
iteration when a soil layer crosses the freezing point is adjusted as fol­lows (Flerchinger and Saxton, 1989):

(15)

(16)

Tr = Tf −

Ea
Ca

(20)

where superscripts j + 1 and j indicate the current and previous timelevels, respectively, k + 1 and k indicate the current and previous iter­
ations for the current time level, respectively, C is the soil water capacity(L− 1), which is defined as the derivative of the water content withrespect to the pressure head, and Δt is the time interval between theprevious and current time levels.

where Tr is the revised temperature, Ea is available energy for freezing.
This adjustment method was found to generate more consistent calcu­
lation results and avoid undesired numerical oscillations.

2.3. Determination of ice content

3. Materials and methods

When estimating the ice content in the iteration procedure, large
errors often appear in calculations of soil water contents and tempera­
tures, and consequently, the numerical convergency is challenging toachieve (Bao et al., 2016; Li et al., 2010). Therefore, in the proposedmodel, the ice content is determined after calculating the new temper­
ature, as shown in Fig. 1. For subfreezing temperatures, liquid water andice coexist in soil pores. Based on the principle of mass conservation, thetotal water content (θt, L3 L− 3) can be defined as the sum of the liquidwater content and the ice content:

3.1. In-situ experimentsρρlθ t = θl + i θi()
Ea = Cp ⋅ Tf − T(21)Experimental data were collected at three field sites: a) the Yu Lin
Desert Ecosystem Research Station (the YD site), b) the Inner Mongolia
Grassland Ecosystem Research Station (the IG site), and c) the Tibetan
Plateau Maqu Grassland (MG) site.
The first in-situ experiment was located at the Yu Lin Desert
Ecosystem Research Station (the YD site), which is a typical area withseasonal freezing in northwestern China (109◦ 41′ 44′′ E and 38◦ 21′ 46′′ N,
Fig. 2). The climate in this region can be characterized as cold and dry(Cheng et al., 2013; Zheng et al., 2020). Observations from the nearby
Yulin meteorological station indicate that the long-term (since 1961)
average air temperature is approximately 6.4 ◦ C, and the minimum airtemperature reaches − 31.4 ◦ C. The average annual rainfall and poten­
tial evaporation are 360 and 2343 mm, respectively.
The soil surface is flat at the study site, and vegetation covers lessthan 5% of it during the freezing period. To determine soil physicalproperties (as listed in Table 2), soil samples were collected at a 20 cm
interval from the surface to the 200 cm depth. Soil water contents andtemperatures were measured at depths of 10, 20, 50, 100, and 200 cm by
the Hydro Probe II sensors (Stevens Inc., USA) every 10 min. Before
being installed, all sensors were calibrated using the original soil sam­ples. The measurement principle of the Hydro Probe II sensors is basedon the frequency domain method (Cheng et al., 2014). Because thepermittivity of liquid water (≈78) is significantly higher than those of
the soil matrix (≈3–5), ice (≈3.4), and air (≈1), observed soil watercontents could be considered to represent unfrozen water contents when

(17)

Then, the liquid water content is updated based on the soil water
freezing curve using the newly calculated soil temperature, and theexistence of ice is determined using the θl < θt criterion. Note that there
exists no ice when the soil temperature is in the interval between thefreezing temperature and zero temperature (Tf < T < 0), as explained
below in detail. There is no ice for the above-zero temperature condi­
tions. The liquid water content needs to be updated at the first iteration
when the soil temperature turns positive from a negative value.
Compared with the previous version of the Freezing module (Hans­
son et al., 2004), differences in liquid water and ice densities are
considered in the updated module. The inclusion of density differences
was expected to decrease simulation errors (Kelleners, 2013).

5

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

Fig. 2. The locations of the Yu Lin Desert (YD) (and its field instrumentation layout in photos on the right), Inner Mongolia Grassland (IG), and Tibetan Plateau Maqu
Grassland (MG) sites. Adapted from Zheng et al. (2020).

the soil is frozen (Sun et al., 2012; Watanabe and Wake, 2009). Micro­
meteorological variables, including air temperature, wind speed, rela­tive humidity, and precipitation, were recorded at a 30 min interval. As
shown in Fig. 3, the total amounts of rainfall and snow were 20.8 and13.6 mm in winters of 2017–2018 and 2018–2019, respectively, mainlyin the form of rains after soil thawing.
The proposed model was further applied to simulate the coupled
movement of water, vapor, and heat at two grassland sites (Fig. 2) toassess the model applicability. The Inner Mongolia Grassland (IG) site is
located at the Inner Mongolia Grassland Ecosystem Research Station(116◦ 42′ 18′′ E and 43◦ 37′ 50′′ N), and this site has a continental climate
with cold winter and warm summer (Zhao et al., 2016). On the other
hand, the Maqu Grassland (MG) site is set up on the north-eastern fringeof the Tibetan Plateau (101◦ 38′ -102◦ 45′ E and 33◦ 30′ -34◦ 15′ N), with theelevation ranging from 3200 to 4200 m above the sea level (Yu et al.,2018). The climate is characterized as wet and cold, with dry winter andrainy summer. Table 3 lists detailed information about the two sites andfield observations.

et al., 2008; Kroener et al., 2014). For heat transport, the calculated
surface temperature and observed temperature at the 200 cm depth
were used as the upper and lower boundary conditions, respectively.
The effects of root water uptake on both water flow and heat transport
were neglected, considering no plant uptake in the winter. The mean
absolute error (MAE), the root mean square error (RMSE), and the co­efficient of determination (R2) were computed to evaluate the model
performance.
To further validate the developed model, two additional datasetscollected at the IG site (the Inner Mongolia Grassland site; from October
14, 2005, to August 31, 2006) and at the MG site (the Tibetan Plateau
Maqu Grassland site; from December 2, 2015, to March 15, 2016) wereused for simulation purposes. Value-specified initial and boundary
conditions were used in the model setup (Yu et al., 2018; Zhao et al.,2016).
4. Results
The Results section is organized as follows. First, model calibrationand validation using measured water contents and temperatures at thethree experimental sites (YD, IG, and MG) are discussed. In the subse­
quent two sections, the more detailed analysis of the modeling results, i.
e., diurnal variations in soil water contents, temperatures, and water
fluxes, is carried out only for the YD site. These sections also include
discussions on the relative contributions of different water fluxes (i.e.,water, vapor, thermal, and isothermal).

3.2. Numerical simulation
The dataset collected at the YD site from November 1, 2017, to
March 31, 2018, was used to calibrate the Hydrus-1D Freezing module.
Soil parameters were initially estimated using the neural network
method based on soil physical properties and then optimized using theLevenberg-Marquardt algorithm, which is directly implemented into theHydrus code (Šimůnek et al., 2008, 2016; Zhao et al., 2013). The cali­brated soil hydraulic parameters are presented in Table 2. Subsequently,
the measured soil data from November 1, 2018, to March 31, 2019, wereused for model validation. The soil profile was 200 cm deep, and a
spatial discretization of 1 cm was applied. Observed soil water content
and temperature data on November 1, 2017, and November 1, 2018,
were used as initial conditions for the two periods, respectively, andvalues between two measurement depths were obtained using linear
interpolation. For water flow, the atmospheric boundary condition was
imposed on the soil surface. The free drainage condition was assumed at
the bottom boundary because the groundwater table fluctuated deep
below the bottom of the soil profile around a depth of 900 cm (Bittelli

4.1. Model calibration and validation
Fig. 3 shows the comparison between measured and simulated un­
frozen water contents and soil temperatures at the YD site during theanalyzed period. The major changes in the unfrozen water content at the
10 and 100 cm depths were distinctive due to soil freeze-thaw cycles (e.
g., varying from 0.074 to 0.009 cm3 cm− 3 at the 10 cm depth during thecalibration period), and the Hydrus-1D Freezing module well captured
these variations. Deviations between simulated and measured values
were lowest in the shallow layer. For instance, the values of MAE and
RMSE for the calibration period were 0.003 and 0.004 cm3 cm− 3 at the6

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

Fig. 3. Comparisons of daily measured and simulated unfrozen water contents (left) and soil temperatures (right) at the Yu Lin Desert (YD) site during the periods:(a) from November 1, 2017, to March 31, 2018, and (b) from November 1, 2018, to March 31, 2019. The MAE and RMSE units are cm3 cm− 3 and ℃ for watercontents and temperatures, respectively.

7

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

periods (the first one for the freezing process from November 30,
2017, to December 2, 2017; the second one for the thawing process fromMarch 6, 2018, to March 8, 2018) were selected to discuss in detailvariations in soil temperatures and water contents. Fig. 6 shows thehourly observed and simulated unfrozen water contents and tempera­
tures at 10, 100, and 200 cm depths during the selected periods. The
average value of RMSE for soil water contents and temperatures are
0.003 cm3 cm− 3 and 0.53 ℃, respectively, which again indicates thatsimulations fitted the measured data well. Generally, diurnal variationsin the soil water content were small, while visible diurnal changes
appeared at the 10 cm depth during the freezing period. This behavior
was caused by diurnal soil freeze-thaw cycles when the unfrozen watercontent in the topsoil layer increased during daytime and decreased atnight.
The simulated spatial–temporal distributions of the unfrozen watercontent and temperature during the selected two periods are shown in
Fig. 7. During the selected freezing period, the shallow layer’s unfrozen
water content was significantly lower than in the deep layer, indicatingthat the soil was frozen. The freezing front dropped from 19 to 28 cm
during this period (which can also be indicated by the zero-temperatureline in the temperature profile). Below the freezing depth, the soil watercontent was decreasing slowly (e.g., from 0.068 to 0.063 cm3 cm− 3 atthe 30 cm depth). This decreasing trend was mainly attributed to up­
ward flow due to the soil matric potential gradient, which caused soil
water to flow from deeper layers towards the frozen layer. In contrast,
soil temperature in the shallow layer (about 40 cm in Fig. 7) usuallyexhibited apparent diurnal variations controlled by solar radiation.
These diurnal variations resulted in the opposite day/night behavior ofthermal liquid water and water vapor flow. Below the 40 cm depth, soil
temperature showed no distinct fluctuations and mainly increased withdepth. For example, the average temperature at depths of 50, 100, 150,
and 200 cm were 1.9, 5.1, 7.8, and 10.3 ℃, respectively. Like the matricpotential gradient, the upward thermal gradient caused soil water flow
towards the frozen layer.
As shown in Fig. 7, the thawing process (as indicated by the zerotemperature line) occurred both at the top and the bottom of thefrozen layer. The downward movement of the upper thawing front (from42 to 54 cm) was faster than the upward movement of the lower thawing
front (from 98 to 93 cm) during the selected thawing period. This phe­nomenon suggested that a significant increase in air temperature inspring had a remarkable effect on the shallow layer’s soil thawing pro­cess, which was also indicated by Cheng et al. (2014). Driven by the soil
matric potential gradient, soil water flowed from both above and below
towards the frozen layer. Continuous upward flow during the entire
freezing-thawing process caused a decreasing trend in the unfrozen
water content in the deep soil (e.g., about 0.114 and 0.099 cm3 cm− 3 at
the 150 cm depth during selected freezing and thawing periods,
respectively). On the other hand, intense solar radiation caused a larger
downward temperature gradient in the topsoil during the daytime. In
comparison, the soil temperature gradient below the freezing depth was
smaller than before.
Fig. 8 shows the simulated vertical profiles of unfrozen water con­tent, ice content, total water content, and temperature at different times
during the freezing period. Comparing the results at 0 and 12 h, there
existed apparent diurnal changes in the topsoil layer. For example, theunfrozen water content at the 2 cm depth increased from 0.012 to 0.046
cm3 cm− 3 during this time interval, while the ice content dropped tozero due to the positive temperature. Below the layer experiencing
freeze-thaw cycles, the ice content gradually increased (e.g., 0.046,
0.061, and 0.065 cm3 cm− 3 at 0, 12, and 72 h at the 15 cm depth,
respectively). Apart from decreasing temperature and corresponding
freezing of unfrozen water, vapor-ice desublimation was another reason
for the ice accumulation (He et al., 2018; Zhang et al., 2016). Mean­
while, the total water content in the frozen layer became larger during
these three days (e.g., 0.066 and 0.075 cm3 cm− 3 at 0 and 72 h at the 24
cm depth, respectively). This data indicated that soil water flowed from

Table 3
Detailed information about the two grassland sites. IG: the Inner Mongolia
Grassland site, MG: the Tibetan Plateau Maqu Grassland site.
Characteristics

IG

MG

Soil texture
Vegetation
Mean annual air
temperature
Groundwaterlevel
Measurement
depth
Soil moistureobservation

Sandy loam
Perennial rhizome grass
0.7 ◦ C

Sandy loam
Alpine meadow
1.3 ◦ C

Below 10 m

8.5–12 m

100 cm160 cmML2x probes (Delta-T
Devices Ltd., UK) at 5, 20,
and 40 cm depths
Pt-100 probes (Omega, UK)
at 2, 5, 20, 40, and 100 cm
depthsMicrometeorologicalobserving system

5TM ECH2O probes (METERGroup, Inc., USA) at 5, 10, 20,
40, 80, and 160 cm depths
5TM ECH2O probes (METERGroup, Inc., USA) at 5, 10, 20,
40, 80, and 160 cm depths
Micrometeorological observing
systemSoil temperatureobservation
Weathervariableobservation

10 cm depth, while they were 0.005 and 0.009 cm3 cm− 3 at the 100 cm
depth, respectively. Compared with the measured data (about 0.055 cm3
cm− 3), the unfrozen water content at the 100 cm depth was under­
estimated by the modeling results (about 0.041 cm3 cm− 3) when the soil
was frozen. This discrepancy was likely caused by a sudden change insoil hydraulic properties above and below the 80 cm depth due todifferent soil textures (Table 2) and the model representation of thephysical process involved in the soil freezing process. Such discrepancy
was also reported by Yu et al. (2020) and Zhao et al. (2016). During thesimulation period, small rains mainly occurred after soil thawing, andthey affected only the topsoil layer. Simulated soil water contents
described well sudden increases due to rainfall events.
During the analyzed period, soil temperatures first decreased be­
tween November and January and then increased in February andMarch. Simulated temperatures at all depths fitted very well measured
data, especially during the calibration period. Correlations between
simulated and measured temperatures were better than for unfrozen
water contents. The average values of MAE and RMSE were 0.48 and0.57 ℃, respectively. These two indices were not calculated for the 200cm depth because the measured temperature was used as the bottom
boundary condition.
The comparisons between observed and simulated unfrozen watercontents and temperatures at the IG and MG sites are shown in Fig. 4.
Due to frequent rainfalls and intense evaporation at the IG site, de­viations between observed and simulated unfrozen water contents at the
5 cm depth were larger than at deeper depths. The errors were slightly
larger than those calculated at the YD site, which may be partly attrib­
uted to the variations in measured unfrozen water contents beingslightly delayed when the phase change occurred. There is a good
agreement between measured and simulated values of both unfrozen
water contents and temperatures at the MG site. For example, theaverage RMSE values for unfrozen water contents and temperatures
were 0.005 cm3 cm− 3 and 0.19 ℃, respectively, which was the lowest
among the three sites.
Fig. 5 shows the scatter plots of the observed and simulated unfrozen
water contents and soil temperatures at the three experimental sites. The
average R2 values, which are used to assess the correspondence between
the observed and simulated unfrozen water contents and temperaturesat the three sites, were 0.92 and 0.98, respectively, reflecting a highcorrelation. In general, the simulated results indicate that the proposedmodel has good applicability in different cold regions. Possible reasonsfor any differences at the three sites are discussed below.
4.2. Diurnal variations in soil moisture and temperature
Based on the simulated results at the YD site, two typical three-day8

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

Fig. 4. Comparison of daily measured and simulated unfrozen water contents (left) and soil temperatures (right) at (a) the Inner Mongolia Grassland (IG) site from
October 14, 2005, to August 31, 2006, and (b) the Tibetan Plateau Maqu Grassland (MG) site from December 2, 2015, to March 15, 2016. The MAE and RMSE units
are cm3 cm− 3 and ℃ for water contents and temperatures, respectively.

9C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

Fig. 5. Scatter plots of the observed and simulated unfrozen water contents (top) and soil temperatures (bottom) at the Yu Lin Desert (YD) site (left), the Inner
Mongolia Grassland (IG) site (middle), and the Tibetan Plateau Maqu Grassland (MG) site (right).

the deep soil layers towards the frozen layer and increased the soil watercontent there, as discussed below.
4.3. Soil water fluxes
Fig. 9 depicts the simulated vertical profiles of soil water fluxes (qlh:
isothermal liquid flux; qlT: thermal liquid flux; qvh: isothermal vapor flux;
qvT: thermal vapor flux) at 0, 12, and 72 h during selected freezing andthawing periods. As indicated by Fig. 9, both liquid and vapor fluxes
were directed towards the frozen layer. The qlh flux dominated waterflow in most depths except in the frozen layer, where it became almost
negligible. In addition to the qlh flux, the effects of the upward qlT and qvTfluxes towards the frozen layer were also observable, indicating thatboth thermal fluxes contributed to the ice formation. For example, theupward qvT fluxes were 6.2e-04, 7.9e-04, and 7.1e-04 cm d− 1 around thefreezing front at 0, 12, and 72 h during the freezing process, representing
about 8.2, 9.6, and 11.5% of the total water flux, respectively.
Fig. 10 shows the diurnal variations of each water flux at different
depths during the freezing and thawing periods. Table 4 then lists con­tributions of each soil water (isothermal and thermal, liquid and vapor)flux to the total water flux during the two selected periods. These ratios
were calculated using the hourly data as follow:∑N
1 |qi |
P = ∑N
(22)
(|q
|
+
|q
lh
lT | + |qvh | + |qvT | )1
where P is the contribution, qi represents each soil water flux (qlh, qlT,
qvh, and qvT), and N is the number of hours during a selected time in­terval (=72 h).
Due to the day/night air temperature changes, there were also
noticeable diurnal variations of water fluxes in the shallow layer. For
example, at the 2 cm depth, the qlh and qlT fluxes were only noticeable atnoon, directed upward and downward due to different driving forces,
respectively. Although the liquid water flux became smaller than thevapor flux during most of the day, the qlh flux still contributed most to
the total flux (about 64%) since it was relatively large (with theFig. 6. Hourly observed and simulated unfrozen water contents (top) and soil
temperatures (bottom) at 10, 100, and 200 cm depths at the Yu Lin Desert (YD)site during the selected freezing (from November 30, 2017, to December 2,2017) and thawing (from March 6, 2018, to March 8, 2018) periods. The RMSE
units are cm3 cm− 3 and ℃ for soil water contents and temperatures, respec­tively. (Ob: Observed value, Si: Simulated value).
10

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

Fig. 7. Simulated hourly variations in unfrozen water content (top) and soil temperature (bottom) profiles at the Yu Lin Desert (YD) site during selected freezing(left, from November 30, 2017, to December 2, 2017) and thawing (right, from March 6, 2018, to March 8, 2018) periods, respectively.

Fig. 8. Vertical profiles of (a) unfrozen water content, (b) ice content, (c) total water content, and (d) soil temperature at 0, 12, 72 h of the selected freezing period(from November 30, 2017, to December 2, 2017) at the Yu Lin Desert (YD) site.

maximum flux of 0.031 cm d− 1) at noon when the soil thawed. Being lessaffected by the daily freeze-thaw cycles, the qvT flux contributed 92% to
the total flux at the 10 cm depth during the freezing period and about

78% at the 80 cm depth during the thawing period, indicating theimportant blocking effect of liquid flow in the frozen layer. The signif­
icant decrease of liquid fluxes (about 1–5 orders of magnitude) was
11

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

Fig. 9. Vertical profiles of soil water fluxes at 0 (left), 12 (middle), and 72 h (right) during selected freezing (top, from November 30, 2017, to December 2, 2017) andthawing (bottom, from March 6, 2018, to March 8, 2018) periods at the Yu Lin Desert (YD) site. The dotted lines represent the locations of freezing/thawing fronts.

attributed to ice, which significantly reduced the soil hydraulic con­ductivity to liquid water. Besides, the qvh flux was negligible at mostdepths except in the shallow frozen layer (representing about 3% of thetotal flux), where the matric potential gradients were relatively large.
Below the freezing front, the qvT flux was smaller due to a decreasing
temperature gradient and usually contributed less than 10% of the total
water flux. In contrast, the qlh flux was the most significant component inmost depths, followed by the qlT flux. Note that the qlT flux (about 2.4e03 cm d− 1) was larger than the qlh flux (about 1.8e-03 cm d− 1) at the
100 cm depth during the freezing period, accounting for about 50% of
the total water flux.

dominated the moisture transfer in the frozen layer due to a sharp
decrease in the liquid flow (e.g., contributing over 90% of the total
moisture flux at the 10 cm depth, as indicated by Table 4). Note that thevapor flux’s dominant contribution to the total soil water flux occurred
only during relatively short periods. For the entire freeze-thaw period(from November to March), the vapor flux contributed only about 7.8%of the total soil water flux at the 10 cm depth despite being larger thanthe liquid water flux most times when the soil was frozen. The overall
total liquid flux was still much larger, mainly due to the higher flux
during the thawing process (including diurnal freeze-thaw cycles andspring thawing). Below the freezing front, vapor flow can still affect thesoil moisture distribution due to its persistence (He et al., 2018; Tenget al., 2018; Yu et al., 2018; Zheng et al., 2020). For example, about 10%of the ice formation in our simulations was due to the upward vapor flux.
Similar results were calculated based on laboratory freezing experiments
by Zhang et al. (2016).
On the other hand, the isothermal vapor flux, which is often
neglected when evaluating soil water flow during nonfreezing periods,
was observable in the shallow frozen layer. The maximum isothermal
vapor flux represented over 15% of the total flux (at the 2 cm depth)
when the soil was frozen and not experiencing freeze-thaw cycles.
However, our simulation results are not consistent with the results re­ported by Yu et al. (2018), who assumed that the isothermal vapor flux
played a dominant role in the total mass flux in the topsoil layers (at the

5. Discussion
5.1. Importance of vapor flow
It has been increasingly recognized that vapor flow has a significant
impact on processes that are important for many agricultural andecological applications, especially in arid areas with large temperature
gradients and high air-filled soil porosity (Bittelli et al., 2008; Zeng et al.,2009b). Although the vapor flux magnitude is relatively low duringfreezing periods due to relatively small temperature gradients, ther­mally driven vapor flow still plays a crucial role in overall soil water
flow. When the soil was completely frozen, thermal vapor flow
12

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

Fig. 10. Temporal variations in soil water fluxes at different depths during selected freezing (top, from November 30, 2017, to December 2, 2017) and thawing(bottom, from March 6, 2018, to March 8, 2018) periods at the Yu Lin Desert (YD) site. Different depths for the two periods were selected based on the location of thefrozen layer. Note different scales on vertical axes.

0.1 cm depth) when the soil was frozen. This discrepancy may be
attributed to different spatial discretizations used in respective numer­
ical models since the soil matric potential gradient varies significantly
over small distances in the top 0.1 cm soil layer. Nevertheless, theseresults indicated that the isothermal vapor flux should be considered instudies of frozen soils. Therefore, in cold and arid regions, only themodel that considers the coupled movement of water, vapor, and heatcan adequately describe the physical mechanisms of the hydrologicalcycle in the vadose zone during freezing-thawing periods.

Table 4Contributions of each soil water flux to the total water flux during selected threeday freezing (from November 30, 2017, to December 2, 2017) and thawing(from March 6, 2018, to March 8, 2018) periods at different depths at the Yu Lin
Desert (YD) site. (qlh and qlT are the isothermal and thermal liquid water fluxes,
respectively, and qvh and qvT are the isothermal and thermal vapor fluxes,
respectively).
Depth
(cm)

Freezing period

Thawing period

qlh
(%)

qlT
(%)

qvh
(%)

qvT
(%)

qlh
(%)

qlT
(%)

qvh
(%)

qvT
(%)

2
10
20
50
80
100
200

64.3
3.8
76.5
69.3
56.9
37.5
79.2

7.1
0.5
9.9
26.1
35.3
49.7
17.2

0.6
3.8
0.3
0
0
0
0

27.9
91.9
13.3
4.6
7.8
12.8
3.6

94.2
94.2
93.7
96.2
16.2
88.2
77.3

0.9
3.5
4.8
3.1
4.4
8.1
16.8

0.2
0
0
0
1.1
0
0

4.7
2.3
1.5
0.7
78.3
3.7
5.9

5.2. Influence of soil texture
Variations in soil water contents, temperatures, and water fluxes
during the freeze-thaw periods were analyzed at three experimentalsites. While sand was the dominant soil texture at the YD site, sandy
loam was prevalent at the IG and MG sites. Comparing the simulation
results between the YD and IG sites, some notable differences due to13

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

different soil textures were observed. First, since the water holding ca­
pacity is higher in sandy loam than in sand, a decrease in unfrozen watercontents due to soil freeze-thaw cycles at the IG site was smaller (e.g.,
about 0.02 cm3 cm− 3 at the 20 cm depth). On the other hand, sandy
loam’s permeability is lower than that of sand (see Table 2), resulting inmuch smaller hydraulic conductivities and liquid fluxes than at the YD
site (He et al., 2020; Teng et al., 2020). Second, due to sand’s small heat
capacity and thermal conductivity, soil temperatures fluctuated moreduring a day at the YD site than at the IG site (He et al., 2020). For
example, the mean temperature amplitude at the 10 cm depth during thefreezing period was about 5.3℃ at the YD site, while it was less than 3℃
at the IG site. Driven by a smaller temperature gradient, the vapor flux
was smaller at the IG site.
In addition to soil hydraulic and thermal properties discussed above,other numerical models’ parameters, such as the enhancement factor η,
the gain factor Gwt, and the impedance factor ψ , are closely related to soil
textures and could also affect the simulation results. For example, theimportance of considering the enhancement factor when calculating
thermal vapor hydraulic conductivities and vapor fluxes has been
demonstrated by Sakai et al. (2009). However, there are ongoing dis­cussions about the enhancement factor’s uncertainty and rationality
because of the absence of direct vapor flow measurements throughporous media (e.g., Du et al., 2017). The impedance factor also played acrucial role in simulating soil water fluxes during freezing periods in our
simulations. When the soil is frozen, liquid water hydraulic conductivity
is significantly reduced due to ice crystals’ impedance. This reduction isaffected by the initial degree of soil saturation, as indicated by soil
freezing experiments (Lundin, 1990; Zhao et al., 2013). As a result, theimpedance factor depends on the soil type, and its accurate estimate is
essential to obtain correct values of the hydraulic conductivity of frozen
soils.

source (e.g., groundwater, precipitation) for vapor flow and/or thecontribution of vapor flow to the surface vegetation.
6. Conclusions
A fully coupled numerical model for simulating the simultaneous
movement of water, vapor, and heat during the freezing-thawing pe­
riods was developed in this study. The proposed Hydrus-1D Freezing
module avoids numerical instability problems when soil temperatures
drop into subfreezing conditions, and the phase change occurs. The
model was evaluated against experimental data collected at the Yu Lin
Desert, Inner Mongolia Grassland, and Tibetan Plateau Maqu Grassland
sites. Simulated soil water contents and temperatures corresponded wellwith the measured data from an in-situ observation site. Results from the
Yu Lin Desert site show that both liquid water and water vapor flows
driven by soil matric potential and temperature gradients were directed
towards the frozen layer, increasing the total water content there. The
isothermal liquid flux was the most significant component in mostdepths, followed by the thermal liquid and vapor fluxes. However, when
the soil was frozen, liquid fluxes became about 1–5 orders of magnitude
smaller than before freezing. Instead, the thermal vapor flux dominatedsoil water flow in the completely frozen layer and contributed about
10% of the total water flux to the ice formation around the freezing
front. Besides, the often-neglected isothermal vapor flux was also
observable in the top frozen layer.
In general, the simulation results suggest that it is crucial to consider
the effects of vapor flow in evaluating the coupled movement of water,
vapor, and heat during freezing-thawing periods. Similar results as thoseobtained in this study could be expected from other areas with similarsoil and climate conditions. The proposed model can be applied in othercold regions as well.
The new module simulating water flow in soils under freezingthawing conditions has been developed for the widely used Hydrus-1Dmodel (Šimůnek et al., 2016). As a result, many current Hydrus users
will be able to operate the new module while using various advanced
Hydrus-1D software features, such as the well-designed user-friendly
GUI, flexible boundary conditions, parameter optimization routines, andmany other features (Šimůnek et al., 2016). The new Freezing module
also provides the current and future Hydrus users higher numericalstability and computational efficiency than the existing module. The
developed model represents a comprehensive tool to numericallyinvestigate many important research problems involving processes inthe vadose zone.
The Hydrus-1D model has been under continuous development forthe past 20 years (Šimůnek et al., 2016). We are currently developing a
module to simulate the fate and transport of soil water isotopes thatconsiders evaporation fractionation. By coupling this additional isotope
module with our Freezing model, future studies will allow us tonumerically investigate many fundamental research problems involving
sourcing and timing of soil water, such as contributions of liquid water
and water vapor to the surface vegetation in cold regions.

5.3. Model evaluation
Due to the highly nonlinear relationships between soil water content,
ice content, and soil temperature caused by substantial latent heat andphase changes, high computational costs and numerical instabilities arecommon problems in frozen soil modeling (Bao et al., 2016). The earlier
version of the Hydrus-1D Freezing module often encountered numericalproblems when simulating complicated conditions. In contrast, thenewly proposed model, which combines the coupled movement ofwater, vapor, and heat with the soil freezing-thawing process, isnumerically more stable and provides a better agreement with theobserved field data.
It should be emphasized that using hourly rather than daily
measured data is preferable. The hourly data are temporally continuous,
resulting in more accurate and stable simulating results, especiallyduring the phase transition periods. In addition to avoiding numericalinstabilities, the main advantage of the proposed Freezing module is thatit can be used for many different applications and situations (e.g., avariety of boundary conditions), since the standard Hydrus-1D model
has been widely applied and well documented for simultaneously
simulating water flow, and solute and heat transport under nonfreezing
conditions (e.g., Saito et al., 2006; Sakai et al., 2009; Scanlon et al.,2003).
Further improvements to this model will focus on the following as­pects. First, only two approaches are currently available in Hydrus-1D to
calculate the thermal conductivity: Campbell (1985) and Chung andHorton (1987). It would thus be beneficial to incorporate othercommonly used formulas (e.g., de Vries method and Johansen method,as summarized by Yu et al. (2020) and Zhao et al. (2019)), which may
provide a better choice for simulating heat transport in various regionswith different soils. Second, it would also be beneficial to implement
into the Freezing module an option to simulate the transport and fate ofstable water isotopes that would allow the identification of water
sources. One could then identify for the cold and arid regions the water

CRediT authorship contribution statement
Ce Zheng: Conceptualization, Methodology, Writing - original draft.
Jiří Šimůnek: Supervision. Ying Zhao: Data curation. Yudong Lu:
Methodology, Supervision. Xiuhua Liu: Methodology, Supervision.
Changchun Shi: Investigation, Data curation. Huanhuan Li: Investi­
gation, Data curation. Lianyu Yu: Data curation. Yijian Zeng: Data
curation. Zhongbo Su: Data curation.
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.
14

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

Acknowledgments

Li, Q., Sun, S.F., Xue, Y.K., 2010. Analysis and development of a hierarchy of frozen soilmodels for cold region study. J. Geophys. Res. 115, D03107. https://doi.org/
10.1029/2009JD012530.
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.
Millington, R.J., Quirk, J.P., 1961. Permeability of porous solids. Trans. Faraday Soc. 57,
1200–1207. https://doi.org/10.1039/tf9615701200.
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.
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.
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.
Niu, G.Y., Yang, Z.L., 2006. Effects of frozen soil on snowmelt runoff and soil water
storage at a continental scale. J. Hydrometeorol. 7 (5), 937–952. https://doi.org/
10.1175/JHM538.1.
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.
Sakai, M., Toride, N., Šimůnek, J., 2009. Water and vapor movement with condensationand evaporation in a sandy column. Soil Sci. Soc. Am. J. 73 (3), 707–717. https://doi.org/10.2136/sssaj2008.0094.
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, p. 315.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J 15 (7), 25.
https://doi.org/10.2136/vzj2016.04.0033.
Slater, A.G., Pitman, A.J., Desborough, C.E., 1998. Simulation of freeze-thaw cycles in a
general circulation model land surface scheme. J. Geophys. Res. 103 (D10),11303–11312. https://doi.org/10.1029/97JD03630.
Sun, Y., Cheng, Q., Xue, X.Z., Fu, L., Chai, J., Meng, F., Lammers, P.S., Jones, S.B., 2012.
Determining in-situ soil freeze-thaw cycle dynamics using an access tube-based
dielectric sensor. Geoderma 189–190, 321–327. https://doi.org/10.1016/j.
geoderma.2012.06.016.
Teng, J., Kou, J., Yan, X., Zhang, S., Sheng, D., 2020. Parameterization of soil freezing
characteristic curve for unsaturated soils. Cold Reg. Sci. Technol. 170, 102928.
https://doi.org/10.1016/j.coldregions.2019.102928.
Teng, J.D., Shan, F., He, Z.Y., Zhang, S., Zhao, G., Sheng, D.C., 2018. Experimental studyof ice accumulation in unsaturated clean sand. Geotechnique 69 (3), 251–259.
https://doi.org/10.1680/jgeot.17.p.208.
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, 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/wrcr.
v53.610.1002/2017WR020451.
Watanabe, K., Wake, T., 2009. Measurement of unfrozen water content and relative
permittivity of frozen unsaturated soil using NMR and TDR. Cold Reg. Sci. Technol.
59 (1), 34–41. https://doi.org/10.1016/j.coldregions.2009.05.011.
Yi, J., Zhao, Y., Shao, M.A., Zhang, J.G., Cui, L.L., Si, B.C., 2014. Soil freezing andthawing processes affected by the different landscape in the middle reaches of the
Heihe River Basin, Gansu. China. J. Hydrol. 519, 1328–1338. https://doi.org/
10.1016/j.jhydrol.2014.08.042.
Yu, L.Y., Zeng, Y.J., Su, Z.B., 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.
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., 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, S., Teng, J.D., He, Z.Y., Sheng, D.C., 2016. Importance of vapor flow inunsaturated freezing soil: a numerical study. Cold Reg. Sci. Technol. 126, 1–9.
https://doi.org/10.1016/j.coldregions.2016.02.011.
Zhang, Y., Cheng, G., Li, X., Han, X., Wang, L., Li, H., Chang, X., Flerchinger, G.N., 2012.
Coupling of a simultaneous heat and water model with a distributed hydrologicalmodel and evaluation of the combined model in a cold region watershed. Hydrol.

This research was supported by the National Natural Science Foun­dation of China (No. 41877179, 41630634, 41977009), the Funda­
mental Research Funds for the Central Universities, CHD (No.
300102298706, 300102290719, 300102291507), the Shaanxi Water
Conservancy Science and Technology plan project (2019slkj-18),
Taishan Scholars Program (201812096) and Fund Project of Shaanxi
Key Laboratory of Land Consolidation. Gratitude is also expressed forthe financial support from the program of China Scholarship Council.
References
Bao, H.Y., 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. 121 (10), 5259–5280. https://doi.org/10.1002/2015JD024451.
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 (3-4), 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.
Campbell, G.S., 1985. Soil physics with BASIC. Elsevier, New York.
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.
Celia, M.A., Bouloutas, E.T., Zarba, R.L., 1990. A general mass-conservative numericalsolution for the unsaturated flow equation. Water Resour. Res. 26 (7), 1483–1496.
https://doi.org/10.1029/wr026i007p01483.
Cheng, D.H., Li, Y., Chen, X.H., Wang, W.K., Hou, G.C., Wang, C.L., 2013. Estimation ofgroundwater evapotranspiration using diurnal water table fluctuations in the Mu Us
Desert, northern China. J. Hydrol. 490, 106–113. https://doi.org/10.1016/j.
jhydrol.2013.03.027.
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.
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.
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.
Du, C.Y., Yu, J.J., Wang, P., Zhang, Y.C., 2017. Analysing the mechanisms of soil waterand vapour transport in the desert vadose zone of the extremely arid region ofnorthern China. J. Hydrol. 558, 592–606. https://doi.org/10.1016/j.
jhydrol.2017.09.054.
Flerchinger, G.N., Saxton, K.E., 1989. Simultaneous heat and water model of a freezingsnow-residue-soil system. I. Theory and development. Trans. ASAE 32 (2), 565–571.
https://doi.org/10.13031/2013.31040.
Guymon, G.L., Luthin, J.N., 1974. A coupled heat and moisture transport model for arctic
soils. Water Resour. Res. 10 (5), 995–1001. https://doi.org/10.1029/
WR010i005p00995.
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.
He, Z., Teng, J., Yang, Z., Liang, L., Li, H., Zhang, S., 2020. An analysis of vapour transfer
in unsaturated freezing soils. Cold Reg. Sci. Technol. 169, 102914. https://doi.org/
10.1016/j.coldregions.2019.102914.
He, Z.Y., Zhang, S., Teng, J.D., Yao, Y.P., Sheng, D.C., 2018. A coupled model for liquidwater-vapor-heat migration in freezing soils. Cold Reg. Sci. Technol. 148, 22–28.
https://doi.org/10.1016/j.coldregions.2018.01.003.
Hillel, D., 1971. Soil and Water: Physical Principles and Processes. Academic Press, NewYork.
Jansson, P.E., 2012. CoupModel: Model use, calibration, and validation. Trans. ASABE
55 (4), 1337. https://doi.org/10.13031/2013.42245.
Kelleners, T.J., 2013. Coupled water flow and heat transport in seasonally frozen soils
with snow accumulation. Vadose Zone J. 12 (4), 108–118. https://doi.org/10.2136/
vzj2012.0162.
Kelleners, T.J., Koonce, J., Shillito, R., Dijkema, J., Berli, M., Young, M.H., Frank, J.M.,
Massman, W.J., 2016. Numerical modeling of coupled water flow and heat transport
in soil and snow. Soil Sci. Soc. Am. J. 80 (2), 247–263. https://doi.org/10.2136/
sssaj2015.07.0279.
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.
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.

15

C. Zheng et al.

Journal of Hydrology 598 (2021) 126250

Process. 27 (25), 3762–3776. https://doi.org/10.1002/hyp.v27.2510.1002/
hyp.9514.
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, Y., Si, B.C., He, H.L., Xu, J.H., Peth, S., Horn, R., 2016. Modeling of coupled waterand heat transfer in freezing and thawing soils. Inner Mongolia. Water 8 (10),424–441. https://doi.org/10.3390/w8100424.

Zhao, Y., Si, B.C., Zhang, Z.H., Li, M., He, H.L., Hill, R.L., 2019. A new thermal
conductivity model for sandy and peat soils. Agr. Forest Meteorol. 274, 95–105.
https://doi.org/10.1016/j.agrformet.2019.04.004.
Zheng, C.e., Šimůnek, J., Lu, Y., Liu, X., Shi, C., Li, H., 2020. Monitoring and modelingthe coupled movement of water, vapor, and energy in arid areas. J. Hydrol. 590,
125528. https://doi.org/10.1016/j.jhydrol.2020.125528.

16