Prediction of soil evaporation measured with weighable lysimeters using the FAO Penman–Monteith method in combination with Richards’ equation

Multiannual data (2016–2018) from 12 weighed lysimeters (four soil types with textures ranging from sandy loam to silt loam, three replicates) of the TERENO SOILCan network were used to evaluate if evaporation (E) rates could be predicted from weather data using the FAO Penman–Monteith (PM) method combined with soil water flow simulations using the Richards equation. Soil hydraulic properties (SHPs) were estimated either from soil texture using the ROSETTA pedotransfer functions, from in situ measured water retention curves, or from soil surface water contents using inverse modeling. In all years, E was water limited and the measured evaporation rates (Em) surprisingly did not vary significantly among the four different soil types. When SHPs derived from pedotransfer functions were used, simulated evaporation rates of the finer textured soils overestimated the measured ones considerably. Better agreement was obtained when simulations were based on in situ measured or inversely estimated SHPs. The SHPs estimated from pedotransfer functions represented unrealistically large characteristic lengths of evaporation (Lc), and Lc was found to be a useful characteristic to constrain estimates of SHPs. Also, when soil evaporation was water limited and Em rates were below Epot (PM evaporation scaled by an empirical coefficient), the diurnal dynamics of Em followed those of Epot. The Richards equation that considers only isothermal liquid water flow did not reproduce these dynamics caused by temperature dependent vapor transport in the soil.


INTRODUCTION
Bare soil evaporation (E) is a critical component of the water cycle at local and global scales. The process of mass transfer across the soil-atmospheric interface is fundamental (Davarzani, Smits, Tolene, & Illangasekare, 2014). Evaporation in arid or semiarid areas can amount to more than half of evapotranspiration (Huxman et al., 2005). Or and Lehmann (2019) estimated that E contributes up to 15% of the total evapotranspiration losses globally and found that this fraction is relatively constant for different climates and biomes.
Processes of E are generally well studied, but it has been a challenge to represent E in simulation models (Seager et al., 2007). E is mediated by a set of processes such as flow of liquid water, and transport of water vapor and heat in the porous medium that are coupled with heat and vapor transport in the atmosphere, also referred to as free flow, and driven by absorbed radiation . The E process can be split up into two phases: an energy-limited E (Phase I) and a water-limited E (Phase II). During Phase I, E is determined by the energy available (radiant energy absorbed at the surface, soil heat flux, sensible and latent heat flux between the soil surface, and the free flow). During Phase II, E is limited by water transfer from the soil to the evaporation surface. There are several approaches of different complexity to model E, and a systematic review of the processes and their simplifications is given by Vanderborght, Fetzer, Mosthaf, Smits, and Helmig (2017). Fetzer et al. (2017) compared these approaches using simulation studies and highlighted differences and similarities between them. To model accumulated evaporation rates (i.e., cumulative evaporation [E c ]) over a longer period, simpler approaches can be used. E rates during Phase I can be estimated using the standardized FAO Penman-Monteith (FAO-PM) approach (Allen, Pereira, Raes, & Smith, 1998). The FAO-PM estimates E rates using the PM equation that solves the surface energy balance for standardized properties of a reference surface representing a wellwatered uniform grass surface, linearizes the relation between saturated vapor pressure and temperature, and approximates the soil heat flux as a function of the absorbed radiation. To simulate the transition between Phase I and II, simpler models use a threshold surface water potential that is fixed and used as a constant potential boundary condition during Phase II. The evaporation rate during Phase II is obtained by solving the Richards equation, which describes isothermal liquid water fluxes in the unsaturated soil, for this constant boundary condition. Despite the fact that this approach decouples the evaporation rate during Phase II from the surface energy balance, heat, and vapor fluxes in the soil, it provided similar estimates of E c as the more comprehensive process models . However, diurnal dynamics of E simulated by the more comprehensive approaches could not be reproduced by the simpler FAO-PM-Richards equation approach due to neglected vapor fluxes and heat flow in the porous medium . Li, Vanderborght, and Smits (2019) compared three different model concepts to determine E: the Richards equation only considering isothermal liquid water flow (Richards' model), a model based on the Richards equation that considered vapor transport (Richards' vapor model), and a nonequilibrium twocomponent two-phase model (Non 2-2 model). They mentioned that predictions of E during Phase II by models based on Richards' equation depend on the accuracy of chosen top boundary conditions like the critical water pressure. Rianna,

