A new top boundary condition for modeling surface diffusive exchange of a generic volatile tracer: theoretical analysis and application to soil evaporation

We describe a new top boundary condition (TBC) for representing the air–soil diffusive exchange of a generic volatile tracer. This new TBC (1) accounts for the multiphase flow of a generic tracer; (2) accounts for effects of soil temperature, pH, solubility, sorption, and desorption processes; (3) enables a smooth transition between wet and dry soil conditions; (4) is compatible with the conductance formulation for modeling air–water volatile tracer exchange; and (5) is applicable to site, regional, and global land models. Based on the new TBC, we developed new formulations for bare-soil resistance and corresponding soil evaporation efficiency. The new soil resistance is predicted as the reciprocal of the harmonic sum of two resistances: (1) gaseous and aqueous molecular diffusion and (2) liquid mass flow resulting from the hydraulic pressure gradient between the soil surface and center of the topsoil control volume. We compared the predicted soil evaporation efficiency with those from several field and laboratory soil evaporation measurements and found good agreement with the typically observed two-stage soil evaporation curves. Comparison with the soil evaporation efficiency equation of Lee and Pielke (1992; hereafter LP92) indicates that their equation can overestimate soil evaporation when the atmospheric resistance is low and underestimate soil evaporation when the soil is dry. Using a synthetic inversion experiment, we demonstrated that using inverted soil resistance data from field measurements to derive empirical soil resistance formulations resulted in large uncertainty because (1) the inverted soil resistance data are always severely impacted by measurement error and (2) the derived empirical equation is very sensitive to the number of data points and the assumed functional form of the resistance. We expect the application of our new TBC in land models will provide a consistent representation for the diffusive tracer exchange at the soil–air interface.

Although the soil diffusive tracer transport problem, particularly for water, has been studied for decades (e.g., Rose, 1963;Hanks et al., 1967;Prat, 2002;Moldrup et al., 2003;Goss and Madliger, 2007), methods on how to compute air-soil diffusive tracer fluxes are still incomplete.Taking soil water evaporation for example, both theoretical and experimental studies indicate there are often two evaporation stages (e.g., Lemon, 1956;Philip, 1957;Gardner and Fireman, 1958;Ritchie, 1972;Pan and Mahrt, 1987;Shaw, 1987;Yamanaka et al., 1997;Saravanapavan and Salvucci, 2000;Suleiman and Ritchie, 2003;Shokri and Or, 2011).During stage-I, water evaporates at the potential evaporation rate, which is determined by atmospheric demand.This stage is maintained mostly by capillary liquid flow, while water vapor diffusion is less significant (Shorki et al., 2009;Shahraeeni and Or, 2010).During stage-II, water evaporates at a lower rate that is determined by how fast water can be supplied from the water source below the soil surface.The contribution of water vapor diffusion in this stage is more significant (Yamanaka et al., 1997;Shokri et al., 2009).A few studies further separate stage-II into two sub-stages (e.g., Lemon, 1956): the first indicates the onset of hydraulic continuity disruption that reduces capillary flow (e.g., Lehmann et al., 2008;Shokri et al., 2008) and the second indicates the dominance of film flow (e.g., Tokunaga, 2009;Shahraeeni andOr, 2010, 2012).
Many experimental studies have provided insights into the soil evaporation problem.Lehmann et al. (2008) showed that stage-I evaporation would last until the first drying front recedes to a depth such that the downward gravity and viscous forces overcome the upward capillary driving force, disrupting hydraulic continuity with the soil surface.Shahraeeni et al. (2012) demonstrated that for a very strong atmospheric demand (characterized by a high wind speed in their wind tunnel experiments), the stage-I evaporation lasted only a short period, or did not exist at all.A few laboratory columnscale experiments (e.g., Yamanaka et al., 1997;Shokri et al., 2009) showed that the stage-II evaporation could be computed accurately as the amount of water vapor diffused from an evaporating front (aka bottom of the dry surface layer (Budyko, 1948) or the second drying front (Shokri et al., 2009)) to the top of the dry surface layer (DSL), which has led to the physics-based DSL approach for computing the stage-II evaporation (see Eq. 7 in Yamanaka et al., 1997, or Eq. 3 in Shokri et al., 2009).However, Yamanaka et al. (1997) pointed out that, under varying incident radiation, equating the surface evaporation with water vapor that diffused from the evaporating front produced significant errors.They asserted that, in their study, such a discrepancy resulted from the water phase change within the topsoil (1 cm) during the course of radiation change.Goss and Madliger (2007) in their field studies also indicated that, with the diurnally cyclic change of environmental variables, the topsoil water was pe-riodically recharged from atmospheric dew formation and deep soil water flux through both capillary and film flow.Because climate and weather forecasting models operate across a wide spectrum of environmental conditions, these findings indicate that a more robust soil evaporation approach should be identified to account for those complicating factors, such as condensation and evaporation within the topsoil.
Other studies have attempted to empirically relate soil evaporation with the topsoil water content (e.g., Sun, 1982;Camillo and Gurney, 1986;Kondo et al., 1990;Lee and Pielke, 1992;Sellers et al., 1992;van de Griend and Owe, 1994;Komatsu, 2003;).They parameterize the soil evaporation as where E (kg water m −2 s −1 ) is the actual evaporation, E p (kg water m −2 s −1 ) is the potential evaporation, β (dimensionless) is the evaporation efficiency defined as the ratio of E with respect to E p , and r a (s m −1 ) is the atmospheric resistance, which is a serial sum of aerodynamic resistance (e.g., Zeng et al., 1998), laminar boundary resistance (Leighly, 1937;Simpson et al., 2012), and possibly litter layer resistance (Schaap and Bouten, 1997;Sakaguchi and Zeng, 2009).We left out the interfacial resistance in this study, as it is only important when the atmospheric tracer concentration is extremely low (Heideger and Boudart, 1962).The soil resistance is denoted by r s (s m −1 ).Equation (1) has led to two different groups of empirical approaches to compute soil evaporation.The first approach empirically parameterizes the soil evaporation efficiency β as a function of soil moisture (Lee and Pielke, 1992 (hereafter LP92); Komatsu, 2003;Merlin et al., 2011), which is achieved by curve fitting with respect to the measurementderived ratio between the actual and potential soil evaporation.The LP92 approach has been applied in many numerical models (e.g., Walko et al., 2000;Jia et al., 2001;Xiu and Pleim, 2001;Lawrence et al., 2011) and has been criticized for a lack of dependence on atmospheric conditions (Komatsu, 2003;Merlin et al., 2011).The second approach empirically parameterizes the soil resistance r s as a function of topsoil moisture, where the topsoil thickness varies from 1 cm to 5 cm (Kondo et al., 1990;Sellers et al., 1992;van de Griend and Owe, 1994).Parameterization of r s is achieved with three steps: (1) use measured soil evaporation and atmospheric variables (e.g., temperature, wind, and humidity profiles) to compute atmospheric resistance, potential evaporation, and soil evaporation efficiency, (2) invert for the soil resistance with the relationship and (3) perform parametric curve fitting against the inverted soil resistance data to obtain the empirical soil resistance  (e.g., Sun, 1982;Camillo and Gurney, 1986;Kondo et al., 1990;Sellers et al., 1992;van de Griend and Owe, 1994).These soil resistance formulations have been extrapolated in many studies to different soils to estimate soil evaporation (e.g., Daamen and Simmonds, 1996;Qiu et al., 1999;Saito et al., 2006;Bittelli et al., 2008;Oleson et al., 2008;Stockli et al., 2008;Smits et al., 2011), although these equations were derived for specific types of soils under specific soil physical and atmospheric conditions and have not been demonstrated to be broadly applicable for all conditions (Bittelli et al., 2008).In addition, the limitations of these empirical equations, including that of the empirical soil evaporation efficiency equation, have not been thoroughly explored.
The lack of consensus in parameterizing soil evaporation has led to ambiguity in formulating consistent top boundary conditions (TBCs) for numerical models that try to model a broad range of volatile tracers that exchange between the soil and atmosphere.This ambiguity has led to many different formulations being applied, such as for soil evaporation (Deardorff, 1978;Camillo and Gurney, 1986;Kondo et al., 1990;Sellers et al., 1992;LP92;Sakaguchi and Zeng, 2009), for water isotope modeling (e.g., Mathieu and Bariac, 1996), and methane modeling (e.g., Walter and Heimann, 2000;Zhuang et al., 2004;Tang et al., 2010;Riley et al., 2011).The use of different TBCs has been shown to substantially impact surface isotope flux estimates (Riley et al., 2002).
Ideally, the problems of specifying a proper TBC to parameterize air-soil tracer exchange could be overcome by modeling the mass transfer in a coupled free (air, Navier-Stokes flow) and porous-medium (soil, Darcy flow) flow system (e.g., Jamet et al., 2009;Shavit, 2009;Mosthaf et al., 2011).However, land models attempt to resolve tracer fluxes at a complex surface that includes plants, litter, and other structures, making such a strategy impractical at the current time.We thus argue that it would be valuable to find a TBC that can consistently describe the air-soil exchange of a generic volatile tracer for numerical models applied at site, regional, and global scales.To be optimally useful, such a TBC should meet the following eight criteria: (1) it should only use parameters associated with soil hydraulic properties (e.g., porosity, conductivity, soil water retention curve) and tracer properties (e.g., diffusivity in different phases, solubility, hydrolysis rate); (2) it can easily incorporate capillary flow and film flow (e.g., Brunauer et al., 1938;Rose, 1963;Tuller and Or, 2001;Goss and Madliger, 2007;Tokunaga, 2009;Lebeau andKonrad, 2010, 2012); (3) it can represent evaporation and condensation (or dissolution for gases other than water vapor) within the soil (e.g., Hanks et al., 1967;Goss and Madliger, 2007); (4) it is able to account for the effects of tracer solubility in water, which is a function of temperature, pH, sorption, desorption, and chemical reactions (e.g., Fang and Moncrieff, 1999;Nassar and Hor- ton, 1999); (5) it can explain the two-stage soil water evaporation behavior that is well documented in various studies (Philip, 1957;Hanks et al., 1967;Salvucci, 1997;Shokri et al., 2009;Shahraeeni et al., 2012); ( 6) it can reveal caveats of empirical soil evaporation formulations; (7) it can be applied across the full range of soil moisture (e.g., Silva and Grifoll, 2007); and (8) it is compatible with the TBC that is used to compute air-water (such as air-sea, air-lake) tracer exchange (e.g., Liss, 1973;Liss and Slater, 1974;Johnson, 2010).
To develop a TBC that meets all these eight criteria, we contend several processes need to be represented (Fig. 1): (1) phase partitioning and its dependence on soil moisture content; (2) representation of the moisture content versus matric potential relationship across a wide range of soil moisture conditions; (3) fluxes associated with the aqueous phase that are not subsumed in the gas phase diffusive calculation in the topsoil control volume (TSCV); and (4) resistances across the laminar sublayer and into the canopy airspace or, in the absence of vegetation, into the atmosphere directly.
The remainder of this paper is structured as follows.Section 2 describes the new TBC conceptual model, relevant mathematics, new formulations for soil resistance and soil evaporation efficiency, applications to soil water evaporation, and evaluation of the caveats of existing empirical soil resistance and soil evaporation efficiency formulations.Section 3 presents the new TBC evaluation results and Sect. 4 presents the caveats and potential applications and improvements of the proposed new TBC.Finally, Sect. 5 concludes the paper with a summary of main findings.
In this section, we first present the conceptual model and theoretical arguments to specify a consistent TBC for modeling volatile tracer exchange between the soil surface and atmosphere (Sect.2.1).We then apply the methodology to soil evaporation to obtain a new formulation of soil resistance and bulk conductance (Sect.2.2).Section 2.3 describes the methods to evaluate the new TBC with existing laboratory experiments and their empirical equations.After that we describe methods to assess the caveats of existing empirical formulations of soil resistance and soil evaporation efficiency (Sect.2.4), as well as the new TBC (Sect.2.5).

