INTRODUCTION

Soils in the arid region are grouped mostly under Entisol and Aridisol orders, characterized by sand domination, diminutive organic matter content, salinity accumulation, carbonate precipitation, and severe low moisture content through an unapologetic horizon (Soil survey staff, 1999; Verheye, 2008). These characteristics substantially affect the physical properties and the soil voids ratio, altering the soil water holding capacity and soil water characteristic curve (SWCC). The soil water characteristic curve or soil water retention curve describes the retained moisture in the soil voids under equilibrium at a given metric potential passing by saturation and unsaturation conditions (Childs, 1940; Hillel, 2003).The soil water characteristic curve is essential in agriculture, environment, hydrology, and geotechnical studies (Ellithy, 2017; Zapata, 1999). However, measurement of SWCC is tedious work that can generate errors, be cost-inefficient, and time-consuming. It can take weeks to determine a curve (Khlosi et al., 2008; Rudiyanto et al., 2021; Shani and Or, 1995; Wesseling, 2009).

Closed-form analytical equations were used to describe the fitting line that linked moisture θ, in the form of volumetric, gravimetric, or degree of saturation, with soil retention ψ, as kpa, bar, or PF (Novák and Hlaváčiková, 2019). In the last five decades, the most recognized classical unimodal functions developed by Brooks and Corey (1964) and by van Genuchten (1980), were validated and verified by many researchers for a wide range of soil types under different conditions (Chen et al., 2016; Du, 2020; Ellithy, 2017; Ellithy et al., 2018; Haghverdi et al., 2020; M. Khlosi et al., 2008; Leong and Rahardjo, 1997; Madi et al., 2018; Morel-Seytoux et al., 1996; Porebska et al., 2006; Seki, 2007).

                The powered law equation presented by Brooks and Corey (1964) described the relation of θ(ψ) as

where θs is volumetric water content at retention ψ, ϕ is porosity, θr is residual water content, λ is pore distribution index, and ψb is a parameter equals the air entry.

van Genuchten (1980) proposed a smooth, closed-form, three-parameter model for the soil-water characteristic curve in the form of

where n and m (m = 1 –n-1) are empirical parameters.

Tow general methods to drive equations 1-3, either by the direct fitting of θ and ψ or by estimating the parameters of the equation using basic soil properties as was named pedotransfer function (PTF) (Guber et al., 2010; Jaiswal et al., 2013; Seki, 2007; Shwetha and Varija, 2013).

Bouma (1989) introduced 'pretransfer function' (PTF). He described PTF as rendering functions for easily and routinely measured raw soil survey data, which was mainly collected over the last four decades, into more useful predictive equations for soil properties with a conscientious accuracy (Odeh and McBratney, 2005; Pachepsky and van Genuchten, 2011; Van Looy et al., 2017). Adapting PTF in soil hydraulic equations was prospered during the last four decades; as a result, tens of equations and models were generated with varying accuracy and reliability depending on sample population and model structure perfection (Du, 2020; Nemes and Rawls, 2006; Ostovari et al., 2015; Rawls and Brakensiek, 1982; Saxton and Rawls, 2006; Vereecken et al., 2010, 1989; Wösten et al., 1999). Botula et al. (2014) reviewed and categorized 35 PTFs collected through 40 years from 35 publications; his finding was that 80% of those PTFs were based on multiple linear regressions and polynomial of the nth order. Recently Abdelbaki (2021) examined 30 PTFs, 11 discrete functions, and 19 continuous functions; both types of PTFs resulted in a different accuracy based on the different soil classes or moisture content levels. Different published PTFs were tested and validated with the local soils and parameters by many researchers. The results quality fluctuated depending on soil parameters, type, and soil potential level (Abbasi et al., 2011; Botula et al., 2012; Cichota et al., 2013; Cornelis et al., 2001; Schaap, 2004).

Guber et al. ( 2009, 2006) investigated the validity of 21 PTFs to estimate the parameters of Brooks and Corey (1964) and van Genuchten's (1980) water retention equation, developed in 2010 to be a computer program to calculate PTFs named (CalcPTF) (Guber et al., 2010). Many researchers have evaluated the performance of the CalcPTF program. Jaiswal et al. (2013) reported, without including the soil physical properties data, adequate accuracy of the program using the equations of Rawls and Brakensiek (1985) and Saxton et al. (1986) in predicting Brooks and Corey equation parameters. On the other hand, Tomasella and Hodnett (1998) gave the best result for the van Genuchten equation. Cassinari et al. (2015) examined the CalcPTF using clayey closed-landfill soil (clay>54%); a general conclusion was that models at a lower suction overestimated the results, but at a higher suction, they underestimated the results. He also concluded that continuous pedotransfer function is the closest to the measured data; on the contrary, the curve by Tomasella and Hodnett (1998) is the worst. The poorly performance of the CalcPTF as a consolidated program or as individual PTF equations were reported frequently by many researchers (Abdelbaki, 2021; Cassinari et al., 2015; Castellini and Iovino, 2019; Dai et al., 2013; Guram and Bashir, 2020; Hewelke et al., 2017; Patil et al., 2016).