Core Ideas
• Lysimeter data is used to evaluate the FAO-PM method and Richards' equation. • Simulated evaporation was very sensitive to soil hydraulic properties (SHPs). • SHPs based on soil texture led to overestimation of evaporation losses. • Hydraulic parameters based on in situ measurements produced better estimates of evaporation (E). • Characteristic length of E can be used to evaluate if SHPs yield realistic predictions of E. • Soil temperature-dependent E dynamics in Phase II cannot be reproduced by Richards' equation.
Reder, and Pagano (2018) tested four models and their ability to predict measured evaporation (E m ) by comparing them with E rates from lysimeters. The models tested did not solve the Richards equation directly to represent soil water movement but used different metamodels of the soil water balance: the Aydin (1998) and Allen et al. (1998), Boesten and Stroosnijder (1986), and Ritchie (1972) (also known as the FAO-56 dual K c approach) models. They stated that these E models with soil hydraulic parameters based on lysimeter data reproduce E losses accurately. In addition to selecting a process model, which fits best to given environmental conditions, soil hydraulic properties (SHPs) must be known in order to predict E. When only liquid flow in the soil is considered, the unsaturated hydraulic conductivity and water retention function must be available. In most studies that estimate soil water fluxes at regional and global scales, water retention and hydraulic conductivity functions are derived from soil texture information using pedotransfer functions (PTFs; Fatichi et al., 2020;Or & Lehmann, 2019). These PTFs have been derived from empirical relations between soil properties (e.g., soil texture) and parameters of functions that describe the water retention and hydraulic conductivity curve (Schaap, Leij, & Van Genuchten, 2001). Parameters of functions that are fitted to retention curves and that are estimated from other soil properties using PTFs are not constrained directly by processes, such as E. This may lead to a large overestimation of evaporation when SHPs with these parameters are used in the Richards equation to estimate evaporation (Lehmann, Bickel, Wei, & Or, 2020). Lehmann et al. (2020) used a characteristic length for evaporation (L c ; Lehmann, Assouline, & Or, 2008), to identify hydraulic parameter sets that would lead to unrealistic high evaporation losses. Especially soils with an intermediate texture were found to be susceptible to overestimation of E. Soil hydraulic properties are also affected by soil structure, which is related to the presence of larger pores and the connectivity of smaller pores between soil aggregates Assouline & Narkins, 2019;Shokri, Lehmann, & Or, 2010). Especially in agricultural soils, the structure of the surface layer during bare soil periods may deviate significantly from that of the deeper soil layers due to tillage practices, which change bulk density and soil structure. The impact of soil structure on soil water fluxes has been associated mainly with (preferential) infiltration. However, also for E, soil structure may play an important role Haghighi & Or, 2015;Schwartz, Baumhardt, & Evett, 2010;Tran et al., 2019).
In vadose zone research, weighable lysimeters are an established method to investigate E and to calibrate and test models (Dijkema et al., 2018;Mutziger, Burt, Howes, & Allen, 2005;Pütz, Fank, & Flury, 2018). By measuring soil water storage changes and seepage water losses with high temporal resolution and precision with weighable lysimeters, diurnal dynamics of E fluxes can be directly calculated giving a closer look into the processes controlling E (Hoffmann, Schwartengräber, Wessolek, & Peters, 2016;Peters et al., 2017;Pütz et al., 2018). Lysimeters are used to observe and compare soil E fluxes under natural conditions for soils with different hydraulic properties, different tillage, or under different climate conditions (Dijkema et al., 2018;Marek et al., 2016;Quinn, Parker, & Rushton, 2018). Dijkema et al. (2018) compared lysimeter data from a desert site with simulations using a model including nonisothermal liquid and vapor transport. They showed that under dry conditions, the simulated E rate was underestimated, which they attributed to underestimations of the vapor transport process through the upper dry soil layer and of the hydraulic conductivity of the dry soil that was estimated using a capillary bundle model. Recent studies emphasize the importance of appropriately choosing the lysimeter's boundary conditions that reflect surrounding (natural) environmental conditions in order to achieve a reliable reproduction of the natural system (Benettin, Queloz, Bensimon, McDonnell, & Rinaldo, 2018;Groh et al., 2018Groh et al., , 2020Heinlein, Biernath, Klein, Thieme, & Priesack, 2017;Kohfahl et al., 2019;Teuling, 2018;Trigo et al., 2018). The mimicking of a natural system also implies the use of undisturbed or monolithic soil cores or blocks in lysimeters. Our overarching objective was to compare lysimetric-based measurements and estimates of E to simulations using an FAO-PM-Richards equation based approach. Subobjectives were to examine (a) how SHPs that are derived from PTFs or from direct in-situ measurements perform in the prediction of E; (b) whether the characteristic length for evaporation (L c ) can be used an indicator of parameter sets that generate unrealistic estimates of E, and (c) the (diurnal) dynamics of evaporation during Phase II, which are influenced by temperature dependent vapor fluxes.