A top boundary condition for generic diffusive tracer exchange between the soil surface and atmosphere
A conceptual description for the exchange of a volatile tracer between the soil surface and atmosphere is presented in Fig. 1.Here and below, we refer to a "volatile tracer" as any substance (or tracer) with an equilibrium vapor pressure that is established rapidly compared to the model time step.We conceptualize the air-soil tracer exchange problem using two assumptions: (1) no tracer accumulates at the air-soil interface, such that the tracer flux from the TSCV to the interface is balanced by the tracer flux from the interface to the atmospheric reference level z a (m); and (2) the equilibrium between gaseous and aqueous phases of the tracer is instantaneous both at the interface and inside the soil column.Assumption (2) implies that the partitioning of a tracer into its different forms (e.g., gaseous, aqueous, sorbed phases) can be modeled by the law of mass action (e.g., Yeh and Tripathi, 1991), and therefore can be computed from the bulk tracer concentration with relevant parameters of solubility, pH dependence, sorption, and desorption (see Tang et al. (2013) for an example of NH 3 and representing its different forms in soil).Some studies have questioned the equilibrium assumption, however no consensus appears to have been achieved (e.g., Novak, 2012;Smits et al., 2012).Nevertheless, these two assumptions lead to where the gaseous flux from the air-soil interface to the atmosphere F a (mol m −2 s −1 ) is With the finite difference approximation, the gaseous flux F g (mol m −2 s −1 ) from the center of the TSCV to the soil-air interface is and the aqueous flux F w (mol m −2 s −1 ) from the center of the TSCV to the soil-air interface is Here D g (m 2 s −1 ) is the effective gaseous phase diffusivity computed as a function of temperature, moisture, and soil properties of the topsoil (see Eqs. A11-A13 for water vapor); z 1 (m) is the thickness of the TSCV; C g,s (mol m −3 air) is the gaseous concentrations of the tracer at the interface; D w (m 2 s −1 ) is the bulk aqueous phase diffusivity, which includes contributions from liquid mass flow (in the form of both capillary flow and film flow) and aqueous molecular diffusivity (Eqs.A1-A2 in Appendix A); and C w,s (mol tracer m −3 water) and C w,1 (mol tracer m −3 water) are, respectively, the aqueous concentrations of the tracer at the soil-atmosphere interface (surface) and at the center of the TSCV.The air filled porosity ε 1 (m 3 air m −3 soil), volumetric soil moisture θ 1 (m 3 water m −3 soil), and tracer concentrations C w,1 and C g,1 are all defined at the center of the TSCV.
Unless stated otherwise, we use the subscript "a" and "1" to indicate variables defined for the atmosphere and the TSCV, respectively.
A few more assumptions are made in obtaining Eqs. ( 4)-(6).First, there exists a linear change in tracer concentration from the soil surface to the atmospheric reference height, and, second, from the center of the TSCV to the soil surface.The reasonableness of the first of these assumptions has been discussed comprehensively in the literature (e.g., Leighly, 1937;Brutsaert, 1965;Deardorff, 1978).The second assumption raises the question of over what z 1 range the tracer concentration gradient can be approximated as linear, which we will revisit in Sect.4.2.2.Third, the law of mass action can be reasonably applied to adjust the aqueous and gaseous phases (and possibly sorbed phase, which can be lumped into the aqueous phase (Tang et al., 2013)) of the tracer, such that possible condensation and evaporation within the TSCV are considered simultaneously.Fourth, the relationship between the air-soil tracer fluxes and tracer concentrations in deeper soil layers can be approximated by the soil matric pressure gradient at the center of the TSCV, though a dynamical relationship is developed by continuously updating the tracer concentration profile in the soil column in a numerical model.In this last assumption, we rely on the parameterization of soil hydraulic properties (Appendix B) to account for the capillary and film flow.
Invoking the equilibrium assumption C w = BC g at both the air-soil interface and center of the TSCV (where B (dimensionless) is the mean TCSV Bunsen solubility coefficient, computed as a function of temperature, pH, sorption, and desorption reactions in the TSCV (Tang et al., 2013)) and Eqs. ( 4)-( 6), one obtains where defines the ratio between atmospheric resistance r a and soil resistance r s , and D 1 (m 2 s −1 ) is the bulk tracer diffusivity defined in the TSCV: From Eq. ( 7), the bulk resistance r T (m s −1 ) for the diffusive tracer exchange between the atmospheric reference level and center of the TSCV is A few comments are relevant to Eq. ( 10).First, the atmospheric resistance r a above the soil surface and the soil resistance r s (calculated as z 1 / [2D 1 (Bθ 1 + ε 1 )]) below the soil surface together determine the bulk resistance r T .Second, when r a is much greater than r s , the bulk resistance r T is effectively determined by r a , indicating that gas exchange is limited by the transport between the soil surface and the atmosphere.Dew formation (Jacobs et al., 2000), stage-I soil evaporation, and water evaporation from standing water are three examples of such conditions.For evaporation from a saturated soil or standing water, r T = r a , and we have the scenario for potential evaporation.Third, when r a is much smaller than r s , the gas exchange is limited by conditions in the TSCV, i.e., transport from the center of the TSCV to the surface, subject to the constraints of soil hydraulic properties and chemical properties of the tracer.Two examples for such cases are evaporation from a very dry soil (Katata et al., 2007) and evaporation subject to a very high atmospheric demand (Shahraeeni et al., 2012).Fourth, when the TSCV is saturated or ponded with standing water, D 1 = D w and Eq. ( 10) is consistent with Eq. ( 6) in Liss (1973), which describes the bulk conductance for air-sea tracer exchange.
For air-sea (Liss, 1973) and air-standing water exchange, D w includes contributions from eddy mixing, molecular diffusion, and possibly dispersion, while the capillary driven mixing is effectively zero.Fifth, for gases of low solubility in water, Eq. ( 10) indicates that the soil resistance, and consequently the bulk resistance, will rapidly become independent of the aqueous diffusivity as the soil dries.This strong dependency on soil moisture explains why field studies that ignored aqueous diffusivity can still accurately compute the air-soil fluxes for some tracers using measured soil-gas concentration profiles (e.g., Davidson et al., 2006).However, Eq. ( 10) should improve air-soil tracer flux estimates under a range of environmental conditions, and consequently enable a better representation of correlations among surface fluxes for different tracers.Finally, when z 1 is set to zero, r T = r a , which states that the bulk conductance is just the atmospheric conductance.
In the following sections, we apply the methods described above to the problem of water vapor exchange between bare soil and the atmosphere.Application of these methods to other generic tracers will be described in future work.

