Simulation of Diurnal Evolution of Evaporation Zone during Soil Drying after Rainfall

: To better understand processes involved in the evaporative drying of the soil, simulations on the dynamics of the evaporation zone, condensation zone, and dry surface layer (DSL) were conducted during a 10-day drying event under diurnal atmospheric conditions. Simulated water contents and soil temperatures matched well with the measured data in the lysimeter. Surface evaporation predominantly occurred during the early period each day, while subsurface evaporation dominated during the remaining part of the day. The evaporation zone presented a distinctly diurnal pattern, moving toward the deeper soil layer during the daytime and back toward the soil surface during the nighttime. The DSL and condensation zone, located immediately above and below the evaporation zone, respectively, also presented diurnal patterns following those of the evaporation zone. As soil drying progressed, both the position of the evaporation zone within the profile and the DSL width exhibited an overall increasing trend, reaching about 4.9 mm by the end of the study period. The occurrence of condensation zones was limited to the daytime when there was a downward surface temperature gradient present. Diurnal patterns observed in both evaporation zones and DSL could potentially be determined by quantifying changes in the near-surface profile’s soil water content, relative humidity, pressure head, and vapor density.


Introduction
Accurate characterization of energy and mass transport in soils during evaporative drying is essential for a better understanding of hydrological processes and many applications in environmental, agricultural, and engineering fields.Bare soil evaporation is a key component both of water and energy balances in semi-arid or raid regions, influencing the development of the atmospheric boundary layer and, therefore, influencing climate change [1][2][3].Understanding the mechanisms underlying evaporation during drying processes holds significant importance for various applications, such as food processing and preservation, the production of ceramics and paper, eye and skin care, and numerous construction activities [4,5].In the field of drying technology, it has been widely recognized that the development and location of drying in a porous medium are dependent on the balance between the liquid water, water vapor, and heat transport processes.
The soil water evaporation is affected by both atmospheric conditions (e.g., air temperature, solar radiation, relative humidity, and wind speed) and soil physical and chemical properties (soil type, pore size, solute composition and content, and thermal and hydraulic conductivities).The rate of evaporation is influenced by complex displacement behaviors involved in mass and energy transfer during soil drying, including liquid water flow, water vapor diffusion, and phase change.Sol water evaporation is often divided into two stages, starting with a relatively high and constant evaporation rate (stage 1) controlled by the Water 2024, 16, 639 2 of 18 land-atmospheric boundary condition, followed by a lower rate (stage 2) controlled by soil hydraulic and thermal properties [6][7][8].During stage 1, water is available at the soil surface, and the evaporation rate is close to the rate of evaporation from open water [9].Most of the surface net radiation is rapidly consumed as latent heat from vaporization during stage 1.Water content at the soil surface decreases with continuous soil drying, and when it reaches a critical value, stage 1 evaporation ends and stage 2 evaporation begins, with the rate of evaporation lower than the potential rate [9].Unlike stage 1, during which evaporation takes place at the soil surface, during stage 2, evaporation occurs at depths just below the soil surface since water at the soil surface is not available anymore [4].This results in the formation of a dry surface layer (DSL), which is easily observed in fields by the color contrast between dry and wet soil.The vaporization of liquid water occurs at the bottom of DSL, which has been confirmed by recent laboratory and numerical studies [10][11][12][13].This shift of the drying front (or evaporation front) from surface to subsurface also has an important impact on the surface energy balance, as the energy used for vaporization must be transferred from the soil surface to the position of subsurface evaporation.The location of drying (or evaporation) is controlled by the balance between the transfer rates of liquid water, water vapor and heat, which has been discussed in a few studies [4,[14][15][16][17].
Several researchers have attempted to better elucidate the dynamics of the evaporation zone during soil drying by using advanced measurement approaches with fine-scale measurement resolution.Heitman et al. [18,19] developed a measurement approach for determining in situ soil water evaporation dynamics by using heat-pulse probes.Following the method developed by Heitman, Deol et al. [16] determined the millimeter-scale subsurface evaporation profiles by using eleven-needle heat pulse probes in the soil column over a drying event.
Since there exists an undetectable zone from the soil surface to the first sensor needle [20], which can cause underestimation of soil evaporation, and a very small (e.g., ≤1.0 mm) sensor spacing near the soil surface is needed to capture vaporization/condensation dynamics within the DSL [21], numerical simulations are often used to evaluate the evolution of evaporative drying processes in the highly dynamic near-surface layer, mainly based on the coupled theory of heat, liquid water, and water transport in the porous media developed by Philp and de Vries (henceforth PDV) and on its modified versions [22][23][24][25][26]. Novak [11] conducted a numerical analysis of evaporation zone dynamics in a bare silt-loam soil under field conditions and observed that the evaporation zone near the soil surface displayed a strong diurnal pattern superimposed on changes associated with the progressive soil drying.The simulation conducted by Novak was based on the PDV model; however, it did not account for the influence of the conventional transfer of soil heat and water on the evaporation process.Therefore, he emphasized the imperative need to continue with a coupled model that incorporates all possible transport mechanisms.Gran [12] conducted evaporation process modeling and observed that the subsurface evaporation rate profile should include the condensation zone located below the evaporation zone, which was not reported in some recent research [11,16,20].
The establishment of the structure and diurnal dynamics of the soil evaporation zone throughout the entire process, from wetting to drying, is essential.This study aims to assess the diurnal dynamics of DSL, evaporation zone, condensation zone, and related variables during soil evaporative drying.To achieve this objective, simulations were conducted on bare field soil under diurnal atmospheric conditions over a 10-day period.The transport mechanisms of heat and water in the soil were described by the extended PDV model, including the movement of liquid water and water vapor driven by both pressure head and soil temperature gradients, the movement of heat by conduction, the convection of sensible heat by liquid water and vapor flow, and the transfer of latent heat by diffusion of water vapor.Simulated values were compared with observed data on lysimeter soil water content and temperature for one year.

Study Site
The observation data used in this study were taken from a lysimeter at the Hanwang Hydrology and Water Resources Experimental Station (HHWES) in the Huaihe River Plain region of China.The strata of the experimental station consist of loose deposits from the Quaternary system, characterized by a surface layer of sandy loam soil.The depth of these deposits generally does not exceed 15 m.Typically, the groundwater level in this area remains within 3 m.The climate at HHWES (34 • 11 ′ N, 117 • 18 ′ E) is a semi-humid temperate continental monsoon with a mean annual temperature and precipitation of 15.5 • C and 750 mm, respectively.Most precipitation (60% to 80%) occurs during the summer from June to September.The surface of the lysimeter was bare during the simulation period, with no need to consider the effect of vegetation.The soil in the lysimeter was homogeneous sandy loam with a dry density of 1.42 g cm −3 and an average particle size of 0.063 mm.
The volumetric water content and matric potential (i.e., pressure head) were monitored by a neutron moisture meter and the tensiometer at depths of 10, 30, 50, 70, and 90 cm, respectively.There was an aluminum alloy tube with a diameter of 44.5 mm in the center of the lysimeter soil column, which was a neutron moisture meter tube.The soil moisture content can be converted based on the readings measured by a neutron meter.A neutron moisture meter employs the neutron scattering technique, which is considered the most accurate technique for measuring the volumetric soil moisture content [27].There was a subsurface observation room around the lysimeter to collect data on surface runoff, groundwater evaporation, and recharge to groundwater for the lysimeter.Within this observation room, the readings from the tensiometer were recorded.To maintain consistency throughout the simulation period, the water table of the lysimeter was regulated by means of a Marriott bottle located inside the observation room, ensuring a constant depth of 100 cm below the soil surface.Regularly observed meteorological variables at HHWES included air temperature, soil temperature, precipitation, solar radiation, wind speed, sunshine hours, and pan evaporation.The soil temperatures were observed at depths of 0, 5, 10, 15, and 20 cm below the soil surface at the same soil type and surface elevation as the lysimeter.In this study, we used these temperature data to represent the temperature condition in the lysimeter.The water content was measured once every 10 days from November to March and every 5 days from April to October during a one-year period, with increased observations before and after rainfall events.Other terms were measured one time, at least once each day, according to the observation rules of the corresponding terms.

Water Flow Equations
Based on the extended PDV theory, the governing equation of the one-dimensional water flow in the soil under variably saturated, non-isothermal conditions can be expressed as [20,22]: where θ L is the volumetric liquid water content (m 3 m −3 ); θ v is the volumetric water vapor content (expressed as an equivalent water content, ρ v θ air /ρ L , m 3 m −3 ), θ air is volumetric air content (m 3 m −3 ), ρ v and ρ L represent the densities of liquid water and water vapor in soil (kg m −3 ) at temperature T ( • C), respectively; q w is the total water flux (m s −1 ); q L and q v 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 flux densities of liquid water (q L ) and water vapor (q v ) driven by pressure head and temperature gradients can be expressed as, respectively: where q TL and q hL are the thermal and isothermal liquid water fluxes (m s −1 ), respectively; q Tv and q hv are the thermal and isothermal water vapor fluxes (m s −1 ), respectively; h is the pressure head (m); D TL (m 2 • C −1 s −1 ) and K (m s −1 ) are the thermal and isothermal hydraulic conductivities for liquid phase fluxes, respectively; and D Tv (m 2 • C −1 s −1 ) and D hv (m s −1 ) are the thermal and isothermal hydraulic conductivities for water vapor fluxes, respectively.The details on four hydraulic conductivities and related parameters can be found in Table 1.The saturated hydraulic conductivity K s is determined during the pre-simulation stage through a calibration process, yielding a value of 20 cm d −1 .The parameters in Table 1 are either known or dependent on temperature and/or water content.The mass continuum Equation ( 1) can be decomposed into two equations, which respectively account for the liquid water and water vapor contents in the soil: ) where e is the depth-dependent evaporation or condensation rate (s −1 ), and E s (=e dz) is the subsurface evaporation rate (m s −1 ) (or the rate of phase change) at each numerical node in the soil column (positive for evaporation and negative for condensation) [10,20].
The soil water retention function, which relates the soil pressure head with the water content, can be described by the van Genuchten model [28]: where θ s is the saturated water content (m 3 m −3 ), θ r is the residual water content (m 3 m −3 ), and α (m −1 ), n, and m (=1 − 1/n) are empirical shape parameters.The van Genuchten function was fitted to the measured water content and pressure head data in the lysimeter, leading to θ r = 0.078 (m 3 m −3 ), θ s = 0.43 (m 3 m −3 ), α = 0.036 (cm −1 ), and n = 1.28.The coefficient of determination is 0.87, showing a relatively good fit.

