Revisiting Salvucci’s Semi-analytical Solution for Bare Soil Evaporation with New Consideration of Vapour Diffusion and Film Flow

Bare soil evaporation is controlled by a combination of capillary flow, vapour diffusion and film flow. Relevant analytical solutions mostly assume horizontal flow conditions and ignore gravitational effects. Salvucci (1997) provided a rare example of a semi-analytical solution for vertical bare soil evaporation. However, they did not explicitly represent vapour diffusion and film flow, which are likely to account for a significant proportion of total flow during vertical evaporation from soils. Vapour diffusion and film flow can be incorporated via Salvucci’s desorptivity parameter, which represents the proportionality constant relating Stage 2 cumulative evaporation to the square root of time under horizontal flow conditions. The objective of this article is to implement vapour diffusion and film flow within Salvucci’s semi-analytical solution and test its performance by comparison with isothermal numerical simulation and relevant experimental data. The following important conclusions are drawn. Analytical solutions that assume horizontal flow conditions are inadequate for understanding vertical evaporation problems because they overestimate evaporation rates and mostly predict vapour diffusion and film flow to be of negligible influence. Salvucci’s semi-analytical solution is effective at predicting the order-of-magnitude reduction in evaporation caused by gravitational effects. However, it is unable to identify the correct importance of vapour diffusion and film flow because these processes can only be represented through its desorptivity parameter.


Introduction
The ability to simulate and understand bare-soil evaporation has been of widespread interest to hydrology researchers for many decades ). In addition to processes associated with Richards' equation, it is common to account for vapour diffusion, film flow and coupled heat transport . The necessary continuum scale equations 1 3 to describe this problem were originally developed by Philip and de Vries (1957). Many associated numerical schemes have since been developed (e.g. Saito et al. 2006;Ciocca et al. 2014;Fetzer et al. 2017;Li et al. 2019;Iden et al. 2021). Nevertheless, an attraction has maintained to develop analytical and semi-analytical solutions with the goal of gaining further physical insight into important processes controlling this phenomenon (Saravanapavan and Salvucci 2000;Landman et al. 2001;Teng et al. 2016;Fetzer et al. 2017;Teng et al. 2019;Mls 2022).
Analytical solutions designed to study bare-soil evaporation in this context mostly assume horizontal flow conditions (i.e., assume that gravitational effects are negligible). The advantage of the horizontal flow assumption is that the associated partial differential equation reduces to a non-linear diffusion equation, which can be further simplified using the Boltzmann transform. Indeed, under isothermal conditions with a fixed surface humidity, horizontal evaporation is described by a constant desorptivity multiplied by the squareroot of time (Gardner 1959;Lisle et al. 1987;Novak 1988;Lockington 1994;Lockington et al. 2003;Fetzer et al. 2017;Mathias and Sander 2021). However, a detailed explanation about why horizontal flow conditions might be reasonable in this context is rarely provided. Teng et al. (2019) argue that the gradient of gravitational head should be negligible compared to the strong matric suction gradient during evaporation. This may be true at the soil surface, but does not rule out the importance of gravitational forces in controlling upward water migration further down in the soil profile (Brutsaert 2014). For semi-infinite domains, numerical simulations have shown that the inclusion of gravity can reduce cumulative evaporation by several orders of magnitude (Novak 2022). Salvucci (1997) provided a rare example of a semi-analytical solution for vertical bare soil evaporation but did not explicitly account for vapour diffusion and film flow. He incorporated gravitational effects by assuming that a vertical moisture profile within the soil preserves geometric similarity and is linearly stretched as the drying front progresses downwards into the soil. His resulting equation is a function of initial hydraulic conductivity and desorptivity. The latter is determined by solving a boundary-value problem for horizontal evaporation (e.g. Mathias and Sander 2021). Salvucci (1997) provides convincing comparisons between his semi-analytical solution and results from numerical solution of the Richards' equation. It is noteworthy that the numerical results were obtained using Brooks and Corey (1966) soil moisture characteristic equations and did not consider the additional effects of vapour diffusion, film flow and coupled heat transfer. Recent numerical analyses suggest that whilst non-isothermal effects are not that important in this context , vapour diffusion and film flow provide a significant contribution to total evaporation (e.g. Walvoord 2002;Li et al. 2019;Iden et al. 2021).
More recently, building on ideas previously considered by van Keulen and Hillel (1974), Fetzer et al. (2017) demonstrated how vapour diffusion can be incorporated within the aforementioned desorptivity term. Given the way in which film flow is typically accounted for by modifying hydraulic conductivity functions (e.g. Scarfone et al. 2020), it should be possible to incorporate film flow into the desorptivity term in the same way. The objective of this article is to implement vapour diffusion and film flow within Salvucci's semi-analytical solution and test its performance by comparison with isothermal numerical simulation and relevant experimental data.
The outline of this article is as follows. The governing equations describing isothermal bare soil evaporation with vapour diffusion and film flow are presented. A numerical solution by method of lines (MoL) is described. An explanation about how to determine desorptivity values using a pseudospectral solution for horizontal evaporation is provided. The semi-analytical solution of Salvucci (1997) is presented. Both the MoL solution and Salvucci's semi-analytical solution are compared to experimental bare soil evaporation data previously presented by Li et al. (2019). A sensitivity analysis is performed to determine where Salvucci's semianalytical solution is likely to perform best. An additional sensitivity analysis is performed to better understand the influence of vapour diffusion and film flow on the desorptivity.

