# Numerical Modeling of Soil Evaporation Process and Its Stages Dividing during a Drying Cycle

Authors: Jiangbo Han; Jin Lin; Yunfeng Dai
Journal: Geofluids
 (2017)

---

## Abstract

The soil water evaporation is a critical component of both the surface energy balance and water balance, affecting the mass and energy exchange between the land and the atmosphere. Evaporation process is involved in the highly complex interactions between media properties, transport processes, and boundary conditions. So, it is difficult to accurately determine these near-surface highly dynamic processes based only on the sparse field data and on the measurement-based methods. The objective of this paper was to obtain a detailed description of the soil water evaporation process and to better understand the evolutions of variables involved in the evaporation process during different stages of evaporation. To do this, a numerical simulation experiment in a bare silt soil was conducted to reproduce the soil drying process during a 20-d period after a 2-cm rainfall event. According to simulation results, the whole 20-d simulation period was divided into two main stages as well as a transient period from stage 1 to stage 2. Diurnal patterns of energy and water balance components, soil moisture, soil temperature, water fluxes, evaporation rate, dry surface layer (DSL), and evaporation zone during this drying process were fully described, which, in turn, could be as the possible indicators for judging the shift of stages of evaporation.


---

## Body

## 1. Introduction

Infiltration and soil water movement play a key element in surface runoff, groundwater recharge, evapotranspiration, soil erosion, and transport of chemicals in surface and subsurface waters [1]. In particular, soil water evaporation and storage during natural wetting/drying cycle exert critical influences on land–atmosphere exchange of water and energy [2, 3], thereby affecting the local or global climate changes. In addition to the important role in the hydrological process, drying and evaporation processes are also of interest for many engineering and industrial applications, such as food processing and preservation, production of ceramics and paper, eye and skin care, and numerous construction activities [4].During drying of soils after a rainfall or irrigation event, shifts between atmospheric demand to soil control evaporation occur with the progressive drying of soils [3]. Soil water evaporation is a dynamic process and this process has frequently been divided into three major stages. Stage 1 evaporation has a relatively high and constant evaporation rate followed by a lower rate (stage 2) [5, 6]. During stage 1 when the soil surface has available water, evaporation takes place at the soil surface and is limited by atmospheric demand. At stage 1, a large portion of net radiation is readily dissipated as latent heat of vaporization. With the soil continuing to dry, evaporation eventually occurs within the subsurface when the soil water content is depleted; and stage 2 evaporation begins, during which the evaporation drops below the potential rate and decreases with time. Stage 3 is characteristic of a very low and relatively constant evaporation rate. When the soil moisture is not available at the soil surface, a dry surface layer (DSL) forms near the soil surface, in which moisture transfer occurs only in the vapor phase [7] with the profile below this layer remaining much wetter [8]. The diffusive vapor transfer within the dry layer is supplied by the liquid water flow from sites below the evaporation zone (where the phase change from liquid water to water vapor is occurring). The transition from the surface to subsurface evaporation also has a significant effect on the surface energy balance and profiles of soil temperature and moisture because the energy required for vaporization of liquid water must be transported from the soil surface to subsurface evaporation zone [4, 9]. Formation of DSL and the shift of the evaporation front from the surface to the subsurface result in a significant change in the surface energy balance as the location of latent heat sink moves from surface to subsurface. The latent heat flux associated with this phase change is important for assessing the surface energy balance during the soil drying [10–12]. The evaporation stage could be an indicator for total water loss, evaporation rate, and soil water content [5], as well as surface energy balance [9]. Therefore, it is important to understand soil evaporation stages for water resource management, environmental conservation, and agricultural production.To evaluate the evaporative drying process in the field condition, especially to determine the timing of stages of evaporation process, there is need to determine the evaporation rate from soils. The approaches for measuring soil water evaporation include the Bowen ratio [13], eddy correlation and automatic weighing lysimeter methods [1]. These techniques are widely used for field measurements to study the water balance dynamics near the soil surface. However, there are some limitations in these measurement methods; for example, it is somewhat time costly and laborious and cannot accurately measure the dynamics of soil water evaporation with time and depth, especially for evaluating the effect of the evaporation process on the surface energy balance and soil temperature and moisture profiles. Although there are some new approaches developed recently to determine the process of soil water evaporation, it is difficult for these measurement-based approaches to capture the highly dynamic near-surface soil water and heat processes involved in soil water evaporation. For instance, some researchers [3, 14] introduced a measurement approach based on sensible heat balance to determine in situ soil water evaporation dynamics by means of heat pulse probes. However, this approach was not available for stage 1 evaporation [3] during which evaporation likely occurred in the days after rainfall or irrigation. And this approach could result in underestimation of total subsurface evaporation due to the existence of undetectable zone [4]. Deol et al. [15] quantify the millimeter-scale subsurface evaporation profiles by using eleven-needle heat pulse probes in the soil column over a drying event. Since these profiles were limited to depth-integrated evaporation rates, the accurate shape and variation of evaporation zone with time and depth remained uncertain, particularly for some cases whereas the width of evaporation zone was less than 1 mm [4, 9]. Thus, because the evaporation process is involved in the highly complex interactions between media properties, transport processes, and boundary conditions, the direct measurement approaches have a limitation to accurately determine the highly dynamic soil water evaporation process. And numerical modeling is an alternative means.The growing computational capacity provides the possibility of numerically reproducing evaporation processes by applying the theory of the coupled heat and moisture transport in variable saturated soils, which is developed in Philip and de Vries [16], henceforth PDV. Currently, there are several numerical codes, for example, SVAT, SHAW, and HYDRUS, developed for simulating coupled water and heat flow in the vadose zone. HYDRUS-1D is a widely used finite-element model for simulating one-dimensional movement of water, heat, and multiple solutes in variably saturated porous media [17] and is applied widely in the study of vadose process. The objective of this paper was to evaluate the process of evaporative drying from a silt soil after a rainfall event by using the HYDRUS-1D model under the field-like atmospheric conditions and to better understand the evolutions of variables involved in the evaporation process.


## 2. Materials and Methods

### 2.1. Water Flow Equations