Heat Transport Equations
The governing equation for the one-dimensional movement of heat in an unsaturated zone is given by [26]: where S h is the heat storage in the soil (J m −3 ), q h is the total heat flux density of the soil (J m −2 s −1 ).S h in the unsaturated soil can be expressed as: where ) are the volumetric heat capacities of dry soil particles, liquid water, and water vapor, respectively; T r is the arbitrary reference temperature ( • C); θ n is the volumetric fraction of solid phase (m 3 m −3 ); and L (=2.495 × 10 9 − 2.247 × 10 6 T) is the volumetric latent heat of vaporization of water (J m −3 ).q h can be described as: where λ is the apparent soil thermal conductivity (W m −1 K −1 ), which can be described as a function of soil water content [29]:

Initial and Boundary Conditions
The simulation was conducted in 1991.To minimize the influence of initial soil moisture conditions on the simulation results, pre-simulation (spin-up) processes were used to obtain initial conditions for the soil profile, which is a common method in land surface models [30].During this process, saturated liquid conductivity K s was determined by minimizing differences between observed and measured water contents at different depths.
The one-dimensional soil domain extended from the soil surface of the lysimeter to the water table, with a vertical length of 100 cm.Boundary conditions for water flow and heat transfer in the soil were determined by surface precipitation, evaporation, runoff, and heat fluxes.For simulating these processes, atmospheric forcing data, including wind speed, air temperature, sunshine hours, precipitation, and air humidity in 1991, were utilized (Figure 1).More specifically, surface boundary conditions for the water and energy equations are respectively given by [31,32]: where q w (0,t) is the total water flux across the soil surface (m s −1 ), P represents the rate of precipitation (m s −1 ), R represents the rate of surface runoff (m s −1 ), and E is the surface evaporation rate (m s −1 ) from the soil.R n 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 ), and G is the surface heat flux density (W m −2 ).R n and G are positive downward, while H and LE are positive upward [26].R n is defined as: where R ns and R nl are net shortwave and longwave radiations (W m −2 ), respectively, a is soil albedo (unitless), S t 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 [33]: where C a is the volumetric heat capacity of air (=1200 J m −3 K −1 ), T s and T a are the temperature ( • C) of the soil surface and the air, respectively, and r h denotes 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 ): where r v is the aerodynamic resistance to water vapor flow (s m −1 ), and r s is the soil surface resistance to water vapor flow that acts as an additional resistance along with the aerodynamic resistance (s m −1 ).
flux density (W m −2 ).Rn and G are positive downward, while H and LE are positive upward [26].Rn is defined as: 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 [33]: where Ca is the volumetric heat capacity of air (=1200 J m −3 K −1 ), Ts and Ta are the temperature ( o C) of the soil surface and the air, respectively, and rh denotes 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 ): 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 ).The details on surface water and energy balance components used in Equations ( 13)-( 15) can be found in Simunek et al. [34].To represent the constant groundwater table, a zero pressure head is set as the lower boundary condition for water.As for the heat transport equation, the lower boundary is expressed as a function of soil temperature, The details on surface water and energy balance components used in Equations ( 13)-( 15) can be found in Simunek et al. [34].To represent the constant groundwater table, a zero pressure head is set as the lower boundary condition for water.As for the heat transport equation, the lower boundary is expressed as a function of soil temperature, which varies with the day of the year.This function is obtained by fitting the measured temperature data at different depths (from 0 to 320 cm) from HHWES over many years (from 1994 to 1997) to an analytic solution for soil temperature at a specific time and depth [35].
The governing equations of water and heat flow subject to initial and boundary conditions were numerically solved by the HYDRUS-1D code using the finite difference method [34].The HYDRUS-1D has been widely used in a variety of studies on heat and water flow in the unsaturated zone [20,26,36,37].The 100 cm soil profile was divided into 501 elements, with thickness increasing exponentially from the soil surface to the water table.The time step was automatically varying between the minimum and maximum time steps of 1 × 10 −6 and 0.05 days, respectively.