Boundary Value Problem
Consider one-dimensional and isothermal water movement in unsaturated porous media. Mass conservation dictates that (Philip and de Vries 1957) where G [ML −3 ] is mass of water per unit volume of soil, t [T] is time, F [ML −2 T −1 ] is the mass flux of water and z [L] is distance from an evaporation face.
The mass of water per unit volume of soil is found from  -] is the angle between the z axis and the direction in which gravity is acting. The relevant initial and boundary conditions are: is the initial pressure head at the evaporating surface, H [L] is the thickness of the porous medium where an impermeable boundary is located, C v0 [ML −3 ] is the absolute humidity of the air at some reference point adjacent to the evaporating surface of the porous medium, C v (z = 0) is the absolute humidity that corresponds to (z = 0) , and r a [L −1 T] is the aerodynamic resistance characterising the mass transfer of vapour from the evaporating surface to the overlying air. Note that when r a = 0 , the boundary condition at z = 0 reduces to The cumulative evaporation from the surface, E 0 [L], is found from

Humidity Equations
The absolute humidity, C v [ML −3 ], is found from the ideal gas law where M w = 18.02 kg kmol −1 is the molecular mass of water, p v [ML −1 T −2 ] is the vapour pressure in the pore-space, R = 8314 J kmol −1 K −1 is the universal ideal gas constant and T [ Θ ] is temperature, specified in Kelvin. Vapour pressure in the pore-space is found from Kelvin's equation (Edlefsen and Anderson 1943) where p vs [ML −1 T −2 ] is the saturated vapour pressure, g [LT −2 ] is gravitational acceleration and [L] is pressure head. Vapour pressure is related to relative humidity (RH), v (−), by The saturated vapour pressure can be found from Teten's empirical equation (Murray 1967) where p vs is calculated in Pa and T is again specified in Kelvin.

Effective Moisture Content Formulation
Assuming w is constant, Eq. (1) reduces to

Soil Moisture Characteristic Equations
Liquid moisture content is assumed to conform to the van Genuchten (1980) moisture content function in conjunction with the Webb (2000) extension is the effective saturation (which is different to the effective water saturation, S w ), c [-] is a critical liquid moisture content, c [L] is a critical pressure head and d [L] is the pressure head at which the porous medium is completely dry.
The effective saturation is found from (van Genuchten 1980) where [L −1 ] and n [-] are empirical parameters and m = 1 − 1 n . The critical liquid moisture content is found from Note that where S s [L −1 ] is the specific storage coefficient and (van Genuchten 1980) The dry pressure head, d , is taken to be −10 5 m (Webb 2000). The critical pressure head, c , can be found by following the procedure described in Appendix A.

Hydraulic Conductivity
The hydraulic conductivity is assumed to conform to the van Genuchten (1980) and Mualem (1976a) functions in conjunction with a film flow extension, attributed to Tokunaga (2009) where K s [LT −1 ] is the saturated hydraulic conductivity, [-] is an empirical parameter, associated with Mualem's model, and f [L] is an empirical critical pressure head below which conductivity attributed to film flow becomes greater than that attributed to more conventional porous medium flow and K f [LT −1 ] is found from Note that Eq. (28) provides an analogous functional response to those previously described by Zhang (2011), Scarfone et al. (2020 and Peters et al. (2021).

Vapour Diffusion Coefficient
The effective diffusion coefficient for vapour diffusion is found from (in m 2 s −1 ) (Moldrup et al. 2000;Li et al. 2019) where D a (in m 2 s −1 ) is the vapour diffusion coefficient in air, found from (Li et al. 2019)

Three Different Model Formulations
Three different model formulations are compared throughout the article. These include:  (1980) and Mualem (1976a) soil moisture content and hydraulic conductivity functions. See Appendix B for a relevant discussion about imposing the upper boundary condition given by Eq. (6) when no vapour diffusion is accounted for.

Solution by Method of Lines
Similar to Tocci et al. (1997), Goudarzi et al. (2016) and Ireson et al. (2023), a numerical solution for Eqs. (19) and (20) can be developed by applying the method of lines. Let us consider a set of N discrete points on the z-axis: z 1 , z 2 , … , z N . The corresponding set of pressure heads can be written as: 1 , 2 , … , N . A corresponding approximation for the equivalent Darcy flux, q, can be found from which corresponds to the alternative set of points: z 1 reduces to the problem solved by Mathias and Sander (2021) using the Boltzmann transform and a Chebyshev differentiation matrix, hereafter referred to as the pseudospectral solution. Mathias and Sander (2021) previously only considered liquid water. However, the method of Mathias and Sander (2021) can be easily modified to account for film flow and vapour diffusion by utilizing the diffusivity function described by Eq. (16).
The cumulative evaporation from the surface, E 0 [L], is found from (assuming r a = 0) where S [LT − 1 2 ] is the desorptivity, which can be calculated using the aforementioned pseudospectral solution (see Mathias and Sander 2021).
An approximate expression for the cumulative evaporation where r a > 0 can be derived using the time-compression method (Salvucci 1997). For horizontal evaporation, this reads as follows: where and C vI = C v ( = I ) , wich is calculated from a combination of Eqs. (10) and (11). Note that the term, e 0 [LT −1 ] is sometimes referred to as the Stage 1 evaporation rate ). Salvucci (1997) derives a simple semi-analytical solution to account for gravitational effects based on an assumption that the initial moisture content distribution is uniform and that the vertical moisture content profile preserves geometric similarity during evaporation (the validity of this assumption is explored further in Sect. 4.2). For the special case where r a = 0 , Salvucci (1997) arrives at

