Improved calculation of soil hydraulic conductivity with the simplified evaporation method

Numerical modeling of soil water dynamics and storage is generally based on the Richards equation. Its solution requires knowledge of the soil hydraulic properties (SHP): the soil water retention function and the hydraulic conductivity function. To determine SHP, laboratory evaporation experiments are particularly popular because they provide data for both SHP functions. The evaluation by the simplified evaporation method (SEM) method, originally proposed by Schindler and subsequently improved by several authors, relies on linearization assumptions that allow for a relatively simple calculation scheme but result in biased conductivity data for some soils. The objective of this study is to propose and test an improved computational scheme for the hydraulic conductivity function. We present the new theory and show that it leads generally to higher accuracy of the conductivity function. The improvement is most pronounced for sandy soils and soil water pressure heads below −100 cm, where the original method provided data with bias.


INTRODUCTION
An accurate determination of soil hydraulic properties (SHP) is required for the prediction of soil water dynamics. A commonly used procedure to estimate the SHP for a specific soil is the simplified evaporation method (SEM), first presented by Schindler (1980) and later improved by Peters and Durner (2008) and Peters et al. (2015). The SEM consists of an experimental laboratory step followed by two computational steps: data evaluation and fitting procedure. The experimental step is carried out on a soil column initially saturated with Abbreviations: SEM, simplified evaporation method; SHP, soil hydraulic properties.
The assumptions of SEM allow SHP to be obtained with relative simplicity, requiring only a simple analytical procedure for data evaluation and a nonlinear fitting procedure to obtain the parameters for the SHP functions. The "current SEM" (considered here as the procedure described in Peters et al., 2015) assumes linear variations of ℎ and θ with depth as well as a linear course of the water flux across the column center between two consecutive time steps in the data evaluation. Generally, the inaccuracies caused by these simplifying assumptions tend to cancel each other out for many soil textures and ranges of suction (Peters & Durner, 2008;Peters et al., 2015). However, the pronounced nonlinearities occurring in well-sorted sandy soils pose a challenge. As a result, the soil hydraulic functions estimated by SEM show systematic errors for sandy soils (Iden et al., 2019;Peters et al., 2015) and differences compared to those obtained by inverse modeling with the Richards equation for soils rich in organic matter (Dettmann et al., 2019). The differences arise primarily from nonlinearities in the suction h as a function of the vertical position in the sample. Bezerra-Coelho et al. (2018) confirmed that using the current SEM gave accurate results for the soil water retention curve, even when estimating parameters of a dual-porosity soil, but problems were encountered in accurately determining the hydraulic conductivity curve for an extremely coarse medium with a very high van Genuchten n value.
In this paper, we introduce a new method for evaluating the data from evaporation experiments to determine point values of the soil hydraulic conductivity function, (ℎ). The new method is composed of three improvements to calculate ( , ℎ) value pairs and relies on modified assumptions, which are better suited to capture the nonlinearities of variably saturated water flow during evaporation, particularly during stage two (Iden, Blöcher et al., 2021;. The objective is to reduce the bias of the (ℎ) function obtained with the current SEM by using the new approach to estimate the hydraulic conductivity of the soil.

Data evaluation with the simplified evaporation method
The procedure of the evaporation experiment for the SEM is characterized by the measurement of the suction ℎ (cm) at two vertical positions ("depths") in a soil sample using tensiometers. At the same time, the average water content in the sample over time is calculated from the total mass of the column and the knowledge of the dry mass of the soil and the mass of all the equipment placed on the scale. Initially, the soil sample is saturated with water, and evaporation occurs from its surface until the sample is relatively dry and tensiometers stop

Core Ideas
• An improved scheme for the calculation of hydraulic conductivity with the simplified evaporation method (SEM) is presented. • Bias in hydraulic conductivity for sandy soils was observed in the original SEM is eliminated. • Improvement is most significant for soil water pressure heads in the medium to dry range, below −100 cm. • The calculation scheme can be easily implemented in existing computer algorithms. functioning. Observations are registered periodically and at the end of the experiment, a time series composed of average water content and tensiometer readings at the two positions in the sample are available. The dry mass of the soil sample is determined after the evaporation experiment by oven-drying. Figure 1 illustrates the experimental setup of an evaporation experiment. The height of the sample is (cm), the average volumetric water content of the soil at a discrete time t i (day) isθ and the vertical positions of the tensiometers are 1 = 0.25 and 2 = 0.75 , respectively. The distance between the tensiometers is therefore 2 − 1 = 0.5 . The vertical axis is positive upwards and = 0 at the bottom of the soil sample. Commonly used soil sample rings for the evaporation experiment have = 4 cm and = 5 cm, where is the radius.
The procedure to calculate points of the water retention curve θ(ℎ) is described in Peters et al. (2015) and is not repeated here because the emphasis lies on improving the calculation scheme for hydraulic conductivity (ℎ). The calculation of (ℎ) is based on the Buckingham-Darcy law. The water flux density across the center of the sample, (cm day −1 ), is assumed to be 50% of the evaporation rate (Schindler, 1980) and is calculated as follows: The gradient in suction is given as: (2) where ℎ 1, and ℎ 2, are the suctions measured by the two tensiometers at time . The original SEM calculates the suction at the center, ℎ ari, (cm), of the soil sample as the arithmetic mean of the four suctions at the two depths and times t i-1 and t i (Peters & Durner, 2008): Based on the suction observations at two subsequent times, the hydraulic conductivity is calculated from the Buckingham-Darcy law as: ( positive upward, ℎ is suction, not pressure head). Peters et al. (2015) proposed improvements to the SEM process to account for the spatial nonlinearities in suction and water content. Specifically, they changed the considerations for the representative mean suction and introduced a mixed mean ℎ mix (mixed between arithmetic and geometric):

Current state of the simplified evaporation method
where ℎ geo, is the geometric mean of the observed suctions at both positions in the sample and the weighting factor ϕ is calculated as the inverse of the gradient in suction ∇ℎ.

New method to remove bias in the hydraulic conductivity function
As shown by Iden et al. (2019), the improved scheme of Peters et al. (2015) still leads to a bias in the calculated conductivity data for light-textured (sandy) soils. This is caused by the fact that the spatial distribution of suction can become highly F I G U R E 2 Illustration of measured suction ℎ at the upper and lower tensiometers over time ( ) and the predicted evolution of the suction at the center of the soil sample by Equation (6)  nonlinear in sandy soils. In general, the errors in the calculated data points of (ℎ) are caused by a wrongly estimated average suction (abscissa of (ℎ) function) and/or the calculated hydraulic conductivity (ordinate). As the latter is calculated from the suction gradient and the water flux across the center of the column, three error sources can be distinguished, corresponding to the calculation of: 1. the average suction to which the conductivity is assigned; 2. the suction gradient at the center of the column; 3. the water flux density across the center.
In the following, we introduce three revised calculation schemes that decrease these errors and their propagation into the (ℎ) function.

2.3.1
Prediction of the suction at the center of the sample The methods described in the previous sections analyze the time series ℎ 1 ( ) and ℎ 2 ( ), where subscript 1 refers to the upper tensiometer position and subscript 2 to the lower one. We propose here an inverted approach, using the reverse functions 1 (ℎ) and 2 (ℎ), as illustrated in Figure 2. These functions return the time of arrival of a suction ℎ at the position of the two tensiometers. Our approach is to calculate the suction at the center of the column as a function of time by applying the reverse function. We select a suction ℎ measured early by the upper tensiometer and later on by the lower one and calculate the average time for this suction to be reached at the center position , by: F I G U R E 3 Time series of suction at different locations in the column. The dashed lines represent measured (upper and lower tensiometer, ℎ 2 and ℎ 1 ) and simulated (true value at center, "ℎ True") values, the solid lines indicate values approximated with the method proposed in this article. The estimated suction at the center of the soil sample profile is ℎ c , the estimated suction at distance Δ * above the center of the soil sample is ℎ + , and the suction at Δ * below the center of the soil sample is ℎ − . The value of Δ * was 0.2 cm.
where (-) is a weighting factor. For = 1, the average becomes arithmetic, and if w → 0, the average converges toward the geometric mean. For > 1, the arrival time c, becomes larger than the arithmetic (and geometric) mean of the two times. In preliminary tests with synthetic data, = 2 resulted in a minimum deviation between the actual time of arrival and the time of arrival approximated by Equation (6). For all missing values of ℎ , simple linear interpolation was applied to fill gaps. Figure 2 illustrates the calculation.

2.3.2
Prediction of the vertical suction gradient at the center of the sample The calculation of the gradient in suction is based on the concept outlined in the preceding section. In general, it allows calculating the position and time of an arbitrary suction ℎ in between the tensiometers. To calculate ∇ℎ at the vertical center of the sample, the values of the suction at a small distance above, ℎ + , and below, ℎ − , the center of the soil sample profile are required. When using the arithmetic mean to calculate c, , which is equivalent to assuming a constant velocity of downward movement of a specific value ℎ, this procedure is intuitive, but for the average considering = 2 in Equation (6) it is less straightforward. Considering Δ * (cm) as the small vertical distance between ℎ and ℎ + (or ℎ − ), to calculate the respective time for ℎ + and ℎ − : where 1 and 2 are the respective values of time for the upper and lower tensiometers to reach a specific value of ℎ, and F I G U R E 4 Process flow for obtaining ℎ , ℎ + , and ℎ − (black, green, and yellow lines, Graph 4) from tensiometer data (orange and blue lines) via their respective reverse functions, , + , and − (black, green, and yellow lines, Graph 3). Graphs (1) and (4) show suction ℎ as a function of time . Graphs (2) and (3) show time as a function of suction ℎ.
Δ( 2 ) eq [d 2 ] is the square of the time it takes a specific value of ℎ to go from = 0.5 + Δ * to = 0.5 − Δ * . Therefore, the respective times for ℎ + and ℎ − can be calculated as follows: Calculating + and − for each specific ℎ value, the respective value pairs (ℎ + , + ) and (ℎ − , − ) are obtained. A simulation was conducted with Rehovot Sand (Mualem, 1976) with a potential evaporation rate of 5 mm day −1 to illustrate the proposed theory. The results can be observed in Figure 3, which shows the shape of ℎ c , ℎ + , and ℎ − as a function of time, together with the interpolated measurements of the upper tensiometer, lower tensiometer, and the true values of ℎ .
Finally, the suction gradient at the vertical center of the soil at each time is calculated as: where 2Δ * is the distance between the vertical positions corresponding to suctions ℎ + and ℎ − at a specific time. Figure 4 illustrates the process of obtaining ℎ c , ℎ + , and ℎ − and the gradient ∇ℎ c, . In steps 1 and 2, the reverse function for the tensiometer data over time is calculated. From steps 2 to3, the values for c (black line), + (green), and − (yellow) are calculated as averages (Equations 6 and 8). Steps 3 and 4 refer to the transformation of (ℎ) into ℎ( ). The gradient is calculated using ℎ + (green) and ℎ − (yellow) in Figure 4 using Equation (9).

2.3.3
Prediction of the water flux across the center of the sample The flux density at the center of the soil sample, c (cm day −1 ), is calculated based on the measured sample weight, the readings of the bottom tensiometer, and the estimate of ℎ c . Two consecutive measurements of ℎ 1 with estimates of ℎ c can be transformed into c as follows. The first step is to estimate the water content at the center of the soil sample, θ c , and the water content related to the lower tensiometer θ 1 . For this, interpolated values from the water retention curve θ(ℎ), that is, data pairs (θ , ℎ mix, ), are used to approximate θ as function of ℎ ("ℎ → θ"). The interpolated values for ℎ are calculated as described in Section 2.3.1. The corresponding water contents are defined as: In the following step, the water storage, wc (cm), is calculated as: considering that = 0.5 ( + 1 ). After calculating wc for both consecutive times, −1 and , the central water flux density c, at the intermediate time 0.5 ( + −1 ) is calculated as: The new c, substitutes the original occurring in Equation (4) for the calculation of . A similar approach to calculate the water flux density at the column center was proposed by Wendroth et al. (1993). Figure 5 shows a schematic illustration of the calculation of c . The θ( 2 , ) and θ( 1 , ) are the true water content profiles at times 2 and 1 , respectively. The respective angles between the profiles with the horizontal at the bottom of the sample are represented by β 1 and β 2 . The areas A and B represent the displaced water depth wc,i − wc,i−1 within the time interval − −1 . A and B are parallelograms in both plots. If the difference between β 1 and β 2 is small enough, the size of A can be considered the same in both images, no matter the angle (parallelogram area rule). This condition is assured with a small time step. Equation (12) considers β 1 = β 2 for the estimative of . The same rule applies for area B.

Data analysis
To evaluate the proposed procedure, synthetic measurement data were created by simulating evaporation experiments using HYDRUS-1D (Šimůnek et al., 2013). The data analyzed in this article were taken from Iden et al. (2019), who simulated evaporation experiments with the coupled liquid water, water vapor, and heat flow model implemented in Hydrus-1D, which was presented by Saito et al. (2006). For the virtual evaporation experiments, the simulations considered water retention and liquid water flow in capillaries and films (Peters, 2013). The SHP were parameterized with the Kosugi-PDI model. Details of the PDI model family (Peters, Durner and Iden model) are given in Peters (2013Peters ( , 2014 and Iden and Durner (2014), and the Kosugi model for the parametrization of the capillary part of the water retention curve is described in Kosugi (1996). Three different sets of SHP, representing a broad range of textures (Rehovot sand, sandy loam, and clay loam) were used for the simulations and for testing the new calculation scheme of the SEM. Table 1 presents the soil hydraulic parameters for each soil. All simulations considered a soil sample height of L = 5 cm, 101 equidistant nodes for the spatial discretization with the finite element method, and a simulation time of 15 days. The initial condition was specified as a hydrostatic pressure head distribution with ℎ = 0 cm at the bottom. At the soil surface, the surface energy balance was solved and the evaporation rate was calculated within Hydrus-1D from the water vapor pressure difference between soil and atmosphere and atmospheric resistance. The details of the simulation and a discussion of the simulated evaporation dynamics are given in Iden et al. (2019). Two levels of white noise were added to the suction data representing the tensiometer readings to investigate the proposed evaluation scheme. In the low noise scenario, an independently, identically normally distributed random error with a standard deviation of σ ℎ = 0.2 cm was added to all F I G U R E 5 Schematic illustration of the simplified evaporation method for the estimate of the water flux density at the center of the soil sample. The areas A and B represent the displaced water depth during a time step. The right border of the light blue area is the water content profile at time 1 , the left border is the water content profile at time 2 ; is the vertical position of the center of the sample, 1 is the vertical position of the tensiometer, and b is the average between c and 1 .
T A B L E 1 Soil hydraulic parameters for the Kosugi-PDI model: ℎ m (cm) is the suction corresponding to the median pore radius, σ * (−) is the standard deviation of the lognormal density function, θ r (−) is the residual water content, θ s (−) is the saturated water content, τ (−) is the tortuosity factor, s (cm day −1 ) is the saturated hydraulic conductivity and ω (−) is the weighting factor for the capillary and film flow.

T A B L E 2
The root-mean-square deviation (RMSD) of the logarithm values (base 10) of the fitted (ℎ) functions compared to the true function using different added noise levels (σ ℎ ) obtained by the modified simplified evaporation method (M-SEM) and current simplified evaporation method (C-SEM) methodology. simulated suctions, and for a high noise scenario, a standard deviation σ ℎ = 1.0 cm was used. For all soils and scenarios, the estimates for (ℎ) points were analyzed and compared considering the new method described in Section 2.3 and labeled "M-SEM" (modified simplified evaporation method). A value of 0.2 cm was used for Δ * . The current scheme by Peters et al. (2015) is labeled "C-SEM" (current simplified evaporation method). For the comparison of the estimated SEM data, a nonlinear fitting procedure was performed to find the best conductivity function parameters for the model considering the noisy data for M-SEM. Then, the root-mean- square deviation (RMSD) between the true and the resulting fitting function for each soil was calculated. For this comparison, values between pF 1.75 and pF 5.00 [pF defined as pF = log 10 (h cm −1 )] were used for RMSD calculation.  Figure 6 shows the true conductivity function together with results obtained by the C-SEM and by the proposed M-SEM for the Rehovot sand soil based on data without noise. Visually, the M-SEM estimate provides a result closer to the true curve. The bias of the C-SEM was discussed in Iden et al. (2019) for sandy soils, and it is caused by nonlinearities of the gradient and the approximation of the flux at the center of the sample. The noisy data in the range pF 1.0-pF 1.8 are due to computational numerical inaccuracies. Such scattering of (ℎ) data is typical for evaporation experiments when the hydraulic conductivity of the soil is much larger than the evaporation rate, because this leads to very small gradients in hydraulic potential which are highly affected by noise (Peters & Durner, 2008). Figure 7 presents the estimated (ℎ) data for Rehovot Sand for the scenarios with added noise. The data in the graph on the left and the right represent the low and high noise levels, respectively. The higher noise level leads to a higher scattering in the estimated (ℎ) data, but the bias appears remarkably small. The data in the gray-shaded tension range area are considered nonidentifiable since the hydraulic gradient is too small (Peters & Durner, 2008). They were thus not considered for the subsequent fitting of the (ℎ) function.

RESULTS AND DISCUSSION
The fitted conductivity functions of the noisy scenarios are shown in Figure 8, together with the true functions that were used to create the synthetic measurement data. All three soils and both scenarios are displayed. For sandy loam and clay loam, the fitted curves are overlapping, showing neither an improvement nor a worsening in the fitting quality when comparing the modified method with the current method. Furthermore, the agreement between the true functions and the identified ones is always good, confirming the results reported by Peters et al. (2015). For the Rehovot sand soil, a considerable improvement is observed when using M-SEM compared with C-SEM. The bias with C-SEM depicted in Figure 6 and described by Iden et al. (2019) disappeared when using the proposed modified method.
The RMSDs between the fitted and true functions quantify the improvement using the M-SEM as compared to C-SEM ( Table 2). The M-SEM provided lower RMSD for all soils and noise levels. For the Rehovot sand soil with a low noise level, the improvement was most pronounced. For the higher noise level ( σ ℎ = 1.0 cm), the difference was smaller, but still significant, with a 26% reduction of the RMSD.

CONCLUSIONS
A new method for estimating the hydraulic conductivity function (ℎ) by the SEM was proposed. Modifications to the current method involved estimating a more representative mean suction, suction gradient, and water flux density across the central plane of the sample. The proposed method was tested using synthetic data with noise added to the suction data. It provided an excellent agreement between the estimated and the true (ℎ) function in all studied scenarios and for both noise levels. Specifically, the bias in (ℎ) for pure sand that was problematic in the former procedure almost vanished with the newly proposed procedure. For the testing, we considered synthetic tension data up to 3000 hPa. This range is beyond the measurement limit of common field tensiometers but can be reached with modern laboratory tensiometers with boiling delay (Schindler, Durner, von Unold, & Müller, 2010). Furthermore, extensions by using the defined air entry point of porous tensiometer cups (Schindler, Durner, von Unold, Mueller, & Wieland, 2010) or hygrometers (Iden, Blöcher et al., 2021; can be included in the evaluation. Thus, our proposed method already lays the foundation for future application of SEM with improved sensors.

SUPPLEMENTARY MATERIAL
Python scripts implementing the new calculation scheme for the SEM presented in this article are available at https:// github.com/infoleon/M.SEM.

AU T H O R C O N T R I B U T I O N S
Leonardo Inforsato: Conceptualization; formal analysis; investigation; methodology; software; supervision; writingoriginal draft; writing-review and editing. Sascha Christian Iden: Methodology; supervision; writing-review and editing. Wolfgang Durner: Project administration; supervision; writing-review and editing. Andre Peters: Methodology; supervision. Quirijn de Jong van Lier: Funding acquisition; project administration; supervision; writing-review and editing.

A C K N O W L E D G M E N T S
This study was supported by São Paulo Research Foundation (FAPESP), grant numbers 2021/10467-0, 2021/10520-8 and 2018/20902-2.