Comparisons between Observed and Simulated Soil Water Contents and Temperatures
In this section, the coupled water and heat transport model was tested by making comparisons between observed and simulated temporal variations in the soil water content and temperature at various depths at HHWES in 1991.The temporal variations in observed and simulated water contents at 10, 30, 50, 70, and 90 cm for the whole year are presented in Figure 2. The simulated soil water contents well matched the observed values at all depths and well captured the important trends, such as rapid increases in the soil water content due to rainfall events and the gradual decrease processes during the subsequent soil drying processes.The simulated and observed water contents both presented more dynamic variations in the upper profile (e.g., at depths of 10 and 30 cm) due to the stronger influence from the atmospheric boundary condition.Conversely, the soil water near the water table (e.g., at depths of 70 and 90 cm) exhibited less variability and remained relatively constant throughout the simulation period because of the impact of the constant water table condition in the lysimeter.
The simulated soil temperatures also closely followed the observed values at depths of 0, 5, 10, 15, and 20 cm during the whole simulation period (Figure 3).Similarly, the soil temperature exhibited more dynamic variations near the soil surface, where it was significantly influenced by atmospheric conditions, displaying a similar trend to that of air temperature (Figure 1).As depth increased, the amplitude of soil temperature decreased due to attenuated energy transport from the soil surface.which varies with the day of the year.This function is obtained by fitting the measured temperature data at different depths (from 0 to 320 cm) from HHWES over many years (from 1994 to 1997) to an analytic solution for soil temperature at a specific time and depth [35].
The governing equations of water and heat flow subject to initial and boundary conditions were numerically solved by the HYDRUS-1D code using the finite difference method [34].The HYDRUS-1D has been widely used in a variety of studies on heat and water flow in the unsaturated zone [20,26,36,37].The 100 cm soil profile was divided into 501 elements, with thickness increasing exponentially from the soil surface to the water table.The time step was automatically varying between the minimum and maximum time steps of 1 × 10 −6 and 0.05 days, respectively.