A new formulation of bare soil evaporation
From Eq. ( 7), the soil evaporation E can be computed as where M w (kg mol −1 ) is the molar mass of water, ρ a (kg m −3 ) is the air density at z a , q a (g water vapor g −1 air) is the atmospheric specific humidity at z a , α is the relative humidity in the TSCV (aka the "alpha factor" that can be computed using the Kelvin equation (e.g., Kondo et al., 1990; LP92)), and q * g,1 (g water vapor g −1 air) is the saturated specific humidity at the center of the TSCV.In obtaining Eq. ( 11), we have used the relationship C a = ρ a q a /M w and C g,1 = ρ a q g,1 /M w .We have also made the assumption that the difference in air density between the atmosphere and the TSCV is small, as has been done in many previous studies (e.g., Kondo et al., 1990;LP92;Oleson et al., 2010).Equation ( 10) is used to calculate the bulk resistance for soil evaporation, with the mean Bunsen solubility coefficient for water vapor in the TSCV defined as where ρ l (kg m −3 ) is the liquid water density and the water vapor density ρ g,1 (kg m −3 ) is calculated as ρ g,1 = ρ a q g,1 .
The new soil resistance term (in Eq. 10) can be further partitioned into the vapor diffusion resistance r g (s m −1 ) and the volatilization resistance r w (s m −1 ), which represents the adjustment to disequilibrium caused by vapor diffusion or liquid water flow, through In general, the aqueous diffusivity D w includes contributions from molecular, Darcy's diffusion, and possibly dispersion (Eq.A1).For water molecules, except its isotopologues, the effective molecular diffusivity D w,m (the subscript "m" indicates the molecular diffusivity of the aqueous tracer) is nil, thus, from Eq. ( A2), one has where we have ignored the gravity and non-isothermal terms (for reasons that will be discussed in Sect.4).K 1 (m s −1 ) and ψ 1 (m) are the TSCV hydraulic conductivity and soil matric potential, respectively.As we will show below, including the liquid mass flow (represented by K 1 ∂ψ 1 ∂θ 1 ) and water vapor adsorption (characterized by B and ψ 1 ) in the model leads to reasonable predictions of the two-stage behavior of soil evaporation.Field data, curve fitting, Narita sand Kondo and Saigusa (1994) Laboratory data, curve fitting, loam Sakaguchi and Zeng (2009) * * * (SZ09) Curve fitting, multiple soils 2+3/b * The soil type is inferred by comparing the Clapp and Hornberger (1978) parameters with those presented in the dataset by Cosby et al. (1984).S 1 is the wetness (WFPS) of the topsoil.* * D 0 is the water vapor diffusivity in the atmosphere (Eq.A11).* * * We corrected an error in the gas diffusivity equation (Moldrup et al., 2004) as it was quoted in Sakaguchi and Zeng (2009).In the relevant calculations, we used w= 5.

Evaluating the new soil resistance formulation
We used two approaches to evaluate the new soil resistance formulation Eq. ( 13) (and consequently the new TBC described in Sect.2.1).First, we evaluated whether using Eq. ( 13) with relevant information about the experiments and documented soil hydraulic property parameterizations can reproduce the behavior of the empirical soil resistance equations that were derived from the inverted soil resistance data.Four different empirical soil resistance formulations derived from four experiments that provided sufficient information for a comparison were used in the evaluation (Table 1).In order to do a comprehensive analysis, we also included the soil resistance formulation proposed by Sakaguchi and Zeng (2009), which was obtained from curve fitting with multiple datasets.Second, we compared the soil evaporation efficiency computed by using different empirical soil resistances and their corresponding theoretical soil resistances from Eq. ( 13).The soil evaporation efficiency is computed as In addition, we compared our predictions with the soil evaporation efficiency formulation proposed by LP92, which was based on a synthesis of multiple datasets: where θ fc (m 3 m −3 ) is the field capacity and θ 1 is the soil moisture content in the TSCV.θ fc is determined by assuming a hydraulic conductivity of 0.1 mm day −1 at the soil water content of θ fc .Equation ( 16) is used in the Community Land Model version 4 (CLM4) to parameterize soil evaporation (Oleson et al., 2010;Lawrence et al., 2011).We compared the difference in the calculated soil evaporation efficiency using Eq. ( 16) and using our new formulation for four very different but typical soils that are parameterized with the Clapp and Hornberg (1978) approach adopted by CLM4 (Table 3), though the results are general for all models that use Eq. ( 16) for parameterizing soil evaporation.We also characterized the gain in prediction accuracy of our new approach compared to the LP92 approach.
Third, we evaluated the importance of phase partitioning in modeling soil evaporation, an analysis made possible by the explicit representation of aqueous and gaseous phases in our new TBC.Combining Eqs. ( 5), (6), and (13), it can be shown that the fraction of surface soil evaporation associated with direct liquid water evaporation (f w ) is and the fraction from water vapor transport (f g ) is We analyzed the partitioning of evaporation into direct liquid water evaporation and water vapor transport for the four different types of soils (Table 3).With these evaluations, we tested whether Eqs. ( 10), ( 13), and our new TBC formulation (Eqs.7-9) are an improvement over existing empirical formulations.
Table 2. Clapp and Hornberger (1978) parameters for the topsoil used in Table 1.Kondo and Saigusa (1994)    * Since no soil hydraulic property data were given in their study, the soil texture of the fine sandy loam soil is assumed as 69 % sand, 11 % clay, and 20 % silt.The algorithm used in CLM4 was then used to obtain the soil hydraulic parameters.

