Studying Unimodal, Bimodal, PDI and Bimodal-PDI Variants of Multiple Soil Water Retention Models: I. Direct Model Fit Using the Extended Evaporation and Dewpoint Methods

This study focuses on the reliable parametrization of the full Soil Water Retention Curve (SWRC) from saturation to oven-dryness using high resolution but limited range measured water retention data by the Hydraulic Property Analyzer (HYPROP) system. We studied the performance of five unimodal water retention models including the Brooks and Corey model (BC model), the Fredlund and Xing model (FX model), the Kosugi model (K model), the van Genuchten constrained model with four free parameters (VG model), and the van Genuchten unconstrained model with five free parameters (VGm model). In addition, eleven alternative expressions including Peters–Durner–Iden (PDI), bimodal, and bimodal-PDI variants of the original models were evaluated. We used a data set consisting of 94 soil samples from Turkey and the United States with high-resolution measured data (a total of 9264 measured water retention data pairs) mainly via the HYPROP system and supplemented for some samples with measured dry-end data using the WP4C instrument. Among unimodal expressions, the FX and the K models with the Mean Absolute Error (MAE) values equal to 0.005 cm3 cm−3 and 0.015 cm3 cm−3 have the highest and the lowest accuracy, respectively. Overall, the alternative variants provided a better fit than the unimodal expressions. The unimodal models, except for the FX model, fail to provide reliable dry-end estimations using HYPROP data (average MAE: 0.041 cm3 cm−3, average r: 0.52). Our results suggested that only models that account for the zero water content at the oven dryness and properly shift from the middle range to dry-end (i.e., the FX model and PDI variants) can adequately represent the full SWRC using typical data obtained via the HYPROP system.