Comparisons between Observed and Simulated Soil Water Contents and Temperatures
In this section, the coupled water and heat transport model was tested by making comparisons between observed and simulated temporal variations in the soil water content and temperature at various depths at HHWES in 1991.The temporal variations in observed and simulated water contents at 10, 30, 50, 70, and 90 cm for the whole year are presented in Figure 2. The simulated soil water contents well matched the observed values at all depths and well captured the important trends, such as rapid increases in the soil water content due to rainfall events and the gradual decrease processes during the subsequent soil drying processes.The simulated and observed water contents both presented more dynamic variations in the upper profile (e.g., at depths of 10 and 30 cm) due to the stronger influence from the atmospheric boundary condition.Conversely, the soil water near the water table (e.g., at depths of 70 and 90 cm) exhibited less variability and remained relatively constant throughout the simulation period because of the impact of the constant water table condition in the lysimeter.The simulated soil temperatures also closely followed the observed values at depths of 0, 5, 10, 15, and 20 cm during the whole simulation period (Figure 3).Similarly, the soil temperature exhibited more dynamic variations near the soil surface, where it was significantly influenced by atmospheric conditions, displaying a similar trend to that of air temperature (Figure 1).As depth increased, the amplitude of soil temperature decreased due to attenuated energy transport from the soil surface.

