Evaluation of Model Concepts to Describe Water Transport in Shallow Subsurface Soil and Across the Soil–Air Interface

Soil water evaporation plays a critical role in mass and energy exchanges across the land–atmosphere interface. Although much is known about this process, there is no agreement on the best modeling approaches to determine soil water evaporation due to the complexity of the numerical modeling scenarios and lack of experimental data available to validate such models. Existing studies show numerical and experimental discrepancies in the evaporation behavior and soil water distribution in soils at various scales, driving us to revisit the key process representation in subsurface soil. Therefore, the goal of this work is to test different mathematical formulations used to estimate evaporation from bare soils to critically evaluate the model formulations, assumptions and surface boundary conditions. This comparison required the development of three numerical models at the REV scale that vary in their complexity in characterizing water flow and evaporation, using the same modeling platform. The performance of the models was evaluated by comparing with experimental data generated from a soil tank/boundary layer wind tunnel experimental apparatus equipped with a sensor network to continuously monitor water–temperature–humidity variables. A series of experiments were performed in which the soil tank was packed with different soil types. Results demonstrate that the approaches vary in their ability to capture different stages of evaporation and no one approach can be deemed most appropriate for every scenario. When a proper top boundary condition and space discretization are defined, the Richards equation-based models (Richards model and Richards vapor model) can generally capture the evaporation behaviors across the entire range of soil saturations, comparing well with the experimental data. The simulation results of the non-equilibrium two-component two-phase model which considers vapor transport as an independent process generally agree well with the observations in terms of evaporation behavior and soil water dynamics. Certain differences in simulation results can be observed between equilibrium and non-equilibrium approaches. Comparisons of the models and the boundary layer formulations highlight the need to revisit key assumptions that influence evaporation behavior, highlighting the need to further understand water and vapor transport processes in soil to improve model accuracy.


Introduction
Soil water evaporation, associated with water movement and heat transfer, plays an important role in the water cycle and energy balance across land-atmosphere interface. Soil water evaporation is a strongly coupled phenomena that involves multiple physical processes (e.g., liquid water flow, phase change, vapor transport, heat transport), and is affected by atmospheric conditions (e.g., wind speed, ambient temperature and relative humidity, radiation), surface conditions (e.g., roughness, surface topography) and soil properties (e.g., hydraulic properties, transport properties and thermal properties). In numerical modeling efforts, soil water evaporation oftentimes serves as the top boundary condition for water flow and is related to soil heat flux at the surface which serves as the top boundary condition of heat flow, thus interacting with the soil moisture and soil temperature and further influencing the carbon (Davidson and Janssens 2006;Koven et al. 2013) and nitrogen cycles (Parton et al. 2001). Therefore, the proper prediction of soil moisture and temperature distributions is necessary for predicting evaporation and ultimately climatic conditions at a variety of spatial and temporal scales. However, the dynamic interactions of heat and mass transfer are oftentimes not considered in numerical models due to the complexity of field scenarios and the lack of field and laboratory data capable of testing such models. Models vary in their level of complexity and oftentimes rely on the use of fitting parameters or key simplifying assumptions without proper understanding of their implications. These models are then tested with limited field data and rarely compared with one another. An alternative approach is to use controlled laboratory-scale experiments to generate data under transient yet controlled conditions. These data can then be used to systematically test key model assumptions and parameters. Use of the same modeling platform with different model assumptions tested against the experimental data then allows for the investigation of the dominant mechanisms and can provide guidance for the improvement of simplified parameterizations and boundary condition development.
Evaporation from initially saturated porous media occurs in two main stages, the liquid water flow dominated stage I and vapor diffusion controlled stage II, with a transition stage in-between (e.g., Lemon 1956;Philip 1957;Idso et al. 1974;Mosthaf et al. 2014). Stage I evaporation is characterized by a high and nearly constant evaporation rate (potential evaporation rate), mainly controlled by atmospheric demand. Stage I evaporation is mainly driven by capillary liquid flow to the soil surface and is maintained until the hydraulic connection between the receding drying front and the soil surface starts to disconnect. At this time, the system enters into a transition stage (Ben Neriah et al. 2014;Mosthaf et al. 2014). The drying front depth at which the soil system transitions from stage I to stage II evaporation is defined by the pore size distribution of the porous media and oftentimes referred to as the characteristic length (Yiotis et al. 2006;Lehmann et al. 2008). When the surface is completely dry and the hydraulic connection between the drying front and the soil surface is completely severed, evaporation enters into the diffusion dominated stage II. During this stage, a dry surface layer (vapor diffusion), film region (capillary liquid flow) and saturated region all exist within the porous media. Phase change between liquid water and water vapor occurs within the film region, and vapor diffuses through the dry layer to the soil surface.
Because of the wide range of scales and purposes of numerical modeling of heat, mass and momentum transport and the exchanges across the land/atmosphere boundary, a variety of different REV scale numerical modeling theories have been developed, varying in complexity and hence assumptions. For a full review of all such modeling concepts, the reader is referred to Vanderborght et al. (2017). Using a fully coupled land-atmosphere model that implicitly couples the soil and the atmosphere with minimal assumptions at the interface is ideal (e.g., Nield 2009;Shavit 2009;Chidyagwai and Rivière 2011;Mosthaf et al. 2011;Davarzani et al. 2014). Davarzani et al. (2014) presented a model based on the coupling between Navier-Stokes free flow with Darcy porous media flow, demonstrating that the concept can predict the different stages of evaporation with great accuracy. However, these models are not often used at the field scale because of high computational costs. Instead, alternative modeling approaches are used that define top boundary conditions (e.g., water and heat fluxes at the soil surface) based on the free flow conditions. It is important to understand how these top boundary conditions affect the estimates of soil water evaporation.
Numerical models based on the widely used Richards equation have been used to describe water transport within unsaturated soil system, and water exchanges across the land-atmosphere interface from the macroscopic (Darcian) scale to the watershed scale (e.g., Nieber and Walter 1981;Vereecken et al. 1991;Schoups et al. 2005;Mortensen et al. 2006). The most widely used form of Richards equation considers only isothermal liquid water flow, assuming that liquid water moves upward through the unsaturated soil region, driven by capillary pressure and vaporizes at the soil surface. This form neglects vapor diffusion and air flow in the soil as well as the influence of temperature gradients. In some Richards equation-based models, vapor transport is incorporated as part of the liquid water transport equations and heat transfer is coupled with water transport (e.g., Philip and De Vries 1957;De Vries 1958;Milly 1982Milly , 1988Saito et al. 2006;Bittelli et al. 2008;Novak 2010;Deb et al. 2011;Zeng et al. 2011b). In these models (i.e., considering vapor transport) a common approach is to assume local thermodynamic equilibrium which the water vapor is always in equilibrium with the liquid water (i.e., equilibrium phase change condition or instantaneous vaporization) at the REV scale. In modeling efforts based on equilibrium phase change, like in the case of Richards equation-based models, the vapor and capillary pressures are coupled through the use of Kelvin's equation. However, studies show that under certain conditions (e.g., low soil water content or soils with large pores), equilibrium is not always established near the land surface boundary (Benet and Jouanna 1982;Udell et al. 1982;Saito et al. 2006;Chammari et al. 2008;Halder et al. 2011;Ouedraogo et al. 2013). Using an equilibrium vapor pressure model, Pruess and Wang (1983) showed severe temporal oscillations in the local vapor pressure when the dry soil region forms. However, if we consider vapor density as an independent variable instead of depending on capillary pressure, this behavior can be avoided, resulting in a smooth transition from partial to zero saturation (Bixler 1985).
Based on the discussion above, three such modeling concepts were selected for this work, including (1) the traditional one-component (water), one-phase (liquid) Richards model, (2) the Richards model with vapor diffusion incorporated (referred to as the Richards vapor model) and (3) a non-equilibrium two-component (water and air) two-phase (liquid and gas) model (referred to as the Non 2-2 model). Heat transfer is considered in all these three models for consistency. Depending on the model concept, the top boundary condition for water flow varies. First, for the Richards model, as long as the pressure at the soil surface has not reached the critical water pressure head, the potential evaporation rate is usually assigned at the soil surface. When the subsurface soil cannot sustain the flux and the surface pressure head reaches the critical water pressure head, the top boundary condition switches to a constant water pressure condition. Theoretically, the critical water pressure is defined as a value corresponds to small hydraulic conductivity and capacity so the simulation results do not change when decreasing water pressure. However, in general, the definition varies between studies resulting in vastly different values in the literature between −1E5 and −1E7 Pa. For example, a critical water pressure of −1E7 Pa is used for a heterogeneous laboratory soil in Bechtold et al. (2012), −1E5 Pa for coarse sand and −1E6 Pa for sandy loam in Assouline et al. (2014). How these values are selected is oftentimes unclear. The influence of the selection of critical water pressure is numerically analyzed using a Richards equation-based model in Fetzer et al. (2017), demonstrating the sensitivity of the numerical model to the critical water pressure varies with soil hydraulic properties. However, there is no comparison with experimental data in the study. Another commonly applied top boundary condition incorporates soil resistance, which accounts for the evaporation rate reduction when the soil surface dries and the vapor pressure drops below the saturated vapor pressure. This approach is widely used in both small-scale and large-scale models (e.g., Saito et al. 2006;Zeng et al. 2011a;Tang and Riley 2013b;Swenson and Lawrence 2014). Difficulty lies in choosing the appropriate soil resistance formulations and the definition of the formulation's parameters. Finally, for Richards vapor model, one other commonly used top boundary condition is evaporation flux which is expressed as the vapor pressure difference between the surface (calculated based on the local equilibrium assumption) and the atmosphere (Novak 2010(Novak , 2016Sakai et al. 2011). The surface vapor pressure can also be obtained based on simulation of non-equilibrium models, and the corresponding evaporation flux has been used as the top boundary condition (e.g., Smits et al. 2011;Trautz et al. 2015).
The objective of this study is to evaluate the performance of three numerical modeling concepts along with multiple top boundary conditions in simulating water and heat transport in subsurface soil, and describing evaporation from soils. These three models vary in complexity in describing mass transport, from multi-phase, multi-component models to single-phase, single-component models. The performance of the three models is evaluated by comparing simulation results with experimental observations (e.g., evaporation, water content dynamics).
This study is organized as follows. In Sect. 2, we provide the model concepts of three numerical models for water and heat transport in subsurface porous media at the continuum scale. Section 3 presents the experimental apparatus and experimental procedure. And finally, in Sect. 4, we provide the evaluation of the models by comparing the simulation results of the three models with the experimental data.

