Studying Unimodal, Bimodal, PDI and Bimodal-PDI Variants of Multiple Soil Water Retention Models: II. Evaluation of Parametric Pedotransfer Functions Against Direct Fits

A high-resolution soil water retention data set (81 repacked soil samples with 7729 observations) measured by the HYPROP system was used to develop and evaluate the performance of regression parametric pedotransfer functions (PTFs). A total of sixteen soil hydraulic models were evaluated including five unimodal water retention expressions of Brooks and Corey (BC model), Fredlund and Xing (FX model), Kosugi (K model), van Genuchten with four free parameters (VG model) and van Genuchten with five free parameters (VGm model). In addition, eleven bimodal, Peters–Durner–Iden (PDI) and bimodal-PDI variants of the original expressions were studied. Six modeling scenarios (S1 to S6) were examined with different combinations of the following input predictors: soil texture (percentages of sand, silt and clay), soil bulk density, organic matter content, percent of stable aggregates and saturated water content (θs). Although a majority of the model parameters showed low correlations with basic soil properties, most of the parametric PTFs provided reasonable water content estimations. The VGm parametric PTF with an RMSE of 0.034 cm3 cm−3 was the best PTF when all input predictors were considered. When averaged across modeling scenarios, the PDI variant of the K model with an RMSE of 0.045 cm3 cm−3 showed the highest performance. The best performance of all models occurred at S6 when θs was considered as an additional input predictor. The second-best performance for 11 out of the 16 models belonged to S1 with soil textural components as the only inputs. Our results do not recommend the development of parametric PTFs using bimodal variants because of their poor performance, which is attributed to their high number of free parameters.


Introduction
The term pedotransfer function (PTF) was first proposed by Bouma and van Lanen [1] to describe regression models formulating relationships between soil properties. Following their proposition, a majority of the published research in this area has focused on estimating soil hydraulic properties, mainly on individual water retention points and the soil water retention curve (SWRC). In many practical applications, PTF-estimated soil hydraulic properties, from readily available basic properties, are used to bridge the gap when adequately measured data are not available. The application of PTFs is particularly important for the larger scales water flow and solute transport modeling studies for which direct measurement of soil hydraulic properties is impractical.
Point pedotransfer functions (Group 1 PTFs, [2]) are developed to estimate the soil water content at one or several predefined soil tensions (typically field capacity and permanent wilting point). Parametric PTFs (Group 2 PTFs, [2]) are the most widely used models for continuous estimation of the soil water content over a wide range of soil tensions. Parametric PTFs estimate parameters of a soil water retention equation, which is subsequently used to predict the SWRC. Haghverdi et al. [3] introduced an alternative non-parametric approach where instead of estimating parameters of soil hydraulic equations for a target soil, the k-nearest neighbor technique is used to find the k most similar soils in the PTF development dataset. The SWRC of the target soil is then predicted as the weighted average of the fitted curves for the k similar soils. Additionally, the pseudo-continuous PTF (PC-PTF, Haghverdi et al., [4,5]) is another alternative approach that uses pF (logarithmic transformation of the soil tension) as an extra input attribute. This extra input allows continuous prediction of the water content by an appropriate data mining method without using any soil water retention equations.
The performance of the parametric PTFs is typically evaluated using data sets with a limited number of soil water retention data per sample, measured via the traditional standard equilibrium approach (i.e., the sandbox apparatus, the sand/kaolin box and the pressure plate extractor). An alternative option for measuring soil water retention in the laboratory is the extended evaporation method (Schindler et al., [6,7]) via the HYPROP system (Hydraulic Property Analyzer, METER Group, Inc., Pullman, WA, USA). The HYPROP system adequately determines the hydraulic properties of the most soils [8,9]. It generates high-resolution data in the wet and intermediate parts of the WRCs, which provides a much better data set than the traditional standard equilibrium approach to determine the accuracy and reliability of the parametric PTFs.
According to Wösten et al. [10], the key factor for developing more accurate and reliable PTFs is to enhance the accuracy of the soil hydraulic properties measurements. A standardized method, however, that guarantees high-quality measurement of soil hydraulic properties has been noticeably lacking [11]. The evidence suggests that the HYPROP system could become the benchmark technique for the laboratory measurement of the soil hydraulic properties and the development of the next generation PTFs. This is suggested given the high-resolution water retention and hydraulic conductivity data generated by the instrument, its relatively fast measurement cycle and its adaptation by researchers from different parts of the world [12][13][14][15][16][17]. In a recent study, Haghverdi et al. [18] showed that artificial neural network-based PC-PTF successfully predicts the SWRC using high-resolution data collected via the HYPROP system.
Various water retention models have been developed over the years to parametrize the SWRC. However, a majority of parametric PTFs [19][20][21][22][23] have been developed using unimodal soil hydraulic equations and, in particular, the van Genuchten soil hydraulic model [24] while alternative models have been rarely explored [11]. The first part of this study [25] focuses on the direct-fit of soil water retention data to a wide range of SWRC models using a data set with high-resolution water retention data mainly obtained via the HYPROP system. This paper aims at evaluating the performance of the respective parametric PTFs and summarizes the overall result. Overall, the results of the companion study showed that the alternative variants provided better fits than the original unimodal models did [25]. A critical unanswered question, however, is whether these alternative variants hold their high performance when model parameters are estimated via parametric PTFs. The accuracy of parametric PTFs is impacted by two sources of error interacting in a complicated manner. The first error occurs when estimating the parameters of the selected soil water retention model using basic soil properties. The second error occurs when determining the SWRC by the chosen model. Consequently, it is crucial to evaluate the direct-fit performance of the soil water retention models and compare it to the performance of the respective parametric PTFs, the main objective of this study.