3.
Temporal in observed and simulated soil temperatures depths of 10, and 20 in 1991.

Surface Energy Balance
In the we this coupled and water transport model further dynamic processes the near-surface layer during drying cycle.

Surface Energy Balance Components
In the subsequent sections, we used this coupled heat and water transport model to further evaluate dynamic processes in the near-surface layer during a drying cycle.Specifically, two rainfall events occurred on days 221 and 223 in 1991, resulting in a cumulative precipitation amount of 2.53 cm.This was followed by a subsequent period of ten consecutive days without any rainfall from days 224 to 234 event with day 224 in the year of study.Figure 4 demonstrates the temporal variation in calculated surface energy fluxes during the 10-day drying event, and Table 2 shows the values of surface energy balance components at noon every day.A decrease in surface albedo resulting from moist soil conditions after rainfall could lead to a reduction in reflected shortwave radiation (according to Equation ( 13)), which eventually resulted in the large value of the net radiation flux during the early period of soil drying (e.g., days 0 and 1).With the drying of soil, net radiation values for days 2 to 10 were primarily influenced by daily sunshine hours (i.e., incoming solar radiation), exhibiting slightly lower values compared to those days 0 and 1, with significantly lower values during days 5 to 7 due to overcast weather conditions.Sensible heat and ground heat fluxes were considerably smaller during the early periods after rainfall since more radiation energy was used to vaporize the liquid water and generally increased with the soil temperature or temperature gradients (Figure 5).The latent heat flux presented a general decreasing trend as the soil progressively dried.

Surface Liquid and Vapor Fluxes
The temporal variations in the calculated surface liquid water flux qL0 and vapor flux qv0 during the 10-day period are illustrated in Figure 6.By observing the temporal variations in qL0 and vapor flux qv0 in this figure, the periods during which evaporation occurred at the soil surface or in the subsurface could be determined.qL0 represents the surface evaporation process in which the liquid water vaporized at the soil surface, while qv0 indicates the subsurface evaporation process in which the liquid water vaporized below the soil surface and the produced water vapor arrived at the soil surface by diffusion.The surface evaporation rate E (the sum of qL0 and qv0) was maximum around noon each day,

Surface Liquid and Vapor Fluxes
The temporal variations in the calculated surface liquid water flux qL0 and vapor flux qv0 during the 10-day period are illustrated in Figure 6.By observing the temporal varia- Simulated soil temperatures at soil surfaces and depths of 1 and 5 cm are shown in Figure 5.The maximum surface temperature occurred around noon (e.g., at 1 p.m.), and the minimum temperature Temperatures at shallower depths capacity and significant evaporation of the soil surface, both of which resulted in relatively lower daily temperature amplitudes for the first few days.The overcast weather conditions for days 5 to 7 also contributed to lower soil temperature amplitudes.