Numerical Models
In this section, three numerical models varying in complexity and assumptions are introduced. These include (1) the Richards model (non-isothermal one-component one-phase model), (2) the Richards vapor model (non-isothermal, Richards model including water vapor transport) and (3) the non-equilibrium phase change two-component two-phase model (non-isothermal, referred to as the Non 2-2 model). The governing equations, parameters and key variables are presented. Although the Richards equation models are mostly employed in a form that only considers isothermal liquid water flow, it should be noted that the heat transfer is coupled with water flow in all these three models in this study to provide a consistent comparison. In addition, even though lateral water flow has been ignored in some models, it has been indicated that the lateral flow has an effect on the fluxes within the soil profiles and also at the soil-atmosphere interface ) and therefore is considered in all the three models.

Governing Equations
The Richards model assumes that there is only liquid water flow in the unsaturated zone and that the pressure of the gas phase is uniform and constant, equal to atmospheric pressure. The equation for liquid water flow in unsaturated soil is described as: where ρ w is the density of water (kg m −3 ), θ w is the volumetric water content (m 3 m −3 ), θ w = φ S w , where φ is the porosity and S w is the liquid water saturation (−), and t is time (s). θ w is related to water pressure ( p w , Pa) through the constitutive relationship-water content and capillary pressure (θ w − p c ). u w is the velocity of water and can be expressed as follows: where K int is the intrinsic permeability which is only related to the configuration of the porous media (m 2 ), k rw is the relative permeability of water (−), μ w is the dynamic viscosity of water as a function of temperature (kg m −1 s −1 ), g is the gravitation acceleration (g = 9.8 m s −2 ), and z is the elevation (m, positive upward). The relations between the water saturation (S w ) and capillary pressure head (H c ) are described by the van Genuchten constitutive relationship (Van Genuchten 1980): where S ew is the effective water saturation and S rw is the residual water saturation (−), S rw = θ rw /φ, where θ rw is the residual water content (m 3 m −3 ). α (m −1 ) and n (−) are van Genuchten soil parameters related to the inverse of the air-entry pressure and the pore size distribution (m = 1 − 1/n), respectively, which are obtained experimentally. H c is the capillary head (m) and is defined as H c = ( p g − p w )/(ρ w g), where p g is the pressure of the gas phase (Pa). It should be noted that this model does not apply to the conditions when water content is less than the residual water content, which may often happen in the field. There are some studies of extending the curve to over-dry conditions. Interested reader can refer to Fayer and Simmons (1995), Webb (2000), Zhang (2011) for extended models. As mentioned before, the pressure of gas phase ( p g ) is assumed to be the atmospheric pressure; thus, the only independent variable in Eq. 1 is the water pressure ( p w ). The relative permeability of water (k rw ) was predicted based on the widely used Mualem-van Genuchten Model (Mualem 1976;Van Genuchten 1980) and as follows: where l is the pore connectivity parameter assumed as 0.5, the average value based on an analysis of 45 soil samples (Mualem 1976). Thermal equilibrium is assumed which indicates that the solid, liquid and gas (air in the Richards-based models) phases have the same temperature within a REV. Therefore, the energy conservation can be described as: where the (ρC) α is the effective heat capacity of solid, liquid and gas and can be expressed as: where C s , C w and C g are volumetric heat capacity of solid, liquid water and gas, respectively (J kg −1 K −1 ), ρ s is the density of the solids (kg m −3 ), and ρ g is the density of the gas (kg m −3 ). The λ T is the effective thermal conductivity (W m −1 K −1 ) (de Vries 1963). The effective thermal conductivity (λ T ) is calculated based on the model proposed by Johansen (1977): where λ sat is the thermal conductivity at full saturation and can be estimated by λ sat = λ where λ s is the thermal conductivity of soil grains related to the quartz content, and equal to the thermal conductivity of quartz (6.5 W m −1 K −1 for this study), λ w is the thermal conductivity of water (0.58 W m −1 K −1 used in this study). λ dry is described by . K e is known as the Kersten number and is described by Côté and Konrad (2005): where κ is an empirical fitting parameter and is fitted as 3.55 for sand.