Materials and Methods
Figure 1 provides an overview of the data collection and analysis steps followed in this study. A detailed explanation of all these steps is provided in the following sections.  Table 1 summarizes the properties of the data set used in this study consisting of 81 soil samples (7729 measured water retention points) collected mainly from the areas surrounding Ankara, Turkey in 2010. The samples were collected from surface soil layer, approximately 0-30 cm of depth, when the soil water content was below the field capacity. The sampling locations were chosen from long term fallow lands and non-agricultural lands to eliminate the effect of agricultural management practices on soil structure and in turn on soil hydraulic properties [5]. The textural distribution of the samples is shown in Figure 2. The soil physical properties were measured in the soil physics laboratory at Ankara University in Turkey. The soil organic matter content (SOM) was estimated using measured soil organic carbon content via the modified method of Walkley and Black [26]. The hydrometer method [27] was used to measure soil texture (percentages of sand, silt and clay). The wet sieving apparatus (Eijkelkamp Agrisearch Equipment, Giesbeek, The Netherlands) was used to determine the percentage of the stable aggregates (SA). The apparatus works based on the principle that when aggregates with different levels of stability are submerged into the water, the unstable aggregates break down more easily. The analysis was done by wet sieving of 1-2 mm air-dried aggregates (using sieves with 60 Mesh screen) as described by Kemper and Rosenau [28]. The undisturbed samples were used to determine the bulk density, BD [29]. More details about this dataset are outlined in Haghverdi et al. [5,18]. The water retention data of the repacked samples were measured using the HYPROP system at the Technical University of Braunschweig in Germany. The companion paper [25] provides detailed information about the HYPROP system and the measurement campaign.   Table 1 summarizes the properties of the data set used in this study consisting of 81 soil samples (7729 measured water retention points) collected mainly from the areas surrounding Ankara, Turkey in 2010. The samples were collected from surface soil layer, approximately 0-30 cm of depth, when the soil water content was below the field capacity. The sampling locations were chosen from long term fallow lands and non-agricultural lands to eliminate the effect of agricultural management practices on soil structure and in turn on soil hydraulic properties [5]. The textural distribution of the samples is shown in Figure 2. The soil physical properties were measured in the soil physics laboratory at Ankara University in Turkey. The soil organic matter content (SOM) was estimated using measured soil organic carbon content via the modified method of Walkley and Black [26]. The hydrometer method [27] was used to measure soil texture (percentages of sand, silt and clay). The wet sieving apparatus (Eijkelkamp Agrisearch Equipment, Giesbeek, The Netherlands) was used to determine the percentage of the stable aggregates (SA). The apparatus works based on the principle that when aggregates with different levels of stability are submerged into the water, the unstable aggregates break down more easily. The analysis was done by wet sieving of 1-2 mm air-dried aggregates (using sieves with 60 Mesh screen) as described by Kemper and Rosenau [28]. The undisturbed samples were used to determine the bulk density, BD [29]. More details about this dataset are outlined in Haghverdi et al. [5,18]. The water retention data of the repacked samples were measured using the HYPROP system at the Technical University of Braunschweig in Germany. The companion paper [25] provides detailed information about the HYPROP system and the measurement campaign. aggregates break down more easily. The analysis was done by wet sieving of 1-2 mm air-dried aggregates (using sieves with 60 Mesh screen) as described by Kemper and Rosenau [28]. The undisturbed samples were used to determine the bulk density, BD [29]. More details about this dataset are outlined in Haghverdi et al. [5,18]. The water retention data of the repacked samples were measured using the HYPROP system at the Technical University of Braunschweig in Germany. The companion paper [25] provides detailed information about the HYPROP system and the measurement campaign.
The bimodal variants formulate the SWRC as the weighted sums of two curves, conceptually representing the textural and structural pore domains of the heterogeneous pore systems. They can provide a better fit to the measured soil water retention data since they have a greater number of free fitting parameters than unimodal expressions. The PDI-variant addresses the undefined dry-end water Water 2020, 12, 896 5 of 29 content issue of the most unimodal expressions by ensuring zero water content at the oven dryness. The PDI and Bimodal-PDI variants were considered for all unimodal expressions except for the BC model. The bimodal variants were considered for the K, the VG and the VG m models. HYPROP-FIT software was used to fit the soil hydraulic models to the measured water retention data. All models are explained in detail, and equations are provided in Appendix A. Table 3 illustrates six modeling scenarios (hereafter referred to as S1, S2, S3, S4, S5 and S6) examined in this study representing different combinations of PTFs' input predictors. The first three scenarios include the soil properties that are typically used to derive SWRC PTFs including soil texture (i.e., sand, silt and clay percentages), BD and SOM. In the remaining scenarios, the saturated water content (θ s ) and SA were considered as additional inputs. Prior to developing PTFs, a preliminary correlation analysis was conducted to check for normality (Shapiro-Wilk W statistic), detect linear relationships between all variables and examine the data structure. The preliminary analysis revealed that logarithmic transformation was necessary for several parameters from different models (response variables), resulting in a greater correlation with input predictors and a more normal distribution. We used a multiple linear regression with interactions and quadratic terms included. A stepwise regression technique was implemented to screen the variables. Various significance levels were examined to enter and retain the variables in the model, and based on the results, 0.15 level was chosen as a suitable criterion to screen the variables for most of the models. The 0.15 level, however, was adjusted whenever no input variables were selected by the stepwise technique. Multiple regression diagnostics were then considered to finalize the list of variables including the first and second moment specification test to check the equal residual variance, the Shapiro-Wilk W statistic to check the normality of the residuals and the condition index to check the collinearity between the variables. To validate the developed equations, the data set was randomly divided into eight, almost equal-sized, subsets (the size of each subset was roughly ten samples). Seven subsets were used to develop the equations, and the remaining subset was utilized to test the PTFs. The process was repeated eight times, each time considering a new subset as the test subset.

