Validating coupled flow theory for bare‐soil evaporation under different boundary conditions

Evaporation from bare soil is an important hydrological process and part of the water and energy balance of terrestrial systems. Modeling bare‐soil evaporation is challenging, mainly due to nonlinear couplings among liquid water, water vapor, and heat fluxes. Model concepts of varying complexity have been proposed for predicting evaporative water and energy fluxes. Our aim was to test a standard model of coupled water, vapor, and heat flow in the soil using data from laboratory evaporation experiments under different boundary conditions. We conducted evaporation experiments with a sand and a silt loam soil and with three different atmospheric boundary conditions: (i) wind, (ii) wind and short‐wave radiation, and (iii) wind and intermittent short‐wave radiation. The packed soil columns were closed at the bottom (no water flux) and instrumented with temperature sensors, tensiometers, and relative humidity probes. We simulated the evaporation experiments with a coupled water, vapor, and heat flow model, which solves the surface energy balance and predicts the evaporation rate. The evaporation dynamics were predicted very well, in particular the onset of stage‐two evaporation and the evaporation rates during the stage. A continuous slow decrease of the measured evaporation rate during stage‐one could not be described with a constant aerodynamic resistance. Adding established soil resistance parametrizations to the model significantly degraded model performance. The use of a boundary‐layer resistance, which takes into account the effect of point sources of moisture, improved the prediction of evaporation rates for the sandy soil, but not for the silt loam.

energy, and phase change (Or et al., 2013;Vanderborght et al., 2017).It is a non-isothermal process and its proper description requires the consideration of the coupling among these fluxes in the soil (Iden, Blöcher, et al., 2021;Novak, 2019) and between soil and atmosphere.Modeling soil evaporation is important for large-and small-scale research and application, covering topics like understanding water and energy fluxes between the ground and the atmosphere (Bittelli et al., 2008), agricultural water management (Allen et al., 1998), and validating laboratory methods for determining soil hydraulic properties (SHPs) (Iden et al., 2019;Iden, Diamantopoulos, et al., 2021).
The main theory of modeling coupled water, vapor, and heat flow in porous media was presented by Philip and de Vries (1957) (PdV hereinafter).In PdV, the total fluxes of liquid water and water vapor are both composed of an isothermal and a thermal flux.Furthermore, liquid water and water vapor are assumed to be in local phase equilibrium (Novak, 2019).In its original form, PdV neglects several processes (Novak, 2010).A good overview about the development of the PdV theory is given in Milly (1982), and extensions have been put forward by Milly (1982) and Nassar and Horton (1997).A comprehensive review of different flow model concepts for evaporation and their assumptions can be found in Vanderborght et al. (2017).They range from coupled equations for non-isothermal two-phase flow with an explicit coupling across the soil-atmosphere interface to the isothermal Richards or Richardson-Richards equation (Richards, 1931;Richardson, 1922), and use of effective descriptions of soil-atmosphere exchange using resistance formulations.
In addition to the selection of an appropriate flow model to simulate evaporation, different approaches exist to choose the right boundary conditions for water and heat flow when modeling the exchange of water vapor across the soil-atmosphere interface.A correct representation of the interface must be able to describe the different stages of bare-soil evaporation with their different evaporation rate characteristics.Although the evaporation rate during stage-one equals the potential rate controlled by the atmospheric demand and is rather constant under constant atmospheric conditions, it drops during stagetwo, in which the rate is controlled by soil hydraulic resistance to upward liquid water flow to the dry, vapor-dominated surface layer (Saravanapavan & Salvucci, 2000).
The surface energy balance is an adequate option to assign dynamic boundary conditions to water and heat flow, because it includes the latent heat flux and allows to solve for the soil surface temperature.Usually, the evaporation rate is calculated from the water vapor pressure difference between the soil surface and the atmosphere at some reference height and an effective vapor flow resistance (Bonan, 2015;Shuttleworth, 2012).There are various resistance formulations, all of which include an aerodynamic resistance that describes the turbulent exchange of water vapor in the lower atmosphere.Although the use of an additional surface resistance for a vegetated surface represents the resistance to vapor flow imposed by the plant canopy (Allen et al., 1998;Bonan, 2015), its use is less easy to justify when modeling bare-soil evaporation.In some model applications, surface or soil resistance is used to describe the water vapor diffusion resistance of the dry topsoil (Fuchs & Tanner, 1967).However, such use is not warranted in process models in which water and vapor flow in the soil