Salvucci's Semi-analytical Solution
where K wI = K w ( I ) , represents the hydraulic conductivity at the depth of the drying front within the soil profile, defined here as the minimum depth at which = I . Note that E 0 (t = 0) = 0. Salvucci (1997) only provides solutions for E 0 for the special cases where K wI t ≪ E 0 and K wI t ≫ E 0 . In fact, it is straightforward to derive a more general analytical solution for Eq. (43). An implicit solution in terms of time, t, takes the form Solving for E 0 then leads to where and W −1 (x) denotes the minus one branch of the Lambert W function. See Appendix C for more information about how to evaluate W −1 (x).
Also of interest is that solving Eq. (43) for E 0 , substituting the result into Eq. (44) and then solving for t leads to (Salvucci 1997) A corresponding approximate expression for the cumulative evaporation where r a > 0 using the time-compression method reads as follows (Salvucci 1997).
where where e 0 is as defined in Eq. (40). Salvucci (1997) did not intentionally account for vapour diffusion and film flow in his semi-analytical solution. However, his semi-analytical solution was designed to more generally approximate the non-linear advection and diffusion equation given by Eq. (14). Analysis in Sect. 2.1.3 shows that the effects of vapour diffusion and film flow are easily incorporated into the associated non-linear diffusion coefficient and desorptivity, associated with Eq. (38). It follows that Eq. (48) should be able to accommodate vapour diffusion and film flow, providing these have been incorporated into the desorptivity term, S.