The general partial differential equation for describing the one-dimensional transient water flow in the soil under variably saturated, nonisothermal conditions can be expressed as [4, 16](1)∂θL+θv∂t=-∂qw∂z=-∂qL∂z-∂qv∂z,where θL is the volumetric liquid water content (m3 m−3), θv is the volumetric water vapor content (expressed as an equivalent water content, ρvθair/ρL, m3 m−3), θair is volumetric air content (m3 m−3), ρL is the density of liquid water (kg m−3) at T (°C),  qw(=qL+qv) is the total water flux (m s−1), qL and qv are the flux densities of liquid water and water vapor (expressed as an equivalent water flux density m s−1), respectively, t is time (s), and z is the spatial coordinate positive upward (m).The PDV model describes liquid water and water vapor movements driven both by the temperature and matric potential gradients. Thus, the flux densities of liquid water,qL, and water vapor, qv are given by, respectively, (2)qL=qTL+qhL=-DTL∂T∂z-K∂h∂z+1(3)qv=qTv+qhv=-DTv∂T∂z-Dhv∂h∂z,where qTL and qhL are the thermal and isothermal liquid water fluxes (m s−1), respectively; qTv and qhv are the thermal and isothermal water vapor fluxes (m s−1), respectively; T is the soil temperature (°C); h is the pressure head (m); DTL (m2 K−1 s−1) and K (m s−1) are the thermal and isothermal hydraulic conductivities for liquid phase fluxes, respectively; and DTv (m2 K−1 s−1) and Dhv (m s−1) are the thermal and isothermal hydraulic conductivities for water vapor fluxes, respectively. The hydraulic conductivities and related parameters are summarized in Table 1. Combining (1) to (3), the governing liquid water and water vapor equation can expressed as(4)∂θL+θv∂t=∂∂zDTL∂T∂z+K∂h∂z+DTv∂T∂z+Dhv∂h∂z+K.Table 1
Hydraulic conductivities for liquid-water and water vapor flux and related parameters.
Hydraulic conductivity
Parameter
Expression/constant
K
=
K
s
S
e
l
1
-
(
1
-
S
e
1
/
m
)
m
2
Pore connectivity coefficient
l
=
0.5
Gain factor
G
w
T
=
7
Surface tension of soil water, g s−2
γ
=
75.6
-
0.1425
T
-
2.38
×
10
-
4
T
2
Surface tension at 25°C, g s−2
γ
0
=
71.89
D
T
L
=
K
h
G
w
T
1
γ
0
∂
γ
∂
T
Vapor diffusivity in soil, m2 s−1
D
=
D
a
Ω
θ
a
i
r
Tortuosity factor in gaseous phase
Ω
=
θ
a
i
r
2
/
3
Vapor diffusivity in air, m2 s−1
D
a
=
2.12
×
10
-
5
T
+
273.15
273.15
D
T
v
=
D
ρ
L
η
H
r
d
ρ
s
v
d
T
+
ρ
s
v
d
H
r
d
T
Density of liquid water, kg m−3
ρ
L
=
1000
-
7.3
×
10
-
3
T
-
4
2
+
3.79
×
10
-
5
(
T
-
4
)
3
Saturation vapor density, kg m−3
ρ
s
v
=
10
-
3
×
exp
⁡
19.84
-
4975.9
T
+
273.15
Relative humidity of soil surface
H
r
=
exp
⁡
h
M
g
R
T
+
273.15
Molecular weight of water, kg mol−1
M
=
0.01805
Gravitation acceleration, m s−2
g
=
9.81
Universal gas constant, J mol−1 K−1
R
=
8.314
D
h
v
=
D
ρ
L
H
r
M
g
R
T
ρ
s
v
Enhancement factor
η
=
9.5
+
3
θ
L
θ
s
-
8.5
exp
⁡
-
1
+
2.6
/
f
c
θ
L
θ
s
4
Mass fraction of clay
0.02The van Genuchten model which described the soil water retention relation was used for the silt [18]:(5)θL=θr+θs-θr1+αhn-m,where θs is the saturated water content (m3 m−3), θr is the residual water content (m3 m−3), and α (m−1), n, and m(=1-1/n) are empirical shape parameters. The parameters employed for the soil are shown in Table 2.Table 2
Soil hydraulic properties used in numerical simulations.
Sample
θ
r (m3 m−3)
θ
s (m3 m−3)
α (cm−1)
n
K
s (cm d−1)
silt
0.034
0.46
0.016
1.37
6


### 2.2. Heat Transport Equations

The governing equation for the one-dimensional movement of energy in a variably saturated vadose zone is given by [19] (6)∂Sh∂t=-∂qh∂z,where Sh is the storage of heat in the soil (J m−3) and qh is the total soil heat flux density (J m−2 s−1). The storage of heat Sh (J m−3) in the soil is (7)Sh=CsT-Trθn+CvT-Trθv+CwT-TrθL+L0θv,where Cs(=1.92 MJ m−3 K−1)  Cw(=4.18 MJ m−3 K−1) and Cv(= 1.8 MJ m−3 K−1) are the volumetric heat capacities of dry soil particles, liquid water, and water vapor, respectively; Tr is the reference temperature (°C); and θn is the volumetric fraction of solid phase (m3 m−3). The soil heat flux density qh (J m−2 s−1), accounting for the sensible heat of the conduction, sensible heat by the convection of liquid water and water vapor, and latent heat by vapor flow, can be described as(8)qh=-λ∂T∂z+CvT-Trqv+CwT-TrqL+L0qv,where λ is the apparent soil thermal conductivity (W m−1 K−1) described by [20](9)λ=b1+b2θL+b3θL0.5,where b1, b2, and b3 are empirical parameters (W m−1 K−1) (for silt: b1=0.243,  b2=0.393, and  b3=1.534); L0(=ρL∗(2.501×106-2369.2T)) is the volumetric latent heat of vaporization of water (J m−3). Combining the continuity equation with (6) to (8) produces the governing equation for the movement of energy in soils:(10)Cp∂T∂t+L0∂θv∂t+CvT-Tr∂θv∂t+CwT-Tr∂θL∂t=∂∂zλ∂T∂z-Cv∂Tqv∂z-Cw∂TqL∂z-L0∂qv∂z,where Cp(=Csθn+Cvθv+CwθL) represents the volumetric heat capacity of the moist soil (J m-3K), and the contribution of air to Cp is considered negligible.


### 2.3. Surface Water and Energy Boundaries

Surface precipitation, evaporation, and heat fluxes are used as boundary conditions for liquid water and water vapor flow and heat transport in field soils. Surface boundary conditions for the water and energy equations are, respectively, given by [21, 22] (11)qL0,t+qv0,t=E(12)Rn-H-LE-G=0,where qL0,t and qv0,t are liquid water and water vapor fluxes through the soil surface, respectively, E is the evaporation rate (m s−1) from the soil, Rn is net radiation (W m−2), H is the sensible heat flux density (W m−2), LE is the latent heat flux density (W m−2), L is the latent heat (J kg−1), and G is the surface heat flux density (W m−2). While Rn and G are positive downward, H and LE are positive upward [19]. Rn is defined as(13)Rn=Rns+Rnl=1-aSt+εsσεaTa4-Ts4,where Rns and Rnl are net shortwave and longwave radiations (W m−2), respectively, a is soil albedo (unitless), St is incoming (global) shortwave solar radiation (W m−2), σ is the Stefan Boltzmann constant (5.6704 × 10−8 W m−2 K−4), εs is soil longwave emissivity (unitless), and εa is the effective atmospheric emissivity (unitless). The sensible heat flux can be defined as [23](14)H=CaTs-Tarh,where Ca is the volumetric heat capacity of air (1200 J m−3 K−1), Ts and Ta are temperature (°C) of the soil surface and the air, respectively, and rh is the aerodynamic resistance to heat transfer (s m−1). The soil surface evaporation rate can be calculated from the difference between the water vapor densities of the air, ρva (kg m−3), and the soil surface ρvs (kg m−3): (15)E=ρvs-ρvarv+rs,where rv is the aerodynamic resistance to water vapor flow (s m−1) and rs is the soil surface resistance to water vapor flow that acts as an additional resistance along with the aerodynamic resistance (s m−1). rv and rh can be defined as [24](16)rv=rh=1uk2lnzref-d+zhzoh+ψhlnzref-d+zmzom+ψm,where u is wind speed (m s−1), k is von Karman constant (which has a value of 0.41), d is the zero plane displacement (m), zh and zm are the height of temperature and wind speed measurement (m), respectively, zoh and zom are the surface roughness length for heat flux and momentum flux (m), respectively, and ψh and ψm are the atmospheric stability correction factors for heat flux and momentum flux (—), respectively. And rs can be expressed as [25, 26](17)rs=10exp⁡35.630.15-θtop,where θtop is the surface water content (m3 m−3). Details on these variables also can be found in Saito et al. [19, 27].


### 2.4. Numerical Modeling Procedure