Surface Liquid and Vapor Fluxes
The temporal variations in the calculated surface liquid water flux q L0 and vapor flux q v0 during the 10-day period are illustrated in Figure 6.By observing the temporal variations in q L0 and vapor flux q v0 in this figure, the periods during which evaporation occurred at the soil surface or in the subsurface could be determined.q L0 represents the surface evaporation process in which the liquid water vaporized at the soil surface, while q v0 indicates the subsurface evaporation process in which the liquid water vaporized below the soil surface and the produced water vapor arrived at the soil surface by diffusion.The surface evaporation rate E (the sum of q L0 and q v0 ) was maximum around noon each day, consistent with the temporal behavior of latent heat flux (Figure 4), and minimum during the nighttime.The evaporation rate E presented a relatively high value in the first few days when the soil was very moist near the soil surface, with its values decreasing with the decreasing surface water content during the soil drying.Before day 1.5, evaporation was dominated by q L0 , suggesting that vaporization of water was occurring only at the soil surface before the subsurface evaporation began.From days 1.5 to 8.5, surface evaporation was dominant during the early period each day during which soil was moist enough to provide water to the soil surface, and subsurface evaporation (q v0 ) dominated during the remaining part of the day, with some overlapping when both q L0 and q v0 occurred simultaneously.The negative value of q v0 during some periods each day represented that the water vapor generated at the soil surface would diffuse back into the soil by the downward temperature gradient [11,20].The rapid increase in q v0 from a small negative to a large positive value essentially matched the decrease in q L0 from a large positive value to a value close to 0, indicating that the transition from surface to subsurface evaporation was occurring.For days 1 to 3, there were apparently small increases in q L0 during nightfall and evening (e.g., on day 1.8 in Figure 6), resulting from the evaporation of the soil by the water from below since the evaporation was reduced by the decreasing radiation condition after the sunset.With the progressive drying of the soil, the water content in the near-surface profile decreased continuously, resulting in more evaporation occurring in the subsurface.q L0 was close to 0 for all hours on the last day, suggesting that all the evaporation occurred in the subsurface.
Water 2024, 16, x FOR PEER REVIEW 12 of 19 a value close to 0, indicating that the transition from surface to subsurface evaporation was occurring.For days 1 to 3, there were apparently small increases in qL0 during nightfall and evening (e.g., on day 1.8 in Figure 6), resulting from the evaporation of the soil by the water from below since the evaporation was reduced by the decreasing radiation condition after the sunset.With the progressive drying of the soil, the water content in the nearsurface profile decreased continuously, resulting in more evaporation occurring in the subsurface.qL0 was close to 0 for all hours on the last day, suggesting that all the evaporation occurred in the subsurface.

Evaporation Zone Dynamics
The subsurface evaporation or condensation rate, Es, represents the phase change rate between liquid water and water vapor within the soil [20,38].Note that the subsurface evaporation rate Es was different from the surface evaporation rate E. The depth interval where the evaporation or condensation occurred would be referred to as the evaporation Diurnal variations in simulated surface liquid water q L0 and vapor fluxes q v0 (as two components of surface evaporation rate E s ) during the soil drying.