Caveats of the empirically derived soil resistance equations
We analyzed the caveats of the empirically derived formulations by mimicking the methods used to obtain them in experiments (Kondo et al., 1990;Yamanaka and Yonetani, 1999) through a synthetic inversion experiment.Specifically, we assumed that using Eq. ( 13) one can compute the theoretically true soil resistance at all soil moisture states provided the soil temperature and soil hydraulic properties are well defined.Next, for simplicity in the analysis, we assumed that the time series of true atmospheric resistance is constant with a value of 50 s m −1 .Using the theoretically true soil resistance and the assumed constant atmospheric resistance, we computed the theoretical soil evaporation efficiency through Eq. ( 15).Then, we imposed a 5 % random noise to both the true constant atmospheric resistance and the true soil evaporation efficiency to mimic a relatively low measurement error.5 % is an arbitrarily chosen value but significant enough to show the effect of measurement error on deriving the empirical soil resistance formulation.For comparison, the error in current eddy flux latent heat measurements is around 5-20 % (Foken, 2008), indicating our analysis is conservative.We then used Eq. ( 2) to invert for the synthetic soil resistance data from the error-convolved soil evaporation efficiency and atmospheric resistance data.Empirical soil resistance equations were then derived through curve fitting with respect to the inverted soil resistance data.We performed this synthetic analysis using the soil hydraulic properties for the loamy soil studied by Kondo and Saigusa (1994).The curve fitting was carried out with the Markov chain Monte Carlo approach (e.g., Green, 1995) by minimizing a cost function that measures the distance between inverted soil resistance data and the calculated soil resistance from using the assumed empirical parameterization.Three sets of curve fitting approaches were designed (Table 4): the CFIT1 used the functional form suggested by Kondo and Saigusa (1994), and CFIT2 and CFIT3 used the functional forms suggested by Sellers et al. (1992), but with different numbers of inverted soil resistance data points.With these synthetic experiments, we evaluated two hypotheses: (1) the inverted soil resistance data is highly affected by measurement error and ( 2) the different functional forms and different number of data points used in the curve fitting produce very different soil resistance formulations.

Caveats of the proposed new top boundary condition
We specifically investigated four issues that would affect predictions with our new TBC formulation.We first investigated the uncertainty of using our new TBC with different soil hydraulic property parameterizations.Both the CH (Clapp and Hornberger, 1978) and VG (van Genuchten, 1980) parameterizations were extended to cover the whole range of soil moisture, from fully saturated to completely dry, using the approach by Silva and Grifoll (2007) (Appendix B; named as CH-SG and VG-SG, respectively).The hydraulic conductivities for the two sets of analyses were parameterized with the default CH and VG approach.The film flow is partially considered in the extended soil water retention curve but not in the hydraulic conductivity curve, which can be done with more complicated parameterizations (Lebeau andKonrad, 2010, 2012), yet we expect the conclusion of our analyses will not change due to this omission.The comparison was done for the silty loam soil (Table 5) documented in Shao and Irannejad (1999), where they provided the standard CH and VG parameterizations that minimized the differences in their modeling experiments of soil mass and heat transfer.Second, we assessed quantitatively whether using a TSCV of 1 cm and 5 cm with our new TBC will produce significant uncertainty for modeling soil evaporation.Results were evaluated by comparing this uncertainty to that from the using different soil water retention curve parameterizations.
Table 5. Silty loam soil hydraulic parameters for different models used for Fig. 7. Parameters are adapted from Shao and Irannejad (1999).
Scale Shape For the third and fourth issue, we analyzed under what conditions ignoring the gravity effect and non-isothermal effect on soil matric pressure can significantly impact the prediction from using our new TBC.We consider these two effects as corrections to the Darcy diffusivity (involved in Eq. 13 for the computation of soil resistance) computed from using the isothermal hydraulic head curve (Eq.A2; e.g., Bear, 1972): where the gravity effect correction is and the non-isothermal effect correction is For Eq. ( 20), we interpret ∂z/∂θ 1 as the change in gravitational hydraulic head per unit change of soil moisture.Since ∂z/∂θ 1 is usually negative, the gravity correction D ψ,grv acts effectively to reduce the liquid mass flow for soil surface evaporation.We again evaluated the impacts for two TSCV thicknesses, 1 cm and 5 cm, to assess if the magnitude of D ψ,grv is small enough to be ignored compared to the isothermal Darcy diffusivity D ψ .For context, most of the climate models involved in the CMIP5 experiment (Taylor et al., 2012 and references therein) used a TSCV of less than 5 cm, while a few used 10 cm.
For the non-isothermal correction, we represent ∂ψ 1 /∂T 1 using the Young-Laplacian equation (e.g., Grant, 2003): where β 0 is a constant dependent on the porous media structure (Grant, 2003).Bachmann et al. (2002) found, for six different soils, that β 0 has a mean value of −457 K, which implies ∂ ln ψ 1 ∂T 1 is generally negative, with a magnitude on the order of a few percent.We also evaluated the non-isothermal effect for two different TSCV thicknesses: 1 cm and 5 cm.During evaporation, the temperature gradient within the TSCV could be either positive or negative, depending on the energy exchange between the atmosphere and soil (see Fig. 1 in Goss and Madliger (2007) and Hanks et al. (1967) for observational examples).Based on our modeling experiences and also on laboratory experiments (Hanks et al., 1967), we imposed a temperature gradient (absolute magnitude) of 400 K m −1 for a 1 cm thick TSCV, and a gradient of 200 K m −1 for a 5 cm thick TSCV, both of which are likely much greater than what could be found under most environmental conditions.We then assessed if the non-isothermal corrections are small enough to be ignored when compared to the isothermal Darcy diffusivity.

Comparison with existing empirical parameterizations
Below we present results from evaluating our new TBC with existing empirical parameterizations of soil evaporation.We describe the results for soil resistance and soil evaporation efficiency.

Soil resistance
Among the four different types of soils represented in the experimental studies, comparison with the only study using  13) for the four soils listed in Table 2.The CH 3 parameterization is used for predictions with Eq. ( 13).For all cases, the vapor diffusivity 4 in air was assumed to be D 0 = 2.4 !10 "5 m 2 s -1 .As specified by the investigators of the 5 original publications, the TSCV thicknesses are 1, 5, 2, and 2 cm for panels (a), (b), (c), 6 and (d), respectively.Note that the soil resistances are plotted on a log10 scale.13) for the four soils listed in Table 2.The CH parameterization is used for predictions with Eq. ( 13).For all cases, the vapor diffusivity in air was assumed to be D 0 = 2.4 × 10 −5 m 2 s −1 .As specified by the investigators of the original publications, the TSCV thicknesses are 1, 5, 2, and 2 cm for panels (a), (b), (c), and (d), respectively.Note that the soil resistances are plotted on a log10 scale.laboratory data (Fig. 2d) shows the best agreement (in terms of root mean square difference) between our calculated soil resistance (using Eq. 13) and the measurement-derived empirical soil resistance.Our calculated soil resistances tend to be higher than the empirically derived (labeled as "Empirical" on figures) soil resistances for the sandy loam soil (Fig. 2a) and the Narita sand soil (Fig. 2c, when the soil is dry).Our calculated soil resistances are lower than empirically derived resistances for the clay loam soil in the study by Sellers et al. (1992) (Fig. 2b).The soil resistance computed using the equation by Sakaguchi and Zeng (2009, SZ09 hereafter) behaves similarly to our new approach, except that the SZ09 calculations are less linear and decrease at a faster rate to zero for wetter soils.Our new approach calculates a very small resistance for saturated soil rather than zero resistance as assumed in SZ09.In addition, for the only empirical soil resistance derived from a laboratory experiment, where conditions could be most tightly controlled, we find that our approach results in a better agreement with the empirically derived soil resistance than with the SZ09 approach, which shows higher resistances when the soil is moderately wet.