The soil profile used in this simulation was considered a 50-cm deep silt horizon. The simulated period lasted for 21 days including a day with a 2-cm daily rainfall following a 20-d soil drying process under clear sky conditions. Boundary conditions at the soil surface for liquid water, water vapor, and heat transport were determined from the surface water and energy balance equations described above. The soil surface was imposed to the field-like atmospheric forcing condition. Thus, the diurnally varying air temperature and air relative humidity could be determined by trigonometric functions with a period of 24 h [4]:(18)Ta=Tmax+Tmin2+Tmax-Tmin2·cos⁡2πt-1324(19)Hr=Hr¯+Hrmax-Hrmin2·cos⁡2πt-124,where Tmax and Tmin are the maximum and minimum air temperatures, respectively, and Hr¯, Hrmax, and Hrmin are the average, maximum, and minimum air relative humidity, respectively. The incoming shortwave solar radiation, St, at any given time can be calculated as follows [19]:(20)St=f×max⁡Stdsin⁡φsin⁡σ+cos⁡φcos⁡σcos⁡2πt-12,0,where Std is the daily solar radiation (W m−2), φ is the site latitude (rad), and σ is the solar declination (rad); f is used to adjust the daily summation of St equal to Std. In this study, we assumed Tmax=30°C and Tmax=20°C,  Hr¯ = 30%,  u=1.2 m s−1, and Std=230 W m−2 for all the days. Figure 1 shows the diurnal variations of Ta,  Hr, and St for a period of 24 h.Figure 1
Diurnal variations of air temperatureTa, relative humidity Hr, and incoming shortwave solar radiation St.Zero pressure head (free drainage) and temperature gradients were used as the bottom boundary conditions of the soil domain. These conditions assume that the water table is located far below the domain of interest and that heat transfer across the lower boundary occurs only by convection of liquid water and water vapor. The initial temperature of 25°C and uniform water content of 0.2 m3 m−3 were applied to the soil.The variably saturated water flow and heat transport equations as well as the surface energy and water balances were solved numerically using HYDURUS-1D software package. The 50-cm long soil profile for simulation was divided into 500 elements, with the mesh thickness increasing linearly from 0.03 cm at the soil surface to 0.27 cm at the bottom of the soil domain. The time step was allowed to vary between the initial and maximum time steps of 1 × 10−5 and 0.01 day, respectively.


## 2.1. Water Flow Equations

The general partial differential equation for describing the one-dimensional transient water flow in the soil under variably saturated, nonisothermal conditions can be expressed as [4, 16](1)∂θL+θv∂t=-∂qw∂z=-∂qL∂z-∂qv∂z,where θL is the volumetric liquid water content (m3 m−3), θv is the volumetric water vapor content (expressed as an equivalent water content, ρvθair/ρL, m3 m−3), θair is volumetric air content (m3 m−3), ρL is the density of liquid water (kg m−3) at T (°C),  qw(=qL+qv) is the total water flux (m s−1), qL and qv are the flux densities of liquid water and water vapor (expressed as an equivalent water flux density m s−1), respectively, t is time (s), and z is the spatial coordinate positive upward (m).The PDV model describes liquid water and water vapor movements driven both by the temperature and matric potential gradients. Thus, the flux densities of liquid water,qL, and water vapor, qv are given by, respectively, (2)qL=qTL+qhL=-DTL∂T∂z-K∂h∂z+1(3)qv=qTv+qhv=-DTv∂T∂z-Dhv∂h∂z,where qTL and qhL are the thermal and isothermal liquid water fluxes (m s−1), respectively; qTv and qhv are the thermal and isothermal water vapor fluxes (m s−1), respectively; T is the soil temperature (°C); h is the pressure head (m); DTL (m2 K−1 s−1) and K (m s−1) are the thermal and isothermal hydraulic conductivities for liquid phase fluxes, respectively; and DTv (m2 K−1 s−1) and Dhv (m s−1) are the thermal and isothermal hydraulic conductivities for water vapor fluxes, respectively. The hydraulic conductivities and related parameters are summarized in Table 1. Combining (1) to (3), the governing liquid water and water vapor equation can expressed as(4)∂θL+θv∂t=∂∂zDTL∂T∂z+K∂h∂z+DTv∂T∂z+Dhv∂h∂z+K.Table 1
Hydraulic conductivities for liquid-water and water vapor flux and related parameters.
Hydraulic conductivity
Parameter
Expression/constant
K
=
K
s
S
e
l
1
-
(
1
-
S
e
1
/
m
)
m
2
Pore connectivity coefficient
l
=
0.5
Gain factor
G
w
T
=
7
Surface tension of soil water, g s−2
γ
=
75.6
-
0.1425
T
-
2.38
×
10
-
4
T
2
Surface tension at 25°C, g s−2
γ
0
=
71.89
D
T
L
=
K
h
G
w
T
1
γ
0
∂
γ
∂
T
Vapor diffusivity in soil, m2 s−1
D
=
D
a
Ω
θ
a
i
r
Tortuosity factor in gaseous phase
Ω
=
θ
a
i
r
2
/
3
Vapor diffusivity in air, m2 s−1
D
a
=
2.12
×
10
-
5
T
+
273.15
273.15
D
T
v
=
D
ρ
L
η
H
r
d
ρ
s
v
d
T
+
ρ
s
v
d
H
r
d
T
Density of liquid water, kg m−3
ρ
L
=
1000
-
7.3
×
10
-
3
T
-
4
2
+
3.79
×
10
-
5
(
T
-
4
)
3
Saturation vapor density, kg m−3
ρ
s
v
=
10
-
3
×
exp
⁡
19.84
-
4975.9
T
+
273.15
Relative humidity of soil surface
H
r
=
exp
⁡
h
M
g
R
T
+
273.15
Molecular weight of water, kg mol−1
M
=
0.01805
Gravitation acceleration, m s−2
g
=
9.81
Universal gas constant, J mol−1 K−1
R
=
8.314
D
h
v
=
D
ρ
L
H
r
M
g
R
T
ρ
s
v
Enhancement factor
η
=
9.5
+
3
θ
L
θ
s
-
8.5
exp
⁡
-
1
+
2.6
/
f
c
θ
L
θ
s
4
Mass fraction of clay
0.02The van Genuchten model which described the soil water retention relation was used for the silt [18]:(5)θL=θr+θs-θr1+αhn-m,where θs is the saturated water content (m3 m−3), θr is the residual water content (m3 m−3), and α (m−1), n, and m(=1-1/n) are empirical shape parameters. The parameters employed for the soil are shown in Table 2.Table 2
Soil hydraulic properties used in numerical simulations.
Sample
θ
r (m3 m−3)
θ
s (m3 m−3)
α (cm−1)
n
K
s (cm d−1)
silt
0.034
0.46
0.016
1.37
6


## 2.2. Heat Transport Equations

The governing equation for the one-dimensional movement of energy in a variably saturated vadose zone is given by [19] (6)∂Sh∂t=-∂qh∂z,where Sh is the storage of heat in the soil (J m−3) and qh is the total soil heat flux density (J m−2 s−1). The storage of heat Sh (J m−3) in the soil is (7)Sh=CsT-Trθn+CvT-Trθv+CwT-TrθL+L0θv,where Cs(=1.92 MJ m−3 K−1)  Cw(=4.18 MJ m−3 K−1) and Cv(= 1.8 MJ m−3 K−1) are the volumetric heat capacities of dry soil particles, liquid water, and water vapor, respectively; Tr is the reference temperature (°C); and θn is the volumetric fraction of solid phase (m3 m−3). The soil heat flux density qh (J m−2 s−1), accounting for the sensible heat of the conduction, sensible heat by the convection of liquid water and water vapor, and latent heat by vapor flow, can be described as(8)qh=-λ∂T∂z+CvT-Trqv+CwT-TrqL+L0qv,where λ is the apparent soil thermal conductivity (W m−1 K−1) described by [20](9)λ=b1+b2θL+b3θL0.5,where b1, b2, and b3 are empirical parameters (W m−1 K−1) (for silt: b1=0.243,  b2=0.393, and  b3=1.534); L0(=ρL∗(2.501×106-2369.2T)) is the volumetric latent heat of vaporization of water (J m−3). Combining the continuity equation with (6) to (8) produces the governing equation for the movement of energy in soils:(10)Cp∂T∂t+L0∂θv∂t+CvT-Tr∂θv∂t+CwT-Tr∂θL∂t=∂∂zλ∂T∂z-Cv∂Tqv∂z-Cw∂TqL∂z-L0∂qv∂z,where Cp(=Csθn+Cvθv+CwθL) represents the volumetric heat capacity of the moist soil (J m-3K), and the contribution of air to Cp is considered negligible.


