Capillary, Film, and Vapor Flow in Transient Bare Soil Evaporation (1): Identifiability Analysis of Hydraulic Conductivity in the Medium to Dry Moisture Range

Evaporation experiments are frequently used to determine soil hydraulic properties. We simulated laboratory evaporation experiments with a coupled water, vapor, and heat flow model which includes the surface energy balance. The simulations are performed with different parametrizations of soil hydraulic properties with a focus on soil hydraulic conductivity in medium to dry soil. In previous studies, conductivity in this moisture range has been shown to be influenced not only by water flow in completely filled capillaries (“capillary flow”) but also by film and corner flow (“film flow”). Our forward simulations highlight the strong influence of an increased conductivity caused by film flow on evaporation rate, cumulative water loss, soil temperature, and soil water pressure head during evaporation. Film flow extends the duration of stage‐1 evaporation and increases the evaporation rate during stage‐2 even if all other physical material properties are the same. The simulated data were used in inverse simulations with the Richards equation to test whether soil hydraulic properties can be identified without bias. This is a priori questionable because the Richards equation is an isothermal flow model and simplifies the true physics considerably, by ignoring thermal liquid and thermal vapor fluxes, as well as temperature effects on the hydraulic properties. Our results show that the identification of the water retention and hydraulic conductivity curves is bias‐free for media with and without film flow. We conclude that the Richards equation can be safely used to identify hydraulic properties from evaporation experiments by inverse modeling.