Developing Parametric Pedotransfer Functions
The 16 soil hydraulic models had a total of 96 coefficients (response variables). In scenario 6, θs was considered as input, which reduced the total number of response variables to 80 (96 − 16 = 80). Therefore, SAS 9.4 software program [37] was used to develop a total of 4480 regression-based PTFs, (5 modeling scenarios × 96 parameters + 1 scenario × 80 parameters) × 8 folds. The raw model parameters estimated by the PTFs were post-processed to ensure the values are in the physically meaningful and feasible range. The coefficients were then utilized to predict the soil water contents at measured soil tensions for all soils.

Evaluation Criteria
Three statistics were calculated including the mean absolute error (MAE, Equation (1)), the root mean square error (RMSE, Equation (2)) and the correlation coefficient (r, Equation (3)) to evaluate the performance of the PTFs. Furthermore, mean bias error (MBE, Equation (4)) was used to assess the degree of model bias (overestimation or underestimation): where M and E are the measured and the PTF-estimated soil water content values (cm 3 cm −3 ), respectively; M and E are the mean measured and the mean estimated water content values, and n is the total number of measured water retention points for all 81 samples (n = 7729). The evaluation statistics were calculated only for the test subset after combining the PTF estimations for the 8 folds resulting in a total of 96 values (16 models × 6 scenarios) for each statistic. When averaged across the scenarios, the K PDI and the VG m bimodal models show the lowest and the highest MAE values of 0.036 cm 3 cm −3 and 0.072 cm 3 cm −3 , respectively. When the MAE is calculated for each sample separately, the VG m bimodal provides the worst fit for 24% to 61% of the samples across scenarios. . Scatter plots of the fitted (dark blue) and the PTF-estimated (white) versus the measured soil water content data (scenario 6). 1 BC: Brooks and Corey [34], FX: Fredlund and Xing [35], K: Kosugi [36], VG and VGm: van Genuchten [24] constrained and unconstrained unimodal models. PDI and b denote Peters-Durner-Iden [30,31] and bimodal variants of the models, respectively. Table 4. Performance evaluation of the unimodal and unimodal-Peters-Durner-Iden (PDI) parametric PTFs for all modeling scenarios (S1-S6).  . Scatter plots of the fitted (dark blue) and the PTF-estimated (white) versus the measured soil water content data (scenario 6). 1 BC: Brooks and Corey [34], FX: Fredlund and Xing [35], K: Kosugi [36], VG and VGm: van Genuchten [24] constrained and unconstrained unimodal models. PDI and b denote Peters-Durner-Iden [30,31] and bimodal variants of the models, respectively.