## 2.3. Surface Water and Energy Boundaries

Surface precipitation, evaporation, and heat fluxes are used as boundary conditions for liquid water and water vapor flow and heat transport in field soils. Surface boundary conditions for the water and energy equations are, respectively, given by [21, 22] (11)qL0,t+qv0,t=E(12)Rn-H-LE-G=0,where qL0,t and qv0,t are liquid water and water vapor fluxes through the soil surface, respectively, E is the evaporation rate (m s−1) from the soil, Rn is net radiation (W m−2), H is the sensible heat flux density (W m−2), LE is the latent heat flux density (W m−2), L is the latent heat (J kg−1), and G is the surface heat flux density (W m−2). While Rn and G are positive downward, H and LE are positive upward [19]. Rn is defined as(13)Rn=Rns+Rnl=1-aSt+εsσεaTa4-Ts4,where Rns and Rnl are net shortwave and longwave radiations (W m−2), respectively, a is soil albedo (unitless), St is incoming (global) shortwave solar radiation (W m−2), σ is the Stefan Boltzmann constant (5.6704 × 10−8 W m−2 K−4), εs is soil longwave emissivity (unitless), and εa is the effective atmospheric emissivity (unitless). The sensible heat flux can be defined as [23](14)H=CaTs-Tarh,where Ca is the volumetric heat capacity of air (1200 J m−3 K−1), Ts and Ta are temperature (°C) of the soil surface and the air, respectively, and rh is the aerodynamic resistance to heat transfer (s m−1). The soil surface evaporation rate can be calculated from the difference between the water vapor densities of the air, ρva (kg m−3), and the soil surface ρvs (kg m−3): (15)E=ρvs-ρvarv+rs,where rv is the aerodynamic resistance to water vapor flow (s m−1) and rs is the soil surface resistance to water vapor flow that acts as an additional resistance along with the aerodynamic resistance (s m−1). rv and rh can be defined as [24](16)rv=rh=1uk2lnzref-d+zhzoh+ψhlnzref-d+zmzom+ψm,where u is wind speed (m s−1), k is von Karman constant (which has a value of 0.41), d is the zero plane displacement (m), zh and zm are the height of temperature and wind speed measurement (m), respectively, zoh and zom are the surface roughness length for heat flux and momentum flux (m), respectively, and ψh and ψm are the atmospheric stability correction factors for heat flux and momentum flux (—), respectively. And rs can be expressed as [25, 26](17)rs=10exp⁡35.630.15-θtop,where θtop is the surface water content (m3 m−3). Details on these variables also can be found in Saito et al. [19, 27].


## 2.4. Numerical Modeling Procedure

The soil profile used in this simulation was considered a 50-cm deep silt horizon. The simulated period lasted for 21 days including a day with a 2-cm daily rainfall following a 20-d soil drying process under clear sky conditions. Boundary conditions at the soil surface for liquid water, water vapor, and heat transport were determined from the surface water and energy balance equations described above. The soil surface was imposed to the field-like atmospheric forcing condition. Thus, the diurnally varying air temperature and air relative humidity could be determined by trigonometric functions with a period of 24 h [4]:(18)Ta=Tmax+Tmin2+Tmax-Tmin2·cos⁡2πt-1324(19)Hr=Hr¯+Hrmax-Hrmin2·cos⁡2πt-124,where Tmax and Tmin are the maximum and minimum air temperatures, respectively, and Hr¯, Hrmax, and Hrmin are the average, maximum, and minimum air relative humidity, respectively. The incoming shortwave solar radiation, St, at any given time can be calculated as follows [19]:(20)St=f×max⁡Stdsin⁡φsin⁡σ+cos⁡φcos⁡σcos⁡2πt-12,0,where Std is the daily solar radiation (W m−2), φ is the site latitude (rad), and σ is the solar declination (rad); f is used to adjust the daily summation of St equal to Std. In this study, we assumed Tmax=30°C and Tmax=20°C,  Hr¯ = 30%,  u=1.2 m s−1, and Std=230 W m−2 for all the days. Figure 1 shows the diurnal variations of Ta,  Hr, and St for a period of 24 h.Figure 1
Diurnal variations of air temperatureTa, relative humidity Hr, and incoming shortwave solar radiation St.Zero pressure head (free drainage) and temperature gradients were used as the bottom boundary conditions of the soil domain. These conditions assume that the water table is located far below the domain of interest and that heat transfer across the lower boundary occurs only by convection of liquid water and water vapor. The initial temperature of 25°C and uniform water content of 0.2 m3 m−3 were applied to the soil.The variably saturated water flow and heat transport equations as well as the surface energy and water balances were solved numerically using HYDURUS-1D software package. The 50-cm long soil profile for simulation was divided into 500 elements, with the mesh thickness increasing linearly from 0.03 cm at the soil surface to 0.27 cm at the bottom of the soil domain. The time step was allowed to vary between the initial and maximum time steps of 1 × 10−5 and 0.01 day, respectively.


## 3. Results and Discussion

During soil drying after a rainfall or irrigation event, the soil continues to dry with the surface water loss by evaporation. This process also affects or is accompanied by changes in the water content, soil temperature, surface energy balance, and soil water balance. In the following sections, the soil water evaporation stages are divided and variables associated with evaporation are discussed at different stages.


### 3.1. Soil Water Evaporation (Drying) Process

As shown in (11), the evaporation rate E (equals the total water flux at the soil surface), actually could be expressed as the sum of surface liquid and vapor water fluxes. Figure 2 presents diurnal variations of the two components of actual surface evaporation rate E. The surface liquid water flux qL(0,t) indicates vaporization of water occurring at the soil surface as long as its values are larger than 0, while the positive surface water vapor flux qv(0,t) indicates that water vaporized within the subsurface and its negative values represent that some water vapor was moving downward from the soil surface to deeper profiles. As shown in Figure 2(a), qL(0,t) totally dominated E until about the noon of day 5, indicating that evaporation was occurring at the soil surface. During this period, the soil was wet enough to supply water to accomplish vaporization of water, and soil evaporation was limited only by the supply of energy. This could also be determined from Figure 3 which shows that the actual evaporation rate E equals to the potential evaporation rate during this period. qv(0,t)<0 (Figure 2(b)) is observed around the noon for some days, suggesting that some water vapor produced at the soil surface was moving downward driven by the large downward temperature gradient at noon. Between days 5 and 8 qL(0,t) was the dominant source early in the day and qv(0,t) took over for the rest of the day with some overlap where both qL(0,t) and qv(0,t) occurred simultaneously. During this period, the actual evaporation rate became less than the potential one (Figure 3) during afternoons. Beginning on day 9 values of qL(0,t) were close zero for the all days, and qv(0,t) became the dominant source of evaporation E (Figure 3), suggesting that vaporization of water occurred in the subsurface at most time of the day during this period. It is also notable that a very small negative qL(0,t) occurred around late afternoon after day 8, which resulted from the condensation of water vapor at the soil surface or in air when atmospheric evaporation demand decreased. From day 9 on, E/E0 presented negative values around the noon, which was caused by the negative qv(0,t) during the same period in the result of the downward temperature gradient. In this case, E actually does not represent the rate of evaporation from the soil but indicates that the flux at the soil surface is moving towards deeper layers as the form of water vapor.Figure 2
The liquid water (red solid line), water vapor (red dashed line) water fluxes, and the total flux (black solid line) at the soil surface. The left vertical dashed lines represent timing point of shift of evaporation process from stages 1 and 2, and periods between the two dashed line indicate the transient stage of evaporation process (same below).
(a)
(b)Figure 3
The temporal variation of ratio of actual evaporation rate (E) and potential (E0) rate E/E0 during the simulation period.Thus, based on these results, the stages of evaporation process for the silt soil could be reasonably determined. Stage 1 evaporation was maintained from day 1 to the noon of day 5, indicated byE/E0=1. Thereafter, stage 2 was maintained until the end of this numerical experiment, with a period as a transient stage in which the shift between qL(0,t) and qv(0,t) occurred alternately every day (from days 5 to 8). Stage 2 was also referred as falling rate stage evaporation because during this stage the actual evaporation rate E fell below the potential evaporation rate E0. There was no clear limit between stage 2 and stage 3 evaporation according to data shown in Figures 2 and 3. The so-called stage 3 evaporation, which was characteristic of the low and constant evaporation rate, was not considered to be very important and often merged into stage 2.