RMSE was used to determine the accuracy of models, but different results were garnered for the equations with the lowest RMSE among those in CalcPTF without consistent justification. Rawls and Brakensiek (1985) then Saxton et al. (1986) equations attainted better results among 16 tested PTFs for Indian tropical soils (Jaiswal et al., 2013). Ghanbarian-Alavijeh and Liaghat (2009) reported that Saxton et al. (1986) estimated the soil water retention curve better than Campbell and Shiozawa (1992). In contrast, Castellini and Iovino (2019), in their work with clay soils, found that Saxton et al. (1986) equations had the highest RMSE compared with Vereecken et al. (1989) and Wösten et al. (1999) equations. Wösten et al. (1999) equation performed poorly for estimating van Genuchten's (1980) equation parameters (Dai et al., 2013; Hewelke et al., 2017; Matula et al., 2007). Patil et al. (2016) reported a good estimation for Tomasella and Hodnett's (1998) equations for predicting van Genuchten's (1980)  equation parameters for Brazilian soils. The previous controversial results and many researchers' findings in different locations and countries confirmed that the application of PTFs to soils different from those used for their development could be erroneous in the estimation process. They recommended using a small local data to develop PTFs than implementing PTFs developed outside the local domain (Castellini and Iovino, 2019; Patil et al., 2016).

The aims of this study were to conduct a comprehensive performance evaluation of CalcPTF (14 PTFs) by using both Brooks-Corey and van Genuchten equations‎, which ‎complied with the available data, using local soil located in a severe arid region. Besides examining the PTFs ‎accuracy and model suitability, a set of statistical measurements will be conducted including correlation coefficient ‎‎(R), root mean square error (RMES), Nash-Sutcliffe efficiency (NSE), the ratio of RMSE to the standard deviation, ‎Percent bias (PB), and Akaike Information Criterion (AIC). The study will be too used the CalcPTF program to calculate the parameters of both Brooks-Corey and van Genuchten equations.‎

MATERIALS AND METHODS

Study area and soil sampling: The study was conducted in the Al-Ahsa region, commonly known as the largest and oldest agricultural and habitation area on the Arabian Peninsula. Al-Ahsa oasis is located about 70 kilometers west of the Arabian Gulf, between latitudes 25 21' and 25 37' N and longitudes 49 33' and 49 46' E (Figure 1). It has a total surface area of 320 square kilometers.  Al-Ahsa is an extremely arid ecosystem, where the average annual precipitation is less than 73 millimeters.

Fig. 1: Al-Ahsa general areal image shows the geographical position and sample location. (image  source: Esri, DigitalGlobe,  GeoEye,  Earthstar  Geographics, CNES/Airbus DS, USDA, USGS, Aerogrid, IGN, and the GISUser Community)

In total, 36 samples were collected from 36 locations in Al Ahsa. With an accuracy of fewer than five meters, handheld GPS devices were used to determine the location of soil samples in the field.  A soil auger with a diameter of 10 cm was used to collect soil samples from disturbed areas at a depth of approximately 30 cm. Samples were ground, air-dried, thoroughly mixed, and sieved through a 2 mm sieve before being stored for physical and hydraulic measurements. The soil particle size (sand%, silt%, and clay%), bulk density ρ, and saturation percentage θ_s were measured in the laboratory using the standard methods published by the Soil Science Society of America (Reynolds et al., 2002).     

Soil water characteristic curve (SWCC):This study measured SWCC using the filter paper method described in ASTM D 5298 (ASTM D6836-02, 2002) and other scientific publications (Al-Khafaf and Hanks, 1974; Bulut and Leong, 2008; Scanlon et al., 2002). The filter paper technique has been extensively investigated and validated (Bulut, 1996; Bulut et al., 2001; Elgabu, 2013; Tripathy et al., 2014) as an indirect method for suction measurement. The filter paper (Whatman No. 42) was sandwiched between two protective filter papers placed between two identical halves of the soil specimen. The soil was packed tightly to ensure perfect contact between soil and filter paper. The cane was sealed to prevent any moisture loss. The soil was packed in moisture canes equal to its bulk density for each chosen moisture content. Canes were kept in an incubator for seven days to ensure equilibrium at a constant temperature of 25°C. At the end of seven days, the filter paper was removed from the soil and weighed with a 0.0001g electronic balance to determine its wet weight. The filter paper was oven-dried at 105 C for 24 hours and weighed again to determine its water content. The metric suction ψ was determined by matching the filter paper moisture content with the calibration curves developed by Al-Khafaf and Hanks (1974) and ASTM D5298-16 (2016). The moisture content values for the field capacity FC and the wilting point WP were calculated using the suction vs. moisture curve at pF values of 2.52 and 4.18, respectively. 