Bare-soil evaporation efficiency
A comparison of the calculated soil evaporation efficiencies using Eq. ( 15) with the same r a but different soil resistance parameterizations (Table 1 and Eq.13) indicates that our new soil resistance formulation (Eq.13) always results in a well-behaved curve of the soil evaporation efficiency versus soil moisture (Fig. 3).Our new approach also accurately  13) and ( 16) for the four soils listed in Table 2.For 3 all cases, we assumed r a = 50 s m -1 and D 0 = 2.4 !10 "5 m 2 s -1 .As specified by the 4 investigators of the original publications, the TSCV thicknesses are 1, 5, 2, and 2 cm for 5  13) and ( 16) for the four soils listed in Table 2.For all cases, we assumed r a = 50 s m −1 and D 0 = 2.4×10 −5 m 2 s −1 .As specified by the investigators of the original publications, the TSCV thicknesses are 1, 5, 2, and 2 cm for panels (a), (b), (c), and (d), respectively.
predicted the typical two-stage behavior of soil evaporation, where in stage-I, the soil evaporates at potential evaporation, and in stage-II, the soil evaporates at lower rate than the potential evaporation.The empirical soil resistance parameterization by Sellers et al. (1992) (Fig. 3b) does not exhibit this typical two-stage soil evaporation behavior, even when the atmospheric demand is low (results not shown), which is contrary to experimental findings (e.g., Shahraeeni et al., 2012).The empirical soil resistance curve derived by Kondo and Saigusa (1994) for Narita sand soil also fails to exhibit the two-stage behavior of soil evaporation, by missing the constant evaporation stage-I (Fig. 3c).One explanation for the lack of an atmospheric controlled evaporation stage in the soil evaporation efficiency computed using the empirical parameterizations by Sellers et al. (1992) and Kondo and Saigusa (1994) could be that the data used to fit their empirical curves might have been affected by the atmospheric resistance.As we will show in Sect.4.1, the high-frequency fluctuations of atmospheric resistance and the uncertainty in measuring actual evaporation and calculating potential evaporation could significantly affect the inversion used to estimate soil resistance (Eq.2).
Compared to the LP92 and SZ09 approach, our approach captures more accurately the onset of stage-II evaporation for the better-controlled experiments (Fig. 3a, d).The soil evaporation efficiency curve computed using the SZ09 soil resistance parameterization always estimates a shorter stage-I evaporation, which would potentially cause underestimation of soil evaporation under many soil moisture conditions if it were used in land surface models.The curve computed  3. The four panels are, 3 respectively, (a) Sand, (b) Loam, (c) Sandy loam, and (d) Organic soil.For all cases, we 4 assumed that D 0 = 2.4 !10 "5 m 2 s -1 and !z 1 = 1.75 cm (as it is used in CLM4).3. The four panels are, respectively, (a) sand, (b) loam, (c) sandy loam, and (d) organic soil.For all cases, we assumed that D 0 = 2.4×10 −5 m 2 s −1 and z 1 = 1.75 cm (as it is used in CLM4).
using the LP92 equation (Eq.16) leads to higher evaporation for all four soils except when the soil moisture is about 25-50 % water-filled pore space (WFPS) (Fig. 3b).The LP92 equation also predicts lower evaporation than our new approach does when the WFPS was less than 10 %.We find, for the four typical soils parameterized with the Clapp and Hornberger (1978) approach adopted in CLM4 (Table 3), that the LP92 approach, which implicitly considers the atmospheric resistance, tends to overestimate soil evaporation compared to our new approach, especially when the atmospheric resistance is low (Fig. 4).Our approach also predicts that a very high atmospheric demand (characterized by low r a ) could result in a very short, or even no, stage-I evaporation (result not shown), in agreement with some experimental studies (Shahraeeni et al., 2012).This feature cannot be captured by the LP92 approach.Therefore, compared to our new approach, implementing the LP92 approach in land surface models will lead to an overestimation of bare soil evaporation for relatively wet soils (a problem that has already been identified in CLM4 (P.Lawrence, personal communication, 2012)), and to an underestimation of bare soil evaporation for dry soils.We address this latter point in another study.
Finally, we note that our new approach predicted the soil evaporation efficiency as a non-monotonic function of WFPS.Soil evaporation efficiency decreased slightly when WFPS increased from a completely desiccated soil, began increasing at a relatively low moisture level, and finally stabilized at a value of 1 as WFPS increased to 100 % (see Fig. 3 and also Figs. 4,6b,7b,8b).This behavior physically reflects the transition from a water vapor transport dominated (a decreasing function of WFPS) to a liquid water direct evapo- vapor transport from below the soil surface for the four soils listed in Table 3.For all 3 cases, we assumed D 0 = 2.4 !10 "5 m 2 s -1 and !z 1 = 1.75 cm (as it is used in CLM4).4 5 Fig. 5. Partitioning of evaporation into direct liquid water vaporization and water vapor transport from below the soil surface for the four soils listed in Table 3.For all cases, we assumed D 0 = 2.4 × 10 −5 m 2 s −1 and z 1 = 1.75 cm (as it is used in CLM4).ration dominated (an increasing function of WFPS) surface evaporation.The measurement based empirical curves were not able to capture such a transition, because, for a very dry soil under natural conditions, (1) dew formation rather than evaporation would be the dominating water flux at the atmosphere and soil interface and ( 2) such a transition (from a decreasing β into a increasing β) is too weak to be accurately measured.

Partitioning of soil evaporation into direct surface liquid water vaporization and below surface water vapor transport
For the four different types of soils analyzed here (Table 3), we find that the direct evaporation from liquid water is about equal to the below surface water vapor transport under moderately dry conditions, depending on the physical properties (e.g., distribution of particle size and particle surface area, which are reflected in the hydraulic parameters used in the model) of the soil (Fig. 5).We also find the differences in partitioning the evaporation into f w and f g for mineral soils is determined mostly by the water sorption capability of the soil.The Clapp and Hornberger (1978) parameterization characterizes the sorption capability to monotonically increase with the shape parameter b.Consequently, for the sandy clay soil that has the greatest water sorption capability, the 50 % evaporation partition point (defined as the WFPS where f w = f g ) occurred at the highest WFPS, followed by the loamy soil (Fig. 5b), and then the sandy soil (Fig. 5c).The partitioning curves for organic soil are quite different from the mineral soils at the same level of volumetric soil water content due to the organic soil's greater porosity (results not shown), though they resembled the sand soil in the scale of relative saturation in our analysis (Fig. 5d).For the four soils, the water vapor transport from below the soil surface dominates the overall evaporation when the WFPS is lower than 25 %.The importance of water vapor flux in the overall soil evaporation for dry soils indicates that it is critical to consider the existence of water vapor to model soil moisture in dry areas, as well as in cold areas, when the liquid water content is low during the frozen period.