### 3.2. Soil Water Content and Temperature

Figure4 presents temporal changes in the volumetric soil water contents at selected depths near the soil surface during a 20-d drying period immediately following a rain event before day 1 (i.e., at day 0). Soil water redistribution and soil drying began at day 1. Water contents near the soil surface decreased rapidly for the first few days, caused by gravity drainage and evaporation imposed by the hot and dry atmospheric conditions (Figure 1). The extent of water loss in the profile declined with depth, with more pronounced decrease occurring at the soil surface (0 cm) where the water content had a significant decrease from the maximum value of 0.43 m3 m−3 at the beginning of day 1 to the minimum value of 0.046 m3 m−3 at the noon of days 5 to 20 during stage 2. It is also notable that the water content near the soil surface (e.g., at 0 and 2.5 cm depths) displayed a distinct diurnal pattern, with it decreasing from late morning through the afternoon and then increasing from afternoon to late morning on the following day. This pattern during evaporation process was consistent with the findings from a number of field measurements [14, 28–31]. This diurnal pattern for surface water content derived from evaporation during daytime and subsequent soil rewetting during nighttime. The consistency between simulated patterns of soil water content in this study and observations in the field conditions supports the reliability of our models on reproducing the evaporative drying process of soils. The amplitude of diurnal pattern decreased with depth, with presenting a small fluctuation at 2.5 cm and no apparent trace for the deeper soil layer, for example, at 5 cm. After day 8, this diurnal pattern of soil water content at the soil surface maintained the same amplitude and magnitude, with maximum values of surface water content maintained at the same value of 0.06 m3 m−3 and the maximum values of surface temperature at around 44.5°C (Table 3).Table 3
Maximum and minimum values of soil water content and the corresponding time points.
Day
Min (m3 m−3)
Day
Max (m3 m−3)
1.6
0.247
1.0
0.43
2.6
0.112
2.2
0.30
3.6
0.058
3.3
0.23
4.6
0.047
4.3
0.15
5.5
0.046
5.3
0.09
6.5
0.046
6.3
0.07
7.5
0.046
7.3
0.07
8.5
0.046
8.2
0.06
9.5
0.046
9.2
0.06
10.5
0.046
10.2
0.06
11.5
0.046
11.2
0.06
12.5
0.046
12.2
0.06
13.5
0.046
13.2
0.06
14.5
0.046
14.2
0.06
15.5
0.046
15.2
0.06
16.5
0.046
16.2
0.06
17.5
0.046
17.2
0.06
18.5
0.046
18.2
0.06
19.5
0.046
19.2
0.06Figure 4
Simulated temporal variations of soil water contents at depths of 0, 2.5, and 5 cm during the soil drying after a rainfall event.Figure5 is the same as Figure 4 but for temporal changes in soil temperatures during the simulation period. It is clear that the soil temperatures at three depths all display a typical diurnal sinusoidal pattern over all the days, with its maximum temperatures occurring around noon (e.g., at 1 pm at 0 cm) and its minimum temperatures occurring just before the sunrise. In the early few days, the relatively uniform and large water content near the soil surface resulted in the uniform and large soil heat capacity (remember that soil heat capacity varies only as function of water content in soil). This resulted in a relatively uniform soil temperature profile for this period, with the daily temperature magnitudes and amplitudes being a little different between the three depths. On the other hand, the initially wet soil conditions combined with high-radiation atmospheric conditions caused a large amount of vaporization of water near the soil surface, causing more incoming solar radiation to be applied to accomplish vaporization rather than heating the soil. These effects caused relatively lower daily temperature magnitudes and amplitudes for the first few days. With the soil drying, both minimum and maximum temperatures in the shallow profile were increasing from day 1 to about day 8 (Table 4) during the first and transient stages since less and less incoming radiation energy was used for evaporation and a growing residual energy was applied to heat this soil during this period. The diurnal amplitude of soil temperature was also increasing during this period resulting from increase in the nonuniformity of soil heat properties in response to an increasing difference between soil water content values at the three depths for the corresponding period (Figure 4). The minimum and maximum temperatures reached around 19 and 44°C at the noon of day 8, respectively, showing larger diurnal temperature difference in the soil surface (23°C) than in air (10°C). Because of attenuation of the heat energy transported from the soil surface, the daily amplitudes of the soil temperature fluctuation also decreased with depth. The diurnal variation of temperature presented the similar magnitude and amplitude after day 8 (Table 4), which also may indicate a timing at which the transient stage was over.Table 4
Maximum and minimum values of soil temperature and the corresponding time points.
Day
Min (°C)
Day
Max (°C)
1.0
12.5
1.5
33.2
2.2
16.1
2.5
33.0
3.2
16.3
3.6
34.9
4.2
16.9
4.6
39.1
5.2
17.5
5.5
41.8
6.2
18.0
6.5
43.0
7.2
18.4
7.5
43.6
8.2
18.6
8.5
44.0
9.2
18.8
9.5
44.3
10.2
18.9
10.5
44.5
11.2
19.0
11.5
44.6
12.2
19.1
12.5
44.6
13.2
19.1
13.5
44.7
14.2
19.2
14.5
44.8
15.2
19.2
15.5
44.8
16.2
19.2
16.5
44.9
17.2
19.2
17.5
44.9
18.2
19.2
18.5
44.9
19.2
19.2
19.5
44.9Figure 5
Simulated temporal variations of air temperature and soil temperatures at depths of 0, 2.5, and 5 cm during the soil drying.


### 3.3. Surface Energy Balance Components

Simulated diurnal variations of surface energy balance components are depicted in Figure6. Rn and its two components (Rns and Rnl in (12)) all reached maximum values at noon and both Rn and Rnl had minimum values just before the sunrise, with Rns being values of 0 during nights. Rns dominated Rn during daytime and Rnl dominated it during nighttime. Net radiation Rn and shortwave Rns presented larger amplitudes for the first few days when the soil was relatively wet as compared than those for the rest of days. The relatively higher soil moisture at the soil surface before about day 5 caused a lower soil albedo a (it had a negative relation with the surface water content in (12)); this, in turn, resulted in a smaller reflection of shortwave radiation which caused the higher shortwave Rns and net radiations Rn. With the drying of the soil, the amplitudes in Rns were decreasing until day 4 after which they remained constant (changes in amplitudes are indicated in Figure 7 which shows diurnal maximum values of surface energy balance components) when diurnal changes of the surface water content and its corresponding soil albedo were kept constant for the remaining days (Figure 4). The longwave radiation Rnl presented increasing amplitude with the increasing soil surface temperature at 1 and transient stages and then approached approximately constant value for the rest of days. Rnl was always downward (negative values), suggesting it would decrease the soil temperature by emitting the surface energy into the atmosphere. Net radiation Rn presented a rapid decrease amplitude at stage 1 and a slow decrease one at the transient stage and then kept a relative constant maximum value for the rest time (Figure 7).Figure 6
Simulated diurnal variations in surface energy balance components (namely, shortwave radiation, long radiation, net radiation, sensible heat flux, latent heat flux, and surface heat flux) during a 20-day drying period following a rainfall event.
(a)
(b)Figure 7
Diurnal maximum values of surface energy balance components depicted in Figure4 (for Rns and LE, negative diurnal maximum values are demonstrated).The net radiation flux at the soil surface must accomplish the partitioning between the latent heat, sensible heat, and ground surface heat fluxes (11), depending on the hydraulic and thermal properties of the near-surface zone. The latent heat flux LE was typically high during the early few days when the soil was very moist and there was adequate water in the soil to maintain the high evaporation rate. As shown in Figure 7, the maximum value of LE approached 316 W m−2 at the noon of day 1, which accounted for 68% of the net radiation Rn at same time. LE decreased rapidly with the deceasing soil moisture at stage 1 and presented a very slow decrease at stage 2 (Figure 7). The maximum value of LE approached 55 W m−2 at the noon of day 5, which accounted for only 16% of Rn at same time, and with a constant ration of LE to Rn being 12% at the last six days when the water stored in surface layer was depleted. The sensible heat flux H, which was more related to the changes of surface soil temperature than the surface moisture, showed an inverse trend with LE and increased from very small values during the first two days to a relatively high and constant amplitude after day 8. The ground heat flux G did not present a distinct variation trend with the changes in the soil moisture or soil temperature during stage 1 and transient stage (Figure 7). That was because G was controlled by both the temperature gradient and moisture-dependent thermal conductivity λ of the soil, both of which had a positive relationship with the value of G. Thus, G reached its lowest value at the noon of day 3 when there were relatively lower values of both soil surface moisture and temperature and a highest value at the noon of day 6 when there were relatively higher values of both soil surface moisture and temperature.G presented a slow decreasing value after transient stage.