Core Ideas
• Evaporation experiments under different atmospheric conditions were modeled with coupled water-vapor-heat flow.• Model predicted evaporation dynamics well, but slightly falling rates during stage-one could not be modeled.• Parameterization of soil hydraulic properties accounts for sorbed water and film flow.• Unjustified empirical soil resistance parameterizations degraded model performance.
• Use of a boundary-layer resistance described the decrease during stage-one for sand, but not for loam.
are already explicitly modeled, that is, models built upon PdV theory, because the effect of the vapor-dominated dry surface layer in the top of the soil is considered twice when the resistance is used (an issue of "double counting, " Novak, 2012;Vanderborght et al., 2017).In agreement with this critique, the inclusion of a soil surface resistance in such models has been reported to cause the onset of stage-two to occur too early and the decrease of the evaporation rate during stage-two to be too rapid (Vanderborght et al., 2017).Examples for such "out-of-context" applications of the surface resistance concept include, but are not limited to, Bittelli et al. (2008) and Saito et al. (2006).Haghighi et al. (2013) derived a boundary-layer resistance parameterization that takes into account the effect that in macroscopically wet soils, three-dimensional vapor diffusion in the boundary layer occurs in the immediate vicinity of individual water-saturated pores.Isolated wet patches on the soil surface can lead to similar effects that limit the evaporation rate by the mass exchange from the soil surface through the boundary layer to the atmosphere and not only by the supply of liquid water from deeper soil layers.The basis for this model concept was established by Suzuki and Maeda (1968) and Schlünder (1988).When the soil dries out, the distance between water-filled areas increases and eventually becomes too large to maintain the supply of water vapor.In contrast to the empirical soil resistances mentioned earlier, adding a boundary-layer resistance that incorporates the effect of an uneven moisture distribution at the soil surface on the evaporation rate is consistent with models describing water and vapor flow in the soil.
Experimental work on bare-soil evaporation often focuses on monitoring soil water content and soil temperature, but this limits the ability to determine the physical transport properties of the soil, particularly its hydraulic properties.This is Vadose Zone Journal especially true for unsaturated soil hydraulic conductivity, which is generally difficult to measure but exerts a strong influence on soil evaporation (Fetzer et al., 2017;Iden, Blöcher, et al., 2021;Vanderborght et al., 2017).Goss and Madliger (2007), Yamanaka et al. (1997Yamanaka et al. ( , 1998)), and Yamanaka and Yonetani (1999) used vapor pressure/relative humidity measurements in the soil to investigate soil water dynamics during evaporation.Iden, Diamantopoulos et al. (2021) used tensiometers and relative humidity sensors to measure soil water pressure head and utilized the resulting data to identify SHP from almost full water saturation to air dryness.Iden, Blöcher et al. (2021) and Iden, Diamantopoulos, et al. (2021) showed that water flow in incompletely filled pores (films, corners) plays an important role for the upward movement of liquid water during evaporation and can extend the duration of stage-one and increase the evaporation rate during stage-two.The role of film and corner flow in the drying of porous media has been highlighted earlier by Prat (2007Prat ( , 2011)).
Both field data and laboratory measurements have been used to validate bare-soil evaporation models and parameterizations (Bittelli et al., 2008;Li et al., 2019;Novak, 2010Novak, , 2016)).Examples range from (i) using field data of baresoil evaporation to compare simulations with different soil surface resistance parametrizations (Bittelli et al., 2008), to (ii) experimental data obtained under controlled atmospheric conditions in a wind tunnel to evaluate model concepts (Li et al., 2019), and (iii) studies using weighable lysimeters (Schneider et al., 2021).Saravanapavan and Salvucci (2000) stated that the adequacy of the PdV theory is still under scrutiny.Although PdV has never really been tested, that is, without calibrating inputs, it has been used to describe the drying of other porous media.Heitman et al. (2008) showed through measurement that the dynamics of the dry surface layer are well described by PdV.Nevertheless, the lack of sufficient experimental data, both in the field and in the laboratory, still hinders model selection and a validation of models (Li et al., 2019;Smits et al., 2011).Common limitations include using destructive sampling to estimate soil water content, which does not allow continuous monitoring, and limited field data that do not allow to properly define SHP.Although there is consensus that theory should be validated under controlled transient conditions, Heitman et al. (2008) criticized that many lab experiments are conducted under a limited range of boundary conditions and often under constant atmospheric forcing.Accurate determination and parametrization of SHP are a prerequisite for adequate modeling of soil evaporation and validation of evaporation models.
In this contribution, we present and evaluate results of evaporation experiments conducted under different boundary conditions.We measured transient evaporation rates, soil temperature, and soil water pressure head for a sand and a silt loam and for three atmospheric boundary conditions: (i) wind, (ii) wind and constant short-wave radiation, and (iii) wind and intermittent short-wave radiation.The SHPs were determined from the experiments without applied short-wave radiation by inverse modeling (Iden, Diamantopoulos, et al., 2021) and include water sorption in dry soil and hydraulic conductivity in partially filled pores (film and corner flow).Apart from the aerodynamic resistance to vapor flow, no further model fitting is involved for predicting the evaporation dynamics of the four experiments with applied short-wave radiation (two soils: constant radiation and intermittent radiation).As the experiments are non-isothermal, the adequate flow model to simulate the experiments is a coupled liquid water, vapor, and heat flow model based on a state-of-the-art formulation of the PdV theory, which solves the surface energy balance and predicts the evaporation rates.The objective of this study is to compare predictions by the coupled water and heat flow model with measurements of transient evaporation rates, soil temperature, and soil water pressure head under different atmospheric boundary conditions.We hypothesize that a constant aerodynamic resistance is suitable for the prediction of the evaporation dynamics.We then test the effect of two frequently applied soil surface resistance parameterizations, which are often incorrectly incorporated within PdV, and investigate the implications.Finally, we test whether the prediction of evaporation rates benefits from accounting for an additional boundary-layer resistance derived by Haghighi et al. (2013) which incorporates the effect of point sources of moisture on the evaporation rate.

Coupled water, vapor, and heat flow
In this section, the equations of the coupled water, vapor, and heat flow model are briefly explained.The flow model of liquid water and vapor is an extended version of the Richards equation, in which isothermal and non-isothermal diffusions of water vapor are included (Novak, 2016;Vanderborght et al., 2017).The governing equation for the flow of liquid water and water vapor in one spatial dimension is where θ tot [-] is the total volumetric water content, the sum of the volumetric liquid water θ l  et al. (2006) but differs from their Equation (5) because it does not contain a thermal liquid flux (TLF) and the respective thermal liquid conductivity.In the original PdV theory (Philip & de Vries, 1957), the flux equation for liquid water used gradients of volumetric water content, temperature, and gravitational potential.As pointed out by Prunty (2009), a TLF need not be included if the water flux is formulated in the gradient of pressure head rather than water content, as is the case in Equation (1).It is noteworthy that Milly (1982) did not include a TLF for capillary water in his h-based formulation of PdV (he did include one for adsorbed water, but this term is controversial, see Prunty, 2009).In fact, he pointed out that the inclusion of the TLF in the h-based formulation is "incorrect" and "erroneously doubles the actual effect of temperature gradients on liquid flow" (Milly, 1982, page 497).
Note that the temperature effect on pressure head was taken into account in our simulations (see Section 2.1.2for details).
The governing equation for heat flow is (de Vries, 1958;Saito et al., 2006) ∂ ∂ where  h ,  w , and  v [J m −3 K −1 ] are volumetric heat capacities of the bulk soil, liquid water, and water vapor, ρ w [kg m −3 ] is the density of liquid water, λ e [J kg −1 ] is latent heat of vaporization,λ h [W m −1 K −1 ] is total soil thermal conductivity, and  l [m s −1 ] and  v [m s −1 ] are the flux densities of liquid water and water vapor, respectively.