Discussion
In this section, we first discuss the uncertainties associated with the new TBC, specifically for soil resistance and phase partitioning in Sect.3. Then we analyze the caveats of existing empirical soil resistance parameterizations, followed by the caveats in formulating our new TBC.Finally, we discuss the limits of using our new soil resistance formulation, measurements that could be used for further evaluation and improvement, and potential extensions to enable a new formulation of surface evaporation and generic tracer exchanges that considers both vegetated and non-vegetated soil surfaces.

Soil resistance
We used four experimental studies to evaluate the soil resistance computed using the empirical approach and our new TBC approach.Our new approach showed best agreement with the empirical approach from the laboratory pan evaporation experiment by Kondo and Saigusa (1994).However, many aspects of those pan soil evaporation experiments have been criticized.One of the criticisms, based on lab experiments (e.g., Lehmann et al., 2008;Shokri, 2009), is that the soil in pan studies is often too shallow to properly account for the capillary liquid flow recharge from deep soil.For example, Kondo and Saigusa (1994) used a depth of 10 cm for their laboratory pan soil evaporation experiment, while the laboratory column experiments, which are typically done with engineered soil particles (e.g., glass beads or sand), have depths ranging from about 30 to 100 cm (Lehmann et al., 2008;Shokri et al., 2009;Smits et al., 2011).In contrast, the field studies we considered involved more complicated soil structure (Campbell and Shiozawa, 1992) and used the whole soil profile, though they provided no information for the depth from soil surface to bedrock, nor for moisture status below the 10 cm soil depth.We were unable to find laboratory experiments with a deeper soil column with sufficient data (such as soil hydraulic parameters, time series of topsoil moisture content and evaporation flux, and relevant atmospheric parameters) for a comprehensive comparison.However, we point out that our approach implicitly assumes that the soil column is sufficiently deep, and that our soil hydraulic property parameterization accurately rep-resents the soil flow regime as the topsoil moisture content evolves.
In addition, laboratory experiments usually use relatively constant atmospheric conditions, whereas field studies always involve the multi-scale variability of all atmospheric variables, including temperature, humidity, wind, turbulence, pressure, and radiation.The strong contrast between environmental conditions in field experiments and those in laboratory experiments makes it difficult to extrapolate laboratory results for land model applications that need to account for large environmental variability.For instance, Yamanak et al. (1998) pointed out significant discrepancies between the measured soil evaporation and that predicted using the DSL approach when radiation was varied in their experiments (see their Fig.11b).Goss and Madliger (2007) also pointed out that their three-month field measurements of the evaporation front never moved deeper than 3 cm into the soil, because the soil water dynamics were completely different from laboratory experiments due to the regular diurnal cycle of atmospheric dynamics, whereas laboratory experiment often observe the evaporation front continuously moving deeper as evaporation continues (Shokri et al., 2009).
Thus, by considering the significant uncertainty in field experiments and the difficulty of extrapolating results from the well-controlled laboratory measurements, we contend that the relatively good agreement with the empirical soil resistances derived from better controlled field chamber experiments (e.g., Fig. 2a) indicates our new top boundary condition can capture the soil resistance well within the range of uncertainty.

Phase partitioning of the soil evaporation
The phase partitioning discussed in Sect.3.2 qualitatively agrees with findings in other studies (e.g., Katata et al., 2007), however we note that our approach simplifies the physics by constraining the surface evaporation to only come from the TSCV.Therefore, no information could be provided on what fraction of the below surface water vapor transport in the TSCV is from deeper soil or from phase change with respect to the liquid water to maintain the pressure equilibrium described by the Kelvin equation.Integrating our new TBC with a more complete multi-phase numerical modeling of soil water and heat transport (Katata et al., 2007;Novak, 2010;Smits et al., 2011) will provide some insights on these two issues.Further, considering the typical strong isotopic gradient in the first few centimeters of an evaporating soil (Miller et al., 1999;Riley et al., 2002), explicit measuring and modeling of isotopic soil water (H 2 18 O and DHO) profiles and fluxes may also help to disentangle the dynamics of f w and f g (e.g., Barnes and Allison, 1988;Yamanaka and Yonetani, 1999).  2 (also used by Kondo and Saigusa (1994)) is 5 used for this analysis.For all cases, we assumed the true r a = 50 s m -1 , D 0 = 2.4 !10 "5 m 2 6 s -1 , and !z 1 = 1.75 cm.There are 51 inverted data points in the moisture range [0.  2 (also used by Kondo and Saigusa, 1994) is used for this analysis.For all cases, we assumed the true r a = 50 s m −1 , D 0 = 2.4 × 10 −5 m 2 s −1 , and z 1 = 1.75 cm.There are 51 inverted data points in the moisture range [0.07, 0.53], of which 35 inverted data points are in the moisture range [0.07, 0.27].

Caveats of the empirical soil resistance formulations
In Sect.2.4, we discussed four types of errors that can affect an empirically derived formulation of r s .The first two are related to the observations that will be used to invert the soil resistance data using Eq. ( 2): (1) the error in estimating the atmospheric resistance r a and ( 2) the error in estimating the soil evaporation efficiency.In experiments, the soil evaporation efficiency is computed as a ratio of the measured actual evaporation and the potential evaporation.The atmospheric resistance is a function of air temperature, wind speed, pan size (for pan evaporation experiment), and water vapor diffusivity in the atmosphere (Kondo et al., 1992).The second two types of errors are associated with the curve fitting process used to derive the empirical equation, which depend on (3) the number of data points being used and (4) the parametric functional forms being used.
We found the inverted soil resistance data (labeled as "Inverted Data") with the first two types of error deviated significantly from the theoretically true soil resistance curve (labeled as "New" in Fig. 6a), indicating that the derivation of soil resistance data from measurements is substantially affected by measurement error.In addition, we found that the inverted soil resistance data at medium to high water content resembled the shape for the Narita sand (Fig. 2c), which implies the soil resistance data inferred by Kondo and Saigusa (1994) could have been impacted by measurement error.Therefore, this finding indicates that the often-adopted assumption that the soil resistance can be reliably inverted using Eq. ( 2) with relevant field measurements from micrometeorological instruments (e.g., Kondo et al., 1990;Yamanak and Yonetani, 1999) probably does not hold under most conditions and that chamber measurements (e.g., van de Griend and Owe, 1994) might be favored due to their tighter con-trol of atmospheric conditions.We thus expect that these two types of error substantially impacted the derivation of the empirical soil resistance equations.
For the empirical equations derived through curve fitting, we found that using the functional form suggested by Kondo and Saigusa (1994) (CFIT1 in Fig. 6a and Table 4) resulted in a better fit to the inferred soil resistance data than using the functional form suggested by Sellers et al. (1992) (CFIT2 in Fig. 6a and Table 4).However, the CFIT1 result is misleading because the systematic drift in the inverted soil resistance data (compared to the theoretical truth labeled as "New") is simply a result from the convolution of the true soil resistance with the first two types of error.The lack of good quality data for dry soil (which is typical in field measurements due to the relatively lower probability of dry soil) also leads to an underestimation in the soil resistance (compared to the theoretically true soil resistance) at low soil moisture by using the functional form suggested by Kondo and Saigusa (1994).CFIT3, which only uses observations in the moderately wet soil moisture range [0.07, 0.27] (the total number of data points used in CFIT3 (confined between points P 1 and P 2 ) is around 70 % of those used in CFIT1 and CFIT2; see Fig. 6a and Table 4), agrees better with the theoretically true soil resistance than that calculated by CFIT2, while they both use the functional form suggested by Sellers et al. (1992).In addition, the resultant shape of the empirical resistance curve from CFIT3 is very similar to the empirical curve in Fig. 2a, which was derived from data that were only collected when the soil had medium and low water content (van de Griend and Owe, 1994).When the three empirical soil resistance curves are used to compute the soil evaporation efficiency (Fig. 6b), we find that CFIT1 leads to a curve similar to the empirical curve in Fig. 3c, CFIT2 leads to a curve similar to that in Fig. 3b, and CFIT3 leads to a curve similar to that in Fig. 3a.These comparisons indicate that the empirical soil resistance equations derived from using field data by Sellers et al. (1992), Kondo and Saigusa (1994), and van de Griend and Owe (1994) are all likely substantially impacted by the four types of errors.Therefore, we suggest that, because of the measurement uncertainty and the uncertainty associated with data coverage across the soil moisture range, one cannot determine which of the two approaches (our new TBC versus those from curve fitting to the observations) is more likely correct.However, the success in explaining the caveats of empirical soil resistance curves indicates that our approach is practically accurate, and more carefully designed experiments are needed to evaluate if our new formulation provides a more accurate parameterization of bare soil evaporation than do those empirical approaches that have been widely used to interpret various observations and perform numerical simulations.5.For all cases, these curves are calculated by assuming r a = 50 s m -1 , 5 D 0 = 2.4 !10 "5 m 2 s -1 , and !z 1 = 1.75 cm.6 7 8 9 10 Fig. 7.The effect of using different parameterizations of soil hydraulic properties (extended Clapp and Hornberger (CH-SG) and extended van Genuchten (VG-SG)) on the computed soil resistance and the corresponding soil evaporation efficiency for the soils listed in Table 5.For all cases, these curves are calculated by assuming r a = 50 s m −1 , D 0 = 2.4 × 10 −5 m 2 s −1 , and z 1 = 1.75 cm.