Brooks & Corey and van Genuchten parameters: Brooks and Corey (1964) eq.(1) and van Genuchten (1980) eq.(3) were calculated from the measured data by using the online program for soil water retention curve.(SWRC.Fit) (https://seki.webmasters.gr.jp/swrc/) developed by Seki in (2007).

CalcPTF model description: CalcPTF, as described by Guber et al. (2010), is a computer program developed to estimate parameters for the Brooks and Corey eq (1) and van Genuchten eq. (3) water retention equations to support the multi modeling approach. Twenty PTFs (Table 1) were derived from a large database (12,625 soils) and categorized into two groups: continuous PTFs, which calculate the parameters of the closed-form equation governing the soil water content and matric potential (1-11), and discrete or point PTFs, which predict the soil water content at multiple matric potentials (12-20). Seven PTFs estimate Brooks and Corey parameters (1-7), four PTFs estimate van Genuchten parameters (8-11), and five models fit the van Genuchten equation to pairs of parameters estimated with PTFs (12-16). Four PTFs calculate the moisture content at field capacity (330 cm) and wilting point (15000 cm).  This code is written in FORTRAN and is invoked from an Excel worksheet.

Statistical analysis: According to Donatelli et al. (2004), limited testing makes it difficult for modelers to verify that the PTFs selected were sufficiently accurate. In general, the more tests conducted in which it cannot be demonstrated that the function is incorrect, the greater the confidence in the function (Donatelli et al., 2004; Schaap, 2004).

The performance of different analytical and PTF models was measured with the benefit of different statistical criteria, including correlation coefficient (R), root mean square error (RMSE), Akaike information criterion (AIC), Nash-Sutcliffe model efficiency coefficient (NSE), and percent bias (PB).

An analysis of the Pearson correlation coefficient was conducted to examine the relationships between laboratory-measured values and the estimated values derived from the PTFs model for both SWCC and the equation 1 and 3 parameters. The significance of the relationships were classified into four levels: no correlation when |R| < 0.28, weak correlation when 0.28 ≤|R|< 0.33, moderate correlation when 0.33 ≤|R|< 0.43, and strong correlation when 0.43 ≤ |R|≤ 1.0 (Addinsoft, 2021). The basic form of the correlation coefficient is shown in equation (4).

The xi and yi are measured and estimated variables, respectively, x̂ and ŷ are the mean.

Root mean square (RMSE) is a commonly used metric for calculating the variance between the estimated and measured values of a model or estimator. Compared with other models, ideal models should have a minimum positive RMSE value (Schaap, 2004). Furthermore, the recommended RMSE value ranged between 50% (lower side) to 30%  (higher side) of the standard deviation (SD) of the measured data (Moriasi et al., 2007; Ouatiki et al., 2020; Singh et al., 2005). The general form of the RMSE equation:

 To provide reliable information about the overall goodness of fit of a model, an NSE(Nash–Sutcliffe efficiency) is recommended as one of the most appropriate objective functions (Legates and Mccabe, 1999; McCuen et al., 2006; Nash and Sutcliffe, 1970; Willmott, 1981).Nash-Sutcliffe efficiency (NSE) is defined as a normalized statistic that indicates the relative magnitude of residual variance ("noise") with respect to estimated data variance (Moriasi et al., 2007).  According to Nash and Sutcliffe (1970), NSE measures how well the plot of measured and estimated data fits a 1:1 curve. The efficiency coefficient ranges from minus infinity to one (-∞ to 1.0), with larger values indicating better agreement. Thus, a zero value indicates that the observed mean is as good a predictor as the model, while negative values indicate the observed mean is a better predictor than the model (Wilcox et al., 1990).  Literatures indicate that the NSE values can be categorized into four groups: < 0.5 is unsatisfactory, 0.50 - 0.70 satisfactory, 0.70 - 0.80 good, and > 0.8 very good (Abdelbaki, 2021; Gupta et al., 1999; Moriasi et al., 2015). The form of the NSE equation:

                Moriasi et al. (2007) provided guidelines for assessing the accuracy of prediction models. He used an equation (5) based on the ratio of RMSE to the standard deviation (SD) of the measured data to characterize the appropriateness of the model fitting as follows: <0.5 (very good), 0.50 < RSR > 0.60 (good), 0.60 < RSR > 0.70 (satisfactory), and >0.7 unsatisfactory. Many researchers used this classification widely (Beharry et al., 2021; Carlos Mendoza et al., 2021; Mekoya, 2019; Pandey et al., 2021).                

n equals the number of samples

Percent bias (PB) reflects the tendency of predictions to be larger or smaller than their measured counterparts. The optimal value is 0.0. Positive values indicate a bias toward underestimation, whereas negative values indicate an overestimation bias (Gupta et al., 1999; Moriasi et al., 2007). The form of the PB equation:

                Akaike Information Criterion (AIC) is a statistical tool used to compare and select the best candidate model among many alternative models. AIC aims to select the model that best explains the variance of the dependent variable with the fewest number of independent variables (parameters). In other words, it facilitates the selection of a simpler model (fewer parameters) over a complex model (many parameters). Selection of AIC reduces the complexity of the model, which can lead to over fitting, as well as, the model is improved by reducing the number of unwanted parameters, which can contribute to additional noise that hampers model fit (Akaike, 1974). The lower the AIC value, the better the fit of the model and the general form of the AIC equation:

SS_e is the sum square of errors, n is the number of observations, and k is the number of parameters.

1-Performance evaluation methods: A total of fourteen PTFs were selected from Table 1. However, PTFs 5, 11, and 17-20 were eliminated for the following reasons:

a-Both PTFs 5 and 11 require organic carbon as an essential parameter. However, the samples collected in the study area do not contain organic carbon.

b-PTF 17-20 predicts moisture only at certain capillary pressures of 330 and 15000 cm, not on a full-scale SWCC.

1-An assessment of the CalcPTF model's performance as associated with equations 1 and 3 in estimating the complete set of soil water characteristics curves (SWCC) in two scenarios.

1. The ClacPTF models were used to predict the entire parameters set of equations 1 and 3 parameters (h_b, θ_r, θ_s, λ, α, and n).
2. Substituting the values of ϕ in equation 1 and θ_s in equation 3 with the value of θ_s estimated from equation (10) (Al-Saeedi, 2022).
3.  

        [10]

• Identify the contribution of the individual parameter to the estimation accuracy of equations 1 and 3 using BCS and VGW. Using BCS and VGW as reference PTFs models, different scenarios were applied and considered scenario one. The second and third scenarios replaced ϕ and θ_s with measured and estimated from equation 10. With measured values, other scenarios from three to six replace ψ_b, θ_r, λ, α, and n,
• Using sample number 35 and equations 1 and 3, a percentage bias ranging between (-30% to +30%) was applied to each of the individual parameters (one by one) of each equation separately to determine the magnitude of the change in the RMSE for SWCC prediction.

To review the accuracy of CalcPTF models and evaluate potential improvements associated with each scenario, the following statistical tests were conducted: R, RMSE, NSE, RSR, PB, and AIC. The use of different criteria for statistical measurements provides researchers with greater flexibility in validating and selecting the best model with fewer complications (Beharry et al., 2021; Donatelli et al., 2004; Golmohammadi et al., 2014; Moriasi et al., 2015, 2007; Ouatiki et al., 2020; Singh et al., 2005; Willmott, 1981).

RESULTS

The physical and hydraulic properties of Soil samples: Analysis of the soil particle content presented in (Fig. 2) revealed that the sand content ranged between 12 and 95.38 %, with a mean of 51.103%, and the silt content ranged between 1.61 to 86% and with mean of 43.538%. The clay content ranged between 1 and 19 %, with a mean of 5.359 (Table 2). Samples ranged between sand, sandy loam, and silty loam. In Table2, the descriptive statistics of the studied soils were presented as ϕ and θ_s that range from 0.198 cm³ cm^-3 to 0.566 cm³ cm^-3, with the mean being 0.398 cm³ cm^-3 and the standard deviation being 0.101 cm³ cm^-3. Bulk density varied from 0.960 g cm^-3 to 1.690 g cm^-3, with a mean of 1.282 g cm^-3 and an SD of 0.208 g cm^-3. According to the logarithm developed by Seki (2007), the equations 1 and 3 parameters are presented in Table 2. Brooks and Corey parameters, BC-θr (θr in Eq.1) ranged from 1×10-10 cm³ cm^-3 to 0.151 cm³ cm^-3, ‎the mean value was 0.044 cm³ cm^-3, and the SD was 0.042 cm³ cm^-3. Air entry hb ranged ‎from 1.452 to 267.700 cm, with a mean of 36.343 cm and an SD of 0.171 cm. pore size ‎index λ exhibited a minimum value of 0.069 and a maximum value of 0.855 with a mean ‎of 0.305 and an SD of 0.171. The van Genuchten equation 3 parameters were also ‎presented in Table 2. VG-θr minimum value was 1×10^-10 cm³ cm^-3, maximum 0.164 cm³ ‎cm^-3, mean 0.072 cm³ cm^-3, and SD 0.045 cm³ cm^-3. α ranged from 0.002 cm^-1 to 0.413 ‎cm^-1 with a mean of 0.078 cm^-1 and a standard deviation of 0.101 cm^-1. The n value ‎ranged between 1.104 and 2.494, with a mean of 1.488 and a standard deviation of ‎‎0.332.‎

BC and VG parameters prediction:

Table 3 presents a statistical evaluation of CalcPTF's accuracy in predicting Brooks and Corey's (eq. 1) parameters. There was a strong correlation between the estimated and the measured porosity ϕ for all models, with R values ranging from 0.863 to 0.896. At the same time, the other parameters (θ_r, ψ_b, and λ) did not show any significant correlation. On the other hand, the other statistical indicators (RMSE, NES, and RSR) for CalcPTFs models prediction for equation 1 parameters (ψ_b, θ_r, and λ) demonstrated an extremely high degree of uncertainty and deviation. They, therefore, were rated as unsatisfactory or invalid models.  As shown in Table 3, the CalcPTF models NES results for predicting equation 1 parameters were very low (< 0.5).

Fig. 2: Texture classes of the soil samples of Al-Hassa (clay (≤2 μm), silt (2–50μm), sand (50–2000 μm)] according to USDA classification.

The RSR values for all of the CalcPTF models were greater than 1.0, which indicated that none of the models could be considered qualified. Table 3 shows that all models were biased toward overestimating PB values for ϕ, with PB values ranging from -23.532 to -30.102. The θ_r was biased toward underestimating PB values for all models, with PB values ranging from 16.092 to 100. H_b was also overestimated for all models with PB values from 41.042 to 92.330. The λ showed underestimation for PB values ranging from 2.297 to 37.628 but overestimation for BCW. Table 4 presents the statistic indicators of CalcPTF models prediction performance for van Genuchten equation (eq. 3) (VGW, VGVA, and VGVE) parameters and discrete model (VGT, VGR1, VGG, VGRA, and VGR2) parameters in the form of van Genuchten equation (eq. 3). The Correlation coefficient (R) showed no significant results for all parameters except θ_s. While estimated θ_s showed a high significant value with measured θ_s for all models ranging between 0.559 and 0.942 (p=0.01), the NES demonstrated an unsatisfactory value of <0.5 for all models and parameters.

Also, RSR analysis revealed an unsatisfactory result with a high value >0.7 for all CalcPTF models and all parameters. PB for all CalcPTF models, θ_s were overestimated by -13.957 to -32.170, except VGW underestimated the values by 2.791. θ_r was overestimated for all CalcPTF models 42.210-100.000 except for VGG (-8.172). α was underestimated by all CalcPTF models with values ranging from 9.096 to 88.014, and VGVE was overestimated by -25.817. CalcPTF models underestimated the value of n, with a PB ranging from 1.765 to 74.611. Both continuous and discrete PTF models showed unsatisfactory and high uncertainty results

SWCC estimation: Table 5 summarizes the statistics performance evaluation values for the examined calcPTF models and the two modelling scenarios for predicting SWCC. 

According to the first scenario, all examined CalcPTFs models (Table 5 and Fig. 3) exhibited highly significant values of correlation coefficients (R), ranging from 0.726 (VGVE) to 0.892 (VGR2) between measured and estimated moisture content. Since correlation alone is not always a valid evaluation criterion when evaluating the validity of a model. NSE results in Table 5 indicate an unsatisfactory rating (< 0.50) for all CalcPTF models BCS, BCC, BCR, BCW, BCO, BCM, VGVE, VGT, VGV, and VGRA, with values of 0.454, 0.424, 0.470, 0.478, 0.500, -0.284, -0.117, 0.194, 0.269, and -0.375, respectively. On the other hand, VGW, VGVA, VGR1, and VGR2 models displayed NSE values greater than 0.5 with satisfactory ratings of 0.612, 0.632, 0.532, and 0.566, respectively. As confirmed by RSR, all of CalcPTF's models were rated unsatisfactory (>0.70) in Brooks and Corey and part of the van Genuchten model. However, the VGW, VGVA, VGR1, and VGR2 models were only rated satisfactory with 0.624, 0.608, 0.686, and 0.660, respectively. Compared to the satisfactory models (Table 5), the PB test results indicated overestimations for VGW, VGR1, and VGR2 and underestimations for VGVA, with values of -5.028, -3.434, -4.674, and 9.425, respectively. VGVA was rated the most satisfactory model based on AIC with a score of -3122.521, followed by VGW, VGR2, and VGR1, with scores of -3090.421, -3021.573, and -2975.013, respectively.

In the second scenario, by substituting ϕ and θs in equations 1 and 3 with θs estimated from equation 10, the statistical performance indicators of SWCC prediction were significantly improved for both Brooks and Corey and van Genuchten equations. Table 5 and Fig 4 show that all CalcPTF models have highly significant correlation coefficients ‎(R) ranging from 0.759 to 0.904, with a marginal improvement over the first scenario ‎of 0.01 to 0.02. ‎In contrast, the NSE and RSR values significantly improved to reduce the number of unsatisfactory models from ten in the first scenario to only five in the second. Additionally, the BCC and VGG models received very good ratings in NSE, equal to 0.777 and 0.750, respectively, and in RSR, equal to 0.480 and 0.500. A total of six models were rated good for NSE: BCR, BCW, VGW, VGT, VGR1, and VGR2, with values of 0.688, 0.683, 0.677, 0.743, 0.652, and 0.685, respectively, as well as good RSR ratings of 0.564, 0.570, 0.509, 0.591, and 0.562. The BCS model achieved satisfactory ratings for NSE (0.689) and RSR (0.559). There were only five unsatisfactory ratings, namely BCO, BCM, VGVA, VGVE, and VGRA, with NSE values of 0.345, -0.353, 0.371, -0.118, and 0.399 respectively, and RSR values of 0.811, 1.166, 0.795, 1.060, and 0.777 respectively. Table 5 and Figure 4 demonstrate that the BCC model underestimated the water content values by 4.068%, while the VGG model overestimated the water content by 3.077%.

The accepted rating models, BCS, BCR, BCW, VGT, VGR1, and VGR2, all underestimated soil water content by 14.488, 12.120, 13.729, 3.673, 17.146, and 15.681%, respectively. Only one model overestimated the water content VGW with a value of -10.058%. The AIC is examined in Table 5, which indicates that the BCC and VGG models have the lowest values of the fourteen models, indicating that they have the lowest prediction complexity levels of -3406.688 -3361.859, respectively. Regarding the other models that have received approval ratings, they are listed according to their lowest AIC values as follows: VGT, BCS, BCR, VGR2, BCW, VGW, and VGR1, having values of -3341.337, -3227.889, -3224.655, -3218.175, -3216.347, -3202.425, and -3157.200, respectively. The AIC value with the lowest value is the easiest and best predictive model. Based on the statistical criteria used in this study, the models BCO, BCM, VGVA, VGVE, and VGRA were rated unsatisfactory.

 ‎

Fig. 3: Estimated versus measured soil water content of all CalcPTFs equations for estimating the parameters of the Brooks and Corey (1964) equation (BC) and van Genuchten equation (VG).

Fig. 4: Estimated versus measured soil water content of all CalcPTFs equations and equation (10) for estimating the parameters of the Brooks and Corey (1964) equation (BC) and van Genuchten equation (VG).

Parameters contribution test: Brooks and Corey equation (Eq. 1), inferring the output of using fully CalcPTF parameters of BCS already listed in Table 5 with high correlation 0.888 while other criteria were unsatisfactory as NSE equal 0.5454 and RSR equal 0.741 (Fig 5). Scenarios two and three showed great improvement with good accuracy for scenario two R (0.909), NSE (0.724), and RSR (0.524), scenario three R (0.890), NSE (0.960), and RSR (0.550). Also, Figure 5 illustrates that despite the high correlation value of all other scenarios, using the measured values of θr, ψb, and λ but did not result in any improvements in the SWCC estimation outputs or the model rating accuracy but instead made it even worse than the first scenario. Fourth scenario resulted R (0.883), NSE (0.405), and RSR (0.770). Scenario five R (0.878), NSE (0.481), and RSR (0.719). Sixth scenario R (0.878), NSE (0.307), and RSR (0.831).

Van Genuchten equation (Eq. 3), Figure 6 illustrated the result of scenario one, which used parameters fully estimated by VGW in the CalcPTF program, as previously listed in Table 5, with satisfactory accuracy for R (0.791), NSE (0.612), and RSR (0.624). The second and third scenarios improved the results for R (0.892 and 0.864), NSE (0.717 and 0.677), and RSR (0.531 and 0.567), respectively. Scenarios four and six showed a similar accuracy criteria parameter for scenario one deuteriations in result quality compared to scenario one R (0.757 and 0.840), NSE (0.516 and 0.610), and RSR (0.695 and 0.623), respectively. Scenario five improves the accuracy with a rating equal to good level R (0.816), NSE (0.662), and RSR (0.581).

Fig.  5

Statistical indicator values for Brooks and Corey equation using Saxton et al., 1986 (BCS) for the different scenarios, correlation coefficient (R), Nash–Sutcliffe efficiency (NES), ratio of RMSE to the standard deviation (RSR). (Sample 36).

Fig.6: Statistical indicator values for van Genuchten equation using Wӧsten et al., 1999 (VGW) for the different scenarios, correlation coefficient (R), Nash–Sutcliffe efficiency (NES), ratio of RMSE to the standard deviation (RSR). (Sample 36).

Estimation error effect of RMSE: Brooks and Corey's equation showed a high sensitivity when varying the value of ϕ as shown in Figure 7. An error in the estimation of -30% or +30% increased the RMSE of SWCC estimation equal to 1383% and 1250%, respectively. λ parameter came as less in influence with a vast difference than ϕ, so the effect of the same range of error will increase the RMSE for SWCC estimation by 518% and 327%, respectively. For the same error range, both θr and ψb had a minor impact compared to previous parameters (117%and 132%) and (188%and 212%), respectively.

An evaluation of the sensitivity of SWCC estimation to parameter error (-30% to +30%) is shown in Figure 8. SWCC estimation RMSE values increased by approximately 7614% to 1730% for n and approximately 2564% to 2500% for θs. The other parameters θr and α showed smaller effects with RMSE (408 to 388) and (479.80 to 3986.68), respectively, for the same standard error range.

DISCUSSION

            When no measured data are available, the output of a PTF may be used as input to other functions. This can positively or negatively affect the degree of uncertainty in the estimation depending upon the level of error propagation and sensitivity of inputs to the PTF outputs (Benke et al., 2018; Gunarathna et al., 2019).Though PTF modeling and data extrapolation are continually improved, they are seldom error-free or completely accurate. The natural variation in soil properties can lead to incorrect results from models (Brown and Heuvelink, 2005; Leenhardt, 1995; Minasny et al., 1999).This introduction is required as a startup to discuss the above results. The means of soil particle size components were shown that samples in this study were within sandy loam (SL) texture and zero percentage of organic carbon (OC). As reported early by Al-Saeedi (2022), the low percentage of clay eliminated any significant effect of clay on the main hydraulic properties. He also showed high sand and silt percentage relations to the main soil properties, i.e., θ_s and ρ. This is in contrast with most of the PTF research and particularly in the CalcPTF program, whereas the clay and OC are the major estimator variables in both continuous or discrete PTFs (Chung, 2021; Guber et al., 2010; Nguyen, 2016; Nguyen et al., 2017; Rawls and Brakensiek, 1982; Vereecken et al., 1989; Wösten et al., 1999; Zhang and Schaap, 2017)

Fig.  7

Effect of parameters percent of bias (PB) on the SWCC estimation accuracy of Brooks and Corey equation.

Fig. 8: Effect of parameters percent of bias (PB) on the SWCC estimation accuracy of van Genuchten equation.

            BC parameters estimation, in all six PTF models (BCS, BCC, BCR, BCW, BCO, and BCM), the estimation was extremely inaccurate, with a very high deviation for all parameters (ϕ,θ_r, ψ_b, and λ).

For ϕ, despite a high correlation coefficient R (0.863 - 0.896), the other statistical criteria revealed an atrocious result. All models displayed negative NES values (<0.5), ranging from -0.403 to -0.10. RSR values were > 0.7, and the PB values ranged between -32.170-2.791; thereby, all suffered from dire performance results annulled the validity of models. As porosity (ϕ) equals saturation (θ_s), BCS used a multiple regression equation with variables sand% and clay%. He used the correlation between groups rather than within groups. This approach increases the potential for variation in the estimation models within the groups themselves. The correlation between groups helps draw directions rather than estimation (Marzban et al., 2013). PTFs (BCC, BCR, and BCO), CalcPTF used the porosity method equation (11) to estimate ϕ, where they assumed ρ_s equaled 2.65g cm^-3, while the value of ρ_s in this study is varied from 1.901 g cm^-3 to 2.412 g cm^-3.

                                [11]

Porosity method equation is reported by many researchers (Khoshkroudi et al., 2013; Mbagwu and Okafor, 1995; Vereecken et al., 1989) as a poor tool for estimating either θ_s or Φ, thus consequently propagating the errors in other related parameters and SWCC estimation PTFs.BCW multiplied Eq. 11 with a factor of 0.93, Williams et al. (1992) used samples dominated by a clay texture where he found  40% of the effect on ϕ and θ_s came from clay, which is not the case in this study (Al-Saeedi, 2022). Mayr and Jarvis's(1999) shown that BCM used an over parameterized multi-linear regression equation, making errors very likely to occur with any small deviation from the mean of the soil texture group, a similar finding described by (Weynants et al., 2009).

θ_r assumed equal zero in (BCS, BCC, BCW, BCO, and BCM). At the same time, in BCR  (Rawls and Brakensiek, 1985), he used a multiple linear regression equation over parameterized with about 12 betas (variables) while clay percentage was the most effective variable at θ_r. This error could be attributed to the effect of clay, not ϕ (Abdelbaki, 2021; Castellini and Iovino, 2019; Karim and Fattah, 2020).

All statistical measurements of h_b showed the invalidity of any of the listed PTFs 1-7 in Table 1. Table 3 demonstrated a high RMSE (62.62 – 69.308), NES negative less than zero, RSR unsatisfactory with values above 1.0, and PB with a high bias under estimation reached 92.330% in the BCO model. These catastrophic results were caused by the approaches used to construct the original PTF models. BCS model estimated ψ_b based on a doubtfulness θ_s value as eq.110 (already discussed). BCC model used the geometric mean particle diameter, geometric standard deviation, and ρ_b. Williams et al. (1992) reported that using geometric techniques in estimating was invalid with his Australian and UK soil samples. BCR model applied an over fitted multi regression equation with 15 parameters including ϕ, Sand%, and clay% with different forms. They used ϕ equal to θ_s from equation 10, which is already discussed as a major source of error. BCW model used θ_s, ρb, clay%, and fine sand%. The model is over fitted. It was built based on Australian soil samples with high clay percent (clay> 40%) (Williams et al., 1992). θ_s was estimated using soils with high clay and a fixed value of ρ_s equal to 2.65 g cm^-3, which could be another source of errors in this model (Dai et al., 2013).With BCO, ψ_b was exerted from the original equation, which was calculated θ as a function of ψ by using ρ_b, clay%, sand%, and D (mean depth). Also, the model assumed the moisture at h_b is near saturation. These factors sabotage the accuracy of the model for other soils (Guber et al., 2009, 2006; Nasta et al., 2021; Oosterveld and Chang, 1980). BCM model, clay soils dominated the sample population, were obtained using backward stepwise multiple regression, including bulk density and organic carbon(Dai et al., 2013; Nasta et al., 2021).

λ pore distribution index (Table3) showed very low and unsatisfactory results for all statistical parameters for all PTF models. BCS used a model built based on the correlation between groups (n=44) with the principal role of clay. For BCC, he used the geometric measurements in his estimation model, which was already criticized by (Williams et al., 1992). For BCR, they again over fitted their model with 12 parameters. BCW, as he used clay soils in his non-linear multi regression equations, he assumed θ_r equals zero at ψ equal 10⁴ bar and θ equal θ_s at ψ_b, also he used ρ_b, clay%, and fine sand% (Williams et al., 1992), so the tow assumptions were not the case in this soil study. BCO used in his model had a fixed value of λ equals 0.190, which is incorrect when applied to all soils with different texture types. BCM model was over fitted with seven parameters, including organic carbon percentage.

VG parameters estimation only three PTF models (VGW, VGVA, and VGVE) related to equation (3). At the same time, the other four PTFs (VGT, VGR1, VGG, VGRA, and VGR2) were originally derived from discrete models. Both model approaches poorly estimated all the parameters (θ_s, θ_r, α, and n) with very high uncertainty.   

Table 4 shows a high correlation coefficient R of θ_s ranging between 0.559 to 0.942 with an RMSE ranging between 0.075 to 0.143 for all models. Despite these lucrative numbers, the other statistical criteria measurements, NES and RSR, exhibited unsatisfactory values, meaning the models were invalid in estimating θ_s. The continues PTFs (VGW, VGVA, and VGVE). VGW, due to the high sand content, the CalcPTF program used tabulated parameters (θ_s, θ_r, α, and n) to represent the average soil hydraulic properties for 11 soil texture classes based on the geometric mean. This approach led to tremendous errors in the estimation models (Abbasi et al., 2011; Dai et al., 2013; Nasta et al., 2021; Weynants et al., 2009). VGVA and VGVE  PTFs were based on the dominant clay content and relative high ρ_b (Esmaeelnejad et al., 2015; Khoshkroudi et al., 2013; Weynants et al., 2009; Xu et al., 2021).CalcPTF program assumed θ_r equal to zero in VGVA, while the other parameters α and n estimation based heavily on clay content and ρ_b, which are the main cause of deviation and errors in the estimation process in this study soils as were proved in early studies (Dai et al., 2013; Tomasella and Hodnett, 2004; Weynants et al., 2009).

Using the parameters of eq. 2 in the discrete PTFs measures the accuracy level of these equations (VGT, VGR1, VGG, VGRA, and VGR2). As shown in Table 4, θ_s, θ_r (at 1500kpa), α, and n the performance of all five PTFs was very poor, reflecting the invalidity and high uncertainty of this PTFs SWCC estimation.

            SWCC estimation, this section included the SWCC estimation accuracy and parameters sensitivity and contribution to the accuracy criteria parameters. Despite the erroneousness of the equations 1-3 parameters for all CalcPTF models but incorporating these parameters together generates sufficient confidence in reproducing adequate estimation similarity or parameters equifinality. This effect could be related to most of these models being affected by their mathematical form rather than by their parameters' physical significance(Du, 2020; Khatami et al., 2019).However, models with more parameters are always preferred in the SWCC estimation models. So, van Genuchten (eq. 3) performed a higher accuracy than Brooks and Corey equation (eq. 1 and 2) in all PTFs models. These results supported by many prior works (Ferreira et al., 2012; Matlan et al., 2014; Weihermüller et al., 2021).Replacing one or more inputs with a measured or supremacy estimated parameter showed a high enhancement to the final result of both SWCC models. It increased the effect and relevance of the physical form. This was endorsed by many SWCC estimation model creators (Rawls and Brakensiek, 1982; Saxton et al., 1986). This study showed the supper role of saturation θ_s input on the quality of SWCC estimation outputs compared with other parameters. BCS and VGW showed high enhancement in NSE and RSR by applying either measured or estimated (eq.10) θ_s (Fig. 8). The significant enhancement in the model’s outputs was attributed to the unique role of θ_s. So, the improvement of θ_s presentation, either by implying measured or well-estimated value, will lead to a magnificent improvement in the SWCC estimation outputs, as already shown in many prior studies(Mohajerani et al., 2021; Rajkai and Varallyay, 1992; Rawls and Brakensiek, 1982; Saxton et al., 1986; Vereecken et al., 2010).The sensitive analysis for both model parameters s revealed a symmetric effect of Brooks and Corey's (1964) equation with high sensitivity for θ_s and insensitive for ψ_b and λ. On the other hand, van Genuchten's (1980) equation parameters exhibited a symmetric effect except with n, where it showed the underestimation resulted in a more severe effect than overestimation. The sensitivity analysis output emphasized the substantial role of having a good measurement or estimation of θ_s over other parameters, as discussed earlier by other articles (Mohajerani et al., 2021; Vereecken et al., 1989).