content (Jury & Horton, 2004). These constitutive relationships will exert an influence on the dynamics of the evaporation rate, soil temperature, soil water content, and soil water pressure head . Conversely, evaporation experiments are frequently used to infer constitutive relationships by a combination of experimentation and (inverse) modeling (Romano & Santini, 1999;Schindler & Müller, 2006;Wendroth et al., 1993).
Evaporation experiments in the laboratory have become a standard method to study the transport of liquid water and water vapor in soils Shokri et al., 2008) and to determine SHP in the laboratory (Bezerra-Coelho et al., 2018). The determination of SHP profits from the simplicity and robustness of the experimental procedure and the reliability of determining the WRC (from full saturation to medium water content) and the HCC (in a more limited moisture range) in high resolution. With suitable instrumentation for matric potential measurement in the medium to dry moisture range, evaporation experiments can also be used to determine SHP towards dryness. However, due to the nonlinear spatiotemporal distribution of soil water state variables during stage-2 Peters et al., 2015), the optimal methodology to determine SHP is an inverse simulation with a flow model capable of describing these nonlinearities.
The current methodology for evaluating evaporation experiments assumes isothermal conditions within the soil column, although the coupling between the water and energy cycle causes bare-soil evaporation to be non-isothermal under most conditions. The resulting temperature gradients in the soil drive thermal fluxes of liquid water and water vapor (Philip & De Vries, 1957) which are often neglected. Moreover, heat fluxes in the soil have an influence on soil surface temperature which affects the saturation water vapor pressure and therefore the evaporation rate. Summing up, heat fluxes influence water fluxes and vice versa. Finally, SHP are known to depend on temperature (Joshi et al., 2019).
The theory of modeling the coupled flow of liquid water, water vapor, and energy is well-established in vadose zone hydrology and its application becomes more widespread, for example, Kelleners et al. (2016). The state-of-the-art model of the coupled flow of liquid water, water vapor, and energy dates back to the classic work of Philip and De Vries (1957) and was later extended by Milly (1982) and Nassar and Horton (1997). A more comprehensive multiphase multicomponent flow model which includes an explicit coupling across the soil-atmosphere interface was presented by Vanderborght et al. (2017). In the same article, the assumptions which are required to derive the coupled soil model of Milly (1982) from the coupled multiphase model are explained. A major difference between these two models is that the multiphase model includes the flow of the gas phase. Sakai et al. (2011) were among the first who applied the coupled water, vapor, and heat flow model to laboratory evaporation experiments. Their primary aim was to analyze the ability to infer subsurface evaporation from the sensible heat balance. Milly (1984) has shown that thermal fluxes of liquid water and water vapor are negligible in field soils for diurnal variations of meteorological variables which drive evaporation from soil and that an isothermal flow model which includes the process of isothermal vapor diffusion can correctly approximate daily evaporation fluxes from bare soil. However, to the best of our knowledge, Milly's fundamental analysis has never been applied to laboratory evaporation experiments, in which unidirectional heat fluxes prevail for extended periods. This is notable, because the Richards equation, which is the de facto standard model to evaluate evaporation experiments (Romano & Santini, 1999;Weber et al., 2017), must, due to its isothermal nature, ignore thermal fluxes of liquid water and water vapor and changes of transport properties caused by temperature variations.
The objective of this paper is to explore the possibility of determining SHP from evaporation experiments by inverse modeling with the Richards equation. We particularly focus on the determination of the HCC in the medium to dry moisture range because film and vapor flow are expected to be most influential in this range. We propose extended instrumentation for matric potential measurement which consists of tensiometers and relative humidity sensors and analyze its information content w.r.t. the determination of hydraulic conductivity. For this purpose, numerical forward simulations with a coupled energy-water transport model which additionally solves the surface energy balance and accounts for the temperature dependence of SHP (Saito et al., 2006) are used to generate data. These data are then evaluated by inverse simulations with the Richards equation and an extended Richards model accounting for isothermal vapor flow.
Specifically, we address the following research questions: (1) Is it possible to match the time series of pressure head obtained with the non-isothermal coupled model with the isothermal flow model (Richards' equation)? (2) If yes, are the SHP identified by inverse modeling correct although the Richards equation provides only a simplified process description? We note that the second point can only be addressed by a process simulation because it ensures that the true material properties are known, which is by principle not possible for experiments on real systems (Peters et al., 2015). In an accompanying article, we apply the methodology tested in this article to real data from evaporation experiments, in which pressure head was measured from almost full water saturation to air dry conditions using tensiometers and relative humidity sensors. The evaluation of the experimental data by inverse modeling in the companion article shows that classic models of soil hydraulic properties cannot describe the observed time series of pressure head and that an adequate match of the observations can only be achieved with a model that contains a "film flow" component in the hydraulic conductivity function.

Materials and Methods
An overview of the workflow of this study is presented in Figure 1. In the next sections, we present (i) the coupled model and (ii) the Richards equation for water flow, the (iii) parametrization of SHP with models of increasing complexity, and (iv) the data generated with the coupled model, and (v) the methodology of the inverse simulations.

Coupled Modeling of Water, Vapor, and Heat Flow
To simulate the coupled flow of liquid water, water vapor, and energy, we used the Hydrus-1D software code (Simunek et al., 2008(Simunek et al., , 2016 which additionally solves the surface energy balance. The code was slightly modified so that the surface energy balance in the lab could be solved and the aerodynamic resistance could be defined by the user. We only provide a short overview of the coupled model in this section and refer the reader to Novak (2016), Saito et al. (2006), Sakai et al. (2011), and the Hydrus manual (Simunek et al., 2013) for more details.
The governing equation for the flow of liquid water and water vapor is (Saito et al., 2006)  The governing equation for heat flow without root water uptake is (Saito et al., 2006)  are the flux densities of liquid water and water vapor, respectively. Since Equations 1 and 2 are coupled, we will refer to the flow model as the "coupled model" in the remainder of this article. The material properties and functions in Equations 1 and 2 are described in Saito et al. (2006), Sakai et al. (2011), andSimunek et al. (2013). A short summary is provided in Table 1 Saito et al. (2006). b Chung and Horton (1987). c Simunek et al. (2013).

Richards Equation Including Isothermal Vapor Flow
If one neglects heat flow in the soil (isothermal water and vapor flow) This equation can be further simplified to the Richards equation if either K vh can be neglected or the effect of the gradient in gravitational potential becomes small relative to the gradient in tension. This occurs when the soil is wet and K lh ≫ K vh , or when the soil is dry and can then be effectively simplified to the Richards equation In Equation 4, the soil water content θ stands for the total water content θ tot to ease notation. The total or effective soil hydraulic conductivity is defined as         lh vh K h K h K h and therefore comprises the isothermal liquid and vapor conductivities. Under the simplifying conditions stated above, the Richards equation, which was originally derived for viscous capillary flow, is extended such that it can additionally model the diffusion of water vapor in the porous medium under isothermal conditions. Note that other authors, for example, Vanderborght et al. (2017), restrict the term "Richards equation" to Equation 4 with a function K(h) excluding the conductivity caused by isothermal vapor flow. Li et al. (2019) refer to the Richards equation including isothermal vapor flow as "Richards vapor model" and formulate the vapor flux by Fick's law without converting it into an equivalent Darcy-Buckingham law with resulting conductivity K vh (h). In this study, the coupled model is given by Equations 1 and 2 is used for the forward simulation of evaporation experiments, and the Richards Equation 4 is used for the determination of the SHP by inverse modeling.

Conceptual Models for Soil Hydraulic Properties
The solution of the partial differential Equations 1 and 4 requires the parametrization of the SHP (Assouline et al., 2013). In this study, we apply four models to parametrize the SHP to answer the two research questions outlined in the introduction. We provide a basic description of the four models in this section. In brief, the models differ in their description of water storage in dry soil ("residual water content" vs. "zero water content at oven-dryness") and the representation of hydraulic conductivity (capillary flow, non-capillary flow, isothermal vapor flow). Table 2 provides a comparison of the basic properties of these models. In the following, hydraulic conductivity is always formulated as a function of the tension h.

M1 -Residual Water Content and Only Capillary Flow
The first model (M1) reflects the current standard approach in the parametrization of SHP for modeling water flow in soils. Water flow is conceptually treated as taking place only in completely filled capillaries, that is, vapor flow and flow in incompletely-filled capillaries (surface films, edges, and corners) are neglected. The WRC curve is expressed by an empirical expression with a "residual" water content. The relative conductivity function is derived from a pore bundle model (Burdine, 1953;Mualem, 1976) which treats the WRC as a representation of effective pore-size distribution. The WRC is parametrized with the van Genuchten model   (van Genuchten, 1980;Mualem, 1976) where K sc [m s −1 ] is saturated hydraulic (capillary) conductivity and τ [-] is a shape parameter, which describes the tortuosity and connectivity of pores (Mualem, 1976). In this standard model, the total hydraulic

M2 -Residual Water Content, Capillary, and Isothermal Vapor Flow
The second model (M2) extends M1 and additionally accounts for isothermal vapor flow. The total hydraulic conductivity is and the hydraulic conductivity due to water vapor diffusion is (Saito et al., 2006)  is the molecular mass of water, g [m s −2 ] is the acceleration due to gravity, R [J K −1 mol −1 ] is the universal gas constant, and T [K] is the absolute temperature. As we deal with packed soil columns in this study, the dependence of the diffusion coefficient on water content is parametrized as (Moldrup et al., 2000)  where D 0 (T) [m 2 s −1 ] is the diffusion coefficient in air, and φ [-] is the porosity of the soil. The relationship between relative humidity (water activity) and tension is given by the Kelvin equation , in which the absolute value of water potential, expressed as head, is assumed to equal the tension h, that is, the osmotic potential is neglected. Equations for the calculation of D 0 (T), ρ vs (T), and ρ w (T) are given in Saito et al. (2006).

M3 -Zero Water Content at Oven Dryness, Capillary, and Isothermal Vapor Flow
Parametric representations of the WRC with a residual water content generally suffer from the drawback that the water content will not become zero at the tension corresponding to oven-dryness,  6.8 0 10 cm h (Schneider & Goss, 2012), if θ r > 0. Note that setting   0 r does not solve this problem in all cases and additionally modifies the HCC, because K(h) approaches zero at θ r . Peters (2013) provides a more comprehensive discussion of why such models are physically inconsistent.
Model M3 solves this problem and ensures that   . Different parametric functions have been proposed in the literature to ensure this condition (Fayer & Simmons, 1995;Rossi & Nimmo, 1994;Webb, 2000). In this study, we use the concept suggested by Iden and Durner (2014) which is an extension to the model introduced by Peters (2013). This so-called PDI ("Peters-Durner-Iden") model describes soil water retention as a superposition of water storage in completely-filled capillaries and non-capillary water. The "residual" water content is kept formally as a model parameter, but now describes the maximum volumetric content of non-capillary ("adsorbed") water. It still represents the water content at which liquid continuity in completely-filled pores breaks and the capillary component of the total hydraulic conductivity becomes zero. This leads to a physically consistent description of the WRC with the same number of model parameters as the original model of the WRC (Iden & Durner, 2014;Rudiyanto et al., 2015).
The WRC of M3 is given by the equation where S c (h) [-] is the capillary saturation function, and S nc (h) [-] is the non-capillary saturation function that describes the sorption of water. The capillary saturation function is obtained by rescaling Equation 5 In the PDI model, different base functions describing U(h) can be used (Peters, 2013). In the following, we apply the van Genuchten (1980) model and refer to the resulting parametrization as "VG-PDI." Rescaling the original van Genuchten function by Equation 12 ensures that the water content becomes zero at 0 h . The function describing storage of adsorbed water is (Iden & Durner, 2014)  is a smoothing parameter which depends on the pore-size distribution (Iden & Durner, 2014) and is calculated with the empirical function Regarding hydraulic conductivity, model M3 accounts for water flow in completely filled capillaries only and neglects film and corner flow. Therefore, the HCC is calculated with Equation 7. The closed-form expression for the hydraulic conductivity in completely filled capillaries is (Peters, 2014)  Note that the original VG-PDI model accounts for film and corner flow which is neglected in M3.
Including these flow, components increase hydraulic conductivity in the mid-moisture range, where most capillaries are no longer fully water-filled and vapor transport is still of limited importance (Peters, 2013;Peters & Durner, 2008b). To simplify terminology, we will refer to the non-capillary component of hydraulic conductivity as a "film" component in the remainder of this article and denote it by K lf [m s −1 ]. This reflects that the PDI model uses the mechanistic film-flow model of Tokunaga (2009) which is based on the dependence of film thickness on tension and a solution of the Navier-Stokes equation. Peters (2013) adopted this model and transferred it into a simplified formulation which ensures the correct slope of K lf in a log-log plot versus tension. In the VG-PDI model, the HCC including film flow is described by the equation Peters' (2013) formulation of Tokunaga's (2009) where K sf [m s −1 ] is the saturated hydraulic conductivity for film flow, , and a = −1.5 [-] is a parameter that defines the slope of K lf (h) in the log-log plot versus tension. Since the hydraulic conductivity is calculated with Equation 17 does not necessarily become zero at h 0 , we rescaled K lf (h) in the same ways as U(h) in Equation 12. In comparison to M3, M4 needs two additional model parameters to describe the film flow component in the HCC, namely the parameters K sf and a.

Comparison of Models
The four model concepts differ systematically in their approach to treating water storage and hydraulic conductivity in the medium to dry moisture range, which is crucial for drying processes in soils. Figure 2 illustrates the changes in the shape of the WRC (left) and the HCC (right). The model parameters used to calculate the SHP are given in Table 3. The parameters were obtained by fitting M4 to the data set "sandy loam" by Pachepsky (1984) which was also analyzed by Peters (2013) and Iden and Durner (2014). For the WRC, the difference between models M1/M2 and models M3/M4 lies in the medium to very dry moisture range, where the PDI model approaches the zero water content gradually and therefore reflects a finite water capacity (a partial derivative of water content w.r.t. tension) until complete dryness. In contrast, the water capacity becomes virtually zero in the traditional parametrization using the residual water content, because there are virtually no pores with a diameter corresponding to tensions greater than 10 5 cm. This is of concern for the numerical stability of algorithms that are used to simulate the drying of soils because a flux of water out of the soil by evaporation leads to an almost arbitrary pressure head change if water capacity is close to zero. In practice, this problem is often solved by limiting the pressure head to some minimum value at the top node of the numerical mesh  but even with this workaround, numerical instabilities have been reported (Li et al., 2019).
With respect to the HCC, model M1 leads to the well-known linear decrease of   10 log ( ) K h with log 10 (h) for pressure heads smaller than the air-entry. Essentially, once capillaries are drained, the liquid hydraulic conductivity in this model concept becomes negligible. The addition of the isothermal vapor diffusion component to hydraulic conductivity (M2) leads to a sharp transition point, from which liquid hydraulic conductivity becomes much smaller than the isothermal vapor conductivity. For well-sorted sands, where most of the pore space drains in quite a narrow pressure head range, this transition point is often located in the tension range of 10 2 cm, that is, relatively close to the air-entry pressure. The HCC of M2 and M3 are almost identical, which illustrates that the change of the WRC by the PDI model, specifically the scaling by Equation 12, hardly alters the HCC if film flow is neglected. However, the consideration of film flow in M4 changes the HCC markedly in the mid moisture range where most capillaries have already been emptied, but isothermal vapor flow is still not playing a dominant role.

Forward Simulation of Laboratory Evaporation Experiments
We simulated bare soil evaporation experiments under laboratory conditions using the software code Hydrus-1D for solving the coupled Equations 1 and 2. The design of the numerical study reflects the experiments presented in the companion paper. Figure 3 provides an overview of the setup. The height of the soil was 10 cm and the domain was discretized into 200 equally long finite elements. Observation nodes were defined at depths 1, 2, 4, and 8 cm. Evaporation was simulated for 16 days. A hydrostatic pressure head distribution with a tension of 0 cm at the bottom was used as the initial condition for the water flow simulation. At the top of the profile, an atmospheric boundary condition was used. The evaporation rate E [m s −1 ] was calculated internally as (Sakai et al., 2011) IDEN ET AL.

10.1029/2020WR028513
where ρ vsoil and ρ vair [kg m −3 ] are the water vapor contents of the soil surface and atmosphere, respectively, ρ w [kg m −3 ] is the density of liquid water, and r a [s m −1 ] is the aerodynamic resistance. In a field situation, r a can be related to wind speed and surface roughness under certain conditions (Bonan, 2015). In this study, r a is an empirical parameter that controls the evaporation rate. We set it to 50 s m −1 to obtain evaporation rates during stage-one (the constant-rate phase) which were similar to the ones observed in the experiments presented in the companion paper. We assumed laboratory conditions, with relative humidity in the atmosphere equal to 50% and a constant temperature of 20°C. A no-flux condition was used as a lower boundary condition.
For the heat flow simulation, a constant soil temperature of 20°C was used as initial condition. . The left plots show the soil water retention curves, the right plots depict the hydraulic conductivity curves. K is total hydraulic conductivity, K lc is hydraulic conductivity due to capillary flow, K vh is hydraulic conductivity due to isothermal vapor flow, and K lf is hydraulic conductivity due to film and corner flow. accounts for the sign convention of the soil model, that is, a positive value of G indicates an upward flux of energy. We assumed that the evaporation experiment is conducted in a dark laboratory and therefore shortwave radiation was neglected in the computation of R n .
To illustrate the influence of film flow on the water and energy fluxes during laboratory evaporation, we performed three forward simulations with the coupled model. The first simulation used M3 for the SHP and therefore ignored film-flow ("M3D", "D" for direct). In the second simulation, film-flow was included and the SHP was parametrized with M4 ("M4D"). Finally, the film flow conductivity K lf was increased by a factor of three ("M4D+") to further illustrate the influence of film flow. The parameters for the SHP are given in Table 3. We kept the WRC identical for the three simulations to highlight the influence of K(h) on the evaporation dynamics. The SHP were passed as external lookup tables to Hydrus-1D. Since the isothermal vapor conductivity K vh is calculated internally, it was not considered in the HCC used as input in Hydrus-1D to avoid double-counting. The forward simulations with the coupled model included the process of vapor flow enhancement in the function K vT . Enhanced vapor diffusion was first discussed by Philip and De Vries (1957) and later analyzed experimentally and theoretically (Cass et al., 1984;Webb & Ho, 1999). However, the existence and relevance of this process are still under debate (Novak, 2016;Shahraeeni & Or, 2012;Shokri et al., 2009). Like previous simulation studies (Saito et al., 2006;Sakai et al., 2011), we used the model of Cass et al. (1984) to parametrize the enhancement factor as a function of water content. The temperature dependence of SHP is included in Hydrus-1D by the temperature effects on surface tension and viscosity (Simunek et al., 2013). The related equations are given in Table 1 Note. M3 is the VG-PDI without film-flow, M4 is the VG-PDI with film-flow, and M4+ is the VG-PDI with an increased amount of film-flow.

Table 3
Parameters of the VG-PDI Model to Parametrize the Soil Hydraulic Properties Figure 3. Setup of the evaporation experiments analyzed in this and the companion article. The five pressure heads h i (t) are measured by tensiometers and the two sensors on the left measure the relative humidity in the soil. The soil column is surrounded by a box containing Styrofoam chips for thermal insulation from the surroundings. The entire setup including data loggers is mounted on a balance to measure the effective water content and the evaporation rate.

Inverse Modeling With the Richards Equation Including Isothermal Vapor Flow
The key question of this study is whether evaporation experiments which are non-isothermal under most circumstances and therefore contain thermal fluxes and temperature effects on SHP, can indeed be evaluated by inverse modeling with the isothermal Richards equation and whether this yields unbiased estimates of the SHP. The simulated pressure head data, which contain the effects of thermal water and vapor fluxes and temperature on the SHP were used in the objective function to identify the SHP by inverse modeling. To convert the simulated data into synthetic measurements, they were disturbed with independent, heteroskedastic Gaussian error (for details, see next paragraph).
We used the Hydrus-1D code for numerically solving the Richards equation for the inverse simulations. The initial and lower boundary conditions for the water flow simulation were identical to those of the forward simulations. At the upper boundary, we specified the actual evaporation rate from the forward simulations as a flux boundary condition ("atmospheric boundary condition" in Hydrus-1D). The maximum tension at the top was set to 2⋅10 6 cm. A constant soil temperature (in space and time) of 20 °C was assumed in the inverse simulations. This temperature is identical to the ambient temperature assumed in the forward simulations. For real evaporation experiments, this approach is straightforward to apply and removes the necessity to measure soil temperature.
The data used in the objective function were the time series of pressure heads at depths 1, 2, 4, and 8 cm and the column-averaged water contents. The latter contains information on the deviation of the simulated cumulative evaporation from the "observed" one, and therefore adds further information on the soil hydraulic properties. Note that the simulated cumulative evaporation can differ from the observed one if the upper boundary condition is switched from a flux to a Dirichlet condition. This is the case if the critical tension at the top of the profile is reached Vanderborght et al., 2017). This methodology is identical to the one used in the study of Weber et al. (2017). The weighted-least-squares objective function was where  p is the vector of unknown parameters (M1 to M4), y i is "observed data" (simulated with the coupled model and disturbed with noise), y(ti,  p) is model-simulated data (Richards equation), and σ i are the standard deviations of the measurement error, which serve as regression weights. Equation 20 assumes independent, normally distributed measurement errors with an expected value of zero. Artificial noise with the same statistical properties was added to the simulated data before they were used in Equation 20. We assumed a standard deviation of 10 −5 for the column-averaged water content. For the tension, we assumed two instruments for the measurement and assigned σ i according to their precision. The first instrument is a mini-tensiometer with boiling delay (Schindler et al., 2010) for which we assumed a measurement range 0 cm < h < 2000 cm and a heteroscedastic error standard deviation    0.2 0.01 tens h [cm]. The second instrument is a relative humidity sensor which measures H r in the soil with a standard deviation   0.0075 hum . As outlined by Goss and Madliger (2007), such sensors can be used to determine the tension by use of the Kelvin equation. The error in h resulting from the error in H r was calculated using Gauβ' law for the propagation of random error and is shown in Figure 4 (black line) together with the Kelvin equation (blue line). Because of the limited accuracy of the pressure head measurement for relative humidity data of h < 10 5 cm, we discarded the corresponding data points from the analysis to adequately reflect the information content of a real measurement campaign. The threshold corresponds to a relative humidity of 0.93 at 20°C and a relative error in the pressure head (coefficient of variation) of approximately 10% (Figure 4).

IDEN ET AL.
10.1029/2020WR028513 Figure 4. Relationship between relative humidity and pressure head at 20°C (Kelvin equation, blue line) and relative error (coefficient of variation) of the pressure head calculated from relative humidity measurement with a standard deviation 0.0075 for H r (black line). It is assumed that the water potential equals the tension (absolute value of pressure head).
The sensors and precisions used in this study are representative of the real experiments presented in the companion paper.
The objective function (Equation 20) was minimized with the shuffled-complex-evolution algorithm (SCE-UA) of Duan et al. (1992). The SCE-UA was run with 8 complexes and 15 points per complex and each optimization was performed five times to ensure the robustness of the results. The run yielding the minimum value of the objective function was used for the analysis. As mentioned above, Moldrup's model was used to describe tortuosity for water vapor flow in the forward simulations. In reality, there is always uncertainty about the exact shape of the gas tortuosity function and different parametrizations for the dependency on water content exist (e.g., Moldrup et al., 2005). To reflect the limited knowledge on the effective gas diffusion coefficient, we included an additional parameter β [-] in the inverse simulations which multiplies the isothermal vapor conductivity K vh defined by Equation 8. This parameter was constrained to the interval 0.1 ≤ β ≤ 10 and estimated jointly with the other parameters describing the soil hydraulic properties. Estimating more parameters is often criticized because it contradicts the principle of parsimony. In this study and the one presented in the companion paper, the additional estimation of β poses an extra challenge to the inverse identification because of a potential compensation between increased vapor diffusion and film flow. Note that the primary objective of this work is not only the closest possible agreement between observations and simulations as in model calibration but the test of whether an unbiased identification of system properties is possible. Therefore, increasing the number of degrees of freedom increases the complexity of the inverse problems and strengthens the analysis.
For model M1, the six parameters θ r , θ s , α, n, K sc , and τ were estimated. For M2, we additionally estimated the vapor flow parameter β. The same seven parameters were estimated for M3. Finally, eight parameters were estimated for M4 because it additionally contains the saturated hydraulic conductivity for film flow K sf . The film flow parameter a was kept at a constant value of −1.5 in M4 as suggested by Peters (2013) based on Tokunaga (2009).
The four models M1-M4 were compared in their ability to match the observed time series of tension in the four different depths. The tension data from relative humidity measurements at depths 1 and 2 cm were included in the objective function. Model performance was assessed by the root-mean-squared-weighted error, computed for each depth independently, and the Akaike Information Criterion, corrected for small sample size AICc as defined in Ye et al. (2008). Parameter interaction (cross-correlation) was quantified by the condition index of Belsley (1991), γ [-], which is defined as the condition number of the column-normalized weighted sensitivity matrix. The uncertainties and cross-correlations of the estimated parameters were quantified by a linear approximation of the parameter covariance matrix and the uncertainties of the SHP were computed by Monte Carlo analysis (Iden & Durner, 2007;Vrugt & Bouten, 2002). The performance of the four models M1-M4 were compared for the data sets obtained from the forward simulations M3D and M4D. A formal confirmation that vapor flow must be included in the inverse simulation would benefit from an additional scenario in which vapor flow is neglected. However, as it is impossible to exclude vapor flow in dry soil during evaporation in reality, we have decided to not include such a scenario in this article.

Forward Simulation of Evaporation Experiments Using Coupled Model
The results of the simulations with the coupled model are shown in Figures 5, 6, and 7. Figure 5 shows the SHP (a, b), the evaporation rate (c), cumulative evaporation Q eva [mm] (d), the tension (log-scale) at the soil surface (e), and the surface temperature (f) for the three simulations M3D, M4D, and M4D+. The results illustrate the effect of the HCC on the water dynamics during an evaporation experiment. Figure 5b shows that liquid hydraulic conductivity K lh differs between the three models in the dry range, with higher conductivity values for the cases where film-flow was included (M4D, M4D+). The WRC and the HCC due to vapor flow K vh are identical for the three models.
For all models, the evaporation rate decreases rapidly at the beginning (Figure 5c), then becomes approximately constant (constant-rate phase, stage-1), and finally decreases during stage-2 (falling-rate phase). The initial decrease of the evaporation rate is caused by a marked decrease in soil surface temperature (Figure 5f), which leads to a decrease in saturation vapor pressure. The temperature decrease is caused by the flux of latent heat at the soil surface which sets in immediately after the start of the evaporation experiment. This illustrates how the process of evaporation couples the fluxes of water and energy. Stage-1 ends when the relative humidity at the soil surface becomes smaller than unity. This corresponds to a tension of approximately 10 4 cm for which the relative humidity is 99% at a temperature of 20 °C according to Equation 11 and Figure 4. At the transition from stage-1 to stage-2, the tension at the soil surface increases quickly and then converges towards a value of 10 6 cm which corresponds to a relative humidity of 50% at 20 °C. This shows that the soil is approaching equilibrium with the surrounding atmosphere. During stage-2, the decrease of the evaporation rate leads to an increase in surface temperature (Figure 5f) because the latent heat loss becomes progressively smaller.
The change of the HCC caused by the inclusion of film-flow leads to an extension of the duration of stage-1 (Figure 5c) for simulations M4D and M4D+. Moreover, the evaporation rate drops less quickly during stage-2 (Figure 5c). In agreement with the differences in the evaporation rates, tension and temperature at the soil surface rise less quickly for the simulations including film-flow (Figures 5e and 5f). The cumulative amount of water evaporating from the soil is increased by film flow (Figure 5d). As a result, the soil column is driest after 16 days for M4D+, followed by M4D and M3d. After 16 days, tension and temperature at the soil surface approach the ambient conditions (20°C,  6 10 cm h ). As the WRC and soil thermal properties are identical for the three simulations, the reason for these effects is the increase in K(h) by film-flow which IDEN ET AL.  . The isothermal vapor conductivity K vh is identical for the three cases. Simulation M3D does not include film flow, simulation M4D includes film flow, and simulation M4D+ includes more film flow than M4D to illustrate the effect of changes in soil hydraulic conductivity on evaporation rate (c), cumulative evaporation (d), surface pressure head (e), and soil surface temperature (f). leads to a higher upward water flux. Overall, the results illustrate that the shape of the HCC affects the duration of stage-1 evaporation and the evaporation rate during stage-2 even if all other physical properties are the same.
In a typical evaporation experiment conducted to determine SHP, the time series of tension at various depths is monitored by tensiometers. Figures 6a and 6b shows the time series of tension simulated with the coupled model for M3D and M4D. The simulation M4D+ is not shown to keep the presentation concise. Accounting for the non-capillary component in the HCC influences the evolution of tension at the five different depths shown. As previously discussed, the tension at the soil surface (z = 0) rises faster for M3D. However, the change in tension at the surface cannot progress as quickly into the soil because of the smaller hydraulic conductivity compared to M4D. The marked differences in the time series of tension at the different depths are a necessary condition to distinguish between different models of the SHP by inverse modeling. This will be analyzed in detail in Section 3.2. The depth profiles of tension are shown in Figure 6c and 6d. During stage-1 (duration ca. 4 days), the vertical distribution of tension is approximately linear (Peters & Durner, 2008a). During stage-2, the vertical distribution of tension becomes nonlinear in both simulations. However, the profiles are very different for M3D and M4D. In simulation M3D, a dry-surface layer develops in which the water flux is dominated by IDEN ET AL. the diffusion of water vapor (Peters, 2013;Shokri et al., 2009). The lower position of this layer is marked by a strong decrease in tension. In contrast, simulation M4D results in smoother depth distribution of tension. The smaller gradients in tension in M4D are a consequence of the higher hydraulic conductivity which overcompensates the higher water fluxes. The water content profiles are shown in Figures 6e and 6f) correspond to the tension profiles closely. The depth distribution is linear at the beginning of the experiment and later becomes nonlinear. A dry surface layer with an almost constant water content with depth becomes evident for M3D but less so for M4D. This illustrates that the formation of a dry surface layer and the sharpness of the transition towards the moist subsurface layer depend on the HCC because the WRCs are identical. The water content profile of M4D shows that the soil becomes drier in M4D as compared to M3D which is in agreement with the cumulative evaporation shown in Figure   We investigated the relative contribution of liquid water and water vapor fluxes, both isothermal and thermal, to the total flux (Figure 7; positive fluxes upwards). During stage-1, liquid fluxes dominate and vapor fluxes are negligible. Therefore, data for stage-1 are not shown. The first profiles correspond to t = 5 days. As the soil surface dries out during the onset of stage-2, vapor fluxes become relevant and as the experiment continues, this "dry front" moves downwards (Figures 7a-7d). However, there is a clear difference between M3D and M4D. In the simulation without film-flow, the liquid flux becomes negligible in the top few centimeters and the flux is dominated by vapor diffusion. In contrast to this, the increase in the HCC (M4D) leads to substantial fluxes of liquid water (15% of the total flux) in the topsoil and although vapor flow is the dominant transport mechanism, it is not the only water flux occurring in the top centimeters of the profile. In M4D, there is always some water that is moving towards the surface in liquid form.
During stage-2 the evaporation rate decreases and there are two driving forces for vapor fluxes. The most important one is the gradient in tension (isothermal vapor flow) which drives water upwards as evaporation continues. However, when the evaporation rate becomes smaller, the soil becomes warmer at the surface as a result of net radiation and the sensible heat flux. This causes a temperature gradient which drives a downward movement of vapor by thermal vapor diffusion (Milly, 1984). The magnitude of this counterflux ranges from 0, at the onset of stage-2, to 20% for M3D, and from 0% to 7% for M4D. This shows that the thermal vapor flux is relatively small in evaporation experiments in which no additional radiation is applied at the top to enhance the evaporation rate. Note that neglect of vapor enhancement (η = 0) will lead to much smaller thermal vapor fluxes and the isothermal fluxes will be even more dominant. Since the existence and magnitude of vapor flow enhancement are still debated and since the parametrization of the enhancement factor by the model of Cass et al. (1984) leads to higher values than those recently reported by Shahraeeni and Or (2012), the analysis presented in this study is conservative in the sense that thermal vapor fluxes are likely overestimated.
The results presented in this section are, of course, conditional on the model assumptions and have to be supported by experimentation. Nevertheless, the theoretical investigation illustrates the processes occurring during an evaporation experiment and allows a quantitative estimation of the relative contribution of water fluxes in liquid and gaseous forms. The coupled model summarizes current knowledge on coupled water, vapor, and heat flow and is a standard that dates back to the theoretical work by Philip and De Vries (1957). Of course, applying the model in one spatial dimension must assume that horizontal flows are negligible and that no energy is exchanged laterally between the soil and the surroundings. We note that the experimental data presented in the companion article, that is, evaporation rates, pressure heads, and soil temperatures, follow the same temporal pattern as the simulated ones in this article. This illustrates that the coupled model is capable of describing the dynamics of water and energy in laboratory evaporation experiments. Figure 8 shows the results of fitting the Richards equation to the data set obtained from simulation M3D in which film flow was not included in the HCC. Table 4 shows the goodness-of-fit measures (GOF), the estimated value of β, and the collinearity index of Belsley γ for the four models. Figure 8a-8d shows the time series of pressure head (simulation M3D, with added noise) and the best-fit obtained by the inversion. The pressure head data have a gap between  2000cm h (upper limit of tensiometer) and  5 10 cm h (lower limit of hygrometer). In the following, we will denote the simulations with the Richards equation by the name of the model which is used to parametrize the SHP (M1-M4).

Data Set Without Film Flow (M3D)
Model M1, the de-facto standard to characterize SHP, was not able to describe the tension data over the full range from  0 h to  6 10 cm h (Figure 8a). The tensiometer data (h < 2,000 cm) are described relatively well, but M1 completely fails to describe the data in the dry range (h > 10 5 cm). This is a result of the steep drop in the HCC caused by the neglect of K vh (Figure 8f). Indeed, the amount of water evaporated in M1 was smaller than the one simulated with the coupled model as indicated by the RMSWE for the column-averaged water content (Table 4). Adding K vh to the HCC improves the fit substantially ( Figure 8b) and M2 describes both types of observations very well. This is supported by the GOF in Table 4. The quality of the fit for M3 and M4 is similar to M2 in Figure 8c-8d, but the values in Table 4 reveal some further improvement by ensuring a water content of zero at oven-dryness by these models. Most importantly, the AICc corroborates the selection of M3 as the best model. The objective function values are identical for M3 and M4, but as M3 has only 7 fitting parameters, the minimum value of the AICc is obtained for M3, the model used for generating the data. This shows that for the data set without film flow in the HCC, the measurements are adequately described by a model in which film flow is not included. The RMSWE are all very close to unity for M3 and M4. The fact that the RMSE at z = 1 cm is the highest indicates that the Richards equation cannot fully describe the time series of tension at this depth. One possible cause for this mismatch is thermal vapor fluxes in the topsoil which are not included in the RE. However, the overall match is still excellent and the correct model of the SHP (M3) is identified. The SHP shown in Figures 8g and 8h

Data Set Including Film Flow (M4D)
The results of the analysis using simulation M4D (including film-flow) are shown in Figure 9. The classic model M1 failed again to describe the data as discussed above. M2 and M3 performed better than M1, but with clear discrepancies between "observed" and simulated data. In fact, M2, M3, and M4 are all able to describe the tensiometer data properly, but only M4 leads to an accurate description of the data in the hygrometer range. This highlights the importance of additional measures of state variables in medium to dry soil and supports the experimental design used in the companion article. As M4 is the only model which is able to describe the full data set, it also has, by far, the best values of the goodness-of-fit criteria (Table 5). Furthermore, M4 is the only model which identifies the correct WRC and HCC (Figures 9e-9h). The values of the collinearity index of Belsley (1991) in Table 5 are similar for M2, M3, and M4. Indeed, the model with the smallest number of parameters (M1) has the highest degree of parameter cross-correlation (smallest value of γ). This shows again that the number of estimated parameters is not a proxy for parameter correlation. Finally, the values of the scaling factor for isothermal vapor flow β are 10 for the models which do not include film-flow in the HCC (M2, M3) but close to unity for M4. In the inverse simulations with M2 and M3, an increased vapor flow is identified to compensate for the absence of film-flow in the HCC. For M4, the value of β is close to unity because film flow is included which eliminates the necessity to increase the vapor flux.

Conclusions
Our conclusions are based on two key results of this work: (1) the sensitivity of the evaporation process on the shape of the SHP in the forward modeling with a fully coupled water-energy-transport model and (2) the identification of SHP by inverse modeling with the Richards Equation as simplified process model from observations of state variables in the soil which are accessible by advanced instrumentation. Our findings with respect to the forward modeling are: 1. Forward simulation with a comprehensive coupled water-energy model shows that the time series of the evaporation rate and the spatiotemporal dynamics of tension and water content are very sensitive to the shape of the HCC in the medium to dry moisture range. 2. The increase in hydraulic conductivity in this moisture range caused by film and corner flow extends the duration of stage-1 evaporation and cause higher evaporation rates in stage-2.
IDEN ET AL. Note. The number of estimated parameters is np, WSSE is the minimal objective function value, AICc is the Akaike criterion adjusted for small sample size, RMSWE is the root-mean-squared weighted error, β is the effective vapor flow parameter, and γ is the collinearity index of Belsley (1991).    Table 5 As Table 3 but for the Inverse Simulation in Which the Data Were Generated With Simulation M4D (With Film-Flow) 3. Altogether, film and corner flow increase the loss of water from the soil by evaporation leading to potential effects on the water and energy balance in arid and semiarid regions. This affects their agricultural potential, the need for irrigation or advanced irrigation techniques, the danger of salinization, and the quantity of groundwater recharge.
Our conclusions with respect to the inverse modeling are: 1. Inverse modeling with the Richards equation showed that the identified WRC and HCC are unbiased although the Richards equation is an isothermal flow model which does neither encompass thermal fluxes of liquid water and water vapor nor temperature effects on SHP. 2. The inclusion of isothermal vapor diffusion in the HCC is an absolute necessity for describing the time series of pressure head and the evaporation rate during evaporation from the soil. 3. In media in which film-flow plays a role, the film-flow component in the conductivity curves is needed to describe the time series of pressure heads. Conversely, in a medium with negligible film flow, the time series can be described without incorporating it in the conductivity curve. This is important because it proves that the identification of an increased hydraulic conductivity in medium to dry soil by inverse modeling, as observed by investigating reals soils in the companion article, is not an artefact caused by using a simplified process model, the Richards equation, for the inversion of experimental data.
We conclude that the SHP can be identified without bias from evaporation experiments using tensiometers and relative humidity data. These findings, which are based on numerical forward modeling, need to be tested with real experiments. In our companion article, we report experimental data of evaporation experiments using advanced instrumentation and the results of the data evaluation by inverse modeling with the Richards equation.