Top Boundary Conditions (TBCs)
Two top boundary conditions are used for the Richards model (Eqs. 9 and 14). The first top boundary condition (TBC 1) (Eq. 9) imposes the potential evaporation flux (Eq. 10) (Neumann-type boundary condition) at the soil surface until the surface water pressure reaches the critical water pressure ( p crit , negative). At this threshold pressure, the water in the subsurface soil cannot sustain the potential evaporation flux and the top boundary condition switches to the constant water pressure boundary condition (Dirichlet boundary condition). As will be discussed later, the definition of the critical water pressure is oftentimes debated.
The potential evaporation flux (E pot , kg m −2 s −1 ) can be deduced from meteorological variables which is defined as: where ρ v,sat(z=0) is the saturated vapor density at the surface and ρ va(z=z ref ) is the vapor density in the air (kg m −3 ) at the reference height (z ref , m) above the surface where wind speed, air temperature and air humidity are measured. The saturated vapor density is a function of temperature which can be calculated by: where the temperature variable (T , K) in this equation is the temperature at the soil surface (T surf ) for saturated vapor density at the surface (ρ v,sat(z=0) ) and the temperature at the reference height (T ref ) for saturated vapor density at the reference height (ρ v,sat(z=z ref ) ). The surface temperature is obtained from experimental measurements in this study.
The vapor density in the air (ρ va(z=z ref ) ) is expressed as: where R H a(z=z ref ) is the measured air relative humidity at the reference height which is 12 cm above the soil surface in this work. The 12 cm is arbitrary as it is limited by the lab environment. A typical reference height of 2 m is usually used in field measurements and large-scale models. r a (m s −1 ) is the aerodynamic resistance parameter which is used to describe the resistance of the vapor transfer in the atmospheric boundary layer, characterized based on the Monin-Obukhov similarity theory (e.g., Camillo and Gurney 1986;Griend and Owe 1994;Yamanaka et al. 1997;Sakai et al. 2011;Zeng et al. 2011a): where d is the zero-plane displacement height (m) (0 for bare soil) and z 0 is the roughness length (0.001 m used in this study). The summarization of z 0 used can be found in Forsythe (2017). k is the von Karman constant 0.41. u a is the wind speed at the reference height (m s −1 ). ψ 1 and ψ 2 are stability correction factors depending on the stability of the atmospheric conditions. According to the calculation, little difference in the value of r a was found under the experimental condition, so the ψ 1 = ψ 2 = 0 which indicates the neutral condition was used in the present study. Readers can refer to Camillo and Gurney (1986), Griend and Owe (1994) for more details about stability factors. It should be noted that although the aerodynamic resistance formulation used here has been widely used in various applications, studies show that vapor transport from the surface to the atmosphere is also affected by the surface water content (Shahraeeni et al. 2012); thus, the aerodynamic resistance formulation needs to be amended along with evaporation. Interested readers can refer to Haghighi et al. (2013), Haghighi and Or (2015) for more discussions about vapor transport through the air boundary layer. An alternative top boundary condition (TBC 2) that is also used with Richards models is by introducing a soil resistance term to account for the reduction of the evaporation rate when the surface soil dries out and vapor pressure at surface is lower than the saturated vapor pressure (Eq. 11) (Kondo et al. 1990;Griend and Owe 1994;Saito et al. 2006) rather than switching to a constant pressure head boundary condition, it is expressed by: where r s is the soil resistance term (s m −1 ) representing the resistance of water vapor diffusion to the soil surface. Because of the significant effect of the soil resistance term on the prediction of evaporation behavior, many efforts have been made to formulate a universal and effective soil resistance, ranging from empirically based (e.g., Camillo and Gurney 1986;Sellers et al. 1992;Griend and Owe 1994;Yamanaka et al. 1997;Swenson and Lawrence 2014) to the physically based (Tang and Riley 2013a). Because of the differences in both soil and experimental setups, the formulas vary largely and only applicable to specific soil types. The soil resistance formula used in the present study is proposed by Griend and Owe (1994), tested by Bittelli et al. (2008) and has been used by many numerical studies (e.g., Saito and Simunek 2009;Smits et al. 2011;Zeng et al. 2011b; Mohanty and Yang 2013): where r sl is the resistance to molecular diffusion across water surface (10 s m −1 ). β is a fitting parameter (35.63 based on the experimental data in Griend and Owe (1994)). θ top is defined as the water content in the top 1 cm layer (Griend and Owe 1994). Defining the thickness of the top surface layer is a problem as the soil resistance parameter is highly dependent on the chosen thickness. Interested readers can refer to Fetzer et al. (2017) for more discussions on the effect of the chosen top layer thickness on the simulation results. Equations 9 and 14 are referred to as Top Boundary Condition 1 and 2 (TBC 1 and TBC 2) in the following analysis. It should be noted that the TBC 1 is employed in most scenarios in Sect. 4 when the simulation results of Richards models are compared with other models unless mentioned in particular. Measured soil surface temperature is used as the TBC for heat transfer. The independent variables in the Richards model are water pressure ( p w ) and temperature (T ).

Governing Equations
The Richards vapor model considers liquid water and water vapor flow. This model assumes the gas phase pressure is constant and uniform, and the vapor pressure is in equilibrium with the liquid phase pressure. It is described as: where ρ v is the vapor density (kg m −3 ), θ v is the water vapor content (m 3 m −3 ) and can be expressed as The vapor density is calculated as: As it is assumed that the vapor is in equilibrium with the liquid water, the relative humidity (RH) can be expressed as Edlefsen and Anderson (1943): where M w is the molecular mass of water (kg mol −1 ), R is the universal gas content (R = 8.314 J mol −1 K −1 ). However, the value of RH calculated based on Eq. 18 remains nearly 100% within the capillary pressure range under most conditions. For example, the calculated R H is 99% when the capillary pressure is 1.5E6 Pa (150 m) which corresponds to the permanent wilting point of plants (Griend and Owe 1994). The R H does not start to decrease considerably until the capillary pressure is over 1E7 Pa (1000 m). D v is the vapor diffusivity in porous media which is expressed as Moldrup et al. (2000): where D a is the gas diffusion coefficient in free air (m 2 s −1 ) and it is a function of temperature D a = 2.92 × 10 −5 × (T /273) 2 , θ g is the gas content (m 3 m −3 ), τ is the tortuosity of the porous media. Several different definitions of τ are reported in the literature. The τ used in this study is characterized in Deepagoda et al. (2016) for the sands used in this study as shown in Table 2.
In addition to TBCs discussed above (Eqs. 9 and 14), another TBC can be used with the Richards vapor model when the vapor pressure at the soil surface can be obtained based on the local equilibrium assumption: where ρ vs(z=0) is the vapor density on the surface calculated based on Eqs. 11 and 18. There is no limit for the water pressure at the soil surface. When the soil surface dries out, the surface water pressure decreases very fast and the vapor pressure on the surface will become unsaturated, which causing the flux estimated by Eq. 20 to decrease. The energy conservation equation is expressed by De Vries (1958), Saito et al. (2006): where L 0 is the volumetric latent heat of vaporization of liquid water (J m −3 ) which is given by L 0 = L w ρ w , and L w is the latent heat of vaporization (J kg −3 ), and it can be described as a function of temperature, L w = 2.501 × 10 6 − 2369.2 × (T − 273.15). Same with Richards model, the surface temperature is used as the top boundary condition for heat transfer. The independent variables (i.e., water pressure, temperature) are the same with that of the Richards model.