Transport coefficients
The soil water retention curve θ(ℎ) and the hydraulic conductivity curve  lh (ℎ) were parameterized with the Peters-Durner-Iden model (Iden & Durner, 2014;Peters, 2013Peters, , 2014)).This model ensures a zero-water content at a pressure head of −10 6.8 cm ("oven dryness") and includes capillary and non-capillary (film and corner) flow components in the hydraulic conductivity function.The water retention curve uses a capillary saturation function  c and a non-capillary saturation function  nc and is parametrized as where θ s [cm 3 cm −3 ] is saturated water content, and θ r [cm 3 cm −3 ] is the maximum volume of non-capillary water per unit volume of soil.The hydraulic conductivity function is given by where  lc (ℎ) is capillary hydraulic conductivity, and  lnc (ℎ) is non-capillary hydraulic conductivity representing film and corner flow (Peters, 2013).
A more detailed explanation of the parametrization of the hydraulic properties can be found in the original publications and in Table S1 in the Supplemental Material (SupMat).For the sand, a unimodal capillary saturation function was used, for the loam a bimodal function was required.The hydraulic properties of both soils are shown in Figure 1a,b.The parameters of the Peters-Durner-Iden model were determined by inverse modeling of the experiments conducted for "WIND" experiments as documented in Iden, Diamantopoulos et al. (2021) and are given in Table S2 in the SupMat.θ s had to be adjusted slightly for the other two experiments to match the initial water content that was determined gravimetrically.This resulted in minor scaling of the water retention curve toward the wet end.
The isothermal vapor conductivity function is (Saito et al., 2006) where  v [m s −1 ] is the diffusion coefficient of water vapor in soil (dependent on temperature and liquid water content), ρ vsat [kg m −3 ] is the saturation vapor density in air (dependent on temperature),  w [kg mol −1 ] is the molecular mass of water,  [m s −2 ] is acceleration due to gravity,  [J mol −1 K −1 ] is the universal gas constant,  [K] is absolute temperature, and  r [-] is relative humidity.The equations for the temperature dependences of  v and ρ vsat can be found in the SupMat.Under equilibrium, relative humidity is related to soil water pressure head by the Kelvin equation if the effect of osmotic potential can be neglected: The isothermal vapor conductivity functions of both soils are shown as dashed lines in Figure 1b.The thermal vapor conductivity function is (Saito et al., 2006) and is shown for both soils in Figure 1c.Note that Equation (7) differs from the respective equation in Milly (1982) but the practical difference in  vT is marginal.Furthermore, Equation ( 7) omits the enhancement factor sometimes used in  vT (Cass et al., 1984), because the related pore-scale processes are very difficult to observe (Shokri et al., 2009) and their magnitude and importance are still controversial (Shahraeeni & Or, 2012;Webb & Ho, 1999).
In addition to the temperature dependence of the two vapor flow conductivities, the SHPs are also temperature dependent.The temperature dependence of soil water pressure head was treated as where σ [kg s −2 ] is the surface tension, and  ref [K] is the reference temperature.Equation ( 8) accounts only for the temperature dependence of the surface tension of pure water (see SupMat for details).In fact, the temperature effect on pressure head in real soils might be higher (Bachmann et al., 2002;Joshi et al., 2019), but simulating this requires a soil-specific parametrization.
The temperature dependence of hydraulic conductivity was modeled as where μ [kg m −1 s −1 ] is the dynamic viscosity of water.The equations for the temperature dependence of σ and μ are given in the SupMat.Soil thermal conductivity as function of liquid water content was parameterized with the Campbell model (Campbell, 1985): with five parameters a 1 , ⋯, a 5 that can be related to basic soil properties.Point values of thermal conductivity were measured for both soils with the method of Markert et al. (2016) at TU Berlin and the five parameters were estimated by fitting Equation ( 10) to the data by nonlinear least squares regression.Parameter values are provided in the SupMat, and the function λ h (θ l ) is shown in Figure 1d for both soils.

Surface energy balance
As vaporization includes the consumption of significant amounts of energy, it is crucial to consider the surface energy balance when simulating the evaporation process.Solving the energy balance assures system-dependent calculations of the top boundary condition, in particular the consideration of the temperature dependence of the saturation water vapor pressure.The surface energy balance is where  n [W m −2 ] is the net radiation (positive implies directed toward the soil),  [W m −2 ] is the total soil heat flux (positive upward),  [W m −2 ] is the sensible heat flux between the soil surface and the atmosphere (positive implies directed away from the soil), λ e ρ w  a [W m −2 ] is the latent heat flux from the soil to the atmosphere (positive means evaporation), and  a is the actual evaporation rate [m s −1 ].The net radiation is the sum of the net short-wave radiation  ns and the net long-wave radiation  nl Net short-wave radiation is calculated as where  [W m −2 ] is incoming short-wave radiation, and α [-] is the albedo of the soil surface.The albedo depends on the surface water content and most of the parameterizations use a piecewise linear function of the surface water content θ surf (Bavel & Hillel, 1976;Novak, 1981Novak, , 2010)).The equation and parameters used in our simulations are provided in the SupMat.
The net long-wave radiation at the soil surface is calculated by where σ sb [W m −2 K −4 ] is the Stefan-Boltzmann constant,  s [K] and  a [K] are the temperatures of the soil surface and laboratory, respectively, and ε a [-] and ε s [-] are the respective emissivities.We assumed an emissivity of one for the surroundings and a dependence of emissivity on volumetric soil water content at the soil surface, θ surf [cm 3 cm −3 ], (Campbell, 1985;Monteith & Unsworth, 2014;Saito et al., 2006) according to The sensible heat flux  is driven by the temperature difference between the soil surface  s and the atmosphere  a and calculated as where  a [J m −3 K −1 ] is the volumetric heat capacity of dry air at constant pressure and  a [s m −1 ] is the aerodynamic resistance which reflects the turbulent airflow in the vicinity of the soil surface.