Effects of different soil hydraulic property parameterizations on soil resistance and soil evaporation efficiency
We found (Fig. 7a) that using the two different soil hydraulic property parameterizations resulted in generally similar soil resistance relationships with soil moisture.Compared to the CH-SG scheme, the VG-SG scheme calculates smaller soil resistance when the soil has 20-60 % WFPS and higher resistance when the soil has higher than 60 % WFPS.Accordingly, the VG-SG scheme results in higher soil evaporation efficiency when the soil has about 20-60 % WFPS (Fig. 7b).
The differences in the calculated soil resistances and soil evaporation efficiencies resulted from the differences in these two parameterization schemes, and we note that by tuning the parameters (e.g., conductivity, shape parameter) for either of the two schemes, one can obtain a good match between the predictions from these two schemes.Thus, given the uncertainties resulting from the different soil hydraulic property parameterizations (or parameterization equifinality (e.g., Beven, 2006;Tang and Zhuang, 2008)), our predicted soil resistance and soil evaporation efficiency values (Sects.5 is used in the calculation.For all cases here, we assumed r a = 50 s m −1 and D 0 = 2.4 × 10 −5 m 2 s −1 .

Choosing the TSCV thickness
We found that using a TSCV of 1 cm and 5 cm resulted in different soil resistance and soil evaporation efficiency dependencies on WFPS (Fig. 8).However, the general variations of soil resistance and evaporation efficiency as a function of topsoil 1 cm or 5 cm mean moisture is maintained.The difference in the predicted onset of stage-II evaporation is, in general, small considering the uncertainty in measuring evaporation during experiments.In addition, the difference in predicted stage-II evaporation efficiency and soil resistance is small compared to the uncertainty derived from using different soil hydraulic property parameterizations (Fig. 7).However, we suggest increasing the TSCV thickness to greater than 5 cm should be done with caution because both measurements (e.g., Hanks et al., 1967;Goss and Madliger, 2007) and high-resolution modeling studies (Grifoll et al., 2005;Novak, 2010) have indicated that the constant soil moisture gradient approximation was acceptable only within the first few centimeters of the soil.Therefore, for numerical applications of the new TBC: (1) we cannot distinguish results from using 1 and 5 cm TSCVs, because the observations are often too poor to resolve such differences; (2) if any differences arise, they can be corrected through model calibration; and (3) we aution against using a TSCV thicker than 5 cm.  3. The circles in panel (a) and (b) indicate the threshold soil moisture where the soil evaporation efficiency equals to 0.90.The ratios are represented in absolute magnitude.The wiggles and discontinuities in the figure resulted from using the Silva and Grifoll (2007) approach to obtain the full range soil water retention curves (see the text for details).

Effect of gravity and non-isothermal soil matric potential on soil resistance
We found the impact due to the gravity correction term D ψ,grv was small even up to a 5 cm thick TSCV for mineral soils (Fig. 9a and b).The organic soil is impacted more by the gravity correction, but this impact can also be considered small (at most 10 % for 5 cm thick TSCV) because it exists mostly during the stage-I evaporation, when the soil evaporates at an almost constant rate.As the stage-II evaporation develops, the capillary pressure gradient rapidly dominates the gravity force and the water vapor contributes more to the evaporating flux, such that D ψ,grv is negligible again.Therefore, D ψ,grv can be safely ignored in computing the Darcy diffusivity and consequently the soil resistance.However, when applying the new TBC in a numerical model, we recommend the time step should be no greater than 1 h, though the exact threshold may depend on the specific application.If one is to calculate soil evaporation at time steps longer than a few hours or even a few days, the gravity term might be important (with its cumulative effect on infiltration) depending on the soil physical properties (e.g., see Salvucci and Entekhabi, 1994;Salvucci, 1997), such that a different mathematical formulation should be considered.
The non-isothermal correction D ψ,T also contributes marginally to D ψ (less than 10 % even for a very strong temperature gradient), and is less important when the TSCV is thicker (Fig. 9c, d), which is in agreement with experimental findings (Hanks et al., 1967).We note that the wiggles and discontinuities at low moisture content in Fig. 9 are a result from using the Silva and Grifoll (2007) approach to obtain a full range extension of the soil water retention curve with the Clapp and Hornberger (1978) parameterization (Appendix B).Those wiggles did not significantly affect the calculation of soil resistance and can be removed through a modified approach to obtain the full range extension of the soil water retention curve.

Usage, limits, and possible extensions of the new top boundary condition
Aside from our goal to develop a TBC to consistently represent the diffusive exchange of a generic volatile tracer at the air-soil interface, an example application of our new approach could be to combine with the force-restore approach (Jacobs et al., 2000) and use the resultant equations to predict the temporal variation of surface soil moisture and evaporation using measured wind speed, soil temperature, and air humidity in semiarid regions.The necessary equations for such applications can be obtained by replacing the soil evaporation formulation in Jacobs et al. (2000) with our new formulation (Eq.9-11).Application of our new formulation to estimate soil evaporation at large spatial scales and long time steps could be enhanced by considering several other factors.First, our development here has not considered the effect of roots on modifying soil evaporation.Theoretically, it is possible to extend our results to include plant roots by homogenizing the horizontal heterogeneity of the root network and assuming the effect of roots on the transport of soil water and relevant tracer can be represented by the capillary bundle model (Frensch and Steudle, 1989).We will leave such an extension to future studies.
Second, because of the short time step used in numerical models (e.g., CLM4, with which we have implemented the new TBC, uses a 30-min time step), infiltration, surface runoff, and evaporation can be computed sequentially.However, application of the new TBC could be problematic for applications (such as estimating soil evaporation at daily time steps) that require longer time steps, when the cumulative effect of gravity in controlling water supply for soil evaporation may be important and D ψ,grv should be included in Eq. ( 14) or a different set of equations should be used (Salvucci, 1997).
Third, our development does not consider the preferential flow in macro-pores, where the flow is dominated by gravity.Such effects can be readily described by adopting the dual permeability and dual porosity model, at the cost of adding additional parameterization uncertainties (Gerke, 2006).Another limit associated with ignoring the gravity effect in the new TBC is that the approach cannot be used to relate the soil evaporation with moisture measured at depths greater than Hydrol.Earth Syst.Sci., 17, 873-893, 2013 www.hydrol-earth-syst-sci.net/17/873/2013/ about 5 cm, for which existing empirical relationships might be more useful (Merlin et al., 2011).However, this last limitation is alleviated if the new top boundary condition is integrated in a land model that prognoses the vertical distribution of water content and fluxes (e.g., CLM4).Fourth, while our approach successfully exposed the caveats of existing empirical soil resistance formulations for soil evaporation, we call for a more comprehensive observational dataset to evaluate the robustness of our developments.Such a dataset should minimally include (1) continuous time series of soil moisture and vapor content in the topsoil, and the corresponding surface water efflux; (2) continuous time series of atmospheric resistances; (3) an accurate division of evaporation into that from direct liquid volatilization and that from water vapor transport; and (4) an accurate characterization of soil hydraulic properties.Other measurements could cover different chemical species other than water, including their surface fluxes and soil concentrations in different phases, which would help evaluate our work's applicability in a broader context.
Finally, we suggest caution when applying our approach to model volatile organic compounds for very dry soils.In our development, we assumed (1) there is always water in the soil, such that the soil moisture is well defined and all other chemicals can dissolve in the water to form solutions and ( 2) there is a linear equilibrium between the aqueous and gaseous phases, where the two are related with the Bunsen solubility coefficient.For very dry soils, the linear partitioning between aqueous and gaseous phase is still a good approximation (Ruiz et al., 1998), provided the VOC concentration is low, which is thus still within the applicability of our theory.However, for such an application, one needs to obtain the linear partitioning parameter (or the solubility coefficient as we used in our development) accurately for different soils, which is experimentally very challenging (Goss, 1993;Ruiz et al., 1998).We call for more collaboration between experimentalists and modelers to solve this problem.

Summary
We developed a new top boundary condition (TBC) to model the diffusive exchange of volatile tracers at the air-soil interface.Application of this TBC to the soil evaporation problem leads to a new formulation for soil resistance and evaporation efficiency.Comparison with measurement-derived empirical equations indicates that our new formulation is a practical candidate to formulate relevant problems in land models.This new formulation also successfully exposed caveats of existing empirical soil resistance formulations through analyses into the impacts of measurement uncertainty, number of data points, and assumed functional forms of soil resistance on the derivations of the empirical soil resistance curve as a function of water-filled pore space in the topsoil.Finally, we found our new approach leads to the largest differences in predicted bare soil evaporation compared with the LP92 approach (which is integrated in many numerical models, such as CLM4) when the soil is relatively dry or the atmospheric conductance is high.(A11)

Figure 2 .
Figure 2. Comparison of calculated soil resistances by the measurement-derived 2 empirical equations and Eq.(13) for the four soils listed in Table2.The CH 3

Fig. 2 .
Fig. 2.Comparison of calculated soil resistances by the measurement-derived empirical equations and Eq.(13) for the four soils listed in Table2.The CH parameterization is used for predictions with Eq. (13).For all cases, the vapor diffusivity in air was assumed to be D 0 = 2.4 × 10 −5 m 2 s −1 .As specified by the investigators of the original publications, the TSCV thicknesses are 1, 5, 2, and 2 cm for panels (a), (b), (c), and (d), respectively.Note that the soil resistances are plotted on a log10 scale.

Figure 3 .
Figure 3.Comparison of predicted soil evaporation efficiency ( ! ) by the measurement- 2 derived empirical equations and Eqs.(13) and (16) for the four soils listed in Table2.For 3

Figure 4 .
Figure 4. Sensitivity of the calculated soil evaporation efficiency to WFPS and the 2 atmospheric conductance for the four different soils listed in Table3.The four panels are, 3

Fig. 4 .
Fig. 4. Sensitivity of the calculated soil evaporation efficiency to WFPS and the atmospheric conductance for the four different soils listed in Table3.The four panels are, respectively, (a) sand, (b) loam, (c) sandy loam, and (d) organic soil.For all cases, we assumed that D 0 = 2.4×10 −5 m 2 s −1 and z 1 = 1.75 cm (as it is used in CLM4).

Figure 5 .
Figure 5. Partitioning of evaporation into direct liquid water vaporization and water 2

Figure 6 .
Figure 6.Sensitivity of inverted soil resistance curves to the different scenarios of 2 empirical curve fitting.The left dashed vertical line indicates the moisture level 0.07 m 3 3 water m -3 soil (P1); and the right dashed vertical line indicates the moisture level 0.27 m 3 4 water m -3 soil (P2).The loam soil in Table2(also used byKondo and Saigusa (1994)) is 5

Fig. 6 .
Fig. 6.Sensitivity of inverted soil resistance curves to the different scenarios of empirical curve fitting.The left dashed vertical line indicates the moisture level 0.07 m 3 water m −3 soil (P 1 ); and the right dashed vertical line indicates the moisture level 0.27 m 3 water m −3 soil (P 2 ).The loam soil in Table2(also used byKondo and Saigusa, 1994) is used for this analysis.For all cases, we assumed the true r a = 50 s m −1 , D 0 = 2.4 × 10 −5 m 2 s −1 , and z 1 = 1.75 cm.There are 51 inverted data points in the moisture range [0.07, 0.53], of which 35 inverted data points are in the moisture range [0.07, 0.27]. 69

Figure 7
Figure 7The effect of using different parameterizations of soil hydraulic properties 2 (extended Clapp and Hornberger (CH-SG) and extended van Genuchten (VG-SG)) on the 3 computed soil resistance and the corresponding soil evaporation efficiency for the soils 4 listed in Table5.For all cases, these curves are calculated by assuming r a = 50 s m -1 , 5

Figure 8 Fig. 8 .
Figure 8 An assessment of the TSCV thickness impacts on (a) predicted soil resistance 2 and (b) evaporation efficiency.The CH hydraulic property parameterization listed in 3Table5is used in the calculation.For all cases here, we assumed r a = 50 s m -1 and 4 D 0 = 2.4 !10 "5 m 2 s -1 .5 6 7 8 9

Figure 9 .Fig. 9 .
Figure 9. Impact of the gravity (a, b) and non-isothermal (c, d) corrections on computing 2 the effective soil diffusivity for the four soils listed in Table 3.The circles in panel (a) 3 and (b) indicate the threshold soil moisture where the soil evaporation efficiency equals 4 to 0.90.The ratios are represented in absolute magnitude.The wiggles and discontinuities 5in the figure resulted from using the Silva andGrifoll (2007) approach to obtain the full 6 range soil water retention curves (see the text for details).7 8 i = 1, • • • , 4 Coefficients of the soil water retention curve None b Clapp-Hornberger shape parameter None f w , f g Fraction of direct liquid evaporation and vapor transport in D w , D g Bulk, aqueous and gaseous tracer diffusivity m 2 s −1 D w,m , D w,ψ Molecular and Darcy diffusivity m 2 s −1 D ψ Total correction to Darcy diffusivity due to gravity and non-isothermal effects m 2 s −1 D ψ,grv Gravity correction to Darcy diffusivity m 2 s −1 D ψ,T Non-isothermal correction to Darcy diffusivity m −2 s −1 E, E p Actual and potential soil evaporation kg H 2 O m −2 s −1 F w , F g , F a Aqueous, gaseous and bulk tracer fluxes mol m −2 s θ r (m 3 m −3 ) is the residual soil water content, κ (m −1 ) is the scaling parameter, and m and n are non-dimensional shape parameters.The molecular diffusivity for water vapor in open air under standard atmospheric pressure (Montgomery, 1947) is D 0 = 2.26 × 10 −3 T 273

Table 3 .
Physical properties of the different typical soils for the comparison of soil evaporation efficiencies computed from different approaches.

Table 4 .
Comparison of different empirical soil resistance curves for the sensitivity study (Fig.6).The loam soil in Table2is used for the comparison.