Governing Equations
In this model, the liquid water flow, gas flow and vapor transport (advection along with air, diffusion) are all included in the mass balance equations. The mass balance for liquid water flow is described as: whereṁ is the phase change rate between the liquid and vapor phases due to evaporation or condensation (kg m −3 s −1 ). The mass balance for the gas phase is as follows: where ρ g is gas density (kg m −3 ). The gas phase is a mixture of dry air and vapor. And the gas density is derived by assuming the dry air and vapor both follow the ideal gas law which can be written as (Lu and Likos 2004): where M a is the molecular weight of dry air (kg mol −1 ). The vapor is part of the gas phase and can be expressed as Bear (2013): where the velocity of the gas phase is: Neglecting gas slippage, relative permeability of air (k rg ) is described by Parker et al. (1987): The energy balance equation in soil could be written as: The phase change rate is proportional to the difference between the equilibrium vapor density and the local vapor density, and the difference between local water content and residual water content (e.g., Benet and Jouanna 1982;Bixler 1985;Lozano et al. 2008;Halder et al. 2011;Nuske et al. 2014): The equilibrium time can be roughly estimated by t eq = L 2 c /D v , in which L c is the characteristic scale of the porous media and D v is the diffusivity. Therefore, the range of b can be estimated based on the properties of the porous media. The top boundary condition for vapor transport is as follows: where ρ v(z=0) is real-time simulated vapor density at the soil surface (kg m −3 ). Because vapor transport is simulated independently, the ρ v(z=0) is no longer the equilibrium vapor density; thus, no soil resistance (r s ) is needed. The top boundary condition for the gas flow is the atmospheric pressure. The independent variables in the Non 2-2 model are water pressure ( p w ), gas pressure ( p g ), vapor density (ρ v ) and temperature (T ). Differing from the system of partial differential governing equations and boundary condition formulations, all the three models in two-dimensional domain (same size with the tank) are implicitly solved using the COMSOL Multiphysics software which is based on the finite element method. The model domain is discretized by triangle elements and the average element size is about 4.5 mm with a total 6800 elements for most scenarios (unless otherwise specified). The model inputs include measured surface soil temperature (T s ), ambient temperature (T a ), ambient relative humidity (RH a ) and wind speed (u a ).