Evaporation rate
The actual evaporation rate was calculated from the differences of the water vapor pressures between soil surface and atmosphere (Li et al., 2019;Novak, 2016;Sakai et al., 2011) and an atmospheric resistance to vapor flow: where ρ vs [kg m −3 ] is the actual vapor content at the soil surface, ρ va [kg m −3 ] is the water vapor content of the atmosphere at reference height, and  a [s m −1 ] is the aerodynamic resistance for water vapor.Actual vapor content at the soil surface ρ vs [kg m −3 ] was calculated as the product of the saturation vapor pressure, which is a function of temperature, and relative humidity, which is related to pressure head by Equation ( 6).Some studies calculate the actual evaporation rate with an additional soil or surface resistance  s [s m −1 ] as follows: However,  s should be set to zero in a coupled water, vapor, and heat flow model applied to bare-soil evaporation (Novak, 2019;Vanderborght et al., 2017).Given that  s is frequently used in combination with PdV, we analyzed the impact of two popular parametrizations for the dependence of surface resistance on water content on the simulation results to illustrate the effect of double counting.The first one is based on Camillo and Gurney (1986) and relates  s linearly to surface water content: The second parameterization is from van de Griend and Owe (1994) and was derived from experiments on a fine sandy loam: where α go = 35.63 is a fitting parameter and θ min = 0.15 is an empirical threshold water content above which the soil is able to deliver vapor at a potential rate, and  sl = 10 s m −1 is Vadose Zone Journal the resistance to molecular diffusion across the water surface itself.These empirical resistance models have been developed in different contexts than PdV theory, that is, for use within the Penman-Monteith combination equation and without detailed soil moisture flow modeling.As a consequence, their application within PdV is conceptually wrong.In contrast, use of the boundary-layer resistance  bl derived by Haghighi et al. (2013) and Haghighi and Or (2015) is consistent with the PdV theory because it describes the effect of point sources of moisture at the surface on the vapor exchange between soil and atmosphere through the boundary layer.This boundary-layer resistance was calculated as where δ [m] is the thickness of the diffusive layer,  p [m] is the median radius of the pores, and  a [m 2 s −1 ] is the temperature-dependent diffusion coefficient in free air.Note that the point source concept was also applied to gas diffusion through stomata by Lehmann and Or (2015).Actual vapor content at the soil surface ρ vs [kg m −3 ] was calculated as product of the saturation vapor pressure, which is a function of temperature, and relative humidity.

Evaporation experiments
We conducted laboratory evaporation experiments with a silt loam from Selhausen (LOAM) and a sand from Schunteraue close to Braunschweig (SAND) in a temperature-controlled lab under different boundary conditions.The packed soil columns had a diameter of 23.2 cm and a height of 10 cm resulting in a volume of the bulk soil of approx.4230 cm 3 .They were saturated with water by upward flow from the bottom controlled by Mariotte bottles and placed in boxes with a Styrofoam lining for thermal insulation.Details of the basic experimental setup and the two soils can be found in Iden, Diamantopoulos et al. (2021).Three types of atmospheric boundary conditions were used in a sequence of experiments: (i) wind forced by a fan in a dark laboratory without short-wave radiation input ("WIND"), (ii) wind plus shortwave radiation through a lamp ( = 200 W m −2 ) ("RAD"), and (iii) wind combined with intermittent short-wave radiation ("RAD-I").To create intermittent short-wave radiation in RAD-I, the lamp was turned on for 18 h and then subsequently switched off for 6 h.This on-off sequence was repeated throughout the duration of the experiment.Vertical airflow directed toward the soil was generated by a CPU fan creating a wind speed of about 0.8 m s −1 .We decided to apply a vertical ventilation to prevent lateral variation The entire setup, including cables and data loggers, was mounted on a balance to determine the effective water content and the evaporation rate.The depth of the sensors can be found in Table S5.
of fluxes and states.The experiments with RAD were performed first, followed by WIND, and the RAD-I experiments were run last.The soil columns were re-saturated between experiments.
Air temperature and relative humidity were measured in the vicinity of the experiments (approx. 1 m distance).Sensors were installed in the columns to measure temperature (PT 100, Testo AG, accuracy 0.2 K, 3-mm diameter), pressure head (mini-tensiometers, UMS T5x, now Meter AG, accuracy 0.2 cm, maximum error 1%, 5 mm diameter) and relative humidity of soil air (Rotronic hygroclip HC2-C04, Rotronic Messgeräte GmbH, accuracy 1.5% rH, 4 mm diameter) at different depths (Table S7).A sketch of the setup is shown in Figure 2.After sensor installation, Styrofoam chips were carefully placed around the soil column and sensors to provide further thermal insulation.The entire setup, including data loggers and all cables, was placed on a scale to continuously monitor the water loss by evaporation.Each of the three consecutive experiments was run for ca.15 days.Pressure heads were calculated from relative humidity data with Equation (6) for ℎ < −10 5 cm to complement the time series of pressure heads measured by mini-tensiometers toward the dry (hygroscopic) range.These pressure head data have a relative error smaller than 10% (Iden, Diamantopoulos et al., 2021).