### 3.4. Soil Water Balance

In this section, variations of soil water balance during the drying process were analyzed to provide an additional way to improve the understanding of the soil water evaporation process. Calculated cumulative infiltration, cumulative evaporation, and change in the net water storage for the whole soil profile during the simulation period are shown in Figure8. The net water storage change was obtained by subtracting the initial water storage value within the 50-cm profile from the water storage values for every node time. Because there was no runoff on the land and no water flowed out from the lower boundary of the soil domain from the result of this simulation, the infiltration water would be totally changed to the storage of water in the soil profile. Therefore, the change in the water storage was due only to the water loss of evaporation from the soil, namely, the water loss in the soil should be equal to the amount of evaporation. Clearly, the water infiltration process began on the rain day (not shown in Figure 8) during which there was a 2-cm precipitation. When this rainfall process terminated at 0 hour of day 1, the amount of cumulative infiltration reached 1.84 cm, which all were contributed to the increased amount of water storage in the profile indicated by the equal values between the cumulative infiltration and net water storage at 0 hour on day 1. At the same time, the cumulative evaporation was equal to 0 before the sunrise. The significant evaporation began before the noon of day 1 and then kept rapidly increasing during each afternoon period until around day 6, which was accompanied by the rapid decrease in the net water storage. At stage 1, it was clearly found that the shapes of “step” on the evaporation or net water storage curves were seen in Figure 8. This step indicated that the evaporation process concentrated on a certain period at each day (e.g., from 10 to 17 hours on day 2) under atmospheric demand conditions, with the little evaporation occurring at the rest of the day. After day 5, cumulated evaporation also kept increasing but with a relatively slow rate until the end of this simulation.Figure 8
Temporal changes in the accumulation infiltration (solid line) and evaporation (dashed line) as well as the net water storage during the simulation period.It is notable that the cumulated infiltration had a small increase from about the day 8 when the transient stage of evaporation terminated. This may be surprising because there was no extra rainfall event to induce the infiltration during these days. Actually, this part of increase in infiltration was due to the condensation of water vapor in atmosphere or at the soil surface. This process occurs when the near-surface soil is extremely dry so that the air humidity exceeds the soil pore humidity and then the water vapor in the atmosphere will diffuse towards soil surface according to Fick’ law and condense there or within the soil. The amount of the increased infiltration due to condensation was about 0.148 cm until the end of day 20. This process may be a relatively important water production mechanism in semiarid or arid regions, especially in the desert areas.


### 3.5. DSL and Evaporation Zone

Once the surface moisture was depleted (i.e., the surface water content dropped too close to the critical value), the soil came in equilibrium with the overlying air and became approximately air-dry. This resulted in formation of DSL at the end of stage 1 evaporation. Figure9 presents the spatiotemporal evolutions of water content, water vapor flux, and liquid water flux in the near-surface zone (the upper 2-cm profile) during the drying period, through which we could find the structure of DSL and its diurnal changes with the soil drying. As shown in Figure 9(a), the critical water content (the minimum value of water content under the drying condition) was about 0.05 m3 m−3, which was also inferred from the temporal change in the water content at 0 cm in Figure 4. There was a distinct line which divided the plot of water content into two zones during the transition and stage 2 evaporation, with the upper zone being the DSL and the lower zone remaining wetter.Figure 9
The spatiotemporal evolutions of (a) water content, (b) water vapor flux, (c) liquid water flux, and (d) evaporation zone during the drying period.
(a)
(b)
(c)
(d)We also could infer the evolution information of the DSL from the diurnal dynamics of vapor and liquid fluxes in the profile; based on that water vapor flux dominated the DSL and liquid water flux dominated the layer below the evaporation zone. As shown in Figures9(b) and 9(c), within the DSL, the liquid water flux was approximately the zero value, while the water vapor flux had a sharp increase within the DSL during daytime hours. From Figures 9(a)–9(c), it was found the distribution of the DSL presented the diurnal pattern, with it being at deeper profiles during daytime hours and shallower profiles during nighttime hours. With the drying of soil, the width of DSL became larger and finally reached a value of about 1.2 cm at the end of simulation.It was also interesting to know the shape and width of the evaporation zone, in which the liquid water changed to the water vapor. This could be found from Figure9(d), in which the distribution of phase change rate was presented. The phase change rate e was calculated from the following expression [4]:(21)e=∂qv∂z+∂θv∂t∗dz,wheredz was the node distance in the numerical simulation. As depicted in this figure, the evaporation zone also displays a diurnal fluctuation with its distribution being more apparent during day time. With the progression of drying, the evaporation zone moved deeper in the soil and its width became narrower.


## 3.1. Soil Water Evaporation (Drying) Process

As shown in (11), the evaporation rate E (equals the total water flux at the soil surface), actually could be expressed as the sum of surface liquid and vapor water fluxes. Figure 2 presents diurnal variations of the two components of actual surface evaporation rate E. The surface liquid water flux qL(0,t) indicates vaporization of water occurring at the soil surface as long as its values are larger than 0, while the positive surface water vapor flux qv(0,t) indicates that water vaporized within the subsurface and its negative values represent that some water vapor was moving downward from the soil surface to deeper profiles. As shown in Figure 2(a), qL(0,t) totally dominated E until about the noon of day 5, indicating that evaporation was occurring at the soil surface. During this period, the soil was wet enough to supply water to accomplish vaporization of water, and soil evaporation was limited only by the supply of energy. This could also be determined from Figure 3 which shows that the actual evaporation rate E equals to the potential evaporation rate during this period. qv(0,t)<0 (Figure 2(b)) is observed around the noon for some days, suggesting that some water vapor produced at the soil surface was moving downward driven by the large downward temperature gradient at noon. Between days 5 and 8 qL(0,t) was the dominant source early in the day and qv(0,t) took over for the rest of the day with some overlap where both qL(0,t) and qv(0,t) occurred simultaneously. During this period, the actual evaporation rate became less than the potential one (Figure 3) during afternoons. Beginning on day 9 values of qL(0,t) were close zero for the all days, and qv(0,t) became the dominant source of evaporation E (Figure 3), suggesting that vaporization of water occurred in the subsurface at most time of the day during this period. It is also notable that a very small negative qL(0,t) occurred around late afternoon after day 8, which resulted from the condensation of water vapor at the soil surface or in air when atmospheric evaporation demand decreased. From day 9 on, E/E0 presented negative values around the noon, which was caused by the negative qv(0,t) during the same period in the result of the downward temperature gradient. In this case, E actually does not represent the rate of evaporation from the soil but indicates that the flux at the soil surface is moving towards deeper layers as the form of water vapor.Figure 2
The liquid water (red solid line), water vapor (red dashed line) water fluxes, and the total flux (black solid line) at the soil surface. The left vertical dashed lines represent timing point of shift of evaporation process from stages 1 and 2, and periods between the two dashed line indicate the transient stage of evaporation process (same below).
(a)
(b)Figure 3
The temporal variation of ratio of actual evaporation rate (E) and potential (E0) rate E/E0 during the simulation period.Thus, based on these results, the stages of evaporation process for the silt soil could be reasonably determined. Stage 1 evaporation was maintained from day 1 to the noon of day 5, indicated byE/E0=1. Thereafter, stage 2 was maintained until the end of this numerical experiment, with a period as a transient stage in which the shift between qL(0,t) and qv(0,t) occurred alternately every day (from days 5 to 8). Stage 2 was also referred as falling rate stage evaporation because during this stage the actual evaporation rate E fell below the potential evaporation rate E0. There was no clear limit between stage 2 and stage 3 evaporation according to data shown in Figures 2 and 3. The so-called stage 3 evaporation, which was characteristic of the low and constant evaporation rate, was not considered to be very important and often merged into stage 2.