Experimental Evaporation Data
Experimental data from Li et al. (2019) are used to verify the applicability of the isothermal evaporation model described above. The complete data set was downloaded from supplementary data provided with an associated PhD dissertation (Li 2020). The experiments involved two, 30 cm deep, soil columns containing a coarse silica sand (#12/20) and a fine silica sand (#50/70). The columns were initially fully saturated with water and encased in impermeable tanks with the top surface of the soil exposed to air in a horizontally orientated wind tunnel. Measurements of RH, temperature and wind speed were recorded within the wind tunnel, 12 cm above the sand surface and maintained at average values of 0.17, 20.8 • C and 1.78 m s −1 (according to the supplementary data provided by Li 2020), respectively.
The average soil surface was found to be 19.9 • C. Both sands were pre-characterised in terms of van Genuchten (1980) and Mualem (1976a) parameters. Cumulative evaporation was measured by monitoring changes in mass within the sand columns. Moisture content was monitored using dielectric soil moisture sensors at 5 cm spacing within the sand columns. A list of model and experimental parameter values are provided in Table 1. The film flow extension point, f , values have been estimated based on considerations described in Sect. 4.1. Figure 1 shows plots of moisture content and hydraulic conductivity against pressure head for the coarse and fine sands studied. The hydraulic conductivity has been further van Genuchten (1980) (1980) and Mualem (1976a) component (the dotted lines). It can be seen that the vapour diffusion component is between 11 and 12 orders of magnitude less than the saturated hydraulic conductivity of the two sands studied. In the absence of film flow, vapour diffusion remains negligible until the pressure head drops below −0.3 m. The film flow curves (the dash-dot lines) have been postulated based on a rough calibration to the cumulative evaporation data, described in more detail in Sect. 4.1. As presented here, it can be seen that film flow is much more effective, as compared to vapour diffusion, in sustaining permeable pathways for pressure heads greater than −100 m. Li et al. (2019) and Li (2020) do not provide information about the film flow characteristics for their coarse and fine sand. To explore the role of film flow in horizontal and vertical evaporation further, we, therefore, also study soil moisture and hydraulic conductivity data from three different soils, with obviously strong film flow components, from a database provided in the supplementary data of Scarfone et al. (2020). These include a sand studied by Nemes (2001), the Shonai sand studied by Mehta et al. (1994) and the Gilat loam studied by Mualem (1976b).

Film Flow Characteristics of Different Soils
Associated moisture content and hydraulic conductivity observations are plotted against pressure head in Fig. 2. Also shown as solid lines are modelled results based on Sects. 2.1.4, 2.1.5 and 2.1.6. Note that the modelled hydraulic conductivity includes vapour diffusion (see Eq. (21)).
The saturated moisture content, s , and saturated hydraulic conductivity, K s , values were taken to be those given by Scarfone et al. (2020). The residual moisture content, r , along with the van Genuchten (1980) parameters, and n, were obtained by calibrating the moisture content function, given by Eq. (23), to the observed moisture content data. The Mualem (1976b) parameter, , and the film flow extension point, f , were then obtained by calibrating the hydraulic conductivity function, given by Eq. (21), to the observed hydraulic conductivity data. Note that the Webb (2000) extension point, c , needs to be uniquely defined for any set of s , r , and n parameters by using the procedure described in Appendix A. Calibration to the observed moisture content data was achieved by minimising the mean absolute error (MAE) between the modelled and observed moisture content data. Calibration to the observed hydraulic conductivity data was achieved by minimising the MAE between the modelled and observed logarithmically transformed hydraulic conductivity data. In both cases, minimisation was achieved using MATLAB's optimisation routine, FMINSEARCH (Lagarias et al. 1998). The resulting set of parameter values is given in Table 2.

Comparison with Experimental Evaporation Data
Results from the MoL solution described in Sect. 2.3 are first compared with two sets of experimental evaporation data from Li et al. (2019) and Li (2020) (see Sect. 3.1). Following Li et al. (2019), a uniform grid spacing with a space-step of 0.25 mm was used for both soil columns. Results from a grid resolution convergence study justifying the use of 0.25 mm are presented in Appendix D. The time-step was selected automatically using the adaptive time-stepping scheme utilised by ODE15s. The model parameters were applied as given in Table 1. Note that these experiments are both vertically orientated so cos = 1.

Fig. 2
Plots of moisture content and hydraulic conductivity against pressure head for the three soil data sets selected from Scarfone et al. (2020). The circular markers are the observed data. The solid lines are due to fitted semi-empirical functions, as described in Sects. 2.1.4, 2.1.5 and 2.1.6. Note that the modelled hydraulic conductivity includes vapour diffusion [see Eq. (21)]. See Table 2 Fig. 3, alongside the observed experimental data of Li et al. (2019) and Li (2020) and results from the semi-analytical solution of Salvucci (1997) Figures 3a and c show results relating to the cumulative evaporation from the two sand columns. The experimental observations exhibit classic Stage 1 and Stage 2 evaporation behaviour . Cumulative evaporation increases linearly with time during Stage 1 due to the constant boundary humidity and constant aerodynamic resistance. Stage 2 evaporation commences when evaporation becomes moisture content limited. Cumulative evaporation during Stage 2 increases non-linearly with time, flattening off as the sands become progressively drier. Similar to as found by Li et al. (2019), the "basic" model (the results shown as thin dash-dot lines) is found to correspond well with experimental data during Stage 1 but then progressively underestimates evaporation during Stage 2. The reason for this is that the "basic" model does not account for evaporation enhancement due to vapour diffusion and film flow .
Again, similar to as found by Li et al. (2019), results from the "+vapour diffusion" model (the results shown as thin dashed black lines) correspond much better with the observed evaporation data during Stage 2 for the coarse sand column (#12/20) (Fig. 3a) but less so for the fine sand column (#50/70) (Fig. 3b). In fact, the vapour diffusion model applied by Li et al. (2019), for the coarse sand column, slightly underestimated the evaporation whereas our model slightly over estimates the evaporation.
The reason why our model evaporates slightly more water is due to the isothermal assumption. Isothermal models assume that heat is perfectly supplied such that temperature is maintained at a constant level despite the presence of evaporation induced cooling. This leads to an overestimation of heat provision and therefore an overestimation of evaporation. Nevertheless, the main point to take from this is that the isothermal assumption is clearly reasonable for investigating the importance of vapour diffusion and film flow in this context.
Our "+vapour diffusion" model can be seen to significantly underestimate evaporation during Stage 2 for the fine sand column (see Fig. 3c). A similar problem was observed by Li et al. (2019). Indeed, Li et al. (2019) attributed the discrepancy with the observed experimental data to ignoring the non-equilibrium between the liquid water and the vapour. They then added a first-order kinetic phase change term to their mass and energy conservation equations, providing them with an additional calibration parameter, enabling them to better match the Stage 2 evaporation response. An alternative approach is to invoke film flow. Recall that we can account for film flow with an alternative single calibration parameter, f , which represents the pressure head below which film flow dominates over conventional porous medium flow. By manually changing f , we found a value of − 0.6 m enabled us to evaporate the correct amount of water during the 40 day observation period. The results from the "comprehensive" model  Fig. 3c. For reference, we found that the maximum value of f that can be sustained before film flow leads to an additional overestimate of evaporation for the coarse sand column was − 0.2 m. The results from the associated "comprehensive" model, for the coarse sand, are shown as a thin solid line in Fig. 3a. Figure 3b shows a comparison of modelled (using the MoL solution) and observed moisture content against time at three different depths within the coarse sand column. All three models correspond well with the moisture content record at 2.5 cm depth and correctly show that evaporation does not affect moisture content at 12.5 cm depth. At 7.5 cm depth, the "basic" model overestimates the moisture content because less water has been evaporated from the column during Stage 2. In contrast, both the "+vapour diffusion" model and the "comprehensive" model underestimate moisture content at this depth because they both slightly overestimate evaporation, probably due to the isothermal assumption. Figure 3d shows a comparison of modelled and observed moisture content against time at three different depths within the fine sand column. Here it can be seen that evaporation affects moisture content as deep as 22.5 cm depth. Furthermore, it can be seen that accounting for film flow leads to significantly improved correspondence between modelled and observed moisture content at 17.5 cm and 22.5 cm depths (compare the thin solid lines and dashed lines).
Results from using the semi-analytical solution of Salvucci (1997) [i.e., Eq. (48)] are shown in Figs. 3a and c, as thick dashed lines. Note that vapour diffusion and film flow are accounted for through the desorptivity term, S. It can be seen that Eq. (48) significantly underestimates evaporation from both of the sands studied. However, despite incorporating vapour diffusion and film flow into the desorptivity term, the results very closely match the evaporation results from the "basic" model for the fine sand column although not for the coarse sand column. Salvucci (1997) derives his equation in a similar way to Green and Ampt (2011). However, Salvucci (1997) points out that his derivation does not explicitly invoke a sharp drying front assumption. Nevertheless, the sharper the soil moisture content relationship, the more likely Salvucci's geometric similarity assumption is to hold. An inspection of Table 1 suggests that the improved correspondence with the fine sand might be because this sand has a higher van Genuchten n parameter as compared to the coarse sand. Indeed, higher n values lead to a sharper soil moisture content relationship.

Numerical Verification of Salvucci's Semi-analytical Solution
A set of numerical simulations was performed to further explore the range of validity associated with Salvucci's semi-analytical solution. These involved repeating the previous "basic" model simulation (ignoring vapour diffusion and film flow) for the coarse sand column but with zero aerodynamic resistance and varied van Genuchten (1980) n parameter. A zero aerodynamic resistance is assumed to avoid additional errors associated with the time-compression procedure (Parlange et al. 2000). The soil moisture content and hydraulic conductivity relationships adopted are plotted in Fig. 4. It can be seen that as n increases the soil moisture content relationship becomes sharper. Figure 5a shows how simulated cumulative evaporation changes with increasing n. The thin solid lines are from the MoL solution, the thick dashed lines are from Salvucci's semi-analytical solution, with sorptivity this time calculated ignoring vapour diffusion and film flow. Note that this evaporation is purely Stage 2 evaporation because r a = 0 . As n increases, the evaporation is reduced. This is because increasing n causes the hydraulic conductivity to reduce faster with decreasing pressure head (see Fig. 4b). As in Fig. 3a, Salvucci's model substantially underestimates the evaporation as compared to the MoL solution when n = 12 . As n is increased, the discrepancy between the two models decreases. However, it seems not possible to achieve perfect correspondence. This is likely due to the Fig. 4 Plots of moisture content and hydraulic conductivity against pressure head, obtained using van Genuchten (1980) and Mualem (1976a) functions with the parameters for the coarse sand (#12/20), listed in Table 1, but with varying van Genuchten n parameter, as indicated in the legend Fig. 5 a Comparison of cumulative evaporation results from the MoL solution (thin solid lines) and Salvucci's semi-analytical solution (thick dashed lines) using the "basic" model for the coarse sand (#12/20) but with the van Genuchten n parameter as stated in the legend and with zero aerodynamic resistance. b Corresponding plots of moisture content against depth from the MoL solution after 1 day (thin dash-dot lines) and 10 days (thin solid lines). The colours relate to the n parameter as stated in the legend for (a). The thin dashed lines were obtained by re-scaling the 1 day profiles using a constant scaling factor on depth, such that the quantity of water evaporated is matched by the 10 day moisture content profiles approximate nature of the limiting assumption made by Salvucci (1997) concerning how the vertical moisture content profile preserves geometric similarity during evaporation.
To explore the validity of this latter assumption further, Fig. 5b shows plots of moisture content against depth at two different times for each of the simulations studied in Fig. 5a. The dash-dot lines and solid lines represent the moisture content profiles after 1 day and 10 days, respectively. The dashed lines were obtained by re-scaling the 1 day profile using a constant scaling factor on depth, such that the quantity of water evaporated is matched by the 10 day moisture content profile. If Salvucci's geometric similarity assumption was perfectly valid, the solid and dashed lines should be close to identical.
In fact, Salvucci's assumption appears to be a reasonably good approximation for all the n values studied. Figure 5b suggests that the validity of the geometry assumption does not provide a strong explanation about why the correspondence between the MoL solution and Salvucci's model improves with increasing n. It is more likely that, contrary to Salvucci's orginal assertion, his equation works better with soils that generate a sharper drying front. The discrepancy between Salvucci's model and the MoL solution increases when vapour diffusion and film flow are included because these lead to additional "tailing" of the hydraulic conductivity function at low pressure heads (see Fig. 4b), which in turn leads to a less sharp drying front.
A question arises as to why Salvucci's model, with the effects of vapour diffusion and film flow included, produces simulation results similar to those from the "basic" model, which ignored vapour diffusion film flow, when considering the fine sand (recall Fig. 3c). Given that vapour diffusion and film flow are only incorporated within Salvucci's model through the desorptivity parameter, S, further insight can be gained by studying the related problem of horizontal evaporation.

Horizontal Evaporation and Desorptivity
A significant advantage of assuming horizontal evaporation when studying the importance of vapour diffusion and film flow is that E 0 (t) = St 1 2 (recall Eq. (38)), where the desorptivity, S, is constant. It follows, assuming a uniform initial effective water saturation, S wI , and a constant boundary RH, v0 , that the relative contributions of vapour diffusion and film flow, to horizontal bare soil evaporation, are constant with time.
Here, we study how the relative contributions of vapour diffusion and film flow vary with initial effective water saturation for the five different soils, collectively described by Figs. 2 and 1, assuming a fixed boundary RH of 17%. This is achieved by calculating desorptivity using the pseudospectral method of Mathias and Sander (2021) It was found that numerically converged results, in terms of desorptivity, for all soil types and all initial saturations could be achieved with 500 Chebyshev nodes. However, obtaining effective water saturation profiles, for the "comprehensive" model formulation, with acceptably low levels of Gibbs phenomenon (for display purposes) required at least 1000 Chebyshev nodes for the coarse sand (#12/20) and fine sand (#50/70); 500 Chebyshev nodes was found to be adequate for the other soils. The extra nodes are needed for the coarse sand (#12/20) and fine sand (#50/70) because of the very high associated values of van Genuchten n parameter. Figure 6a shows plots of desorptivity against initial effective water saturation for each of the five soil types (indicated by the colours in the legend for Fig. 6b) using the "comprehensive", "+vapour diffusion" and "basic" model formulations (recall Sect. 2.2), shown as solid, dashed and dash-dot lines, respectively. For reference, corresponding self-similar effective water saturation profiles for the "comprehensive" model formulation with S wI = 0.98 are shown as plots of effective water saturation, S w , against a dimensionless similarity transform variable, [-], defined by The similarity transform, , has two interesting qualities (see Fig. 7). For a given value of initial water saturation, S wI , plots of water saturation, S w , against will converge on the same asymptotic line for large values of , irrespective of the value applied for the surface

Fig. 6
Desorptivity study for five different soils (see Tables 1 and 2 for model parameters). a Plot of desorptivity, S, against effective initial water saturation, S wI . for different soils, as indicated by the colour legend in (b), and model formulation, as indicated by the line-type legend in (a). b Plots of effective water saturation, S w , against dimensionless similarity transform variable, , assuming S wI = 0.98 . c Plots of % contribution of vapour diffusion to desorptivity against initial effective water saturation. d Plots of % contribution of film flow to desorptivity against initial effective water saturation boundary relative humidity, v0 . For a given value of surface boundary relative humidity, v0 , plots of water saturation, S w , against will converge on the same asymptotic line for small values of , irrespective of the value applied for the initial water saturation, S wI .
It can be seen that there is no difference in desorptivity derived using the three different model formulations for the coarse sand (#12/20), down to an initial effective water saturation of 0.1 (see Fig. 6a). The difference between the three model formulations, at moderate effective water saturations, becomes progressively more pronounced for the fine sand (#50/70), the sand, the Shonai sand and the Gilat loam, respectively. The order is largely due to the residual moisture content values associated with those soils, which are 0.012, 0.030, 0.040, 0.051 and 0.126, respectively (see Tables 1 and 2). The larger the residual saturation, the greater the proportion of pore-water that can only be removed by vapour diffusion and film flow. Figure 6c and d show plots of % contribution of vapour diffusion and film flow to desorptivity, respectively. The % contribution of vapour diffusion is obtained by determining the difference between desorptivities calculated using the "+vapour diffusion" model and "basic model" and dividing by the desorptivity from the "comprehensive" model. The % contribution of film is obtained by determining the difference between desorptivities calculated using the "comprehensive" model and the "+vapour diffusion" model and dividing by the desorptivity from the "comprehensive" model.
With the exception of the Gilat loam, vapour diffusion and film flow represents less than 5% of evaporation for initial water saturations greater than 0.3. In contrast, for the Gilat loam, film flow represents 60% of evaporation where the initial water saturation is 0.4 and 2.5% where the initial water saturation is 0.98. The Gilat loam performs differently in this way because of its much larger residual moisture content, compared to the other soils, combined with a relative large film flow extension point, f (see Table 2).

Fig. 7
Plots of effective water saturation, S w , against dimensionless similarity transform variable, , for different initial effective water saturations, S wI , and boundary RH, 0 , for the Gilat loam and utilizing the "comprehensive" model formulation The reason why vapour diffusion and film flow have such little relevance to horizontal evaporation when the initial saturation is high, can be further understood by studying Figs. 8a and b, which show plots of effective water saturation, S w , against the dimensionless similarity transform, , with S wI set to 0.98, for the Gilat loam and the fine sand (#50/70), respectively.
The blue, green and red lines were derived using the pseudospectral solution with the "comprehensive", "+vapour diffusion" and "basic" model formulations, respectively. The pseudospectral solutions were evaluated with 2000 Chebyshev nodes (extra nodes are needed to adequately capture the sharper fronts associated with the "+vapour" diffusion model formulation for the fine sand (#50/70)). The green lines in Fig. 8 are from MoL solutions for the "comprehensive" model formulation, shown for model verification purposes. These were achieved using 400 logarithmically distributed spatial grid points spanning from zero to 1000 m, with the smallest space-step starting at 1 m. The profiles were evaluated after one day of evaporation. There is excellent correspondence between the pseudospectral solution and the MoL solution, verifying the equivalence of both schemes.
The variable is linearly proportional to the distance from the fixed RH boundary, and is shown on a logarithmic scale. The area above the red lines (the "basic" model) (in Fig. 8) represent the amount of water that will be evaporated in the absence of vapour diffusion and film flow. The differences between the red and the green (the "+vapour diffusion" model) lines represent the additional evaporation achieved due to vapour diffusion. The differences between the green and the red (the "comprehensive" model) lines represent the additional evaporation achieved due to film flow.
For the Gilat loam, it can be seen (Fig. 8a) that the vapour diffusion contribution exists mostly in the region, < 10 −4 , and the film flow contribution exists mostly in the region, 10 −4 < < 10 −3 . In contrast, the non-vapour diffusion and non-film flow contribution of evaporation exists mostly in the region, < 10 −1 , making it almost two orders of magnitude larger than the film flow contribution. For the fine sand (#50/70), the vapour diffusion Fig. 8 Plots of effective water saturation, S w , against dimensionless similarity transform variable, , assuming an initial effective water saturation, S wI = 0.98 , and a boundary RH, 0 = 0.17 . The thick green lines are from MoL solutions with the "comprehensive" model formulation. The black, blue and red lines are from pseudospectral solutions with the "comprehensive", "+vapour diffusion" and "basic" model formulations, respectively. (a) shows results for the Gilat loam. (b) shows results for the fine sand (#50/70) and film flow contribution are mostly derived from the region, < 10 −6 , explaining why these processes represent close to zero contribution to desorptivity at this initial saturation (recall Figs. 6c and d).
The results show that the vapour diffusion and film flow fronts only need to travel a very short distance into the porous medium to hydraulically link the evaporation face with water behind the main drying front (beyond which the moisture content is still at its initial value). This explains why desorptivity is largely insensitive to vapour diffusion and film flow at high initial water saturations, and also explains why Salvucci's model, for the fine sand (#50/70), produced close to identical results as compared to the "basic" model despite incorporating vapour diffusion and film flow through the desorptivity parameter (recall Fig. 3c).

Comparison of Horizontal and Vertical Evaporation
The results in Sect. 4.3 show that vapour diffusion and film flow are almost negligible during horizontal evaporation from soil with a relatively high initial saturation. In contrast, the results in Sect. 4.1 show vapour diffusion and film flow represent a significant proportion of total evaporation during vertical evaporation due to the presence of gravity. Figure 9, shows results from simulating horizontal and vertical evaporation in the fine sand (#50/70) using the MoL solution with the "comprehensive" model formulation. The parameters used for both models were identical to those listed in Table 1 except that the aerodynamic resistance, r a = 0 , and the boundary depth was set to 1000 m for the horizontal evaporation simulation. The vertical evaporation simulation used the same 0.25 mm grid spacing, previously employed. The horizontal evaporation simulation used the logarithmically growing grid described in the previous section. The need for a logarithmically growing grid and a 1000 m deep boundary for the horizontal evaporation problem is due to horizontal evaporation being much more efficient, leading to a drying front that develops much further way from the boundary surface, as compared with vertical evaporation scenario.
Figures 9a) and b) compare cumulative evaporation from the two soil columns. As previously shown by Novak (2022), accounting for gravitational effects leads to a significant decrease in evaporation. After 40 days, the horizontal column has produced 274 cm whereas the vertical column has produced 6.26 cm. Note that the reduced evaporation amount for vertical flow conditions is not due to the smaller column length, because it is observed that the drying front does not reach the model boundary for this scenario. The results for the vertical column are similar to those previously shown in Fig. 3c except that the initial evaporation rate is much higher, on account of the zero aerodynamic resistance.
The dash-dot and dashed lines in Fig. 9b) are from the "basic" model and Salvucci's semi-analytical solution, respectively. As discussed in the previous section, the effect of film flow and vapour diffusion are negligible for horizontal evaporation. However, for vertical evaporation, the difference between the "comprehensive" and "basic" model formulations is substantial, suggesting that, after 40 days, film flow and vapour diffusion account for 28% of total evaporation. Also note that the results from Salvucci's semi-analytical solution are closest to the "basic" model despite the desorptivity accounting for vapour diffusion and film flow. This is because vapour diffusion and film flow provide a negligible contribution to horizontal evaporation in this context.
Recall that the depth of the drying front was previously defined as the minimum depth at which = I . Unfortunately, for the horizontal evaporation case, this definition of drying front occurs where the similarity transform, = ∞ . Let us therefore consider an approximate location of the drying front, defined as the maximum depth at which the effective moisture content is less than 99% of its initial value. The approximate location of the drying front after 40 days is 96.9 m depth for the horizontal evaporation case (see Fig. 9c) and 0.319 m depth for the vertical evaporation case (see Fig. 9d). This substantial difference explains the massive variation we see in cumulative evaporation plotted in Figs. 9a and b. The movement of the drying front is slowed down by gravity.
The vapour diffusion and film flow fronts move into the sand more slowly under horizontal flow conditions as compared to under vertical flow conditions (compare Figs. 9e and f). After 40 days, the film flow front is at 1 mm depth for the horizontal evaporation problem and 100 mm depth for the vertical evaporation problem. Hence, vapour diffusion and film flow are found to represent a much larger contribution to evaporation when flow is vertical.
The main reason for the differences observed between horizontal and vertical evaporation can be described as follows. In the absence of gravity, capillary flow is highly efficient at bringing water to the boundary surface where evaporation is taking place. Consequently, the boundary surface remains insufficiently dry for vapour diffusion and film flow to become dominant processes. For vertical evaporation, gravitational forces inhibit upward capillary flow of water such that the boundary surface becomes much dryer. In this dryer environment, vapour diffusion and film flow become much more important. Unfortunately, Salvucci's semi-analytical solution is unable to accommodate this point because vapour diffusion and film flow can only be represented through the desorptivity parameter, which is based on horizontal evaporation.

Summary and Conclusions
The objective of this study was to implement vapour diffusion and film flow within Salvucci's semi-analytical solution and test its performance by comparison with isothermal numerical simulation and relevant experimental data.
In a first set of examples, we compared our numerical model with results from Salvucci's semi-analytical solution (Salvucci 1997) and experimental data from Li et al. (2019). Salvucci's semi-analytical solution was modified to account for vapour diffusion and film flow through its desorptivity term. Excellent correspondence between the numerical model results and experimental data were observed when vapour diffusion and film flow were accounted for. However, Salvucci's semi-analytical solution was found to significantly underestimate Stage 2 evaporation and yielded similar results to those of the numerical solution when vapour diffusion and film flow were ignored. A further set of comparison simulations revealed that Salvucci's approximation is less effective for situations where hydraulic conductivity curves exhibit long "tailing" behaviour at high capillary pressures. This is particularly the case where vapour diffusion and film flow are important.
The desorptivity parameter in Salvucci's semi-analytical solution represents the proportionality constant relating Stage 2 cumulative evaporation to the square root of time under horizontal flow conditions. A set of horizontal evaporation simulations were therefore developed to explore the relative contributions of vapour diffusion and film flow to desorptivity as a function of initial saturation.
For four out of the five soil types studied, vapour diffusion and film flow were found to account for less than 5% of the total evaporation for initial saturations greater than 30%. This is because, under horizontal flow conditions, the vapour diffusion and film flow fronts only need to travel a very short distance into the porous medium ( < 1 mm) to hydraulically link the boundary surface with water behind the main drying front (beyond which the moisture content is still at its initial value). Desorptivity is therefore largely insensitive to vapour diffusion and film flow at high initial water saturations. This explains why Salvucci's approximation yields similar results to those from numerical solutions where vapour diffusion and film flow were ignored.
A set of simulations was repeated for a high initial saturation case but assuming vertical flow conditions (i.e., with gravitational effects included). Gravitational effects were found to significantly inhibit the upward flux of water from deep within the soil column. Hence total evaporation was greatly reduced by the presence of gravity. However, pores at the near surface became much drier, as compared to when gravity was ignored, because it was harder for liquid water from depth to replenish evaporated water at the surface. A consequence of this drier region at the near surface was that vapour diffusion and film flow played a more important role in the transport of water when gravitational effects were accounted for. In the presence of gravitational effects, vapour diffusion and film flow were found to account for 28% of total evaporation despite having a negligible effect when gravity was ignored.
The following important conclusions can be drawn. Analytical solutions that assume horizontal flow conditions are inadequate for understanding vertical evaporation problems because they overestimate evaporation rates and mostly predict vapour diffusion and film flow to be of negligible influence. Salvucci's semi-analytical solution is effective at predicting the order-of-magnitude reduction in evaporation caused by gravitational effects. However, it is unable to identify the correct importance of vapour diffusion and film flow because these processes can only be represented through its desorptivity parameter. Webb (2000) tells us that the critical pressure head, c , in Eq. (23), should be chosen to ensure that the gradient, w , is continuous at = c . Given Eq. (26), it follows that c is found by iterative solution of the following equation

Finding Ã c for Webb's Extension
For this study we used MATLAB's solver FMINSEARCH Lagarias et al. (1998), which requires a seed value for c , we will denote as c0 . When r = 0 , it is found that c0 = −10∕ works well. When r > 0 it is better to use an asymptotic result due to Kelly and Mathias (2018), namely that where and W 0 (x) is the zero branch of the Lambert W function and The Lambert W function can be approximated in this case using the asymptotic expansion (Corless et al. 1996) (53) where L 1 = ln(x) and L 2 = ln(L 1 ).

Note About the Top Boundary Condition for the "Basic" Model
The "basic" model described in Sect. 2.2 ignores the presence of water vapour in the soil and only solves for liquid water. Previous researchers chose not to apply the boundary condition specified in Eq. (6) in this context, due to associated uncertainty about the value of C v (z = 0) . Two alternative approaches are generally considered instead. The first takes the form Li et al. 2019) where crit [L] is an additional empirical parameter. In this way, Stage 1 evaporation occurs when (z = 0) > crit and Stage 2 evaporation occurs when (z = 0) ≤ crit . Unfortunately there is a great deal of uncertainty about how to quantify crit Fetzer et al. 2017;Li et al. 2019). Furthermore, the magnitude of predicted Stage 2 evaporation is strongly dependent on its value (consider Fig. 5 in Li et al. 2019). An alternative approach is to assume Li et al. 2019) where r s [L −1 ] is a soil resistance term representing the resistance to vapour diffusion between the soil surface and the point within the soil where p v = p vs . The r s term is thought to continuously increase with decreasing soil moisture content in the near surface of the soil, therefore negating the need for the discontinuous transition to Stage 2 evaporation imposed by Eq. (58). Several empirical equations, relating r s to the average moisture content in the top layer of the soil, have been proposed in the literature (see Li et al. 2019, for a list of relevant references) and these vary widely depending on soil-type and experimental setup (Li et al. 2019).
Interested readers should compare results from the "basic" model in Figs. 3a and c with those presented in Fig. 6 of Li et al. (2019) to confirm that Eqs. (58) and (59) provide very little added value.

Evaluating the Lambert W Function for Salvucci's Equation
The Lambert W function in Eq. (45) is particularly challenging due to the negative exponential argument in Eq. (46). The numerical challenge of evaluating the Lambert W function with negative exponential argument (i.e. W −1 (−e z ) ) was dealt with previously by Mathias and Roberts (2013). The solution is to consider the asymptotic expansion of W −1 (x) for x < 0 and x → 0 − (Chapeau-Blondeau and Monir 2002) where L 1 = ln(−x) and L 2 = ln(−L 1 ). The exponential term in Eq. (46) is eliminated by substituting Eq. (46) into the definition of L 1 to obtain Grid Resolution Convergence Study Following Li et al. (2019), we chose to use a space-step of 0.25 mm for the finite difference discretistion in the MoL solutions presented in Sect. 4.1. Figure 7 in Li et al. (2019) shows results from a grid resolution convergence study confirming that 0.25 mm should be adequate in this context. Nevertheless, we chose to perform an additional grid-convergence study using the "comprehensive" model for the coarse sand column (#12/20) described in Sect. 4.1. The results are shown in Fig. 10. Simulations were performed with space steps of 1, 0.5, 0.25, 0.1 and 0.05 mm.
Figures 10a, b and c show results in terms of moisture content after 40 days, cumulative evaporation and evaporation rate, respectively, using space steps of 1, 0.25 and 0.05 mm. A minor difference can be observed between the 0.05 mm and the 0.25 mm simulations. However, at the scale of plotting, the 0.25 mm and 0.05 mm simulations are indistinguishable. Figure 10d shows plots of percentage error for the cumulative evaporation value at 40 days. Percentage error was determined in reference to the value obtained from the finest resolution simulation, which had a space step of 0.05 mm. It can be seen that the simulation accuracy increases almost linearly with decreasing space-step. Nevertheless, it is also confirmed that the error associated with the 0.25 mm simulation is less than 1%, which is arguably sufficient for graphical inspection purposes as applied in this article. Funding The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.
Data availability Spreadsheets containing the experimental data from Li (2020) are available at the associated URL.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 10
Results from the grid resolution convergence study of the "comprehensive" model for the coarse sand column (#12/20) described in Sect. 4.1 with space steps, Δz , as indicated in the legends. a Moisture content plotted against depth after 40 days. b Cumulative evaporation plotted against time. c Evaporation rate plotted against time. d Percentage error relating final cumulative evaporation from the different models to that from the 0.05 mm model, plotted against space step, Δz