Calculation of the resistance to vapor flow during stage-one evaporation
A resistance to vapor flow was calculated from the measured evaporation rates and Equation ( 17) as a diagnostic variable to quantify the observed resistance.Note that this resistance was not used in the numerical modeling.As the relative humidity in the soil is above 99% as long as the surface pressure head is greater than −10 4 cm (Equation 6), it is justified to assume that the actual vapor content at the soil surface equals the saturation vapor pressure, which is only a function of soil surface temperature  0 [K].The resistance to vapor flow was calculated after rearranging Equation ( 17) and setting ρ vs = ρ vsat (  ) [kg m −3 ]: where   [s m −1 ] is the apparent resistance to water vapor flow that replaces   in Equation ( 17), and   [m s −1 ] is the measured actual evaporation rate that is calculated from the change in mass of the soil column.The resistance calculated from Equation ( 22) can be interpreted as an "atmospheric resistance" during stage-one (100% relative humidity at soil surface) but includes transport processes within the soil during stage-two.We used the soil temperature measured at depth 0.8 cm as a proxy for the soil surface temperature.As temperature gradients in the soil during stage-one were relatively small, the error in the calculated resistance to vapor flow resulting from this simplifying assumption is small.The vapor content in the air was calculated from the measured ambient temperature and relative humidity data.