Models
When averaged across the models, scenarios 3 and 6 have the highest and lowest MAE values of 0.055 cm 3 cm −3 and 0.037 cm 3 cm −3 , respectively. The minimum MAE values for all models belong to S6. The second-best performance for 11 out of the 16 models belongs to S1. The maximum MAE values belong to multiple scenarios for different models with S2 and S5 having the highest MAE values for a total of ten and four models, respectively.
For S6, the first to fifth best models are the VG m (MAE = 0.024 cm 3  . When the MAE is calculated for each sample separately, the VG m bimodal followed by the VG bimodal and the VG bimodal-PDI provide the worst fit for 56%, 13%, and 13% of the samples, respectively. The VG PDI, the K PDI, and the VG provide the best fit for 16%, 13% and 11% of the samples, respectively. Table 4. Performance evaluation of the unimodal and unimodal-Peters-Durner-Iden (PDI) parametric PTFs for all modeling scenarios (S1-S6).   [35], K: Kosugi [36], VG and VG m : van Genuchten [24] constrained and unconstrained unimodal models. PDI and b denote Peters-Durner-Iden [30,31] and bimodal variants of the models, respectively; RMSE: root mean square error (cm 3 cm −3 ), MAE: mean absolute error (cm 3 cm −3 ), r: correlation coefficient, MBE: mean bias error (cm 3 cm −3 ). Table 6 shows the correlation coefficient values between the PTF-estimated and fitted model parameters for all models and the six modeling scenarios. The average of positive correlation coefficients is consistently low across the scenarios ranging between 0.45 to 0.54. The lowest positive r values were close to zero for all scenarios, while the highest values varied between 0.89 to 0.95. On average across scenarios, strong correlations (r > 0.70) between the PTF-estimated and fitted model parameters are only observed for 31% of the parameters from which only 7% were parameters besides θ s and θ r .     [30,31] and bimodal variants of the models, respectively; S1 to S6 denote modeling scenarios 1 to 6, respectively. Figure 4 depicts the performance of the models across four soil textural classes for the six modeling scenarios. Information is also provided regarding the total number of the measured water retention pairs for each soil texture. respectively. Figure 4 depicts the performance of the models across four soil textural classes for the six modeling scenarios. Information is also provided regarding the total number of the measured water retention pairs for each soil texture. Figure 4. Performance of the 16 models for 4 textures and the six modeling scenarios (S1-S6). BC: Brooks and Corey [34], FX: Fredlund and Xing [35], K: Kosugi [36], VG and VGm: van Genuchten [24] constrained and unconstrained unimodal models. PDI and b denote Peters-Durner-Iden [30,31] and bimodal variants of the models, respectively. L: Loam (1087 data pairs), CL: Clay Loam (1277 data pairs), SL: Silt Loam (1348 data pairs), C: Clay (3708 data pairs). MBE: mean bias error (cm 3 cm −3 ).  [34], FX: Fredlund and Xing [35], K: Kosugi [36], VG and VG m : van Genuchten [24] constrained and unconstrained unimodal models. PDI and b denote Peters-Durner-Iden [30,31] and bimodal variants of the models, respectively. L: Loam (1087 data pairs), CL: Clay Loam (1277 data pairs), SL: Silt Loam (1348 data pairs), C: Clay (3708 data pairs). MBE: mean bias error (cm 3 cm −3 ).

Performance across Textural and Tension Classes
The data set contained only a limited number of samples with loamy sand, sandy clay loam and silt loam textures; therefore, these textures were excluded from the analysis. Some high MAE values are observed for the loam soils across multiple models and modeling scenarios. Average MAE values (data not shown here) across models are highest for the loam soils for four out of six modeling scenarios. For S6, the MAE values of the models for the sandy loam soils show relatively lesser fluctuations compared to that for the other textures. Overall, the best models (i.e., unimodal and PDI variants) performed superior for all textures while models with lower performance (i.e., bimodal and bimodal-PDI variants) performed poorly for all textures. Figure 5 illustrates the performance of the models at different parts of the curve divided into a total of 11 tension classes: 10 classes from pF 0 to 3 (100 cm increments) and one class for pF values larger than 3. Overall, the relative performance of the models across tension classes was not sensitive to the modeling scenarios. The MAE values show relatively low fluctuations across tension classes for the original unimodal expressions ranging roughly between 0.02 cm 3 cm −3 to 0.06 cm 3 cm −3 . A similar trend is observed for the PDI variants. However, the bimodal and bimodal-PDI variants, in particular variants of the K, the VG and the VG m models, show relatively larger MAE values as well as greater degrees of fluctuations across tensions classes. For the VG m bimodal, the K bimodal and the VG m bimodal-PDI models, MAE values peak at the wet and the dry parts of the PTF. For the VG PDI and the VG bimodal-PDI, MAE values increase from the wet-end up to the intermediate tension range but level off toward the dry part. larger than 3. Overall, the relative performance of the models across tension classes was not sensitive to the modeling scenarios. The MAE values show relatively low fluctuations across tension classes for the original unimodal expressions ranging roughly between 0.02 cm 3 cm −3 to 0.06 cm 3 cm −3 . A similar trend is observed for the PDI variants. However, the bimodal and bimodal-PDI variants, in particular variants of the K, the VG and the VGm models, show relatively larger MAE values as well as greater degrees of fluctuations across tensions classes. For the VGm bimodal, the K bimodal and the VGm bimodal-PDI models, MAE values peak at the wet and the dry parts of the PTF. For the VG PDI and the VG bimodal-PDI, MAE values increase from the wet-end up to the intermediate tension range but level off toward the dry part. Figure 5. Performance of the 16 models across tension classes and the six modeling scenarios (S1-S6). BC: Brooks and Corey [34], FX: Fredlund and Xing [35], K: Kosugi [36], VG and VGm: van Genuchten [24] constrained and unconstrained unimodal models. PDI and b denote Peters-Durner-Iden [30,31] and bimodal variants of the models, respectively. pF classes (cm of water): C1 (<100), C2 (100-200), C3 (200-300), C4 (300-400), C5 (400-500), C6 (500-600), C7 (600-700), C8 (700-800), C9 (800-900), C10 (900-1000) and C11 (>1000). MBE: mean bias error (cm 3 cm −3 ). Figure 5. Performance of the 16 models across tension classes and the six modeling scenarios (S1-S6). BC: Brooks and Corey [34], FX: Fredlund and Xing [35], K: Kosugi [36], VG and VGm: van Genuchten [24] constrained and unconstrained unimodal models. PDI and b denote Peters-Durner-Iden [30,31] and bimodal variants of the models, respectively. pF classes (cm of water): C1 (<100), C2 (100-200), C3 (200-300), C4 (300-400), C5 (400-500), C6 (500-600), C7 (600-700), C8 (700-800), C9 (800-900), C10 (900-1000) and C11 (>1000). MBE: mean bias error (cm 3 cm −3 ).

Overall Performance of the Parametric PTFs
The RMSE values reported for parametric PTFs in the literature typically range from 0.034 cm 3 cm −3 to 0.085 cm 3 cm −3 [11]. Therefore, PTFs developed using the HYPROP data, both in this study and by Haghverdi et al. [18], rank high among previously published parametric PTFs. Our results (Table 6), however, show that a majority of the SWRC model parameters cannot be estimated with high accuracy as a function of basic soil properties. The low correlations and differences between parametric PTFs estimated versus fitted soil hydraulic parameters were also reported by other researchers [20,38]. Tomasella et al. [21] argued that the relationships between the basic soil properties and the model parameters are often complicated and also vary at different parts of the curve, which makes them difficult to estimate using parametric PTFs. Schaap et al. [20] mentioned that macroscopic variables typically used as inputs for parametric PTFs could not predict information contained by the SWRC.
In our study, not highly accurate PTF-estimated parameters (average positive r ranged between 0.45 to 0.54 among scenarios) resulted in reasonable soil water content estimations by the water retention models (average r ranged between 0.86 to 0.93 among scenarios), a result also reported by Schaap et al. [20]. The VG m parametric PTF with an RMSE of 0.034 cm 3 cm −3 was the best PTF when all input predictors were considered. When averaged across modeling scenarios with different combinations of input predictors, the K-PDI model with an RMSE of 0.045 cm 3 cm −3 showed the highest performance. This performance is also comparable to the performance of neural network PC-PTFs developed using similar soils with RMSE values of 0.033 cm 3 cm −3 and 0.045 cm 3 cm −3 with the HYPROP and sandbox/pressure plates data, respectively [5,18].

Performance Across Soil Textures and Tension Classes
Overall, the performance of the PTFs across soil textures was consistent such that highly ranked models (i.e., unimodal and PDI variants) performed superior for almost all soil textures. Haghverdi et al. [5] showed that the accuracy of PC-PTFs developed using equilibrium data for each textural class was affected by the percentage share of the textural class in the PTF training dataset. Haghverdi et al. [18] noticed the same trend when the HYPROP data was used to develop PC-PTFs and suggested to include at least 13% of data for each textural class in the PTF training data set. However, we did not see a clear relationship between the number of soil samples per textural class and the performance of the parametric PTFs for each texture. This is because PC-PTFs developed by Haghverdi et al. [5,18] are data-driven models, which means they need a large enough training data set to learn the shape of the WRC for different textures. For parametric PTFs, however, the soil hydraulic models control the shape of the WRCs across all textural classes. Since our data set is mainly dominated by four soil textures and given the limited and uneven number of samples across soil textural classes, further studies are needed with global data sets with a more comprehensive soil textural distribution to confirm our findings.
Cornelis et al. [2] reported a reasonable performance for most parametric PTFs they tested at wet-end (saturation) and dry part of the SWRC (near PWP), while the lowest performance occurred around the field capacity. We only observed high fluctuations in MAE values at different parts of the SWRC for PTFs with low performance (i.e., bimodal and bimodal-PDI PTFs). The PTFs with the highest performance in our study (i.e., PDI variants and unimodal expressions) showed a good stable performance from the wet to the dry boundaries of the PTF.

Importance of the Input Parameters
We could not detect any consistent improvement in the performance of the PTFs when BD (S2: average RMSE: 0.065 cm 3 cm −3 ) and SOM (S3: average RMSE: 0.068 cm 3 cm −3 ) were considered as input predictors in addition to soil texture (S1: average RMSE: 0.061 cm 3 cm −3 ). When averaged across the models, the first five scenarios (S1-S5) showed a relatively similar performance with average RMSE values varying between 0.061 to 0.068 cm 3 cm −3 . The best performance of all models occurred at S6 (average RMSE: 0.052 cm 3 cm −3 ) when θ s was added as an additional input predictor. The low effect of SOM and the high impact of θ s as input predictors were also reported for PC-PTFs developed using similar soils by Haghverdi et al. [18]. The low impact of SOM is attributed to the fact that a majority of the soils used in the two studies were collected from the dry Central Anatolia Region of Turkey with low organic matter content. Haghverdi et al. [4] also reported only a slight improvement in the performance of their PTFs when OC and BD were added as extra predictors. Cornelis et al. [2] studied multiple parametric PTFs and observed no correlation between SOM, BD and the performance of the PTFs. Vereecken et al. [11] reviewed several studies on van Genuchten based parametric PTFs and noticed improvements in the performance of parametric PTFs when the water content data were considered as input in addition to the basic soil properties.
When the HYPROP data was used to develop PC-PTF by Haghverdi et al. [18], excluding θ s as an input predictor remarkably decreased the accuracy of the most PTFs down to unacceptable performance levels (RMSE values ranged from 0.061 to 0.159 cm 3 cm −3 ). PTFs without θ s as input performed much better in this study (average RMSE values ranged from 0.061 to 0.068 cm 3 cm −3 ) compared to Haghverdi et al. [18], which is a promising result given the practical difficulties to measure θs directly. We attribute this to the differences between the structure of the PC-PTFs developed by Haghverdi et al. [18] and the parametric PTFs derived in this study. The PC-PTF relies on neural networks to relate basic soil properties to the whole SWRC, whereas in parametric PTFs, soil hydraulic models govern the shape and continuity of the estimated SWRC. Consequently, including θ s as a measure of the soil water content upper boundary seems to be much more critical when PC-PTFs are developed using the HYPROP data.   [34], FX: Fredlund and Xing [35], K: Kosugi [36], VG and VGm: van Genuchten [24] constrained and unconstrained unimodal models. PDI and b denote Peters-Durner-Iden [30,31] and bimodal variants of the models, respectively. MBE: mean bias error (cm 3 cm −3 ).

Conclusions
We developed and evaluated a set of parametric PTFs for a total of 16 soil water retention models, including five original unimodal models and their PDI, bimodal and bimodal-PDI variants. Table A1 provides regression equations (parametric PTFs) for estimating the parameters of the selected soil water retention models for the six modeling scenarios (S1-S6). Our results showed that considering saturated water content as an additional input predictor enhanced the performance of the parametric PTFs for all models. However, parametric PTFs that were developed only using soil texture as input also showed promising performances. When all input predictors were considered, the parametric PTF of the van Genuchten model [24] with five free parameters showed the highest performance (RMSE of 0.034 cm 3 cm −3 ). The parametric PTF of the PDI variant [30,31] of the Kosugi model [36] showed the highest average performance across six modeling scenarios with different combinations of input predictors (0.045 cm 3 cm −3 ). The high number of model parameters of bimodal  [34], FX: Fredlund and Xing [35], K: Kosugi [36], VG and VG m : van Genuchten [24] constrained and unconstrained unimodal models. PDI and b denote Peters-Durner-Iden [30,31] and bimodal variants of the models, respectively. MBE: mean bias error (cm 3 cm −3 ).
When MAE values were correlated to the number of model parameters, a positive strong correlation was observed (r = 0.71) for the fitted models while the correlation was negative (r = −0.79) for the parametric PTFs. This is because a greater number of free parameters make models very robust and flexible when directly fitted to the high-resolution HYPROP data. On the other hand, a high number of response variables to estimate when parametric PTFs are developed means more potential sources for error which, as our results show, substantially diminishes the performance of the bimodal parametric PTFs. We should highlight that in this paper, bimodal models were fitted to repacked samples to investigate how the number of free model parameters impacts the direct fit performance of the models as well as the efficacy of the respective parametric PTFs. In practice, bimodal SWRC models should only be fitted to soils with evident soil structure and matrix effects. Since the structure of the samples was impacted during repacking, the increase in the performance of the bimodal variants is mainly attributed to the rise in the number of fitting parameters rather than bimodality.
The original unimodal expressions (except for the FX model) ranked lowest when fitted to the observations (average MAE: 0.011 cm 3 cm −3 ), but their parametric PTFs on average performed as high as the PDI variants (S6, average MAE: 0.029 cm 3 cm −3 ). In the first part of this study [25], we showed that the PDI variants are suggested for full SWRC parametrization using HYPROP data due to their superior dry-end performance (average MAE: 0.017 cm 3 cm −3 ) compared to the traditional unimodal expressions (average MAE: 0.041 cm 3 cm −3 ). Therefore, the combined results of both parts of this study suggest the application of PDI variants for full SWRC parametrization (direct fit and parametric PTFs) using HYPROP data. More investigation is needed to determine whether these findings hold for different data sets, parametric PTFs developed using data mining techniques and the nonparametric k-nearest PTF of Haghverdi et al. [3].

Conclusions
We developed and evaluated a set of parametric PTFs for a total of 16 soil water retention models, including five original unimodal models and their PDI, bimodal and bimodal-PDI variants. Table A1 (in Appendix B) provides regression equations (parametric PTFs) for estimating the parameters of the selected soil water retention models for the six modeling scenarios (S1-S6). Our results showed that considering saturated water content as an additional input predictor enhanced the performance of the parametric PTFs for all models. However, parametric PTFs that were developed only using soil texture as input also showed promising performances. When all input predictors were considered, the parametric PTF of the van Genuchten model [24] with five free parameters showed the highest performance (RMSE of 0.034 cm 3 cm −3 ). The parametric PTF of the PDI variant [30,31] of the Kosugi model [36] showed the highest average performance across six modeling scenarios with different combinations of input predictors (0.045 cm 3 cm −3 ). The high number of model parameters of bimodal variants negatively impacted the performance of their parametric PTFs. PDI variants, on the other hand, have the same number of free parameters as original unimodal and showed an excellent fit to the complete SWRC [25], and their parametric PTFs highly ranked among all models. Consequently, combined results of both parts of this study recommend the application of PDI variants and FX model for estimation and parametrization of the complete SWRC (i.e., from saturation to oven dryness) using typical data obtained by the HYPROP system.  The BC model [34] 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, and θ r and θ s are the residual and saturated soil water contents, respectively. The FX model [35] is given by: and where α, n, and m are the SWRC shape parameters, h r is the tension corresponding to θ r , h 0 is the soil tension at zero water content, and e is the Euler number. The K model [36] is expressed as: 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 [24] is represented as: where α and n are the curve shape parameters and other parameters as previously defined. The VG m model [24] with m as an additional shape parameter has five free parameters:

Appendix A.2 PDI Expressions:
The general form of the PDI model [30,31] 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 (A10) The water adsorption saturation function is given as [30]: 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 (A15)

Appendix A.3 Bimodal Expressions:
The bimodal expressions for each unimodal model are 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).

Appendix A.4 Bimodal-PDI Expressions:
The bimodal-PDI expression for each model is 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.

Appendix B
Note that the goal of this study is only to develop PTFs for soils with identical conditions and within the range of the soils used in this study. We realize that a much bigger data set is needed to develop PTFs applicable to other regions in the world with different soil conditions. Table A1. List of regression equations (parametric PTFs) for estimating the parameters of the selected soil water retention models for the six modeling scenarios, S1-S6 (bold text).