Lysimeter installation
Two lysimeter clusters (SE3 and SE4), managed as arable land, were used in the present study. At each cluster, six lysimeters were arranged in a hexagonal pattern ( Figure 1). Each lysimeter cylinder, made of stainless steel with a surface area of 1.0 m 2 and a depth of 1.5 m , is placed in a porous concrete ring. The gap between lysimeter and the concrete housing was closed with a synthetic resin collar, which was be filled with surrounding soil. The small (0.01 m) remaining gap between the collar and the lysimeter was closed with silicone rubber sheet to ensure that soil temperature dynamics reflect those of the surrounding. Water content is measured by time domain reflectometry probes (TDR100, Campbell Scientific) at 10-, 30-, and 50-cm depth. The weight of the lysimeters is recorded in 1-min intervals. The matric potential is measured with matric potential sensors in 0.1-m depth (MPS-1, Decagon Devices); at other depths, tensiometers are installed (TS1, METER Group). The lysimeters' weighing systems have an accuracy of 0.01 mm (Hannes et al., 2015). The seepage water of a lysimeter is collected in a weighable tank with a resolution of 1 g, which corresponds to 0.001 mm . The matric potential near the bottom of the lysimeter (1.4-m depth) is adapted to the measured matric potential in the surrounding field at the same depth by water pumping out or in through the suction rake, which consists of five parallel suction cups, at the bottom of the lysimeter. If the matric potential in the lysimeter increases compared to the field, water is pumped out of the lysimeter into the seepage water tank and vice versa when the matric potential decreases in the lysimeter. This allows downward  and upward vertical water fluxes across the lysimeter bottom and ensures that the water potential dynamics at the bottom of the lysimeter are as observed in the surrounding field (Karimov, Šimůnek, Hanjra, Avliyakulov, & Forkutsa, 2014;Luo & Sophocleous, 2010). The adaption of the matrix potential to the surrounding conditions is based on the space-for-time concept of the TERENO SOILCan project. A weather station (WXT510, Vaisala Oyj) next to each lysimeter station monitors meteorological parameters: wind speed, air temperature, relative humidity, barometric pressure at 2-m height, and precipitation at 1-m height (OTT Pluvio 2 , OTT HydroMet). Net radiation sensors (LP Net07, Delta OHM) above one lysimeter of each cluster and cameras document the state of the lysimeters .
The 12 lysimeters were filled with undisturbed soil originating from the four sites within the SOILCan network in Germany: Bad Lauchstädt (BL), Dedelow (DD), Sauerbach (SB), and Selhausen (SE) (for information related to these stations, see Groh, Vanderborght, Pütz, & Vereecken, 2016 and more detailed profiles are described in Annex 1 (supplemental material). Soil surrounding the lysimeters was treated in the same manner as the lysimeters. Soil evaporation was monitored in the arable land lysimeters for three consecutive years (2016-2018) between harvest and sowing when there was no or a little vegetation and E was dominated by bare soil characteristics. During each bare soil period, the soil was tilled twice. In 2017, herbicides were sprayed due to the growth of weeds.  Table 2.
The observation period in 2018 was characterized by dry conditions with only 119 mm of precipitation during the bare

Data processing
The distinction of precipitation and E from lysimeter mass changes needs a suitable data preprocessing and postprocessing to reduce the effect of errors and noises on the water fluxes. Lysimeter and leachate tank weights were processed in the following manner: (a) the raw data were corrected by adjusting systematic errors and outliers, (b) random errors and noise were smoothed, (c) E and precipitation fluxes were calculated, and (d) data were checked for plausibility and data gaps were filled in a post-processing step (Schrader et al., 2013).
Steps i and ii: Error correction and smoothing. A distinction was made between measurement noise (e.g., induced by external sources such as wind) and outliers (defined as singular disturbances or technical malfunctions) (Peters et al., 2017). In some cases, outliers and noise could be removed by automated plausibility controls or suitable filter schemes; in other cases, manual corrections were required (Schrader et al., 2013).
To smooth the time series of mass measurements obtained with 1-min resolution and avoid smoothing errors on the estimation of precipitation and evapotranspiration, the adaptive window and adaptive threshold filter method (AWAT filter) developed by Peters, Nehls, Schonsky, and Wessolek (2014); Peters, Nehls, and Wessolek (2016); and Peters et al. (2017) was used. To manage time-variable noise levels, the filter uses an adaptive threshold (δ) and adaptive averaging window width (ω) . For our calculations, d varied between 0.02 (0.02 kg) and 0.2 mm (0.2 kg), and w between 3 and 31 min. An intermediate output of the filter algorithm is a step function of lysimeter weights. This step function is subsequently interpolated using spline functions to produce a time series with 1-min temporal resolution of noise-removed and interpolated lysimeter weights. Further information and an extensive validation of this method can be found in Hannes et al. (2015) and Peters et al. (2014Peters et al. ( , 2016Peters et al. ( , 2017. Step iii: Calculation of evaporation and precipitation. E and precipitation were calculated from time series of lysimeter (WL) and seepage tank weights (WT) obtained in Step ii as in Schrader et al. (2013), assuming no simultaneous E and precipitation during a 1-min time interval: where P i (kg) and ET i (kg) are the precipitation and evapotranspiration amounts during the time interval i. Time series of P and E were subsequently aggregated to hourly values.
The obtained E rates were tested for plausibility. For times with missing data, three cases were considered to fill the data gaps. In case data were missing in one or two of the three replicate lysimeters (2.9% of the dataset), the mean E rate of the other lysimeters with the same soil profile was taken to fill the data gap. In case data were missing in the three lysimeters of the same soil type for a period shorter than 6 h (0.31% of the dataset), fluxes were assumed to be constant during this period and estimated from the weights at the end and beginning of the time period. If data gaps exceeded 6 h (0.05% of the data set), E fluxes were set equal to E rates of the other lysimeters.

Statistical analysis
A multivariate nested ANOVA was carried out to evaluate the variability in E flux among the different soil types and within lysimeter replicates (i.e., filled with a given soil type). For this, 33 d with gap-filled data were excluded from the analysis. The considered factors in the ANOVA analysis were soil type, time, and lysimeter. The lysimeter factor was treated as a random factor nested within each soil type. The daily E rate was the investigated dependent variable and expressed as where E i,j,k is the daily E rate (mm d −1 ) at day i for soil j and lysimeter k, E . . . is the overall average of the daily E rate over all days, soils and lysimeters (the dotted subscript index means that E is averaged over that index), E i.. is the average effect of day i on the E (the mean value of E rate at day i minus the overall E rate mean value), E .j. is the average effect of soil j on the E (the E rate mean value for soil j minus the overall E rate mean value), E ij. is the time-soil interaction effect (the difference between the average E in the three lysimeters of soil j at day i and the overall average plus the mean time and soil effects), E .(j)k is the mean effect of lysimeter k of soil j (mean difference between the E in lysimeter k and the mean E of all three lysimeters of soil j), and ε is the error. Since only one measurement was done in one lysimeter for a given time, the statistics of the error term had to be derived from the lysimeter-time interaction effect (the difference between the interaction between the lysimeter and time effect and the average interaction effect of the soil and the time). The ANOVA analysis was carried out with the MAT-LAB (Mathworks, 2019) anovan function. The same ANOVA model was used to evaluate the soil moisture content measurements in the different lysimeters and soils.

Computation of potential evaporation
Reference evapotranspiration, ET 0 (mm d −1 ) was estimated with the PM method (Allen et al., 1998): where λ (MJ kg −1 ) is the latent heat of evaporation, R n is the net radiation (MJ m −2 d −1 ), G is the soil heat flux (MJ m −2 d −1 ), (e se a ) is the vapor pressure deficit (kPa), ρ a is the air density (kg m −3 ) at a constant pressure, and c p is the specific heat of the air (MJ kg −1˚C−1 ). The slope of the saturation vapor pressure temperature relationship is expressed as Δ (kPa˚C −1 ), whereas γ is the psychometric constant (kPå C −1 ). Aerodynamic and surface bulk resistance (s m −1 ) are denoted by r a and r s , respectively. R n was measured over one bare soil lysimeter and used as a representative measure for all lysimeters. Differences in soil color or darkness between different lysimeters could therefore not be considered. The aerodynamic resistance is derived from the wind speed considering a logarithmic wind profile. For calculating the evapotranspiration of a reference grass surface, the surface resistance represents the stomatal resistance and is defined to be r s = 50 s m −1 for daytime and r s = 200 s m −1 for nighttime periods (Allen, Pruitt, Raes, Smith, & Pereira, 2005). The soil heat flux G (MJ m −2 d −1 ) is calculated from the net radiation and put equal to 0.1R n during daytime and 0.5R n during nighttime (Allen et al., 1998). To obtain the potential bare soil evaporation E pot (mm d −1 ), the calculated reference ET 0 (mm d −1 ) was scaled by an empirical coefficient :

Simulation of soil water flow
E pot values reflect situations when the soil surface is water saturated, and E is energy controlled (Phase I; Idso, Reginato, Jackson, Kimball, & Nakayama, 1974). When the E m rate becomes less than E pot , E is soil controlled (i.e., movement of water to the soil surface is a limiting factor [Phase II]). In order to simulate the periods when the E takes place at a potential rate (Phase I), a soil water flow model was used that simulates the temporal evolution of the surface soil water content so that time periods when the surface is sufficiently wet to maintain the E pot rate can be identified. In a second instance, the model was used to simulate the evaporation rate (E sim ) during Phase II when the flux of soil water to the surface is limiting the E rate.
To describe the dynamics and vertical distribution of the soil water content and water fluxes in the soil profile, the Richards equation was used. The Richards equation simplifies the two-phase (liquid and gas) two-component (dry air and water) flow and transport, which are coupled to the heat transport, to an isothermal (no coupling to heat transport in the soil) one-component (only water) and one-phase (only liquid) equation . It thereby assumes that the effect of the air phase on the liquid flow is minor, that flow due to thermal gradients can be neglected, and that transport of water vapor in the gas phase can be neglected compared with the transport of water by the liquid flow (Li et al., 2019). The one-component, one-phase Richards' model for one-dimensional vertical flow in a soil profile without plants given by Šimůnek, Šejna, and van Genuchten (2013): where h is the water pressure head (m), θ the volumetric water content (m 3 m −3 ), t the time (s), x the vertical spatial coordinate (m, positive upward), and K(h,θ) the hydraulic conductivity function (m s −1 ), which depends on pressure head or water content. The soil water retention and hydraulic conductivity function complete Equation 5. To solve Equation 5, boundary conditions at the soil surface must be defined. A common approach is to define so-called "atmospheric" or mixed boundary conditions (Šimůnek et al., 2013;Vanderborght et al., 2017). During Phase I and during precipitation events, the E pot flux (Equation 4) and the precipitation rate define the flux boundary condition. These flux boundary conditions are used as long as the water pressure head at the surface stays above or below a critical threshold pressure head for evaporation or infiltration conditions, respectively. When the threshold pressure head is reached, the boundary condition is shifted to a pressure head boundary condition. For the minimum allowed pressure head at the soil surface, a critical pressure head of −1,000 m was used. The periods when the pressure head boundary condition is used correspond with Phase II periods and the calculated water fluxes at the soil surface with the Phase II E rates. This approach to estimate Phase II E rates is based on the idea that these rates are predominantly limited by liquid water fluxes in deeper soil layers (Salvucci, 1997). Assuming that thermal and vapor fluxes have a low impact, the so calculated fluxes were found by Assouline et al. (2013); Assouline, Narkis, Gherabli, Lefort, and Prat (2014); and Assouline and Narkis (2019) to reproduce quite well the Phase II E rates. However, simulations by Fetzer et al. (2017) demonstrated that these assumptions may not hold true during daytime, when water transport through the dried out and heated up surface layer is mainly due to vapor transport. The higher temperatures during the day increase the vapor pressure in the soil resulting in higher vapor concentration gradients and vapor transport. This incomplete description of the subsurface vapor diffusion and evaporation dynamics is also described by Shahraeeni and Or (2010, 2012a, 2012b; Shokri, Lehmann, and Or (2009);and Or, Lehmann, Shahraeeni, and Shokri (2013).

HYDRUS-1D model setup, soil hydraulic properties, and initial conditions
To solve the Richards equation, we used the HYDRUS-1D code. HYDRUS-1D is an established model that has been applied successfully for the simulation of various hydrologic processes and problems in the vadose zone (Šimůnek, Šejna, Saito, Sakai, & van Genuchten, 2018). For estimating the unsaturated hydraulic conductivity function in terms of soil water retention parameters, the description by van Genuchten-Mualem (vGM) was chosen (van Genuchten, 1980): where α (m −1 ) is the inverse of the air entry value, n the pore size distribution, l the pore connectivity, and m is defined as m = 1 − 1/n, where n > 1. S e is the effective saturation, and θ s (m 3 m −3 ) and θ r (m 3 m −3 ) are the saturated and residual water contents, respectively. Two sets of hydraulic parameters were used in this study: a first set for the vGM model were estimated from soil texture  using the ROSETTA PTF implemented in HYDRUS (Schaap et al., 2001) (Figure 2). Values for θ s were derived from measured maximal soil water contents in 2016, 2017, and 2018 (I). A second set of SHP parameters was derived by fitting the water retention curve (Equation 6) to measurements of water content and water potentials in the lysimeters at the 10-, 30-, and 50-cm depth (Table 4, II). Comparing texture dependent SHP (I), there is a difference between the two different soil types, silty loam (BL, SB, and SE), and sandy loam (DD) (Figure 2). The fitted SHPs (II) are very similar for all soils and do not show differences that can be related to difference in texture ( Figure 2).
For one soil, SE, a third set of hydraulic parameters was considered. These parameters were derived from L-band brightness temperatures, which are related to the soil moisture content of the surface soil layer and which were measured using a tower-based L-band radiometer that looked on three plots with different soil surface structure: a compacted, harrowed and tilled soil plot (Dimitrov, 2015). Parameters of the surface soil layer were estimated using inverse modeling with a radiative transfer model that was coupled with the Richards equation. In order to represent the effect of the soil structure changes due to soil surface management, the bimodal model of Durner (DBM) was used to represent the SHP of the tilled surface layer (Durner, 1994;Priesack & Durner, 2006). The tilled surface layer was assumed to be 5 cm thick, which corresponds with the common tillage depth of the lysimeters. The SHP parameters of the soil below the tilled surface layer were derived from the soil texture using the ROSETTA PTF. The effect of tillage on the SHP was assumed to be nondynamic or static. In reality, tillage is done each growing season, and rainfall events after tillage might change SHP. Further information about the estimation of the hydraulic parameters is provided by Dimitrov et al. (2014).
The tilled surface layer had a larger saturated hydraulic conductivity (K s ) than the undisturbed soil ( Figure 2). For more negative h, K tilled is smaller than that of the nondisturbed soil representing the impact of the empty interaggregate pores and the hydraulic disconnection of the aggregates.
We assumed a vertically homogeneous soil texture profile with one single set of parameters (Table 4). This assumption was made since the soil texture did not vary considerably with F I G U R E 2 Hydraulic conductivity (K) plotted against the pressure head (h, top) and the water retention curve (bottom) for the four different soil types: Selhausen (SE), Bad Lauchstädt (BL), Sauerbach (SB), and Dedelow (DD). As a reference, curves for clay and sand defined by the ROSETTA database were added. Pedotransfer function (PTF) soil hydraulic property (SHP) presents texture-based van Genuchten-Mualem parameters, and the fitted SHPs are derived from in situ measurements of water contents and pressure heads. Parameters for the tilled surface (DBM) are given by radiometer measurements as described in Dimitrov et al. (2015). θ is the volumetric water content depth and since simulation results that used hydraulic properties derived from soil layer and lysimeter specific soil texture were similar to simulations that assumed vertically uniform soil texture (results not shown).
Initial conditions were set up by interpolation of soil moisture contents measured at 10-, 30-, and 50-cm depth at the beginning of each observation period. There were neither sensors at the surface nor bottom of the lysimeter. Hence, the value of the shallowest sensor at 10 cm in the lysimeter was considered instead, whereas saturation was considered to prevail at 1.50-m depth as initial condition. "Atmospheric" boundary conditions were defined at the upper boundary, and a variable pressure head boundary condition was determined at the bottom boundary for 2017 and 2018, whereas a flux boundary condition was applied in 2016.

Characteristic length of evaporation
We defined the characteristic length of evaporation, L c , as the depth to which a wet soil dries out during Phase I evaporation. In order to define this depth, we calculated the amount of water that can be evaporated during Phase I, E cum,I (m). This amount was divided by the change in water content in the dried out soil layer, which was approximated by (θ i − θ r )/2: where θ i is the initial water content and θ sur is the water content at the soil surface E cum,I was derived from the soil "desorptivity" (S evap , m d −0.5 ) and the potential evaporation rate (E pot ), which was T A B L E 4 The hydraulic parameters of the van Genuchten-Mualem model (vGM) for the different soils estimated from soil texture using ROSETTA pedotransfer functions implemented in HYDRUS-1D Schaap et al., 2001)  Note. The saturated water content (θ s ) is taken from maximum field values. Also, the soil hydraulic parameters estimated from fits of the retention curve to water content and pressure head measurements in the lysimeters (fitted SHP) are shown. For the Selhausen soil, parameters of the Durner bimodal hydraulic functions that were obtained by inverse modelling of surface soil moisture and account for the impact of the changed structure of the surface soil layer of a tilled plot (Dimitrov, 2015) are also shown. θ r is the residual water content, θ s is the saturated water content, α and n are shape parameters, K s is the saturated hydraulic conductivity, w represents the weight of the pore domain to total pore volume, and τ is the tortuosity parameter in the conductivity function. θ s values were derived from field measurements.
taken to be 4 mm d −1 , using the time compression analysis : 2 evap was derived from the hydraulic conductivity and the water retention curve according to Parlange, Vauclin, Haverkamp, and Lisle (1985): where h i is the initial pressure head and h sur is the pressure head at the soil surface. θ sur was set equal to the residual water content, θ r and h sur to −15,000 cm [at this pressure head θ(h) ≈ θ r ]. The initial water content was defined as the water content at h i = −10 cm.

Evaluation of model performance
Deviations of predictions and observations were evaluated by using the RMSD. The RMSD can be divided into the unsystematic difference RMSD u and the systematic difference RMSD s , which are defined as (Wilmott, 1982;Wilmott et al., 1985) where x i refers to E m and y i is E SHP/PTF/tilled , and̂is calculated from the ordinary least-squares (OLS) linear regression between E SHP/PTF/tilled and E m .

Potential and measured evaporation: Total amounts, monthly and daily averages
In Figure 3, E pot , E m from the different soils, and the precipitation (P), are shown for the monitoring periods in the three different years. For all years, E m < E pot , E pot > P, and E m does not differ considerably among the soil types. The difference between E pot and E m was the largest in 2018 and the smallest in 2017, whereas precipitation was the largest in 2017 and the smallest in 2018. Despite the fact that E pot > P, E m was  smaller than the precipitation in all years so that there was an infiltration surplus in all years that wetted up the soil profile. The infiltration surplus was largest in 2016. Since the timing and duration of the monitoring period in 2016 is different from that in the other two years, a direct comparison of the values for this year with the other two years is difficult. According to the ANOVA analysis, the difference in timeaveraged E m among soil types was not significant ( Table 5). The time-averaged E rates of the 12 different lysimeters were significantly different independently from the soil type. Their rates varied too much within a given soil type to observe significant soil effects on the time-averaged E, and consequently cumulative E. The standard deviation of the soil effect was 0.014 mm d −1 , whereas the standard deviation of the time-averaged E rates from the different lysimeters was 0.027 mm d −1 . In other words, the soil effects (Table 6) were too small compared to the variation of the E between the different lysimeters. The soil-time interaction effects were significant. This means that despite the nonsignificant differences between the time-averaged E rates of the different soils, the temporal variations of the E rates varied significantly between the different soils. It should be noted that the analysis was carried out using daily averaged E rates. Therefore, the analysis does not provide information about differences in evaporation rates between soils and individual lysimeters at a higher temporal resolution (e.g., diurnal variation in E). The day-to-day "unexplained variation" or "model error" was 0.09 mm d −1 (

F I G U R E 4
Monthly mean diurnal variation of potential and actual evaporation rates in each soil for the three different observation periods (2016, 2017, and 2018). BL, SE, SB, and DD are the different soil types: Bad Lauchstädt, Selhausen, Sauerbach, and Dedelow, respectively. E pot describes the potential evaporation. The vertical bars represent the standard deviation of the day-to-day hourly evaporation rate at a given hour of the difference between the E rate at a given day of a single lysimeter, which is corrected for its average lysimeter effect, and the average rate of lysimeters from the same soil at that day. The time effects were strongly significant (i.e., average E rates in the different lysimeters of the different soils varied significantly over time compared with the "unexplained" variation). Figure 4 presents the monthly mean diurnal variation of the E pot and E m in the four soils. In all cases, the peak of the mean diurnal E m was less than that of E pot . In months when E m was significantly less than E pot (August 2016 and 2018, September and October 2018), the peak of the mean diurnal E m flattened off. The mean diurnal evaporation values show that waterlimited or Phase II evaporation occurred during all months. E m values during Phase II depend on the SHPs (Mosthaf, Helmig, & Or, 2014). We therefore expected to observe differences in the E rates between the different soils in these months, especially between the DD soil with a sandier texture and the other soils. But the mean diurnal E m variation was very similar for the different soil types.

Water content
Temporal dynamics in soil water content (θ) were stronger in shallow layers (10-and 30-cm depths) than at 50 cm, regardless of the investigated bare soil period ( Figure 5). The A systematic difference in water content at a given soil depth among the different soil types was not detected in the ANOVA analysis (i.e., it was systematically smaller than the difference in water content at the considered depth among lysimeters of the same soil type). The standard deviation of the lysimeter effects on soil moisture measurements in a certain soil was between 0.023 m 3 m −3 at 10-and 30-cm depth and 0.03 m 3 m −3 at 50-cm depth, whereas the standard deviation of the soil effect on soil moisture was only 0.011 m 3 m −3 at 10-cm depth, 0.019 m 3 m −3 at 30-cm depth, and 0.021 m 3 m −3 at 50-cm depth. The unexplained error in the ANOVA model was 0.023, 0.024, and 0.019 m 3 m −3 at 10-, 30-, and 50-cm depth, respectively.

3.4
Model results

Cumulative evaporation rates
Simulations using SHPs that were derived from soil texture with PTFs highly overestimated measured evaporation.
Although there was barely any difference in E PTF between the silty loam soils (BL, SB, and SE), the sandy loam soil (DD) E PTF rates were less in comparison with BL, SB and SE in all years. Although E PTF >> E m , differences between textures can be seen in simulations contrary to F I G U R E 6 Comparison of cumulative precipitation (P), potential (E pot ), measured (E m ), and simulated evaporation for the different soils in different years. Simulations using hydraulic parameters obtained from pedotransfer functions (E PTF ) and from in situ measurements of water content and pressure heads (E SHP ) are shown. For the SE soil, simulations using parameters obtained from surface soil moisture measurements and inverse modeling are also shown (E tilled ). Shaded areas show the range in evaporation of lysimeters filled with the same soil. BL, SE, SB, and DD are the different soil types Bad Lauchstädt, Selhausen, Sauerbach, and Dedelow, respectively observations. The highest E PTF rates were simulated for 2016 with E sim ≈ E pot for the silty loam soils. The lowest overestimation by E PTF was in 2018 for the DD soil (+35.5 mm). Using SHP parameters obtained from fits to in situ measurements of water content and pressure heads led to better estimates of E (E SHP ) compared with E PTF (Figure 6). For all years and soils, E SHP < E PTF and E SHP corresponded better than E PTF with E m . E SHP was smaller for the silty loam BL and SE soils than for the silty loam SB and the sandy loam DD. For the latter two soils, E SHP was still larger than E m , especially in 2016.
Further simulations with HYDRUS-1D were carried out for the SE soil using SHPs that were estimated from surface soil moisture measurement of a tilled soil plot, E tilled . Due to tillage, the hydraulic properties of the upper layer showed similarities to those of a coarse textured layer . A coarser layer on top of the soil profile accelerates the decrease in hydraulic connectivity of the soil surface with the soil profile, which reduces the duration of Phase I and causes a rapid transition towards Phase II and much smaller E rates during Phase II (Assouline et al., 2014). This was also observed indirectly in the tilled plots that were monitored by Dimitrov (2015) at the SE site. Observed surface temperature differences between the tilled and compacted plots corroborated the differences in simulated evaporation fluxes between both plots that resulted from the different SHPs of the plots due to the different structure of the surface layer. Implementing the tilled layer in HYDRUS-1D resulted in E tilled < E pot , in contrast with the simulation for a homogeneous soil profile using the PTF-derived SHPs, which predicted E m ≈ E pot for the silty loam (BL,SB, and SE) in 2016 (191 mm). The tilled surface layer limited E tilled to 107 mm, which is closer to the observations (58 mm). A reduced E tilled compared with E PTF was also simulated in 2018 (99 mm), and in 2017, E tilled = E m (Figure 6).
For a more detailed analysis of this system and the influence of a tilled layer on E rate dynamics, more field information is needed (e.g., detailed measurements above 10-cm depth, nearsurface atmospheric conditions, tillage, etc.). Note again that T A B L E 7 Characteristic lengths of evaporation calculated from soil hydraulic parameters (SHPs) that were derived from soil texture using pedotransfer functions (PTFs) (L cPTF ), from in situ measured soil water contents and pressure heads (L cSHP ), and from soil surface water contents in a tilled soil using inverse modeling (L c,tilled  a tilled layer of 5-cm thickness is an assumption, and using a thicker or thinner surface layer may have an impact on the E tilled rates.

Characteristic lengths of evaporation
The soil water retention curves that are obtained for the different soils using PTFs, using in situ measurements of soil water contents and pressure heads, or using inverse modeling show important differences (Figure 2). For instance, the water retention curves derived using PTFs suggest a narrower pore size distribution than the retention curves that were measured in situ, which is reflected in the smaller n values of the in situ measured curves. These differences are translated in the soil hydraulic conductivity curves. For the same saturated hydraulic conductivity, the Mualem model predicts considerably smaller unsaturated hydraulic conductivities when the n parameter decreases. The hydraulic conductivity curves that are based on the in situ measured soil water retention curves resemble the unsaturated hydraulic conductivity of a clay soil. Also, the unsaturated hydraulic conductivity of the tilled surface layer is considerably smaller than the unsaturated hydraulic conductivity that was derived from PTFs. How these differences in soil hydraulic properties translate into differences in simulated soil evaporation can be characterized using the characteristic length of evaporation (L c , Table 7). Soil hydraulic parameters that are derived from soil texture using PTFs corresponded for the silty loam soils with very high characteristic lengths (i.e., around 300 cm). The large overestimation of E m by E PTF demonstrates that large characteristic lengths do not represent the evaporation properties of these soils. This confirms the argument of Lehmann et al. (2020) that estimates of SHPs should best be constrained by characteristic lengths of evaporation. Soil hydraulic parameters that were derived from fitting the water retention curve to in situ measured water content and pressure head measurements corresponded with smaller L c (10-50 cm) and bet-ter predicted E m . Additionally, the tilled surface layer had a smaller L c . It should be noted that the ranking of L c reproduced well the ranking of the simulated E in the different soils and different years. L c can therefore be considered as a good characteristic that translates SHPs to evaporation properties of a soil. Since the measured cumulative evaporation rates were reproduced the best by the simulations that used hydraulic parameters derived from in situ measured retention curves, E SHP , we evaluated the simulated daily evaporation rates only for E SHP . The overall difference between daily E m and E SHP is quantified by the root mean square difference [RMSE = (RMSD 2 s + RMSD 2 u ) 0.5 ].
In Figure 7, the scatter plots for E SHP show us frequent under-and overestimation of evaporation. In 2016, the overestimation of the cumulative E m in the DD and SB soils by E SHP is reflected in the large RMSD s . The RMSD s is small for 2018 when cumulative E m was predicted well in all soils. The nonstructural error between the daily E m and E SHP varied between 0.41 mm d −1 in 2018 and 0.69 mm d −1 in 2017. This error represents the day-to-day error between the model and measurements that cannot be removed after a linear transformation of the simulation results. It is considerably larger than the day-to-day "unexplained variation" of the daily evaporation rates among lysimeters of the same soil, (0.09 mm d −1 , Table 5). This indicates that the nonstructural error must be attributed to model uncertainty or errors rather than to measurement errors. In 2017, the good fit of E cum is achieved due to under-and overestimations of E, which are canceling out each other. Compared with 2016 and 2018, 2017 has more days where E SHP < E m . A closer inspection revealed that the period in 2017 when E SHP < E m was before the pesticide application on 1 September, and during this period, the soil surface was not completely bare and the measured water loss from the lysimeter could represent transpiration by weeds.
We can identify three main cases for daily E rates: (a) E SHP > E m during Phase II as seen in 2018 (Figure 8a), (b) E SHP < E m during Phase II (August 2018) (Figure 8b), and (c) E pot = E m = E SHP , which is mostly the case in early Phase I (Figure 8d). This leads to the following conjectures: at an intermediate soil wetness (θ = 0.15-0.3), E m decreases to <E pot faster than E SHP does (Figure 8a). This faster decrease is due to an overestimation of the unsaturated conductivity for an intermediate wetness or due to E pot overestimating Phase I evaporation. Under dry conditions, vapor transport and liquid film flow play important roles but are not taken into account in the hydraulic conductivity function used Figure 8b), which may lead to an underestimation of E m . The role of vapor transport during dry conditions was also apparent from the course of the diurnal E m rate, which followed the typical course of the E pot but at a lower rate (Figure 8a, b). Another indication for vapor flow is F I G U R E 7 Scatter plots of measured daily (E m ) versus simulated (E SHP ) using soil hydraulic parameters estimated from in situ measured retention curves. The systemic (RMSD s ) and unsystemic (RMSD u ) RMSD values are mean values for the Bad Lauchstädt (BL), Dedelow (DD) Sauerbach (SB), and Selhausen (SE) soil in each year; individual values for each soil can be found in Annex 4 (supplemental material). In 2017, dates of high underestimation are not considered for calculating RMSD s , RMSD u , and R 2 , as completely bare soil conditions are not given in this timeframe F I G U R E 8 Examples highlighting three differences between measured and simulated evaporation: (a) overestimation of evaporation by the Richards equation during Phase II, (b) underestimation of measured evaporation rates during Phase II, (c) measured evaporation rate underestimated by the model during nighttime due to an underestimation of the potential evaporation, and (d) measured and simulate evaporation equal to the potential evaporation. Shaded areas show the range in evaporation rates of lysimeter replicates the increase of daily E m when E pot increases during Phase II (Figure 8b). The diurnal course of E m and the increase of E m with E pot during Phase II cannot be reproduced by a model that only considers liquid flow, uses a pressure head threshold upper boundary condition, and does not account for an increase of vapor transport due to higher vapor concentrations and concentration gradients in the surface soil layer when the soil surface heats up .
Having a closer look at diurnal dynamics, we can also identify a fourth case relevant for hourly E rates: (d) E m > E pot and E SHP ≤ E m by night (Figure 8a, c), which is especially relevant for wet soils as in 2017. Large differences were especially visible during nighttime, when E m is mainly driven by wind speed (Figure 8c); Groh, Pütz, Gerke, Vanderborght, & Vereecken, 2019). Recent investigations highlight that water losses during night are of relevance at the global scale (Padrón, Gudmundsson, Michel, & Seneviratne, 2020). Also, non-rainfall water such as dew might play a role for larger E m , which can make up 2-7 mm per month of E on grassland lysimeters (Brunke, Groh, Vanderborght, Vereecken, & Pütz, 2019).

DISCUSSION AND CONCLUSION
During all three observation periods, daily potential evaporation (E pot ) was never reached by daily measured evaporation (E m ), even under wet conditions as in 2017. Even though evaporation was water limited, no differences in E m were observed between different soils. This does not agree with findings of Merlin et al. (2016) and Tolk, Evett, and Schwartz (2015), who found clear effects of soil texture on evaporation fluxes and of the soil type on evaporation rate during Phase II. We can assume that a lack of differences in E is due to small differences in soil texture and soil hydraulic properties in our lysimeters in combination with the impact of a disturbed surface layer. Additionally, it must be considered that E pot is much smaller in this study than experienced by Tolk et al. (2015). Three model approaches with different SHPs were compared with each other. Using SHPs that were estimated from soil texture using the ROSETTA PTFs resulted in a considerable overestimation of the measured evaporation, especially for the silty loam soils (BL, SB, and SE). This suggests that the hydraulic conductivity of the fine soils was overestimated.
For one soil, SHPs were derived from surface measurements of soil moisture using inverse modeling (Dimitrov, 2015), and it was observed that the disturbance of a surface layer by, for instance, tillage affected the hydraulic properties of the surface layer. Loosening the soil surface disrupts the connection of the soil surface via fine capillaries with soil water in deeper soil layers. This reduces the amount of water that can be evaporated from a wet soil at a potential rate. Under unsaturated conditions, the hydraulic conductiv-ity of a loosened aggregated soil layer can be much smaller than that of a less structured soil so that the loosened aggregated soil behaves as a soil with a coarser texture. Not considering this effect of soil loosening or soil structure on unsaturated soil properties could lead to an overestimation of E from finer textured soils, as was also demonstrated in a simulation study by Schlüter et al. (2012). In addition, the thickness of the loosened layer could have an effect on the evaporation, but this was not investigated further in our study, since we did not have information about this surface layer thickness. Assouline et al. (2014); Li, Vanderborght, and Smits (2020); and Shokri et al. (2010) discussed the effect of soil layering on evaporation. They considered layers of different texture and quantified the effect of soil layering in terms of the unsaturated hydraulic properties of the soil layers, (i.e., the characteristic length; Lehmann et al., 2008). These concepts could also be applied to layers with a different soil structure. Implementing a surface layer with different SHPs because of a change in soil structure due to soil loosening improved the simulation results, but except for 2017, E tilled was larger than E m .
A second set of SHP parameters were derived from in situ measurements of water content and soil water pressure heads. Even though the saturated hydraulic conductivity was still estimated from soil texture, using in situ measured water retention curves improved the estimation of the evaporation rates considerably. The in situ measured retention curves suggested a wider pore size distribution than the retention curves that were derived from PTFs. This resulted in lower estimates of the unsaturated hydraulic conductivity using the capillary model of Mualem (van Genuchten, 1980) and consequently in a lower simulated evaporation and a better agreement with the measured evaporation.
The simulated evaporation rates were very sensitive to the SHP used. In analogy to Lehmann et al. (2008), we derived a characteristic length of evaporation (L c ) from the soil hydraulic soil properties and found that L c can be used to compare evaporation that is predicted for different SHPs. We also confirmed the conclusion by Lehmann et al. (2020) that SHPs derived from soil texture using PTFs may have unrealistically large L c values for soils with a medium soil texture and consequently lead to overestimations of evaporation losses from these soils.
The Richards equation considers only liquid flow in the soil and uses a pressure head boundary at the soil surface to simulate the evaporation when the water pressure head reaches a threshold value. The Richards equation simulated bare soil lysimeter evaporation losses quite well when appropriate SHPs were used, including how they were reduced compared with the potential evaporation losses, estimated with the FAO-PM method scaled by an empirical coefficient . The use of a threshold boundary condition in the Richards equation decouples simulated evaporation dynamics from the dynamics of evaporative forcing and heat fluxes during Phase II and forces the simulated evaporation rate to decrease monotonously over time when the soil dries out . As was already shown using numerical simulations (Assouline et al., 2013;Fetzer et al., 2017) and is now further confirmed by measurements, this approach it is not able to reproduce diurnal dynamics when vapor transport in the porous medium is more important than liquid flow and vapor concentration in the soil increases during a day due to an increase in soil temperature. Neither can an increase in daily evaporation during Phase II, due to for instance a higher radiative forcing, be reproduced by this model. Other boundary conditions (e.g., using resistance terms that depend on surface soil moisture; Vandegriend & Owe, 1994) could be implemented to represent diurnal evaporation dynamics. Especially in loosened surface layers with large interaggregate pores, advective air flow may enhance vapor transport (Farrell, Greacen, & Gurr, 1966;Kimball & Lemon, 1971;Ishihara, Shimojima, & Harada, 1992;Maier et al., 2012;Scotter & Raats, 1969). Further research is needed to evaluate and quantify the impact of soil surface management practices, which influence soil structure and surface roughness and which aim at reducing soil evaporation losses.
Since the soil evaporation was in general smaller than E pot , except during some nights, we could not test the FAO-PM approach to calculate the potential evaporation of a wet soil surface. An important assumption in the PM is that the soil heat flux at the soil surface can be approximated as a fraction of the net radiation. When fluxes are averaged or accumulated over a longer time, this approximation may provide quite accurate results. However, over shorter times, when soil heat fluxes can make an important contribution in the surface energy balance, a more accurate estimation of the soil heat flux may be required-for instance, to obtain better estimates of the variation in daily evaporation. In this case, if soil heat flux is of concern and a better estimate for diurnal variations in potential evaporation is needed, models as developed by Evett et al. (2012) can give estimates that are more accurate.

A C K N O W L E D G M E N T S
This research was funded by the German Science Foundation (DFG) Project VA 351/15-1 and supported by funding from the Helmholtz Association and the Federal Ministry of Education and Research (BMBF) in the framework of TERENO (TERrestrial Environmental Observatories). Thanks to Werner Küpper for the technical support and providing the data.

D A T A AVA I L A B I L I T Y S T A T E M E N T
Data of the corresponding lysimeter at Selhausen are available at TERENO data portal (http://teodoor.icg.kfajuelich.de/ddp/index.jsp; Lysimeter: SE_Y_031-SE_Y_036 and SE_Y_041 -SE_Y_046).