Numerical simulation setup
The Hydrus-1D software code (Simunek et al., 2008(Simunek et al., , 2013(Simunek et al., , 2016) was used for the numerical simulations.The changes made to the code are summarized in the SupMat.The initial condition for heat flow was a constant temperature set equal to room temperature in the entire soil.The surface energy balance was solved at the top boundary by the Hydrus code.We extrapolated the temperature data using a second-order polynomial at each measurement time to the bottom of the column at 10 cm and assigned the resulting temperature time series as a transient Dirichlet condition at the bottom.
For water flow, the initial condition was hydrostatic.A noflux boundary condition was used at the bottom.At the top, an atmospheric boundary condition with a critical pressure head of −10 7 cm was applied.This value was small enough to prevent the surface boundary condition to be switched to a Dirichlet boundary condition (Li et al., 2019).The actual evaporation rate was calculated internally in Hydrus-1D with Equation ( 17) or ( 18).The required input data air temperature and relative humidity of the laboratory air are shown in Figure 3.Note that Figure 3b shows the saturation deficit (1 −  r ) because it is more directly related to the evaporation rate than the humidity  r .The atmospheric variables were measured in close proximity of the experiments and therefore differ between experiment types.
We performed three sets of simulations.In the first set, the aerodynamic resistance  a was kept constant throughout the simulation of the experiments.Its value had to be adjusted such that the observed evaporation rates in the beginning of stage-one were well matched by the simulations.The reasoning behind this is that the soil surface is fully saturated with water vapor during stage-one and the aerodynamic resistance is the only adjustable parameter controlling the evaporation rate in Equation ( 17).The aerodynamic resistance was not further adjusted during stage-two, reflecting the constant wind speed throughout the experiments.The values are given in Table 1, which are identical for both soils and differ only slightly between boundary conditions.
In the second set of simulations, the soil surface resistance parametrizations of Camillo and Gurney (1986) and van de Griend and Owe (1994) were added.Finally, the boundarylayer resistance  bl of Haghighi et al. (2013) and Haghighi and Or (2015) was used.The median radius of the pores was calculated from the water retention curves and was 65 and 4.8 μm for sand and silt loam, respectively.The thickness of the diffusive layer was set to 1.6 mm for all six experiments as that led to best match with the measures evaporation rates during stage-one.The applied theory of the boundary-layer resistance was developed for lateral flow whereas air was blown vertically toward the soil surface in the experiments.Therefore, the application of the model of Haghighi et al. (2013) to our experiments is not per se warranted.However, the theory is able to describe a decreasing evaporation rate during stage-one, an important feature of our experimental data as shown later, and the airflow becomes horizontal when it hits the soil surface.

Evaporation experiments under different boundary conditions
In the following, we present evaporation rates, soil temperature, and soil water pressure heads observed in the experiments.The ambient conditions near the soil columns, air temperature, and water vapor saturation deficit are shown in Figure 3.For the experiment WIND, the air temperature (Figure 3a) was constant at 20˚C.The air temperature was also constant for RAD but due to the heating by the lamp the air was almost 3˚C warmer.For RAD-I, the temperature reflects the transient input of short-wave radiation and varies systematically between 20 and 21˚C.Figure 3b shows that the saturation deficit fluctuated, which influenced the evaporation rates.
Figure 4 shows the measurement results for the sand column.The top row shows the evaporation rates, the second row shows cumulative evaporation, the third row the soil temperatures, and the bottom row the pressure head (shown as suction, log scale), as a function of time.For the experiment without applied short-wave radiation (WIND, left), three evaporation stages can be distinguished (Iden, Diamantopoulos, et al., 2021) based on the evaporation rate and soil temperature.The first phase (stage-0) shows a falling evaporation rate that is caused by the initial cooling of the soil due to the loss of latent heat at the surface.The initial drop is sudden and pronounced and leads to a quick temperature decrease of about 7˚C.The second phase is the commonly described stage-one during which the evaporation rate and soil temperature are relatively constant.This phase lasts for less than 2.5 days.The evaporation rate fluctuates during stage-one but does not show a decreasing trend.A comparison with the saturation deficit shown in Figure 3b highlights the close relationship between the driving force and the evaporation rate.Stageone is followed by the falling-rate phase (stage-two) with a continuous decrease of the evaporation rate during which the soil temperature increases due to the decreasing loss of latent energy.Finally, the evaporation rate approaches zero and the temperature approaches the ambient air temperature.
In the experiment using constant radiation (RAD, middle column), the soil temperature decreases until the end of stageone after approximately 2.2 days.The reason is the increasing evaporation rate, which reflects the rising saturation deficit in the air (Figure 3b).The evaporation rate during stage-one is higher than in the experiments with WIND because of the higher soil temperature which increases the vapor content at the soil surface.At the onset of stage-two, the evaporation rate falls, similar to the WIND experiment, the soil temperatures start to increase and approach a constant value that reflects a new equilibrium defined by the surface energy balance (with a negligible latent heat flux).
For the experiments with intermittent radiation (RAD-I, right), the evaporation rate drops rapidly in the beginning, which is similar to the experiment using WIND.Overall, the dynamics is similar to RAD but it shows the superimposed effect of the intermittent short-wave radiation.The direct influence of the short-wave radiation on the evaporation rate and soil temperature becomes evident.As soon as the lamp was switched off, the soil temperature dropped markedly, which led to a significant decrease in the evaporation rate.The maximum temperature of RAD-I is a few degrees colder than in RAD because the energy input was not constant and the surrounding temperature was cooler (Figure 3a).In the case of RAD-I, the saturation deficit of the air was almost constant during stage-one (Figure 3b), which resulted in the smoother course of evaporation rate and soil temperatures.
The cumulative evaporation data shown in the second row show a similar temporal pattern for the different boundary conditions.The differences in the total evaporation are caused by different energy inputs and different initial saturations caused by the rewetting of the soil columns between experiments.The time series of pressure head for the sand are very similar for the three boundary conditions and are typical for evaporation experiments on sandy soils (Iden, Diamantopoulos, et al., 2021).Notably, the measured pressure head data reached values as low as −10 −6 cm due to the use of the humidity sensors.In the beginning, the vertical gradients in the soil are small because the evaporation rate is much smaller than the hydraulic conductivity.When the hydraulic conductivity decreases, it becomes much smaller than the evaporation rate after a certain time, and accordingly the pressure head gradient builds up.The similarity in the time series shows that the SHP remained largely unchanged during the three successive experiments.
Figure 5 shows the measurement results for the loam column.Overall, there is a close similarity to the data in Figure 4 and thus there is no need to repeat the detailed discussion.A striking difference is that the evaporation rate during stage-two decreases less for the loam compared to the sand because the former has a higher hydraulic conductivity in medium to dry soil.This is also the reason for the smaller vertical pressure head gradients in Figure 5 (bottom) compared to Figure 4.

Simulations with the coupled water, vapor, and heat flow model
Simulating the evaporation experiments under different boundary conditions with the coupled water, vapor, and heat flow model and a constant aerodynamic resistance overall leads to an excellent description of the observed dynamics of the surface water flux and the soil temperatures (Figures 4 and 5).It is important to recall that the experimental data of the WIND experiments were used to identify the SHP by inverse modeling (Iden, Blöcher, et al., 2021).Therefore, these experiments should not be regarded as part of the validation-only the experiments RAD and RAD-I are validation experiments.In all cases, the time of the onset of stage-two is matched remarkably well and this indicates that the transition from liquid-dominated to vapor-dominated flow in the top of the soil is adequately described by the model (Figures 4 and 6 in Iden, Diamantopoulos, et al., 2021 give a more detailed picture of this).Notably, this holds true for all boundary conditions and both soil textures.In addition, the evaporation rates during stage-two are represented remarkably well across all scenarios and soil textures.For RAD-I, the ups and downs of evaporation rate and soil temperature that are induced by the intermittent application of short-wave radiation are nicely represented.This shows that the model, which calculates the system-dependent evaporation rate internally with Equation ( 17), can reproduce the water flux across the top surface if the radiative input is altered.The slight mismatch between measured and simulated temperatures indicates that the energy balance cannot be properly described by the model.We hypothesize that the soil exchanged some energy with the insulating box and that this lateral exchange cannot be captured by a flow model in one spatial dimension.
In addition to the surface water flux and the soil temperatures, the soil water pressure heads shown in the bottom rows of Figures 4 and 5 are also described quite well.This shows that not only the water balance of the entire system is correctly described by the model, but also the local energy status of soil water.There are, of course, still some systematic discrepancies between measured and simulated values, but the temporal dynamics and the development of the pressure head gradient in both soils are described well for all boundary conditions.At this point, it is important to recall that a parametrization of the SHP accounting for non-capillary storage and conductivity was used.Because film and corner flow are accounted for in the hydraulic conductivity function, even the pressure head data in the hygroscopic range (ℎ < −10 5 cm), monitored with the humidity sensors, are reasonably matched (Iden, Diamantopoulos, et al., 2021).For RAD-I, the interruptions in the short-wave radiation become evident both in the observed and simulated pressure head data except for the initial phase.Switching off the lamp leads to a cooling of the soil, a decrease in the evaporation rate and a vertical redistribution of moisture because the gradient in pressure head adapts to the changed water flux density.In addition to this, the cooling and reheating induce a temperature effect on pressure head, which is accounted for in the model by the temperature dependence of soil surface tension.
Despite the overall excellent description of the data by the model, a systematic overprediction of the evaporation rate at the end of stage-one occurs for both soils and all boundary conditions.This overestimation resulted in a stronger surface cooling in all simulations.The increase of the simulated evaporation rate for WIND and RAD was caused by a decrease in ambient relative humidity, which increased the vapor deficit (Figure 3b).The model captures this change in the driving force correctly but the resulting increased evaporation is clearly not in agreement with the observations.For RAD-I, the laboratory air remained at relatively constant humidity during stage-one (Figure 3b) and the simulated evaporation rates thus do not show an upward trend.In contrast, the observed evaporation rates show a decreasing trend during stage-one.Summing up, the decrease in the evaporation rate during stage-one which was observed for both soils and all boundary conditions cannot be explained by the applied model which uses a constant resistance to vapor exchange,   , for calculating the evaporation rate.This issue deserves a closer look and will be analyzed in more detail in the next section.

Analysis of evaporation rate and resistance to vapor flow during stage-one
The observed decrease of the evaporation rate during stageone cannot be explained by a change in the ambient conditions.Indeed, the fluctuations in the saturation deficit in the air should cause an increase in the evaporation rate for RAD and WIND (only final part of stage-one) and an almost constant evaporation rate for RAD-I (Figure 3b).This is exactly what is predicted by the water, vapor, and heat flow model (Figures 4 and 5,top row).According to Equation (17), a lower observed evaporation rate is either caused by a smaller vapor pressure difference between soil and air, or a higher resistance to vapor flow.The first option is not supported by the saturation deficit in the air (Figure 3b).A lower gradient could also be caused by an increase in soil temperature during stage-one but this was not observed, except for WIND.Of course, soil temperature and evaporation rate are mutually coupled through the surface energy balance and the temperature dependence of saturation vapor pressure.Therefore, we are confronted with a "chicken or the egg" type of problem.An increase in the aerodynamic resistance used in Equation ( 17) can be excluded because the air was blown from the top toward the soil with a constant velocity and the soil temperatures were lower than air temperature throughout stage-one leading to stable atmospheric conditions for all six experiments.As a consequence, an additional resistance to vapor flow must exist.
Figure 6 shows the apparent resistance to vapor flow   calculated by Equation ( 22) for all six experiments.The horizontal axes in Figure 6a are the elapsed time for each experiment, normalized by the duration of stage-one.Interestingly, an increase of the resistance with time becomes evident during stage-one and the increase is of similar magnitude and structure across all soils and boundary conditions.Note that the assumption behind Equation ( 23) that the soil surface is saturated with water vapor is no longer valid during stagetwo.The steep increase of the calculated apparent resistance at the onset of stage-two (rel.time greater than 1) is simply a consequence of this assumption.As the soil surface dries out, the relative humidity becomes smaller than one and the total water flux (liquid and vapor) toward the soil surface becomes the limiting factor for evaporation.The increase of the resistance during stage-one is not in agreement with Equation (17) because the latter uses a constant   as the only resistance to vapor flow.Therefore, the coupled model must overestimate the evaporation rate as time progresses during stage-one.Interestingly, although the resistance used to model the evaporation rate is smaller than the calculated effective resistance, we only observe an overestimation of the evaporation rate at the end of stage-one.This is because the underestimated temperature mitigates the effect of a lower resistance.Overestimations of evaporation rates by coupled models have also been reported in Bittelli et al. (2008) but the reported over-estimation reached more than 100% and was not related to a particular evaporation stage.
Figure 6b shows the total resistance as function of soil water content.Because the water content at the surface was not measured, and to avoid the use of modeled data, we selected the water content which was calculated from the measured pressure head at  = −0.8cm and the soil water retention curve.The apparent resistance to vapor flow clusters nicely for the two soils.This implies that for a given soil, the relationship may be invariant with respect to the boundary conditions.Interestingly, the total resistance for the loam increases during stage-one although the water content is almost constant, suggesting that a water content-dependent parameterization of the resistance may be insufficient.
In the literature, a soil surface resistance has been proposed as an additional factor to be included in evaporation models for bare soils to decrease the evaporation rate during stage-two.The results of using the soil resistance formulations of Camillo and Gurney (1986) and van de Griend and Owe (1994) are shown in Figure 7.It becomes evident that the use of these models causes the transition to stage-two to occur too early and leads to a mismatch with the observed evaporation rates during stage-two, in particular for the sand.This is in accordance with theoretical arguments that using this type of resistance in a coupled water, vapor, and heat model is double counting (Novak, 2016;Vanderborght et al., 2017).The increased evaporation rate of the models with added soil resistance during stage-two in some cases of Figure 7 is caused by the decrease of the evaporation rate during earlier times.If the water loss is decreased at early times, the soil water content becomes larger and this implies a higher soil water pressure head and soil hydraulic conductivity.This causes a more intense upward movement of liquid water and therefore a higher evaporation rate.
The results of the simulations with the boundary-layer resistance of Haghighi et al. (2013) are shown in Figure 7, too.Note that this resistance is able to predict a decreasing evaporation rate during stage-one and improves the match between model and observations for the sandy soil (Figure 7).However, for the loamy soil, the effect of the resistance is minimal and does not improve the prediction of the evaporation rate.The reason for the difference between soils is that the resistance model predicts a stronger increase with decreasing water content for soils with a greater mean pore size.This assumption is contradictory to the calculated effective resistance shown in Figure 6b.Although the boundary-layer resistance leads to promising results, its dependence on mean pore size requires a revision, at least for finer soil textures.
Other potential sources for the observed inability of the coupled model to describe the evaporation rates are nonequilibria that are not accounted for in the model.Phase equilibrium between the pressure head ℎ and the relative humidity is assumed by the Kelvin equation (Equation 6).This equilibrium assumption has recently been challenged (Li et al., 2019) but was supported by Novak ( 2019) by a theoretical analysis, at least for soils that are not too coarse and evaporation rates which are not too high-both are warranted for the soils investigated in this study.Another potential process that may take place as a consequence of quick changes in the upper boundary condition is dynamic nonequilibrium of soil water flow, that is, a nonequilibrium between soil water content and pressure head (Diamantopoulos & Durner, 2012;Diamantopoulos et al., 2012) that has been observed and modeled under a variety of boundary conditions (Diamantopoulos et al., 2015).It is beyond the scope of this study to analyze this issue in more depth.Finally, the coupled model assumes thermal equilibrium between the phases of the porous media, that is, the temperatures of pore water, soil air, and solid matrix are identical.As the phase transition occurs in the water phase, the water is likely cooler than the surrounding soil particles that would decrease the saturation vapor pressure and hence the evaporation rate.Although local thermal nonequilibrium has been shown to be relevant for infiltration (Heinze & Blöcher, 2019), we can only speculate about its magnitude and relevance for evaporation.

CONCLUSIONS
The validation of coupled water, vapor, and heat flow models needs carefully designed experiments under controlled and differing boundary conditions.We have presented the results of laboratory evaporation experiments for two different soils under three different boundary conditions.The experiments were numerically simulated with a coupled water, vapor, and heat flow.One may criticize the use of one of the experiments for the identification of SHP, but an accurate parametrization of these is crucial for validating flow models.The coupled flow model is very sensitive to SHP (Fetzer et al., 2017;Vanderborght et al., 2017) and any error in their parametrization may lead to the wrong conclusion that the flow model is incorrect, although the dominant error lies in the description of the transport properties.
The overall temporal patterns of evaporation rate, soil temperature, and soil water pressure head observed in the experiments were represented very well with a constant aerodynamic resistance   .A prerequisite for this is that an adequate parametrization of SHP is used (Iden, Blöcher, et al., 2021;Iden, Diamantopoulos, et al., 2021).The timing of the transition from stage-one to stage-two evaporation was matched well.However, the observed evaporation rate decreased in the course of stage-one and this could not be represented by the simulation.As a consequence, a too pronounced cooling of the soil was simulated.These discrepancies were observed for both soils and for all boundary conditions.
Apparent resistances to vapor flow calculated from the measured actual evaporation rates and the vapor pressure difference between soil and air increase during stageone.This is not in agreement with most of the currently applied models and contradicts the concept of a constant evaporation rate during stage-one.Simulations with empirical surface resistance parametrizations led to large errors in evaporation rates and should be discarded in coupled watervapor-heat models because of their intrinsic double counting.As the calculated apparent resistance hinges on the assumption of a full saturation of the complete soil surface with water vapor, a deviation from this assumption is a possible explanation for the observed mismatch of the evaporation rates during stage-one.A possible cause for this is an uneven spatial distribution of water at the soil surface caused by the simultaneous occurrence of water-filled and already emptied pores or wet and dry patches.We showed that the use of the boundary-layer resistance derived by Haghighi et al. (2013) improved the predictions of the evaporation rate for the sandy soil but fell short for the loamy soil.This finding indicates that the boundary-layer resistance is a promising approach but requires further study and refinement, in particular for soils with wider pore-size distributions.Other causes for the observed decrease during stage-one might be thermal nonequilibria and dynamic nonequilibrium in soil water flow.The setup of our experiments does not enable to identify the exact cause for the observed behavior and to discriminate between possible causes.

AU T H O R C O N T R I B U T I O N S
Johanna R. Blöcher: Conceptualization; data curation; formal analysis; investigation; methodology; software; visualization; writing-original draft; writing-review and editing.Efstathios Diamantopoulos: Data curation; investigation; methodology; software; writing-review and editing.Wolfgang Durner: Funding acquisition; writing-review and editing.Sascha C. Iden: Conceptualization; data curation; formal analysis; investigation; methodology; project administration; software; supervision; visualization; writing-original draft; writing-review and editing.

A C K N O W L E D G M E N T S
This study was financially supported within the DFG Research Group FOR 1083 "Multi-Scale Interfaces in Unsaturated Soil" (MUSIS), grant DU283/10-1.We thank Jirka Simunek for sharing the source code of Hydrus-1D and Michael Facklam (TU Berlin) for measuring soil thermal conductivity.Michael D. Novak gave us the hint on the inappropriate use of the thermal liquid flux in potential-based formulations of the PdV theory.

C O N F L I C T O F I N T E R E S T S T A T E M E N T
The authors declare no conflicts of interest.

D A T A AVA I L A B I L I T Y S T A T E M E N T
The measurement data presented in this article is accessible from LeoPARD, the publication and research data reposi-

F
Soil water retention curves (a), soil hydraulic conductivity curves (b), thermal vapor flow conductivity (c), and thermal conductivity curves (d) for the two investigated soils.

F
Setup of the evaporation experiments.The pressure heads ℎ  () are measured by tensiometers (right).The two sensors on the top left measure the relative humidity   () in the soil and the black sensors measure soil temperature   ().All sensors were connected to data loggers.The soil column was placed in a box with Styrofoam lining and Styrofoam chips for thermal insulation.The top of the soil column was open to exchange water and energy with the surrounding.

F
Temperature (a) and water vapor saturation deficit of the air (b) in the laboratory, measured in the vicinity of the soil columns.Blue indicates "fan only" (WIND), red indicates fan and lamp (RAD), and orange shows fan and intermittent radiation (RAD-I).

F
I G U R E 4 Time series of actual evaporation rates (top row), cumulative evaporation (second row) temperatures (third row), and pressure heads (shown as common log of suction, bottom row) for the sample SAND.The dots indicate observed data, lines are time series simulated with the coupled water, vapor, and heat flow model.The different colors refer to the boundary conditions; blue is fan only (WIND, left), purple is fan plus short-wave radiation (RAD, middle), and orange is fan and intermittent short-wave radiation (RAD-I, right).
The same as Figure4but for sample LOAM.

F
Effective resistance to vapor flow as function of relative time in stage-1 (a) and water content 0.8 cm below the surface (b).The different colors refer to the boundary conditions (WIND, RAD, RAD-I).

F
I G U R E 7 Evaporation rates with different resistance parametrizations.The black line is the default for the coupled model (r s = 0), dark blue lines ("GO94") indicate the soil resistance model of van de Griend and Owe (1994), light blue lines (CG86) show the soil resistance parameterization of Camillo and Gurney (1986) (CG), and orange lines (H13) indicate the boundary-layer resistance parameterization of Haghighi et al. (2013).Observations are shown as grey circles.Soil samples from top to bottom, experiments with different forcing from left to right.
lh [m s −1 ] is the liquid hydraulic conductivity,  vh [m s −1 ] is the effective hydraulic conductivity for the isothermal diffusion of water vapor,  vT [m 2 K −1 s −1 ] is the thermal vapor flow conductivity,  [m] is the vertical coordinate, positive upward, and  [s] is time.Equation (1) is taken from Saito Aerodynamic resistances used for the simulations.
T A B L E 1