Evaporation Zone Dynamics
The subsurface evaporation or condensation rate, E s , represents the phase change rate between liquid water and water vapor within the soil [20,38].Note that the subsurface evaporation was different from the depth interval where be Equation ( 5), E from (cm d −1 ): where the subscripts i and j represent the indexes of depth and time steps, respectively.The temporal and spatial variations in the subsurface evaporation rate E s , liquid water flux (q L ), and water vapor flux (q v ) in the near-surface profile (the 1 cm depth) during the 10-day period are illustrated in Figure 7, generated using MATLAB software (version number: 2012a).Note that the positive values for q L and q v indicate the upward movements of liquid water and water vapor fluxes, and vice versa.The evaporation zone can be clearly found in Figure 7a (light blue area), as indicated by the dense isolines of the phase change rate.The evaporation zone divided the profile of soil water flux into two different parts (Figure 7b,c).The liquid water flux q L below the evaporation zone presented a relatively large value, while the value of q L above it was very close to 0, as indicated by the dense and sparse isolines for q L below and above the evaporation zone, respectively (Figure 7b).In contrast, the vapor flux showed a relatively large value above the evaporation zone and a value close to 0 below it, as indicated by the dense and sparse isolines for q v above and below the evaporation zone, respectively (Figure 7c).The liquid water in the deeper profile by q L moved upward to the evaporation zone, where it totally changed to water vapor.The water vapor produced within the evaporation zone by q v diffused upward from the evaporation zone to the soil surface, from which it eventually entered the atmosphere by turbulent diffusion.The diurnal evolution of subsurface evaporation rate E s in the profile corresponded to the temporal trends of q L and q v for all days.The evaporation zone initially formed just below the soil surface on day 1.5 (Figure 7a), when subsurface evaporation just began (Figure 6).Between days 1 and 9, the evaporation zone moved toward the deeper layer during the daytime, and when it reached the deepest location in the profile around the sunset, it gradually moved back toward the soil surface during the nighttime.The location of the evaporation zone from days 1 to 9 presented a continuous increase with the increasing drying of the soil, except for days 5 to 7 under overcast weather conditions.The temporal snapshots of the 1 cm profile of the simulated phase change rate E s during a one-day cycle (from day 8.3 to 9.2) are presented in Figure 8, in which positive and negative values indicate rates of evaporation and condensation, respectively.The evaporation zone was typically narrow, showing a nearly symmetric shape around its peak value.The width of the evaporation zone was on the millimeter scale, with a width of around 2 mm for days 8.3 to 9.2.During the daytime, the peak value of the subsurface evaporation rate at different times increased with the increasing drying of the soil, and the location of the peak also became deeper and deeper in the profile (Figure 8a).However, this trend was totally reversed during the nighttime (Figure 8b), with a decreasing peak value and its depth with time.At day 8.3, the vaporization of water occurred at the soil surface, as indicated by the peak value of the evaporation zone at the soil surface.The temporal evolution of the peak of the subsurface evaporation rate and its location for the 10-day period are depicted in Figure 9.It was clear that the peak and its corresponding position in the profile both displayed a strong diurnal pattern, with the peak decreasing and its depth generally increasing with the progressive soil drying.
The condensation rate was not easily found in Figure 7a because of the smaller value compared with the rate of evaporation.However, it was clearly found in Figure 8a.The condensation zone is always located immediately below the evaporation zone during the daytime from day 8.3 to 8.7 (hours 7.2 to 16.8 on day 8), with its width presenting much wider than it is for the evaporation zone.The magnitude of the condensation rate (absolute value) at a time point decreased with depth from its peak.For instance, the condensation rate presented a peak value of −0.013 cm d −1 at a depth of 0.21 cm on day 8.4, with a value) at a depth 0 at during the daytime, also increasing with the drying of soil.It could be inferred from Figure 8a,b that the condensation process occurred mainly during the daytime.Actually, the vaporization of liquid water within the evaporation zone could produce a peak of water vapor density in the profile.Some water vapor would move downward, driven by the downward temperate gradient, and eventually change to liquid water when the relative humidity value in the soil pore approached the saturation value of 1 (Figure 10b).This process resulted in the condensation process occurring mainly during the daytime when the temperature gradient was downward in the near-surface layer, while it nearly did not occur during the nighttime when the temperature gradient was upward (Figure 8b).The condensation rate was not easily found in Figure 7a because of the smaller value compared with the rate of evaporation.However, it was clearly found in Figure 8a.The condensation zone is always located immediately below the evaporation zone during the daytime from day 8.3 to 8.7 (hours 7.2 to 16.8 on day 8), with its width presenting much wider than it is for the evaporation zone.The magnitude of the condensation rate (absolute value) at a time point decreased with depth from its peak.For instance, the conden-  DSL formed above the evaporation zone 9 shows variations in the peak value of the subsurface evaporation rate and its location during the period of soil drying.The peak depth of the subsurface evaporation rate can be approximated as the width of the DSL.The diurnal pattern of DSL followed the temporal of the evaporation zone (Figure 7a), with its width becoming wider with the increasing drying of the soil during the daytime and thinner with the rewetting of the soil during the nighttime.The width of DSL reached about 4.9 mm on the last day of the simulation.This relatively small width of DSL suggested that the downward development of DSL was eventually restricted by the continuous supply of water from the shallow water table in our study (the depth of 100 cm below the soil surface).This also indicated the role of the shallow water table in the control of the development of soil drying, and it was reasonably expected that the width of DSL would become wider when the water table became deeper.
Figure 10 shows diurnal variations in soil water content, relative humidity, pressure head (always negative), water vapor density, and temperature in the near-surface profile during the simulation period, all of which were closely associated with the evaporative drying process.The evolutions of soil water content, relative humidity, pressure head, and water vapor density in the profile were apparently consistent with diurnal patterns of the evaporation zone and DSL.Water content presented relatively small and constant values (e.g., 0.1 m 3 m −3 ) within the DSL and a relatively high value below the evaporation zone (e.g., 0.2 m 3 m −3 ) during the daytime and showed a sharp gradient within the evaporation zone (Figure 10a).Relative humidity (RH) within the soil pore during the daytime remained significantly lower than unity within the DSL (Figure 10b), whereas it approached the saturation value of unity below the evaporation zone.The pressure head displayed a diurnal pattern similar to relative humidity, with substantial negative values within the DSL during the daytime and smaller negative values below the evaporation zone (Figure 10c).The density of water vapor above the evaporation zone was relatively low (e.g., 0.02 kg m −3 ), but notably higher below it (e.g., 0.06 kg m −3 ) during daytime hours (Figure 10d).These findings suggest that evaluating diurnal patterns of soil water content, relative humidity, pressure head, and water vapor density using numerical or measurement methods could potentially provide insights into the dynamics of the evaporation zone.We also examined variations in soil temperature in the near-surface layer (Figure 10e), and it was found that the temperature distribution in the profile was not clear enough to determine the drying processes in the highly dynamic near-surface layer.