Various models have been developed over the years to parametrize the SWRC. These models provide a continuous representation of the water content and soil tension relationship over the wet, intermediate and dry parts of the curve where typically sparse or even no measured data are available. Moreover, the soil hydraulic equations are often implemented in hydrologic models to characterize saturated/unsaturated water flow of different soil types and soil horizons and to describe the spatiotemporal heterogeneity of soil hydraulic properties across the scales from subfields, to fields, up to the landscape and beyond. However, the most widely used SWRC models such as models developed by Brooks and Corey [7] and van Genuchten [15] may not adequately describe the complete SWRC for all soils since they are unimodal and assume a residual water content for high dry-end tensions.
Undisturbed soils may have pore systems that are different from the approximately normally distributed unimodal type, which makes widely used unimodal water retention models inadequate for these heterogeneous pore systems. For such soils, multimodal retention functions can provide a better description of the SWRC. To address this shortcoming, Durner [27] suggested superimposing unimodal SWRCs of the van Genuchten model [15]. Later, Romano et al. [28] proposed a bimodal water retention function by conceptually dividing the pore domain into textural and structural parts. They assumed a lognormal pore size distribution function for each component and formulated the bimodal function as the weighted sums of the two components.
The widely used SWRC models typically consider residual water content for very high tensions. The residual water content has been assumed as the water held by adsorptive forces [29], the air dry water content [30], the water content at which water movement stops [31] and the water content at a very large tension head such as the permanent wilting point [15]. In practice, generally, the residual water content simply acts as a fitting parameter at which the slope of the SWRC becomes zero toward the SWRC dry-end. However, this assumption is incorrect because for the very dry part of the SWRC water content eventually becomes zero [14]. To overcome this issue, the Peters-Durner-Iden (PDI) SWRC model variant [10,14] was proposed which reaches zero water content at oven dryness through a linear water content reduction in the dry range of the SWRC against the logarithmic transformation of the soil tension, pF [10]. The main advantages of PDI variants are that they do not require more parameters than the original SWRCs and showed a great match to experimental data [32].
The parameters of SWRC's models are obtained by fitting the equations to the measured water retention data or are estimated via PTFs. The performance evaluation of the soil hydraulic models is typically done using data sets with a limited number of measured data points per sample and often sparse or no data on the SWRC dry-end [33,34] due to the difficulties in measuring the SWRC at matric suctions larger than 1500 kPa. This shortcoming exists mainly because the traditional standard equilibrium approach (i.e., the sandbox apparatus, the sand/kaolin box, and the pressure plate extractor) typically produces a few discrete measured water retention data pairs per soil sample. Consequently, much attention has been paid to extending the predictive capability of the existing SWRC models to the dry region from limited measurements in the wet region [35,36].
An alternative laboratory approach for measurement of the soil water retention is the extended evaporation method [2,3]. The HYPROP system (Hydraulic Property Analyzer, METER Group, Inc., Pullman, WA, USA) is a commercial laboratory instrument that works based on the extended evaporation method. HYPROP generates high-resolution data typically close to 100 water retention data points in the 0 to 100 kPa soil tension range plus an extra point close to the wilting point by considering the air-entry value as an additional measurement [37]. In recent years, the HYPROP system has been increasingly used and evaluated around the world [4,[38][39][40][41][42][43][44] and promising results have been reported. For instance, Bezerra-Coelho et al. [38] evaluated the HYPROP system as applied to the van Genuchten [15] soil hydraulic functions for a wide range of soil textures using the HYDRUS-1D software package [45]. They reported an extremely well HYPROP performance, especially for the SWRC. In another study, Zhuang et al. [44] reported an excellent agreement between the soil hydraulic properties measured by multistep flux and the hanging water column experiments against the HYPROP system.
While this paper focuses on the direct-fit of soil water retention data to SWRC models, the second part of this study [17] concentrates on the performance of the respective parametric PTFs and summarizes the overall result. HYPROP data, although high resolution, are almost entirely in the wet-end and the middle part of the SWRC since HYPROP cannot generate dry-end data. It is crucial for the soil hydraulic models to provide an accurate estimation of soil water content for the entire range of the SWRC from saturation to over dryness in particular in arid regions. Accurate estimation of water content on the dry part of the curve is essential for understanding many critical processes related to soil including wind erosion, plant growth, methane oxidation, water conservation, biological and microbial activities, and nitrogen mineralization [8,[46][47][48][49]. Consequently, it is crucial to determine how HYPROP limited range of measurement impacts the performance of the soil water retention models describing the entire SWRC from saturation to over-dryness, the main objective of this study. The specific objectives of this study were to: (i) investigate the performance of 16 unimodal, bimodal, PDI-unimodal and Bimodal-PDI SWRC models at different parts of the SWRC and (ii) investigate the extrapolation ability of the models beyond the HYPROP measurement range using additional dry-end measured data via the WP4C instrument (METER Group, Inc., Pullman, WA, USA). Figure 1 illustrates the main parts of the HYPROP system (accuracy: 0.15 kPa (0 to 82 kPa); resolution: 0.001 kPa) as well as four typical phases of an optimum measurement campaign.

Extended Evaporation Approach Using the HYPROP System
Water 2020, 12, x FOR PEER REVIEW 3 of 19 between the soil hydraulic properties measured by multistep flux and the hanging water column experiments against the HYPROP system. While this paper focuses on the direct-fit of soil water retention data to SWRC models, the second part of this study [17] concentrates on the performance of the respective parametric PTFs and summarizes the overall result. HYPROP data, although high resolution, are almost entirely in the wet-end and the middle part of the SWRC since HYPROP cannot generate dry-end data. It is crucial for the soil hydraulic models to provide an accurate estimation of soil water content for the entire range of the SWRC from saturation to over dryness in particular in arid regions. Accurate estimation of water content on the dry part of the curve is essential for understanding many critical processes related to soil including wind erosion, plant growth, methane oxidation, water conservation, biological and microbial activities, and nitrogen mineralization [8,[46][47][48][49]. Consequently, it is crucial to determine how HYPROP limited range of measurement impacts the performance of the soil water retention models describing the entire SWRC from saturation to over-dryness, the main objective of this study. The specific objectives of this study were to: (i) investigate the performance of 16 unimodal, bimodal, PDI-unimodal and Bimodal-PDI SWRC models at different parts of the SWRC and (ii) investigate the extrapolation ability of the models beyond the HYPROP measurement range using additional dry-end measured data via the WP4C instrument (METER Group, Inc., Pullman, WA, USA).  Hydraulic Property Analyzer (HYPROP) system components (right) and the four typical phases (P1 to P4) of a measurement campaign (left) for the optimal measuring curve including the regular measurement range phase (P1), the boiling delay phase (P2), the cavitation phase (P3) and the air entry phase (P4), adapted from HYPROP user manual [50].

Extended Evaporation Approach Using the HYPROP System
During phase one, the tensiometers are in the regular measurement range wherein the tension values increase steadily until the boiling point of the water is reached. If the system is filled completely air-free (an ideal case) the boiling delay phase (phase two) will occur in which the tension increases to above the ambient air pressure (up to the boiling delay area). During phases one and two, the measured soil tension reflects the real soil matric potential of the soil surrounding tensiometers. Phase three, the "Cavitation phase", happens when water vapor is generated in the tensiometers, which causes a sharp reduction in the tension value down to the boiling point. During this phase, the tensiometers no longer provide readings that are representative of the real soil tension. After this phase, the tension value decreases only slightly until a second sharp reduction in tension Hydraulic Property Analyzer (HYPROP) system components (right) and the four typical phases (P1 to P4) of a measurement campaign (left) for the optimal measuring curve including the regular measurement range phase (P1), the boiling delay phase (P2), the cavitation phase (P3) and the air entry phase (P4), adapted from HYPROP user manual [50].
During phase one, the tensiometers are in the regular measurement range wherein the tension values increase steadily until the boiling point of the water is reached. If the system is filled completely air-free (an ideal case) the boiling delay phase (phase two) will occur in which the tension increases to above the ambient air pressure (up to the boiling delay area). During phases one and two, the measured soil tension reflects the real soil matric potential of the soil surrounding tensiometers. Phase three, the "Cavitation phase", happens when water vapor is generated in the tensiometers, which causes a sharp reduction in the tension value down to the boiling point. During this phase, the tensiometers no longer provide readings that are representative of the real soil tension. After this phase, the tension value decreases only slightly until a second sharp reduction in tension values occurs, this time down to zero, as air enters the ceramic tip of the tensiometer because the tension in the soil surrounding tensiometers becomes greater than the air-entry pressure of the ceramic material (phase four). The air entry point depends on the material characteristic of the ceramic tip and equals to approximately 880 kPa for the HYPROP system's tensiometers. HYPROP works based on the extended evaporation method, meaning the initiation of the air-entry phase can be used as an extra soil matric potential measurement point which further extends the range of the readings to the dry part of the SWRC.

Dew Point Technique Using WP4C System
WP4C uses the chilled-mirror dew-point technique to measure the tension of a soil sample. Figure 2 illustrates the main parts of the WP4C system. Inside the chamber, a 15 mL sample cup is sealed against a sensor block containing an internal fan, a dew-point sensor, a temperature sensor, and an infrared thermometer. The air within the sample chamber is circulated by the fan, which helps to speed up the equilibrium time and to also control the conductance of the dew-point sensor's boundary layer. The dew-point and the infrared thermometer sensors simultaneously measure the dew-point temperature of the air and the sample temperature, respectively. The sample block temperature is monitored and stabilized using an internal thermo-electrical module. The equilibrium of the liquid phase water of the sample with the vapor phase water in the headspace of the sealed chamber is reached when the water potential of the sample becomes the same as the water potential of the air. WP4C directs a beam of light onto the mirror and uses the photoelectric cell to precisely determine the change in reflectance when condensation occurs on the mirror. At the condensation point, the temperature is recorded using the thermocouple attached to the mirror. At equilibrium, the vapor pressure of the air in the headspace is determined as the saturation vapor pressure at dew-point temperature and the saturation vapor pressure is computed from sample temperature. Finally, the water potential of the sample is calculated using the headspace and saturation vapor pressure values. values occurs, this time down to zero, as air enters the ceramic tip of the tensiometer because the tension in the soil surrounding tensiometers becomes greater than the air-entry pressure of the ceramic material (phase four). The air entry point depends on the material characteristic of the ceramic tip and equals to approximately 880 kPa for the HYPROP system's tensiometers. HYPROP works based on the extended evaporation method, meaning the initiation of the air-entry phase can be used as an extra soil matric potential measurement point which further extends the range of the readings to the dry part of the SWRC.

Dew Point Technique Using WP4C System
WP4C uses the chilled-mirror dew-point technique to measure the tension of a soil sample. Figure 2 illustrates the main parts of the WP4C system. Inside the chamber, a 15 mL sample cup is sealed against a sensor block containing an internal fan, a dew-point sensor, a temperature sensor, and an infrared thermometer. The air within the sample chamber is circulated by the fan, which helps to speed up the equilibrium time and to also control the conductance of the dew-point sensor's boundary layer. The dew-point and the infrared thermometer sensors simultaneously measure the dew-point temperature of the air and the sample temperature, respectively. The sample block temperature is monitored and stabilized using an internal thermo-electrical module. The equilibrium of the liquid phase water of the sample with the vapor phase water in the headspace of the sealed chamber is reached when the water potential of the sample becomes the same as the water potential of the air. WP4C directs a beam of light onto the mirror and uses the photoelectric cell to precisely determine the change in reflectance when condensation occurs on the mirror. At the condensation point, the temperature is recorded using the thermocouple attached to the mirror. At equilibrium, the vapor pressure of the air in the headspace is determined as the saturation vapor pressure at dewpoint temperature and the saturation vapor pressure is computed from sample temperature. Finally, the water potential of the sample is calculated using the headspace and saturation vapor pressure values.

Data Collection and Laboratory Analyses
A total of 94 soil samples were used in this study including 15 samples collected from multiple irrigation experimental sites in southern and central California, USA and 79 samples collected from the areas surrounding cities of Ankara and Anamur located in central and southern Turkey, respectively. HYPROP measurement campaigns were completed at the Technical University of

Data Collection and Laboratory Analyses
A total of 94 soil samples were used in this study including 15 samples collected from multiple irrigation experimental sites in southern and central California, USA and 79 samples collected from the areas surrounding cities of Ankara and Anamur located in central and southern Turkey, respectively. HYPROP measurement campaigns were completed at the Technical University of Braunschweig in Germany and the University of California, Riverside for the soils from Turkey and the United States, respectively. The soil texture was determined by the hydrometer method [52] and the PARIO soil particle analyzer (METER Group, Inc., Pullman, WA, USA) for the Turkish and Unites States samples, respectively. More details about the soil data set and the laboratory procedures are outlined in Haghverdi et al. [1,16]. Table 1 summarizes the characteristics of the soils used in this study. For the United States samples, 250 cm 3 stainless-steel cylinders (inside diameter: 8 cm, height: 5 cm) were used to collect undisturbed samples. All samples were brought up to saturation and then two vertically aligned tensiometers were installed in two holes made via a small auger, such that the tensiometers' tips were positioned 1.25 cm below and above the center of the cylinder (2.5 cm). For the Turkish samples, the air-dried disturbed sieved samples (10 mesh sieve; <2.0 mm) were packed and brought up to saturation under vacuum in the same size cylinders to dry BDs close to the field condition (repacking was done in multiple steps to ensure consistent BD of the repacked samples).
Only the upper part of the samples was open to the atmosphere for evaporation. For each sample, depending on the soil type, the measurement campaign lasted between 4 and 14 days (on average, ten days per sample for the entire data set). Throughout the measurement campaign, the soil tensions were automatically monitored using HYPROP-VIEW software at the two depths using the tensiometers. The integral of the water content distribution over the entire column divided by the height of the sample was considered by HYPROP as the mean water content of the sample. The measurement frequencies of the tensiometers were set to be higher at the beginning of the experiment (approximately once a minute for the first hour) and further apart afterward (every ten minutes). The weights of the United States samples were continuously recorded, while the Turkish samples were weighed only twice a day following the single balance approach as outlined in the HYPROP operation manual [50].
Immediately after the HYPROP measurement campaign was finished, the WP4C Dew Point PotentioMeter instrument (METER Group, Inc., Pullman, WA, USA) was used to get data at the SWRC dry-end for the US samples. A total of 4-5 subsamples (each roughly 7 cm 3 ) were sliced off the original undisturbed sample and stored in caped small sample cups (15 cm 3 capacity). The subsamples were taken from the top, middle and bottom of the original sample with different water contents. The soil tension of each subsample was measured using WP4C and the weight of the sample was immediately recorded afterward. Finally, the subsamples and the remaining soil material from HYPROP were put in the oven to get the oven-dry weight of the sample. The water retention data measured by the HYPROP and WP4C are shown in Figure 3. The HYPROP measurement campaign yielded high-resolution data (on average, 98 data points per soil sample) between saturation and up close to the wilting point with an average pF of 1.76. On average WP4C measurements had a pF value of 4.68.

Soil Hydraulic Models
A total of 16 models were evaluated in this study. In the following, the models are grouped into four categories including original unimodal expressions, PDI variants, bimodal variants, and bimodal-PDI variants. The unimodal variants were fitted to all the data. The bimodal variants were only fitted to the 15 samples from the United States as these were undisturbed cores and therefore retained structural and textural pores while disturbed samples from Turkey likely mainly reflect textural pores.

Original Unimodal Expressions
Five original unimodal water retention models were evaluated in this study including the Brooks and Corey model (BC model; [7]), the Fredlund and Xing model (FX model; [9]), the Kosugi model (K model; [12]), the van Genuchten constrained model with four free parameters (VG model; [15]), and the van Genuchten unconstrained model with five free parameters (VGm model; [15]).
The BC model is one of the first models developed to parametrize the SWRC consisting of a constant straight line and a power function joining around the bubbling pressure. The BC model is represented as: where θ(h) is the volumetric water content at soil tension h, α (1/cm) is the inverse of the air entry value (the bubbling pressure), λ (−) is the pore size distribution index (which affects the slope of the curve) and θr and θs are the residual and saturated soil water contents, respectively. The FX model assumes that the shape of the SWRC is dependent upon the pore size distribution of the soil and is given by: and ℎ = 1 − ln 1 + ℎ ℎ ⁄ ln 1 + ℎ ℎ ⁄ (4)

Soil Hydraulic Models
A total of 16 models were evaluated in this study. In the following, the models are grouped into four categories including original unimodal expressions, PDI variants, bimodal variants, and bimodal-PDI variants. The unimodal variants were fitted to all the data. The bimodal variants were only fitted to the 15 samples from the United States as these were undisturbed cores and therefore retained structural and textural pores while disturbed samples from Turkey likely mainly reflect textural pores.

Original Unimodal Expressions
Five original unimodal water retention models were evaluated in this study including the Brooks and Corey model (BC model; [7]), the Fredlund and Xing model (FX model; [9]), the Kosugi model (K model; [12]), the van Genuchten constrained model with four free parameters (VG model; [15]), and the van Genuchten unconstrained model with five free parameters (VG m model; [15]).
The BC model is one of the first models developed to parametrize the SWRC consisting of a constant straight line and a power function joining around the bubbling pressure. The BC model is represented as: where θ(h) is the volumetric water content at soil tension h, α (1/cm) is the inverse of the air entry value (the bubbling pressure), λ (−) is the pore size distribution index (which affects the slope of the curve) and θ r and θ s are the residual and saturated soil water contents, respectively. The FX model assumes that the shape of the SWRC is dependent upon the pore size distribution of the soil and is given by: with and where α, n, and m are the SWRC shape parameters, h r is the tension corresponding to θ r and h 0 is the soil tension at zero water content (in this study was set to 10 6.8 cm which is the suction at oven dryness for 105 • C, [53]) and e is the Euler number.
The K model is a modified version of the Kosugi [11] model that is derived by applying a three-parameter lognormal distribution law to the soil pore radius distribution function: where "erfc" is the complementary error function, σ is the standard deviation of the log-transformed pore-size distribution density function, and h m is the suction corresponding to the median pore radius. The VG model is considered to be the most widely used parameterization of the SWRC. It is represented as: where α and n are the curve shape parameters and other parameters as previously defined. The VG m model with m as an additional shape parameter has five free parameters:

PDI Expressions
The PDI variants were derived for all the unimodal original expressions except for the BC model. The general form of the PDI model [10,14] consists of a capillary retention term, θ cap (h) and an adsorptive retention term, θ ad (h), and is given as: where S cap and S ad are the capillary and the water adsorption saturation functions and θ r is the maximum water content for the water adsorption.
To guarantee that the water content reaches zero at h = h 0 , the S cap is substituted by scaled versions of the original functions: where Γ(h) represents basic saturation functions and Γ 0 is the basic function at h = h 0 . The basic classic saturation functions for the abovementioned unimodal expressions are: Γ(h) = ln e + (αh) n −m for the FX model (10) for the VG m model (13) The water adsorption saturation function is given as [10]: where x a and x 0 are pF values at suctions equal to h a and h 0 , respectively, h a is the suction at air entry for the adsorptive retention and b is the shape parameter which is given by: for the FX, the VG, and the VG m models (15)

Bimodal Expressions
The bimodal expressions for each unimodal model were the unscaled weighted sum of the two unimodal subfunctions without adsorption (S ad = 1): where w i is the weighting factor for the subfunction i, subject to 0 < w i < 1 and Σ w i = 1. The Γ(h) is calculated using Equations (10)-(13).

Bimodal-PDI Expressions
The bimodal-PDI expression for each model was the scaled weighted sum of the two unimodal subfunctions with adsorption considered: The Γ is calculated using Equations (10)- (13). The S ad is calculated using Equation (14) for which the shape parameter b is calculated using Equations (15) and (16). The shape parameter is calculated only for the "coarsest" subfunction which is the subfunction with the lowest h m value for the K model or the highest α value for the FX, the VG and the VG m models.

Model Parametrization
HYPROP-FIT software was used to fit the soil hydraulic models to the measured water retention data. HYPROP-FIT works based on a revised version of SHYPFIT2.0 (Soil Hydraulic Properties Fitting, [54]). The fitting was simultaneously accomplished for both water retention and hydraulic conductivity models through minimizing the sum of weighted squared residuals between model estimations and measured data pairs (this study only focuses on the SWRC). The shuffled complex evolution algorithm [55] is used for parameters estimation, which guarantees that the best parameter combination (i.e., global minimum in the multidimensional parameter space) for the appropriate model combination is found. Table 2 summarizes the default bounds for the model parameters imposed during the fitting process by the HYPROP-FIT software. The bounds are chosen to provide a high degree of flexibility as well as physical consistency.   [7], FX: Fredlund and Xing [9], K: Kosugi [12], VG and VG m : van Genuchten [15] constrained and unconstrained unimodal soil water retention models. PDI and b denote Peters-Durner-Iden [10,14] and bimodal variants of the models, respectively.

Evaluation Statistics
The mean absolute error (MAE, Equation (19)), the root mean square error (RMSE, Equation (20)) and the correlation coefficient (r, Equation (21)) were used as the main statistics to assess the performance of the models. In addition, Mean Bias Error (MBE, Equation (22)) was calculated to determine the degree of systematic model bias, with positive and negative values indicating net overand under-estimate of the data by the fitted model.
where E and M are the fitted and the measured soil water content values (cm 3 cm −3 ), E and M are the mean fitted and the mean measured water content values and n is the number of the measured water retention points for all samples considered in each analysis. Table 3 summarizes the overall performance of the unimodal and PDI models for the entire data set. Figure 4 shows the scatter plots of measured versus fitted water retention data for all the models and both datasets. The correlation coefficient values are close to 1 for all models, while the RMSE values vary between 0.007 cm 3 cm −3 and 0.020 cm 3 cm −3 . The MBE values are close to zero indicating systematic bias for none of the models. Among unimodal expressions, the FX and the K models with the MAE values equal to 0.005 cm 3 cm −3 and 0.015 cm 3 cm −3 have the highest and the lowest accuracy, respectively. The PDI variants outperformed the unimodal expressions for VG, VG m and K models. This improvement was more pronounced for clay and clay loam samples compared to loam and sandy loam samples when performance statistics were separately calculated (data not shown here) for the textures with highest number of samples (i.e., clay, clay loam, loam, sandy loam). When MAE for each sample is separately calculated (data not shown here), the FX, FX-PDI and VG m -PDI model exhibit the best fit to the measured data for 47%, 23% and 22% of our soil samples, respectively. The K model followed by the BC model show the worst fit (highest MAE values) for 73% and 26% of the samples, respectively. The lower performance of the K and BC models were observed across clay, clay loam, loam, and sandy loam textures (data not shown here).   [7], FX: Fredlund and Xing [9], K: Kosugi [12], VG and VGm: van Genuchten [15] constrained and unconstrained unimodal soil water retention models. PDI and B denote Peters-Durner-Iden [10,14] and bimodal variants of the models, respectively. Table 5 shows the dry-end performance of the 16 models for the US samples with and without the additional WP4C data included in the fitting process. With the dry-end WP4C data included, the MAE range between 0.005 cm 3 cm −3 and 0.023 cm 3 cm −3 and r values vary between 0.57 and 0.97. Without the dry-end data, the MAE values vary from 0.008 cm 3 cm −3 to 0.078 cm 3 cm −3 and the correlation coefficients range between 0.42 and 0.96. The MBE values range from −0.001 to 0.018 with WP4C data and from 0.001 to 0.078 without WP4C data included in the fitting process.

Extrapolation Capability beyond HYPROP Range
Without the dry-end data, there is a noticeable increase in the MAE and MBE values for the unimodal expressions except FX (more pronounced for the K, the VG and the VGm models), indicating a low performance mainly due to overestimation of the dry-end data. The FX model shows comparable performance with and without the dry-end data included in the fitting process with the MAE values equal to 0.009 cm 3 cm −3 and 0.011 cm 3 cm −3 , respectively. All variants of the FX model show similar performance. Without dry-end data, for the K, the VG and the VGm models, the MAE values of the PDI, the bimodal, and the bimodal-PDI variants were substantially lower than that of the unimodal expressions. The performance improvement of these alternative variants over the original models is less pronounced with the dry-end data included in the fitting process. Except for the FX model, the none-PDI variants of the other models (i.e., the original expression as well as the bimodal variants) have an average r value of 0.56 showing a lower performance than the PDI variants with an average r value of 0.92.  [7], FX: Fredlund and Xing [9], K: Kosugi [12], VG and VG m : van Genuchten [15] constrained and unconstrained unimodal soil water retention models. PDI and B denote Peters-Durner-Iden [10,14] and bimodal variants of the models, respectively. Table 4 summarizes the overall performance of all 16 models for the USA soils. The bimodal and the bimodal-PDI variants show similar statistics. For all models, the ranking of the variants from the highest to the lowest performance is: (i) bimodal-PDI, (ii) bimodal, (iii) PDI and (iv) unimodal. No systematic bias is observed since the MBE values are close to zero for all models. Among all 16 models, the K model shows the lowest accuracy with the MAE and RMSE values of 0.009 cm 3 cm −3 and 0.012 cm 3 cm −3 , respectively. When MAE values for individual US soil samples are considered (data not shown here), the VG m -b-PDI and K model exhibit the best and worst fit for 47% and 93% of the samples, respectively.  [7], FX: Fredlund and Xing [9], K: Kosugi [12], VG and VG m : van Genuchten [15] constrained and unconstrained unimodal soil water retention models. PDI and b denote Peters-Durner-Iden [10,14] and bimodal variants of the models, respectively; 2 MBE: mean bias error (cm 3 cm −3 ), RMSE: root mean square error (cm 3 cm −3 ), MAE: mean absolute error (cm 3 cm −3 ), r: correlation coefficient. Table 5 shows the dry-end performance of the 16 models for the US samples with and without the additional WP4C data included in the fitting process. With the dry-end WP4C data included, the MAE range between 0.005 cm 3 cm −3 and 0.023 cm 3 cm −3 and r values vary between 0.57 and 0.97. Without the dry-end data, the MAE values vary from 0.008 cm 3 cm −3 to 0.078 cm 3 cm −3 and the correlation coefficients range between 0.42 and 0.96. The MBE values range from −0.001 to 0.018 with WP4C data and from 0.001 to 0.078 without WP4C data included in the fitting process. Without the dry-end data, there is a noticeable increase in the MAE and MBE values for the unimodal expressions except FX (more pronounced for the K, the VG and the VG m models), indicating a low performance mainly due to overestimation of the dry-end data. The FX model shows comparable performance with and without the dry-end data included in the fitting process with the MAE values equal to 0.009 cm 3 cm −3 and 0.011 cm 3 cm −3 , respectively. All variants of the FX model show similar performance. Without dry-end data, for the K, the VG and the VG m models, the MAE values of the PDI, the bimodal, and the bimodal-PDI variants were substantially lower than that of the unimodal expressions. The performance improvement of these alternative variants over the original models is less pronounced with the dry-end data included in the fitting process. Except for the FX model, the none-PDI variants of the other models (i.e., the original expression as well as the bimodal variants) have an average r value of 0.56 showing a lower performance than the PDI variants with an average r value of 0.92.  [7], FX: Fredlund and Xing [9], K: Kosugi [12], VG and VGm: van Genuchten [15] constrained and unconstrained unimodal models. PDI and b denote Peters-Durner-Iden [10,14] and bimodal variants, respectively.

Unimodal Expressions
Cornelis et al. [33] evaluated the performance of ten closed-form unimodal SWRC models including BC, VG, VG m , and K models and equations developed by Tani [56], Russo [57], Rossi and Nimmo [58], Kosugi [12,13], and Assouline et al. [6]. They found that the VG m model with a mean RMSE value of 0.008 cm 3 cm −3 showed the highest match with the measured water retention data providing the best fit for 67% of the soil samples. The BC, the VG and the K models in their study showed lower performance with mean RMSE values of 0.014 cm 3  The higher error near saturation for the BC model (MAE: 0.012 cm 3 cm −3 for pF < 2) is attributed to its discontinuous form, which causes a sharp corner where its two parts meet near saturation. The poor match of the BC model near-saturation was also reported by Cornelis et al. [33]. Our results also showed relatively high wet-end MAE values for the K model (MAE: 0.016 cm 3 cm −3 for pF < 2), the VG model (MAE: 0.012 cm 3 cm −3 for pF < 2) and less pronounced for the VG m model (MAE: 0.008 cm 3 cm −3 for pF < 2), a trend that was not reported by Cornelis et al. [33]. The differences between the results of the two studies might be related to the limited number of measured data near saturation in the data set used by Cornelis et al. [33] as opposed to the high-resolution wet-end data that was utilized in this study.
We realized when no dry-end data is included in the fitting process, except for the FX model, the unimodal models fail to provide reliable dry-end estimations (average MAE: 0.041 cm 3 cm −3 , average r: 0.52), a drawback also highlighted by Cornelis et al. [33], hence their extrapolation beyond the driest measured water retention point is not suggested ( Table 5). Our results show that when additional dry-end measurements are included in the fitting process, the unimodal models perform better in the SWRC dry-end (average dry-end MAE improved from 0.041 cm 3 cm −3 to 0.014 cm 3 cm −3 ). We tested the feasibility of setting the residual water content to zero in the model data fitting scheme as a possible strategy to improve the dry-end estimation by the unimodal VG and VG m models. This strategy was ineffective, however, and caused a decrease in the overall performance of both models. For instance, the MAE values increased from 0.007 to 0.009 cm 3 cm −3 for the VG model and from 0.005 to 0.006 cm 3 cm −3 for the VG m model when the residual water content was set to zero. The FX model was the only unimodal expression that accounted for the zero water content at the oven dryness and accurately transitioned from the middle range to dry-end. A similar result was reported by Lu et al. [36] who studied the extrapolative capability of the FX model to estimate the complete SWRC for soils with various textures. They reported that when measured data pairs in the tension range of 0 to 1500 kPa are available, the fitted FX model provided a good agreement with experimental data from the complete SWRC. They, therefore, concluded that the need for the dry-end measurement can be eliminated with the FX model.

Bimodal and PDI Variants
The PDI and bimodal-PDI variants showed the best fit to the samples collected from Turkey (repacked) and the USA (undisturbed), respectively. The bimodal-PDI variants maintained high accuracy across a wide range from the saturation to the oven dryness with little to no fluctuations in their performances ( Figure 5). The higher performance of bimodal variants compared to original unimodal models for undisturbed USA samples can be partially attributed to their ability to account for both structural and textural pore domains. This improvement is also attributed to the fact that free fitting parameters of the bimodal variants are twice more than unimodal expressions.
The low correlation coefficients of K-b, VG-b, and VG m -b models (Table 5) reveal that despite their high number of free parameters, bimodal variants fail to accurately estimate the dry-end data using HYPROP data. The PDI variants, on the other hand, outperformed unimodal and bimodal variants in estimating the SWRC dry-end with average dry-end MAE values of 0.008 cm 3 cm −3 (WP4C included in the fitting process) and 0.014 cm 3 cm −3 (WP4C not included in the fitting process). Lu et al. [17] evaluated the performance of three models (developed by Fayer and Simmons, [59]; Webb [60]; Khlosi et al., [61]) describing the complete SWRC using data points in the pressure plate measurement range (i.e., tension less than 1500 kPa) on eight soils. The results of their study illustrated that when measured data pairs are available in the 0 to 1500-kPa range, all three models could provide a reliable fit to the complete SWRC range. Our results suggested that typical HYPROP data (i.e., high-resolution data from saturation to 100 kPa soil tension range and an additional point close to the wilting point) are sufficient for an accurate representation of the full SWRC (from saturation to oven-dryness) when PDI variants are fitted to the experimental data. Since all of our undisturbed samples were coarse-textured and the number of soil samples was relatively limited, we recommend further investigations to determine whether the same result holds for other soil types.

Conclusions
This study focused on the full parametrization of the SWRC using high resolution but limited range measured soil water retention data generated by the HYPROP system. A total of sixteen soil water retention models were evaluated in this study using a high-resolution data set (a total of 9264 measured water retention data pairs) obtained via the extended evaporation method (the HYPROP system) and supplemental WP4C data covering the SWRC dry-end (available for some samples). The data set contained soil samples from Turkey and the United States, representing multiple soil textural classes. The models studied included five original widely used unimodal expressions and eleven additional PDI, bimodal, and bimodal-PDI variants. Table A1 summarizes the average values of the fitted model parameters for the 16 models across textural classes. Overall, the alternative variants (i.e., PDI and bimodal expressions) outperformed the original unimodal expressions. In the absence of additional dry-end measured data by the WP4C instrument, most of the original unimodal expressions provided undefined and unreliable water content estimations at the SWRC dry-end hence should not be used for parametrization of the full SWRC using HYPROP system. Our result only recommends the application of the PDI variants and FX model for complete SWRC estimation using the HYPROP system since these models account for zero water content at oven dryness and accurately model the SWRC dry-end.