## 3.2. Soil Water Content and Temperature

Figure4 presents temporal changes in the volumetric soil water contents at selected depths near the soil surface during a 20-d drying period immediately following a rain event before day 1 (i.e., at day 0). Soil water redistribution and soil drying began at day 1. Water contents near the soil surface decreased rapidly for the first few days, caused by gravity drainage and evaporation imposed by the hot and dry atmospheric conditions (Figure 1). The extent of water loss in the profile declined with depth, with more pronounced decrease occurring at the soil surface (0 cm) where the water content had a significant decrease from the maximum value of 0.43 m3 m−3 at the beginning of day 1 to the minimum value of 0.046 m3 m−3 at the noon of days 5 to 20 during stage 2. It is also notable that the water content near the soil surface (e.g., at 0 and 2.5 cm depths) displayed a distinct diurnal pattern, with it decreasing from late morning through the afternoon and then increasing from afternoon to late morning on the following day. This pattern during evaporation process was consistent with the findings from a number of field measurements [14, 28–31]. This diurnal pattern for surface water content derived from evaporation during daytime and subsequent soil rewetting during nighttime. The consistency between simulated patterns of soil water content in this study and observations in the field conditions supports the reliability of our models on reproducing the evaporative drying process of soils. The amplitude of diurnal pattern decreased with depth, with presenting a small fluctuation at 2.5 cm and no apparent trace for the deeper soil layer, for example, at 5 cm. After day 8, this diurnal pattern of soil water content at the soil surface maintained the same amplitude and magnitude, with maximum values of surface water content maintained at the same value of 0.06 m3 m−3 and the maximum values of surface temperature at around 44.5°C (Table 3).Table 3
Maximum and minimum values of soil water content and the corresponding time points.
Day
Min (m3 m−3)
Day
Max (m3 m−3)
1.6
0.247
1.0
0.43
2.6
0.112
2.2
0.30
3.6
0.058
3.3
0.23
4.6
0.047
4.3
0.15
5.5
0.046
5.3
0.09
6.5
0.046
6.3
0.07
7.5
0.046
7.3
0.07
8.5
0.046
8.2
0.06
9.5
0.046
9.2
0.06
10.5
0.046
10.2
0.06
11.5
0.046
11.2
0.06
12.5
0.046
12.2
0.06
13.5
0.046
13.2
0.06
14.5
0.046
14.2
0.06
15.5
0.046
15.2
0.06
16.5
0.046
16.2
0.06
17.5
0.046
17.2
0.06
18.5
0.046
18.2
0.06
19.5
0.046
19.2
0.06Figure 4
Simulated temporal variations of soil water contents at depths of 0, 2.5, and 5 cm during the soil drying after a rainfall event.Figure5 is the same as Figure 4 but for temporal changes in soil temperatures during the simulation period. It is clear that the soil temperatures at three depths all display a typical diurnal sinusoidal pattern over all the days, with its maximum temperatures occurring around noon (e.g., at 1 pm at 0 cm) and its minimum temperatures occurring just before the sunrise. In the early few days, the relatively uniform and large water content near the soil surface resulted in the uniform and large soil heat capacity (remember that soil heat capacity varies only as function of water content in soil). This resulted in a relatively uniform soil temperature profile for this period, with the daily temperature magnitudes and amplitudes being a little different between the three depths. On the other hand, the initially wet soil conditions combined with high-radiation atmospheric conditions caused a large amount of vaporization of water near the soil surface, causing more incoming solar radiation to be applied to accomplish vaporization rather than heating the soil. These effects caused relatively lower daily temperature magnitudes and amplitudes for the first few days. With the soil drying, both minimum and maximum temperatures in the shallow profile were increasing from day 1 to about day 8 (Table 4) during the first and transient stages since less and less incoming radiation energy was used for evaporation and a growing residual energy was applied to heat this soil during this period. The diurnal amplitude of soil temperature was also increasing during this period resulting from increase in the nonuniformity of soil heat properties in response to an increasing difference between soil water content values at the three depths for the corresponding period (Figure 4). The minimum and maximum temperatures reached around 19 and 44°C at the noon of day 8, respectively, showing larger diurnal temperature difference in the soil surface (23°C) than in air (10°C). Because of attenuation of the heat energy transported from the soil surface, the daily amplitudes of the soil temperature fluctuation also decreased with depth. The diurnal variation of temperature presented the similar magnitude and amplitude after day 8 (Table 4), which also may indicate a timing at which the transient stage was over.Table 4
Maximum and minimum values of soil temperature and the corresponding time points.
Day
Min (°C)
Day
Max (°C)
1.0
12.5
1.5
33.2
2.2
16.1
2.5
33.0
3.2
16.3
3.6
34.9
4.2
16.9
4.6
39.1
5.2
17.5
5.5
41.8
6.2
18.0
6.5
43.0
7.2
18.4
7.5
43.6
8.2
18.6
8.5
44.0
9.2
18.8
9.5
44.3
10.2
18.9
10.5
44.5
11.2
19.0
11.5
44.6
12.2
19.1
12.5
44.6
13.2
19.1
13.5
44.7
14.2
19.2
14.5
44.8
15.2
19.2
15.5
44.8
16.2
19.2
16.5
44.9
17.2
19.2
17.5
44.9
18.2
19.2
18.5
44.9
19.2
19.2
19.5
44.9Figure 5
Simulated temporal variations of air temperature and soil temperatures at depths of 0, 2.5, and 5 cm during the soil drying.


## 3.3. Surface Energy Balance Components