Conclusions
The process of soil drying by evaporation in field conditions involves complex mass and energy transfer processes, including the phase change between liquid water and water vapor and the associated energy change with latent heat.To gain a better understanding of the dynamics of this process, we conducted a simulation for a 10-day drying event under atmospheric boundary conditions using a coupled model for heat, liquid water, and water vapor.Our simulated values for soil water content and temperature matched well with observed values from HHWES in 1991.We then used our validated model to simulate the dynamics of the evaporation zone, condensation zone, and DSL during this period.All three zones presented distinct diurnal patterns as the soil progressively dried out.Surface and subsurface evaporation alternated throughout each day; surface evaporation was dominant early on, while subsurface evaporation dominated later in the day.During one-day cycles, we observed that the location of the evaporation zone trended towards deeper profiles during daytime hours before returning to the soil surface at night.The DSL and condensation zones are located above and below the evaporation zone, respectively.As soil continued to dry out over time, the location of the evaporation zone in the profile and the width of DSL both presented the general increase trend.The width of DSL reached only about 4.9 mm on the last day of the simulation, restricted by the continuous supply of water from the shallow water table in our study.Under the same atmospheric conditions, the width of DSL should be related to the depth of the groundwater table.The condensation zone occurred only during the daytime when the soil temperature gradient was downward.The dynamics of the evaporation zone and DSL in the near-surface layer were potentially

Figure 1 .
Figure 1.The observed daily atmospheric data at HHEMS in 1991 included maximum and minimum air temperatures, air relative humidity, wind speed, sunshine hours, and rainfall.

Figure 1 .
Figure 1.The observed daily atmospheric data at HHEMS in 1991 included maximum and minimum air temperatures, air relative humidity, wind speed, sunshine hours, and rainfall.

Figure 2 .
Figure 2. Temporal variations in observed and simulated soil water contents at depths of 10, 30, 50, 70, and 90 cm in the 100 cm profile above the water table in 1991.Water 2024, 16, x FOR PEER REVIEW 9 of 19

Figure 4 .
Figure 4. Diurnal variations in simulated surface net radiation, sensible heat, latent heat, and ground heat fluxes during the 10-day soil drying.

Figure 4 .
Figure 4. Diurnal variations in simulated surface net radiation, sensible heat, latent heat, and ground heat fluxes during the 10-day soil drying.

Figure 4 .
Figure 4. Diurnal variations in simulated surface net radiation, sensible heat, latent heat, and ground heat fluxes during the 10-day soil drying.

Figure 6 .
Figure 6.Diurnal variations in simulated surface liquid water qL0 and vapor fluxes qv0 (as two components of surface evaporation rate Es) during the soil drying.

Figure 6 .
Figure6.Diurnal variations in simulated surface liquid water q L0 and vapor fluxes q v0 (as two components of surface evaporation rate E s ) during the soil drying.

Figure 7 .
Figure 7. Temporal and spatial evolutions in the subsurface phase change rate (E s ) (a), liquid water flux (q L ) (b), and vapor flux (q v ) (c) in the near-surface soil profile during the soil drying.

Figure 8 . 19 Figure
Figure 8.The estimated phase change rate profiles at different times during the daytime (a) and the nighttime (b) from day 8.3 to 9.3.Water 2024, 16, x FOR PEER REVIEW 15 of 19

Figure 9 .
Figure 9. Temporal variations in the peak value of the subsurface evaporation rate profile and the width of the DSL occurred during the soil drying.

Figure 10 .
Figure 10.Temporal and spatial evolutions in soil water content (a), relative humidity (b), pressure head (c), water vapor density (d), and temperature (e) in the near-surface layer during the soil drying.

Table 1 .
Hydraulic conductivities and related parameters.

Table 2 .
Values of surface energy balance components at noon every day (W m −2 ).