Experimental Setup
To evaluate the three theoretical models presented above, a soil tank/boundary layer wind tunnel experimental apparatus was employed, which was equipped with a sensor network to continuously monitor moisture-temperature-humidity variables in both the free flow and porous media domain. Two laboratory experiments were conducted with coarse (#12/20) and fine sand (#50/70), as shown in Table 1. Experiments were run simultaneously with the same open channel wind tunnel so that the ambient conditions were the same (i.e., air temperature, humidity and wind speed). Precision data for wind speed, ambient temperature and humidity at the reference height were generated as inputs to the models. Also, soil moisture and evaporation flux data were generated to compare with the simulation results for the purpose of evaluating and validating the models.

Soil Material
Two types of uniform specialty silica sand #12/20 (coarse) and #50/70 (fine) (Unimin Corp., Ottawa, MN), which are identified by the effective sieve number, were used in this study. Both sands are rounded with a composition of 99.8% quartz, and the grain density is 2.65 g cm −3 . The uniform coefficient is 1.23 for the #12/20 sand and 1.20 for the #50/70 sand. The measured θ w − p c relationship and the respective fitting curves by van Genuchten model (Eq. 3) are shown in Fig. 1. The fitting parameters (i.e., van Genuchten α and n) and selected properties of the sands are summarized in Table 2. Interested readers can refer to the following for detailed description of measurements (e.g., Sakaki and Illangasekare 2007;Smits et al. 2010;Deepagoda et al. 2016). Figure 2 shows the relationship between effective water saturation (S ew ), relative permeability (k rw ) and the capacity (dθ w /dp w ) as a function of the absolute water pressure (| p w |) for the #12/20 and #50/70 sands based on their respective hydraulic properties. The relative permeability (k rw ) is estimated based on the Mualem-van Genuchten model (Eq. 4). It can be seen that the effective water saturation (S ew ), relative permeability (k rw ) and the capacity (dθ w /dp w ) decreases dramatically when the absolute water pressure (| p w |) exceeds some value (i.e., 1000 Pa for #12/20 sand and 4000 Pa for #50/70 sand). The values shown in the figure are used for analysis in Sect. 4.2.  Fig. 2 The effective water saturation (S ew ) and relative permeability (k rw ) as a function of the absolute water pressure (| p w |) for the #12/20 sand (a), and the #50/70 sand (c); the capacity (dθ w /dp w ) as a function of absolute water pressure (| p w |) for the #12/20 sand (b), and the #50/70 sand (d)

Experimental Apparatus
Experiments were conducted in a rectangular soil tank constructed with plexiglass (45 cm long, 30 cm tall and 9 cm wide) connected to an open channel low-speed wind tunnel, as shown in Fig. 3. The wind tunnel was used to control the surface boundary condition (i.e., wind speed, air temperature). Both the tunnel and soil tank were equipped with a network of sensors to continuously monitor environmental changes in the free flow and porous media.
In the soil tank, moisture and temperature were continuously monitored using dielectric soil moisture sensors (ECH2O EC-5, Decagon Devices Inc., accuracy ±3%) and thermistors (EC-T, Decagon Devices Inc., accuracy ±0.5 • C. The sensors were installed from 2.5 cm below the ground surface with increments of 5 cm depth. Four relative humidity and temperature sensors (EHT relative humidity & temperature sensor, Decagon Devices Inc., accuracy ±2% from 5 to 90% RH, ±3% from 90-100% RH) were placed in good contact with the soil surface. Temperature, moisture and humidity data were collected using five-channel data loggers (ECH2O EM-50, Decagon Devices Inc.). The tank was placed on a scale (65 kg, Sartorius Corporation, resolution ±1 g) to continuously record the weight of the tank. Therefore, the cumulative evaporation and evaporation rate could be accurately calculated by the mass loss. Two tanks with the same size were used in the experiment to run experiments in parallel as mentioned above and in Table 1. The wind tunnel was constructed out of galvanized steel ductwork. A 15.2-cm-diameter duct fan (Suncourt, Inc.) along with a speed controller (Suncourt, Inc.) was installed in the downstream side of the wind tunnel to control the wind speed across the soil surface.
The wind speed in the free flow was monitored using a pitot-static tube (Model 167-12, Dwyer Instruments Inc). The pitot tube was connected to a differential pressure transmitter (Model PX 653, OMEGRA Engineering Inc.), which then was connected to the data acquisition system (Model USB-6218, National Instruments Corp.). The free flow relative humidity (RH a ), temperature (T a ), air pressure ( p a ) and vapor pressure were continuously monitored (VP-4, Decagon Devices Inc., accuracy of relative humidity: from ± 2% to ± 5% for temperature between 20 • C and 40 • C, accuracy for temperature: ± 0.2 • C, accuracy for air pressure: ± 0.4 kPa, accuracy of vapor pressure: from ± 0.06 to ± 0.49 kPa for temperature between 20 and 40 • C).

Experimental Procedure
The soil tank was wet-packed using deionized water and sand. The weight of the two types of sand used in each experiment was first determined based on the prescribed porosity. The deionized water was poured into the tank in small 3-5 cm increments, followed by 1 cm thick sand increments while always keeping at least 1 cm water over the sand to avoid trapped air. The sand in the tank was first tamped, and then the tank was tapped at the walls to achieve uniform bulk density and a tight compaction in accordance with the procedure outlined in Sakaki and Illangasekare (2007). The process was repeated until all the measured sand filled the tank. Then a plastic sheet was temporarily used to cover the tank surface to avoid evaporation; the plastic sheet was removed at the start of the experiment. The valve at the bottom of the tank was closed throughout the experiment so water loss only occurred from the soil surface. Above the soil surface, plastic panels were used adjacent and above the ductwork to channel the air properly across the soil surface and achieve laminar flow conditions. The pitot tube and VP-4 (i.e., humidity, temperature, air pressure and vapor pressure) sensors were installed in the free flow 12 cm above the soil surface. The tank weight, water content, temperature, relative humidity and wind speed were continuously monitored every 1 hour. The experiments were started from a fully saturated condition and continued for 20 days. Because of the change in the free flow conditions (i.e., temperature and relative humidity in the air), the experiments could not be directly repeated/replicated and were therefore only conducted once. However, the experimental procedure and setup are well documented (e.g., Smits et al. 2011Smits et al. , 2012Davarzani et al. 2014;Moradi et al. 2015;Trautz et al. 2015). In addition, although not shown below, the experimental results were compared to similar experiments conducted with the same experimental apparatus and showed consistent trends in moisture, temperature and evaporation behavior. Therefore, we believe that the well specified experimental procedure and well-characterized experimental material can assure the reliability of the experimental measurements. In addition, multiple sensors at the same depth are installed in the tank and free flow for comparison and repeatability of measurements.

Results and Discussion
This section first presents a brief summary of key experimental results. This is followed by a comparison of the three models by evaluating their performance against experimental measurements. Figure 4 shows the measured free flow, soil surface and evaporation behavior for the two experiments (#12/20 and #50/70). This includes the temperature and relative humidity in the free flow and at the soil surface as well as cumulative evaporation, allowing for comparison of the experimental results between the two experiments. As these two experiments were conducted in the same wind tunnel simultaneously, the temperature and relative humidity in the ambient air are identical. As can be seen from the figure, the ambient air temperature (T a ) remains relatively stable, while the ambient relative humidity (RH a ) fluctuates and is always lower than the relative humidity on the soil surface (RH s ).

Experimental Observations and Discussion
The evaporation is strongly associated with the dynamics of temperature and relative humidity at the soil surface (T s and RH s ) and the two soils present similar patterns. For each Fig. 4 Observed cumulative evaporation (cumulative E), soil surface (T s ) and ambient air (T a ) temperature, surface (R H s ) and ambient air (R H a ) relative humidity for the #12/20 sand and #50/70 sand experiments soil, the stage I, transition and stage II evaporation can be clearly seen from the cumulative evaporation (cumulative E) curve. The temperature on the soil surface (T s ) is lower than the temperature in the ambient air (T a ) (approximately 1-5 • C), especially at the beginning during stage I which illustrates the energy consumption effect during evaporation (Shahraeeni and Or 2010). Over time, the temperature difference between the soil surface and air decreases; this is then followed by a stable soil surface temperature for the remainder of the experiment (approximately 10-15 days). Initially, the temperature difference between the soil surface and ambient air is 4-5 • C which suggests that the assumption of using ambient air temperature as the surface temperature is not appropriate when the evaporation rate is high. On the soil surface, the relative humidity is relatively stable during stage I evaporation followed by a steep decrease during the transition and stage II evaporation for the two soils.
The experimental results of the two soils display both similarities and differences. The cumulative evaporation curves of the two soils initially coincide because they share the same ambient condition. The curves deviate over time as the soil properties rather than the ambient condition play a dominant role during the transition and stage II evaporation. The #50/70 sand has a longer stage I duration and a more gradual transition compared with the #12/20 sand as shown in the figure and can be attributed to the differences in the characteristic length, or the maximum drying front depth of the two sands. Through computation, the characteristic length for #12/20 (coarse-textured) and #50/70 (fine textured) sand is 3.4 and 10.8 cm, respectively. The similar pattern can be seen from the curves of temperature and relative humidity on the soil surface (T s and RH s ). This section only presents a few key experimental results. Additional results and comparisons such as the soil moisture and pressure behavior will be further presented in subsequent sections.

One-Component One-Phase, Richards Model
This section evaluates the application of the two top boundary conditions (TBCs) in the Richards model compared with the experimental data from both the coarse and fine sand experiments.

Evaluation of TBC 1
As explained above, TBC 1 (Eq. 9) is the threshold boundary condition in which a flux boundary condition is initially used until the water pressure decreases to the critical water pressure. At that threshold, the flux boundary condition switches to a constant water pressure condition. The selected critical water pressure is oftentimes debated. In theory, it should correspond to a small hydraulic conductivity and capacity (dθ w /dp w ) so the simulation results remain relatively constant when the water pressure continues to decrease . In practice, different values varying from (−1E5 to −1E7 Pa) are used in the literature (e.g., Bechtold et al. 2012;Assouline et al. 2014;Fetzer et al. 2017).
Based on the hydraulic properties and capacity shown in Fig. 2a and b, a water pressure of −1400 Pa, for example, corresponds to an effective water saturation of 2.0E−2, relative permeability of 2.2E−5 and capacity (dθ w /dp w ) of 4.5E−5 (Pa −1 ) as marked in the figure.
The water pressure of −1400 Pa was chosen as the starting critical water pressure for the simulation based on the knowledge of the soil water retention behavior. We then selected additional critical water pressures from −1400 to −2000 Pa to show the impact of the critical water pressure on the evaporation behavior in Fig. 5a. For each selected critical water pressure, the corresponding S ew , k rw and (dθ w /dp w ) is shown in Fig. 2a and b. To better observe the change of simulated cumulative evaporation along with the chosen critical water pressures, Fig. 5b presents a snapshot of the cumulative evaporation at 20 days based on the data in Fig. 5a as a function of the critical water pressure. The cumulative evaporation firstly increases with the absolute critical water pressure (| p crit |) and then becomes stable. With the increase in the absolute critical water pressure (| p crit |), the Richards model predicts a slightly longer duration of stage I and a larger magnitude of cumulative evaporation. This relation between the duration of the stage I evaporation and the selected critical water pressure is closely related to the hydraulic property of the soil (k rw and dθ w /dp w shown in Fig. 2). In the Richards model, the evaporation rate only depends on the ability to transport liquid water to the soil surface, which relies on the pressure gradient and unsaturated hydraulic conductivity near the surface. As the surface dries, the water pressure grows more negative which causes competing effects between a larger pressure gradient and lower hydraulic conductivity near the surface (Salvucci 1997). The stage I evaporation ends when the soil dries to a point where the effects of the corresponding reduction in unsaturated hydraulic conductivity exceed the effects of the increase in pressure gradient.
Compared with the observed data, the Richards model with relatively low critical water pressures (−1650 Pa to −2000 Pa) agrees well with the observed evaporation across the entire range of water saturation, even though the model consistently underestimates the magnitude of stage II evaporation. Therefore, the critical water pressure (| p crit |) should be defined from an analysis of the behavior of the cumulative evaporation with changing | p crit | and use the | p crit | that gives the highest cumulative evaporation. It should be noted that uniform sands with large van Genuchten parameters, n, are used both in the experiments and simulations. For field soils such as silts, sandy loam, more negative critical water pressures are appropriate as shown in Fetzer et al. (2017).
The underestimation of the stage II evaporation rate using the Richards model is primarily due to the models' inability to simulate vapor diffusion. When the vapor diffusion is considered in the Richards model, a higher stage II evaporation rate is observed in the coarse sand scenario shown in Sect. 4.5. As can be seen in Fig. 1, the van Genuchten model fit diverges at low water contents as the fitted soil moisture remains unchanged while the observed soil moisture continues to decrease slightly. Correspondingly, the Mualem-van Genuchten underestimates the unsaturated hydraulic conductivity. In addition, the Mualemvan Genuchten model assumes cylindrical pores and disregards films and corner flow, which has been demonstrated to underestimate the unsaturated hydraulic conductivity especially in the low water content region Or 2001, 2005). Follow-on work investigating the impact of film and corner flow on the simulation results is needed to fully estimate the performance of the Richards model.
An equivalent analysis as presented above was performed for the #50/70 fine sand to demonstrate the influence of the soil hydraulic properties on the system behavior. Although the results are not shown here, the optimal critical water pressure is between −6000 Pa to −7000 Pa. Similar conclusions can be drawn based on the simulation results of the #50/70 sand compared with the #12/20 sand.

Evaluation of TBC 2
Differing from TBC 1 (Eq. 9), TBC 2 (Eq. 14) introduces a soil resistance parameter (r s ) to account for the evaporation reduction as the soil surface dries and the vapor pressure at the surface is no longer saturated. The soil resistance parameter (r s ) increases during the drying process, which results in a decreasing flux (Eq. 14), rather than the constant potential evaporation flux as in TBC 1. Simulation results comparing TBC 2 (r s < 1235 s m −1 and r s < 665 s m −1 for the coarse and fine sands, respectively, r a = 213 s m −1 for both scenarios) with TBC 1 ( p crit = −2000 Pa and −7000 Pa for the coarse and fine sands, respectively) can be seen in Fig. 6. In general, TBC 2 leads to an underestimation of the cumulative evaporation for both sands under the current conditions mainly due to a shorter stage I duration. Comparing the TBC 1 and 2 results, TBC 2 predicts a more gradual transition stage which agrees better with the observed transition. It is important to note that the prediction capability of TBC 2 is associated with the soil properties, soil resistance formulation and magnitude of the aerodynamic resistance. There are abundant discussions about the formulation of r s and how

Richards Vapor Model
TBC 1, 2 and 3 are all applicable to the Richards vapor model. As TBC 1 and 2 were compared above, this section will focus on the evaluation of TBC 1 and 3.
The relative humidity remains close to 100% (saturated vapor pressure) under the normal capillary pressure range for sands. If a capillary pressure within the normal range is imposed at the surface as the critical water pressure for TBC 1, there will be no vapor gradient within the porous media domain. Therefore, a series of very low critical water pressure (−1E6 to −5E7 Pa) for TBC 1 is used in the simulation so the vapor pressure calculated based on the equilibrium approach is smaller than the saturated vapor pressure. When such very low critical water pressures (−1E6 to −5E7 Pa) are employed in TBC 1, the large gradient that emerges near the surface can cause numerical artifacts. Therefore, different grid sizes are employed in the simulation with varying critical water pressure of TBC 1 to avoid the numerical artifacts. Figure 7a presents the effect of the grid sizes on the simulation results with a critical water pressure of −1E7 Pa (TBC 1) as an example. Although not shown here, the simulated cumulative evaporation results with other several critical water pressures (−1E6, −5E6 and −5E7 Pa) show the similar trend, and all the simulation results stabilize when the element size is approximately 0.25 mm.
TBC 3 (Eq. 20) is based on the vapor density gradient between the soil surface and the air. The surface vapor density is a function of the surface capillary pressure and starts to decrease rapidly when the surface water pressure decreases to approximately −1E7 Pa (R H = 0.93). The decrease in vapor density results in a decreasing flux boundary condition similar to TBC 2. However, in TBC 2, the decreasing flux depends on the decreasing water content of the top surface layer rather than the water pressure. This is an important distinction as sometimes the change in water content and water pressure are not simultaneous, especially when the water content approaches the residual value. Figure 7b shows the simulation results with TBC 3 (Eq. 20) with different grid sizes; the simulation results also stabilize when  Figure 7a, b indicates that the stage I evaporation is not sensitive to the space discretization, while the discretization influences the transition and stage II evaporation prediction. It should be noted that such analysis has also been conducted in the Richards model, and the simulation results do not change with a decreasing grid size (4.5 mm is used in the Richards model and is chosen as a starting grid size in this section).
With the same discretization (grid size = 0.25 mm), Fig. 8 presents the comparison of the simulated cumulative evaporation by TBC 1 with different critical water pressures (−1E6, −5E6 and −5E7 Pa) and TBC 3 for #12/20 coarse and #50/70 fine sands. It can be seen that the simulated cumulative evaporation increases as the critical water pressure decreases from −1E6 to −5E7 Pa. And the simulation results of TBC 1 with a critical water pressure of -5E7 (corresponds to an R H of 0.7) is identical to the simulation results of TBC 3. It should be noted that for TBC 3, the lowest water pressure emerged at the soil surface is approximately −2.2E8 Pa. This value corresponds to an R H of 0.2 which is close to the averaged air R H (R H a = 0.17). Similarly, in the numerical modeling software Fig. 9 Using the Richards vapor model with TBC 3 for #12/20 sand: a simulated surface flux as a function of time; b simulated corresponding near-surface pressure gradient (between surface and 2 mm depth) and average relative permeability (between surface and 2 mm depth) Hydrus, the minimum water pressure head allowed at the surface is set to a value, which is calculated based on the air humidity using the equilibrium approach between soil water and water vapor (Simunek et al. 2005). This indicates that the critical water pressure can be approximately chosen as the value which corresponds to the air humidity based on the equilibrium approach. However, as the air humidity changes over time, this minimum value is also dynamic. Figure 9a shows the simulated surface flux by the Richards vapor model with TBC 3. It can be seen that the liquid flux decreases rapidly to zero, while the vapor flux increases and dominates the evaporation in the transition and stage II evaporation. The rapid decrease in liquid flux can be explained with Fig. 9b, which shows the near-surface water pressure gradient and the relative permeability. The relative permeability decreases sharply to 1E−16 which leads to a nearly zero liquid water flux and an transition to stage II evaporation.
In summary, TBC 1 with a very negative critical water pressure or TBC 3, both with fine grids, can be used in the Richards vapor model. As the use of fine grids brings great challenges to the computational cost, a pre-analysis of the optimum meshing should be done at first. In addition, the convergence of the numerical model cannot always be assured even though with very fine grids. A free triangle mesh was found to produce numerical oscillations in the simulations of this section, so structured quadrilateral mesh was used. In addition, the time-stepping method in COMSOL allows the solver to adjust the timestep size to boost efficiency as well as satisfying the tolerance. However, how to balance speed, robustness and accuracy when solving Richards equation numerically is beyond the scope of this study. Interesting readers can refer to Farthing and Ogden (2017) for a review of the advances and challenges in solving Richards equation numerically.

Non-equilibrium Phase Change (Non 2-2) Model
The Non 2-2 model considers the phase change from liquid water to water vapor that occurs within the porous media; the local equilibrium assumption is not used. Instead, the vapor pressure is treated as an independent variable in the simulation. The sensitivity of the phase change rate and the surface flux are discussed in this section. Only one top boundary condition (Eq. 30) is introduced and discussed in this study. Interested readers can refer to Smits et al. (2012) for more top boundary conditions and their effects.

Sensitivity of the Simulated Evaporation to the Phase Change Rate in the Non-equilibrium Model
The phase change rate coefficient in Eq. 29 determines how quickly the vapor concentration in the air will reach equilibrium. Figure 10 shows the observed and simulated cumulative evaporation using different phase change rates, adjusted via the b parameter in Eq. 29. It can be seen that the Non 2-2 model can generally capture the full range of evaporation. However, the b parameter in phase change rate formulation has a marked effect on the simulation results, mainly during stage I evaporation. An increase in the phase change rate coefficient results in an increase in the slope of the cumulative evaporation curve at early times (i.e., during stage I) but has little effect of the evaporation during stage II. This corresponds with the results of previous work (Smits et al. 2011). Based on the equilibrium time analysis which is related to the pore size distribution, the range of b is between 10 −3 and 10 −2 (s m −2 ) for coarse #12/20 sand, and at the magnitude of 10 −1 (s m −2 ) for fine #50/70 sand (Halder et al. 2011;Deepagoda et al. 2016). However, for the #50/70 sand, the calibrated b (b = 0.005) based on the experimental data is not within the estimated range from theoretical analysis, which indicates the difficulty of defining the proper phase change rate formulations when the experimental data is not available. Theoretically, equilibrium is reached when the mass transfer coefficient (or b coefficient) is infinitely large. But numerically, b cannot be infinite because the phase change rate, which is proportional to b, serves as the source/sink for the liquid/gas phase equation and this messes up the math. This has also been confirmed numerically, the convergence of the models cannot be maintained when continue increasing b. Further work focusing on developing phase change rate formulation based on limited fluid and porous media information is needed to strengthen the application of this model. For a more thorough review of the effect of the b parameter on evaporation and non-equilibrium phase change formulations, the reader is referred to Smits et al. (2011) andTrautz et al. (2015).

Surface Flux
The contribution of liquid water and water vapor flux over time using the Non 2-2 model (b = 0.01) is shown in Fig. 11 for the coarse sand scenario. During stage I, the evaporation is mainly driven by the liquid water flux at or close to the potential evaporation rate. However, when the surface starts to dry and the liquid water flux decreases, vapor flux becomes the dominant contributor to the total flux. The increase in the vapor flux at t = 0-5 days is because of the increase in the diffusivity as shown in Fig. 12a. At approximately day 5, the surface dries out, reaching the residual water saturation (S rw ) as indicated in Fig. 12b, and the liquid water flux goes to zero as shown in Fig. 11. After day 5, the drying front begins to migrate or develop below the soil surface, and the dry surface layer (DSL) thickens. For example, the thickness of the dry surface layer at t = 10 and 15 days is 1 and 2 cm, respectively. Although not shown, this is consistent with the experimental observed data. After day 5, even though the diffusivity stays at a high value, the vapor flux decreases (Fig. 12a) because of the increase in the thickness of the dry surface layer (DSL). And the thicker DSL results in the increase in the diffusive path for the vapor.

Overall Model Comparison with Observed Data
In this section, the three models are evaluated together with the experimental data in terms of cumulative evaporation and state variables (i.e., water pressure, water saturation, vapor density). The simulation results of the Richards and Richards vapor models are influenced by the top boundary conditions and the associated critical water pressure, and their impacts on the simulation results have been discussed above. Therefore, for consistency in the comparisons, the setting with the best fit with the observed cumulative evaporation for each model is adopted in this section for comparison unless otherwise specified. That is to say, the TBC 1 with the optimized critical water pressure (−2000 Pa for #12/20 coarse sand and −7000 Pa for #50/70 fine sand) are used in the Richards model. The simulated cumulative evaporation of TBC 1 with a critical water pressure of −5E7 Pa and TBC 3 for the Richards vapor model are identical and fit the observation well, so both TBCs with a fine meshing are included in this section . For the Non 2-2 model, the calibrated phase change rate from above is used for the model comparison (0.01 for #12/20 coarse sand and 0.005 for #50/70 fine sand). Figure 13 shows the observed and simulated cumulative evaporation using the three numerical models (Richards, Richards vapor and Non 2-2 models) for the #12/20 coarse and the #50/70 fine sands with the best setting for each model. It can be seen that all the approaches predict similar trends with small differences. All the models predict a similar and longer stage I duration for both sands. In terms of stage II, the models vary in their ability to properly capture the slope for the #12/20 coarse sand. Including the vapor transport, the Richards vapor model predicts a larger stage II evaporation compared with the Richards model for the #12/20 coarse sand scenario. However, based on the simulation results of the fine #50/70 sand, the three models give the same stage II evaporation, which indicates that the Richards model can represent the stage II evaporation (with an optimized critical water pressure) under some conditions. All the three models underestimate the stage II evaporation for the #50/70 sand scenarios, implying that there are still limitations in the models' ability to capture all the physics. Again this work is not aimed at perfectly fitting observed data but rather to com-  14 Simulated absolute water pressure on the surface of #12/20 sand as a function of time using the three modeling approaches. The dotted line marks the time of 2.4 days which is also the time when the stage I evaporation ends (Fig. 13a) pare the different model concepts and top boundary formulations, so additional empirical parameters or parameter calibration processes are not employed.
One main difference between the three model approaches can be seen in the simulated water pressure at the soil surface over time. As shown in Sect. 4.4.2, for the non-equilibrium model, the vapor flux plays a role in the total flux even at early stages, contributing to the evaporation flux. In the equilibrium models (i.e., Richards vapor model), the noticeable vapor gradient emerges only when the water pressure is extremely low (Fig. 9). To illustrate this, Fig. 14 shows the simulated absolute water pressure (| p w |) on the surface of #12/20 sand as a function of time by different modeling approaches. As shown in the figure, there is a large and unrealistic jump at 2.4 days which is also the time when the transition stage happens. The vapor gradient difference before and after this pressure jump is shown in the next figure. This pressure jump results from the conflict between the high evaporation flux demanding (i.e., TBCs) and the low water capacity (dθ w /d p w ) as well as the low ability of the soil to transport water to the surface (k rw ) when the water content is at residual. Based on the measured relative humidity data at (or near) the soil surface (vapor density in Fig. 16a), there is no sharp decrease over time which implies that the large pressure jump of the Richards vapor model under the equilibrium assumption is not realistic. Figure 15 shows the vapor density profile of #12/20 sand near the surface (5 cm depth) at two time points (2 and 3 days), simulated by the Richards vapor model (TBC 3) and Non 2-2 model (b = 0.01). These two time points are chosen because the large pressure jump happens between them (2.4 days shown in Fig. 14). At day 2, limited by the equilibrium (liquid water and water vapor) approach, the vapor gradient is small and contributes little to the evaporation flux, which also can be seen from Fig. 9a. After the pressure jump at the surface at day 2.4 ( Fig. 14), the relative humidity at the surface (Eq. 18) is less than 1, creating a larger vapor gradient near the surface (Richards vapor model t = 3 days). By contrast, the Non 2-2 model predicts a relatively larger vapor gradient near the surface at day 2 and the vapor flux starts to contribute to the evaporation flux.
The simulated vapor density at the soil surface and 0.5 cm depth from Fig. 16 also help to illustrate the vapor gradient differences above. As can be seen from Fig. 16, the vapor density difference between surface and 0.5 cm depth predicted by the Richards vapor model grows largely after 2.4 days for #12/20 sand (Fig. 16a) and 5 days for #50/70 sand (Fig. 16b). Figure 16 also presents the measured vapor density which is obtained from the surface temperature and surface relative humidity measurements (Eq. 17). Compared with the measured surface vapor density, the simulated surface vapor density using both the models deviate from Fig. 15 Vapor density profile of #12/20 sand (t = 2 and 3days), simulated by the Richards vapor (TBC 3) and Non 2-2 (b = 0.01) models

Fig. 16
Measured and simulated vapor density at the soil surface and 0.5 cm below the soil surface for the Richards vapor (RV) and Non 2-2 models. In the legend, "Saturated surface" means the calculated saturated vapor density at the surface based on measured the surface temperature (Eq. 11), "surface" means the vapor density is at the surface and "0.5 cm" means vapor density is at the 0.5 cm depth the observation. The Richards vapor model predicts a nearly saturated vapor density at the surface from 0 to 2.4 days for #12/20 sand scenario and 0-5 days for #50/70 sand scenario, which agree well with the observation, followed by a rapid decrease and a steady vapor density after that. On the contrary, the vapor density at the surface predicted by the Non 2-2 model shows a slow decline with time. However, these two models have similar predictions of the surface vapor density for both the scenarios after some time (2.4 days for #12/20 scenario and 5 days for #50/70 scenario), although there is a larger distinction in prediction of surface water pressure (Fig. 14). On the other hand, the figure shows that the simulated vapor density at 0.5 cm depth by the Richards vapor model and Non 2-2 model fit the experimental soil surface data better than the simulation results at the surface (z = 0 cm), which also has been reported by Davarzani et al. (2014). In the experiment, the relative humidity sensor requires good contact with the soil to ensure that the readings are not influenced by the surrounding free flow air. As argued in Davarzani et al. (2014), the relative humidity sensor is reading the relative humidity below the soil surface rather than directly on the soil surface. Figure 17 shows a comparison of the observed and simulated water saturation dynamics using three modeling approaches (Richards, Richard vapor and Non 2-2 models with their Fig. 17 Comparison of the observed and simulated water saturation dynamics using the three modeling approaches at different depths of a #12/20 sand and b #50/70 sand, respectively. The simulated curves of water saturation at 12.5 cm depth by different modeling approaches coincide with the observation best setting same as the setting in Fig. 13) at three different depths (2.5 cm, 7.5 cm and 12.5 cm) for the #12/20 and #50/70 sand experiments. It can be seen that for the #50/70 sand experiment, the simulation results of the three models are similar and generally agree well with the observed saturation as a function of time (MIA values range from 0.967 to 0.996). The over-and underestimation of the saturation at different times and depths is possibly due to the limitation of the van Genuchten model in fully representing the measured soil water retention curve, and the changes of the retention curve caused by local density difference in the soil tank (Assouline 2006). However, for the #12/20 sand scenario, three models vary in their ability in predicting the water saturation at 2.5 cm and 7.5 cm depth. The Richards and Richards vapor models predict identical water saturation of 2.5 cm depth of the whole range and at the beginning (approximately from 0 to 3 days) of 7.5 cm depth. However, the water saturation curves at 7.5 cm predicted by the Richards and Richards vapor model deviate, in which the Richards vapor model can better capture the observed water saturation at 7.5 cm. This is also statistically confirmed by the MIA values, which is 0.971 for the Richards vapor model and 0.816 for the Richards model. The Non 2-2 model has a slight better fitting at 2.5 cm depth compared with the other two models, but underestimates the water saturation at 7.5 cm depth after 5 days (MIA is 0.884).
Based on the simulated water saturation results at different depths between different modeling approaches shown in Fig. 17, it implies that the different cumulative evaporation predictions (Fig. 13) result from the models' ability to transport water from deeper depth to the soil surface when the surface dries, instead of the near-surface water evaporation. It can be better illustrated by Fig. 18, which shows the comparison of the simulated and measured water saturation profiles predicted by the three models and limited data (at 2.5 cm and 7.5 cm depths) for #12/20 sand at two times (1 day and 20 days). The predicted water saturation profiles show little difference between the three models and agree with the experimental data at the two depths at day 1. At 20 days, the simulated near-surface (approximately 0-3 cm depth) water saturation by three models is the same, but the differences between water saturation profiles can be seen at approximately from 3 to 10 cm depth, which is closely associated with the cumulative evaporation. Even though the three curves seem close, the water saturation difference at the same depth can be large (a horizontal line can help present the difference).

Fig. 18
Water saturation profile predicted by three models (Richards, Richards vapor and Non 2-2 models) at a t = 1 day and b t = 20 days. Measured water saturation at 2.5 cm and 7.5 cm depth is also shown in the figure

Conclusion
In this study, three numerical models (Richards, Richards vapor and Non 2-2 models) which are used to describe the water flow and heat transfer in subsurface soil and across the soil-air interface were presented. The models were critically evaluated by comparing with precision data from laboratory experiments performed in a well-controlled wind tunnel equipped with a sensor network capable of monitoring variables in both the free flow and porous media region.
Results show that simulation results of both the Richards and Richard vapor model are affected by the TBCs, and their ability to capture the stage II evaporation vary with scenarios. For the Richards-based model, the critical water pressure in TBC 1 should be defined as the water pressure which gives the maximum cumulative evaporation after sensitivity analysis of the critical water pressure and discretization to the evaporation behavior. TBC 1 with a very negative critical water pressure and TBC 3 are prone to numerical errors when applied in the simulation of uniform sand and special attention is needed in selecting the critical water pressure and discretization. A fine meshing approach helps to enhance the reliability of the simulations but introduces high computational cost. When the critical water pressure is properly selected, in spite of the incapability of Richards model to simulate vapor diffusion, the Richards model can still predict all the evaporation stages which generally fits the observed evaporation for the fine #50/70 sand experiment. However, for the coarse #12/20 sand scenario, the underestimation of stage II evaporation predicted by the Richards model can be clearly seen, which indicates the inapplicability of using limited liquid flow to express stage II evaporation for all soil types. The contribution of vapor transport in the Richards vapor model can be seen under very dry conditions. TBC 1 with a very low critical water pressure and TBC 3 have been successfully applied in Richards vapor model when a fine meshing is used, and the simulation results agree well with the observations. Compared with the Richards-based models, the non-equilibrium two-component two-phase (Non 2-2) model is able to generally capture the evaporation of all stages when a fitting parameter within the phase change rate term is properly chosen. The application of the coarse grid indicates that the non-equilibrium approach helps to reduce numerical oscillations which has also been reported (Bixler 1985).
Besides the evaporation rate, the dynamics of the state variables (e.g., water pressure, vapor density, water saturation) were simulated by the models with the best settings and compared with available experimental data. The Richard vapor model (with TBC 3) and Non 2-2 model both give a better prediction in terms of cumulative evaporation and water saturation dynamics compared with the Richards model in #12/20 scenario, indicating that it is important to include vapor transport to reproduce soil water dynamics in the shallow subsurface. On the other hand, although similarity in modeling predictions between equilibrium and nonequilibrium approach (i.e., Richards vapor and Non 2-2 model) can be seen, the surface water pressure and vapor density profile present large differences, calling for more experimental observations to evaluate the different approaches.
This study evaluates the different modeling approaches in many aspects. All the modeling approaches and corresponding TBC formulations have various limitations and no one can be deemed most appropriate for every scenario. The choice depends on the research focus, the requirement of the accuracy and efficiency and availability of the data set. Future work focusing on the characterization of the complex interactions at the soil-atmosphere interface by the combination of numerical modeling and detailed experimental measurements, incorporating all kinds of soil textures and atmospheric conditions, is needed to improve the predictions of models across scales.