Simulated diurnal variations of surface energy balance components are depicted in Figure6. Rn and its two components (Rns and Rnl in (12)) all reached maximum values at noon and both Rn and Rnl had minimum values just before the sunrise, with Rns being values of 0 during nights. Rns dominated Rn during daytime and Rnl dominated it during nighttime. Net radiation Rn and shortwave Rns presented larger amplitudes for the first few days when the soil was relatively wet as compared than those for the rest of days. The relatively higher soil moisture at the soil surface before about day 5 caused a lower soil albedo a (it had a negative relation with the surface water content in (12)); this, in turn, resulted in a smaller reflection of shortwave radiation which caused the higher shortwave Rns and net radiations Rn. With the drying of the soil, the amplitudes in Rns were decreasing until day 4 after which they remained constant (changes in amplitudes are indicated in Figure 7 which shows diurnal maximum values of surface energy balance components) when diurnal changes of the surface water content and its corresponding soil albedo were kept constant for the remaining days (Figure 4). The longwave radiation Rnl presented increasing amplitude with the increasing soil surface temperature at 1 and transient stages and then approached approximately constant value for the rest of days. Rnl was always downward (negative values), suggesting it would decrease the soil temperature by emitting the surface energy into the atmosphere. Net radiation Rn presented a rapid decrease amplitude at stage 1 and a slow decrease one at the transient stage and then kept a relative constant maximum value for the rest time (Figure 7).Figure 6
Simulated diurnal variations in surface energy balance components (namely, shortwave radiation, long radiation, net radiation, sensible heat flux, latent heat flux, and surface heat flux) during a 20-day drying period following a rainfall event.
(a)
(b)Figure 7
Diurnal maximum values of surface energy balance components depicted in Figure4 (for Rns and LE, negative diurnal maximum values are demonstrated).The net radiation flux at the soil surface must accomplish the partitioning between the latent heat, sensible heat, and ground surface heat fluxes (11), depending on the hydraulic and thermal properties of the near-surface zone. The latent heat flux LE was typically high during the early few days when the soil was very moist and there was adequate water in the soil to maintain the high evaporation rate. As shown in Figure 7, the maximum value of LE approached 316 W m−2 at the noon of day 1, which accounted for 68% of the net radiation Rn at same time. LE decreased rapidly with the deceasing soil moisture at stage 1 and presented a very slow decrease at stage 2 (Figure 7). The maximum value of LE approached 55 W m−2 at the noon of day 5, which accounted for only 16% of Rn at same time, and with a constant ration of LE to Rn being 12% at the last six days when the water stored in surface layer was depleted. The sensible heat flux H, which was more related to the changes of surface soil temperature than the surface moisture, showed an inverse trend with LE and increased from very small values during the first two days to a relatively high and constant amplitude after day 8. The ground heat flux G did not present a distinct variation trend with the changes in the soil moisture or soil temperature during stage 1 and transient stage (Figure 7). That was because G was controlled by both the temperature gradient and moisture-dependent thermal conductivity λ of the soil, both of which had a positive relationship with the value of G. Thus, G reached its lowest value at the noon of day 3 when there were relatively lower values of both soil surface moisture and temperature and a highest value at the noon of day 6 when there were relatively higher values of both soil surface moisture and temperature.G presented a slow decreasing value after transient stage.


## 3.4. Soil Water Balance

In this section, variations of soil water balance during the drying process were analyzed to provide an additional way to improve the understanding of the soil water evaporation process. Calculated cumulative infiltration, cumulative evaporation, and change in the net water storage for the whole soil profile during the simulation period are shown in Figure8. The net water storage change was obtained by subtracting the initial water storage value within the 50-cm profile from the water storage values for every node time. Because there was no runoff on the land and no water flowed out from the lower boundary of the soil domain from the result of this simulation, the infiltration water would be totally changed to the storage of water in the soil profile. Therefore, the change in the water storage was due only to the water loss of evaporation from the soil, namely, the water loss in the soil should be equal to the amount of evaporation. Clearly, the water infiltration process began on the rain day (not shown in Figure 8) during which there was a 2-cm precipitation. When this rainfall process terminated at 0 hour of day 1, the amount of cumulative infiltration reached 1.84 cm, which all were contributed to the increased amount of water storage in the profile indicated by the equal values between the cumulative infiltration and net water storage at 0 hour on day 1. At the same time, the cumulative evaporation was equal to 0 before the sunrise. The significant evaporation began before the noon of day 1 and then kept rapidly increasing during each afternoon period until around day 6, which was accompanied by the rapid decrease in the net water storage. At stage 1, it was clearly found that the shapes of “step” on the evaporation or net water storage curves were seen in Figure 8. This step indicated that the evaporation process concentrated on a certain period at each day (e.g., from 10 to 17 hours on day 2) under atmospheric demand conditions, with the little evaporation occurring at the rest of the day. After day 5, cumulated evaporation also kept increasing but with a relatively slow rate until the end of this simulation.Figure 8
Temporal changes in the accumulation infiltration (solid line) and evaporation (dashed line) as well as the net water storage during the simulation period.It is notable that the cumulated infiltration had a small increase from about the day 8 when the transient stage of evaporation terminated. This may be surprising because there was no extra rainfall event to induce the infiltration during these days. Actually, this part of increase in infiltration was due to the condensation of water vapor in atmosphere or at the soil surface. This process occurs when the near-surface soil is extremely dry so that the air humidity exceeds the soil pore humidity and then the water vapor in the atmosphere will diffuse towards soil surface according to Fick’ law and condense there or within the soil. The amount of the increased infiltration due to condensation was about 0.148 cm until the end of day 20. This process may be a relatively important water production mechanism in semiarid or arid regions, especially in the desert areas.


## 3.5. DSL and Evaporation Zone

Once the surface moisture was depleted (i.e., the surface water content dropped too close to the critical value), the soil came in equilibrium with the overlying air and became approximately air-dry. This resulted in formation of DSL at the end of stage 1 evaporation. Figure9 presents the spatiotemporal evolutions of water content, water vapor flux, and liquid water flux in the near-surface zone (the upper 2-cm profile) during the drying period, through which we could find the structure of DSL and its diurnal changes with the soil drying. As shown in Figure 9(a), the critical water content (the minimum value of water content under the drying condition) was about 0.05 m3 m−3, which was also inferred from the temporal change in the water content at 0 cm in Figure 4. There was a distinct line which divided the plot of water content into two zones during the transition and stage 2 evaporation, with the upper zone being the DSL and the lower zone remaining wetter.Figure 9
The spatiotemporal evolutions of (a) water content, (b) water vapor flux, (c) liquid water flux, and (d) evaporation zone during the drying period.
(a)
(b)
(c)
(d)We also could infer the evolution information of the DSL from the diurnal dynamics of vapor and liquid fluxes in the profile; based on that water vapor flux dominated the DSL and liquid water flux dominated the layer below the evaporation zone. As shown in Figures9(b) and 9(c), within the DSL, the liquid water flux was approximately the zero value, while the water vapor flux had a sharp increase within the DSL during daytime hours. From Figures 9(a)–9(c), it was found the distribution of the DSL presented the diurnal pattern, with it being at deeper profiles during daytime hours and shallower profiles during nighttime hours. With the drying of soil, the width of DSL became larger and finally reached a value of about 1.2 cm at the end of simulation.It was also interesting to know the shape and width of the evaporation zone, in which the liquid water changed to the water vapor. This could be found from Figure9(d), in which the distribution of phase change rate was presented. The phase change rate e was calculated from the following expression [4]:(21)e=∂qv∂z+∂θv∂t∗dz,wheredz was the node distance in the numerical simulation. As depicted in this figure, the evaporation zone also displays a diurnal fluctuation with its distribution being more apparent during day time. With the progression of drying, the evaporation zone moved deeper in the soil and its width became narrower.


## 4. Conclusions

In this paper, we conducted a numerical simulation on the continuously evaporative drying process with 20-d clear sky conditions immediately following a 2-cm rainfall using the HYDRUS-1D software package. Our simulations provided a range of soil conditions from very wet after rainfall to thorough dry. Various stages of evaporation process could be determined form results of simulation. At stage 1, the actual evaporation rate was same as the potential evaporation rate before the noon of day 5 during which the liquid water flux at the soil surface was the sole source of evaporation. During the transient stage between days 6 and 9, the evaporation of the soil came from both the liquid water and vapor fluxes at the soil surface accompanying the ratio of actual to potential evaporation ratesE/E0 less than 1. And after day 9, the water vapor flux became the dominant source of evaporation for stage 2. There was a distinct rewetting of soils during the soil cooling period at stage 2, either because of the liquid water moving upwards from below or condensation of water vapor from the atmosphere. When the soil moisture was depleted, the location of evaporation zone moved into the subsurface and DSL formed immediately above the evaporation zone. Both the DSL and evaporation zone displayed a diurnal pattern during the transition and stage 2 evaporation. The width of DSL became wider with the progressive drying of the soil, while the width of evaporation zone behaved in the inverse trend. The corresponding changes in the soil moisture, soil temperature, surface energy balance, and soil water balance during different stages of evaporation were fully discussed in this paper, which, in turn, may be as some indicators for the division of stages for the evaporation process.



---

*Source: 5892867-2017-11-14.xml*