2.4. Statistics and data evaluation
The HYDRUS-1D outputs were used to analyse the fluxes of the soil water balance, i.e. actual evaporation (Ea ), actual evapotranspiration (ETa ), drainage (D ), and runoff (R ). For the data of the 13 different PTFs the cumulative flux was selected and the arithmetic mean, 0.05 and 0.95 percentiles (which spans the 90 percent tolerance interval), and the 0.15 and 0.85 percentiles (which spans the 70 percent tolerance interval), were calculated for the last data point of the cumulative fluxes (tend = day 10988) for each soil textural class. In a next step, outliers based on the percentiles were calculated and flagged according to Eq. [8] and [9]:
\(\text{If}\ \text{value}\ \text{of}\ a\ \text{PTF}\ \text{for}\ \text{flux}\left(x\right)>0.95\ \text{percentile},\ \text{value}=1,\ \text{else}\ \text{value}=0\)[8]
\(\text{If}\ \text{value}\ \text{of}\ a\ \text{PTF}\ \text{for}\ \text{flux}\left(x\right)<0.05\ \text{percentile},\ \text{value}=-1,\ \text{else}\ \text{value}=0\)[9]
where flux (x ) is the cumulated flux (e.g., actual evapotranspiration) at tend for each individual PTF used to simulate the soil class for each model scenario. In other words, if the flux value, flux (x ), attend exceeds (is less than) the percentile span calculated, the simulation was flagged with a 1 (-1), whereas if the flux value, flux (x ), at tend lies within the given percentile, the simulation was flagged with a 0. The same procedure as for the 90 percent tolerance interval was repeated for the 70 percent tolerance interval. In order to present variability of the hydraulic parameters or simulated fluxes, the coefficient of variation (CV) was calculated, relating the standard deviation to the mean value. The CV was expressed as a percentage.
For the analysis of differences in soil hydraulic properties estimated by the PTFs, the comparison of two group means was performed with the non-parametric Wilcoxon rank sum test using Matlab® ranksumfunction using the probability p = 0.05. Principal component analysis (PCA) was performed, also with Matlab®.
To help interpret the differences between model runs with regards to the water fluxes, the matric flux potential, MFP , the characteristic length, LC , the sorptivity, S , and characteristic time tgrav was calculated as explained next.
The matric flux potential MFP [cm2d−1] is a convenient bulk soil hydraulic property that is often used in soil water movement studies (e.g., Raats, 1977; Pullan, 1990; Grant & Groenevelt, 2015), which is defined as the integral of the hydraulic conductivityK (h ) [cm d−1] over the pressure head h [cm] starting at an arbitrary reference pressure head href :
\(\text{MFP}=\ \int_{h_{\text{ref}}}^{h}{K_{(h)}\text{dh}}\) [10]
where href was set to permanent wilting point ath = -15,000 cm according to Pinheiro et al. (2018) and hset to full saturation (h = 0).
As a second soil physical feature, the characteristic length of bare soil evaporation, LC [cm], was calculated. According to Lehmann et al. (2008 and 2018), LCis the maximal extent of the capillary flow region to supply water to evaporating surface. LC is determined by the range of capillary pressure between large and small pores driving the capillary flux against gravity (expressed as head difference, denoted as gravity length LG ) and the hydraulic conductivityKeff of the supply region. Formally,LC is defined via:
\(L_{C}=\frac{L_{G}}{1+\frac{E_{0}}{K_{\text{eff}}}}=\frac{\frac{\left(1-m\right)}{\alpha}\left(1+\frac{1}{m}\right)^{(1+m)}}{1+\frac{E_{0}}{4K(h_{\text{crit}})}}\)[11]
where α and m are the van Genuchten parameters used in Eq. [4] and E0 is potential evaporation rate. To calculate effective conductivity, Keff was estimated as 4Kcrit (Haghighi et al., 2013; Lehmann et al., 2018), whereby Kcrit is the unsaturated hydraulic conductivity at critical water content and capillary pressure when capillary pathways start to disconnect (Lehmann et al., 2008). The critical capillary pressure and gravity length are determined based on linearization of the soil water retention curve. For the van Genuchten formulation used in Eq. [11], the linearized retention curve consists of the tangent to the inflection point, andLG and hcrit can be expressed analytically. For Brooks and Corey, the values are determined numerically with a line passing through the air-entry value (Se = 1, h = hb ) and a particular point on the retention curve that is closest to (Se = 0, h = hb ).
The soil sorptivity, S [cm d0.5], which is defined as a measure of the capacity of a porous medium to absorb or desorb liquid by capillarity, was calculated assuming a soil column with uniform initial water content and infinite length, following the approach of Parlange (1975) as described in Moret-Fernández et al. (2017), Lotorre et al. (2018), and Rahmati et al. (2019):
\(S^{2}(\theta_{s},\theta_{i})=\int_{\theta_{i}}^{\theta_{s}}{D(\theta)\left[\theta_{s}+\theta-2\theta_{i}\right]d\theta}\)[12]
where D (θ) [cm2 d-1] is the diffusivity defined by Klute (1952) as:
\(D\left(\theta\right)=K(\theta)\frac{\text{dh}}{\text{dθ}}\)[13]
For the initial water content, θi , the maximum reported residual water content (θr ) for all MvG parameters used in this study (0.192 cm3cm-3) was used for all PTFs based on MvG and 0.12 cm3 cm-3 for all PTFs based on BC.
Finally, the so-called characteristic time (tgrav ) was calculated according to Philip (1957), which determines the ‘time’ t where gravitational forces become dominant (ttgrav ), while forttgrav capillary forces remain dominant over gravitational forces (Rahmati et al., 2020):
\(t_{\text{grav}}=\left(\frac{S}{K_{s}}\right)^{2}\) [14]
Two final related soil properties that were calculated were the characteristic time for the attainment of field capacity FC , τFC (d), and the elapsed time required for attainment of FC , denoted as tFC (d]), see Assouline and Or (2014). To do so, the effective soil saturation at field capacity SFC [Eq. 15] and the unsaturated hydraulic conductivity at field capacityK (SFC ) [Eq. 16] were calculated by:
\(S_{\text{FC}}=\left[1+\left\{\left(\frac{n-1}{n}\right)^{(1-2n)}\right\}\right]^{\left(\frac{1-n}{n}\right)}\)[15]
\(K_{\text{FC}}={K_{s}S_{\text{FC}}^{0.5}\left[1-\left(1-{S_{\text{FC}}}^{(1/m}\right)^{m}\right]}^{2}\)[16]
Moreover, the quantity of drainable water from the soil profile to depthz at field capacity FC , QFC , expressed as equivalent water depth (dimensions of length), has to be calculated as:
\(Q_{\text{FC}}=z(\theta_{s}-\theta_{r})\left(1-S_{\text{FC}}\right)\)[17]
Here, z was set to 30 cm to match the maximal rooting depth for convenience, as z only scales with totalQFC .
Consequently, a characteristic time [d] for the attainment ofFC , τFC , can be deduced from the ratio of these two quantities, QFC andKFC :
\(\tau_{\text{FC}}=\frac{Q_{\text{FC}}}{K_{\text{FC}}}\) [18]
The drainage dynamics can now be linked to the elapsed time [d] required for attainment of FC , denoted astFC :
\(t_{\text{FC}}=\ \frac{z\left(\theta_{s}-\theta_{r}\right)}{K_{m}}\ln\left(\frac{K(S_{\text{FC}})}{\text{Ks}}\right)\)[19]
where Km is the effective hydraulic conductivity that represents the mean value of K (SFC ) weighted by SFC (representing the available relative cross section of flow):
\(K_{m}=\frac{\int_{0}^{1}{S_{\text{FC}}K(S_{\text{FC}})dS_{\text{FC}}}}{\int_{0}^{1}{S_{\text{FC}}d}S_{\text{FC}}}\)[20]
Results and Discussion
Predicted hydraulic parameters and hydraulic functions
In a first step, the retention and hydraulic conductivity curves for the 13 PTFs and the 12 USDA soil classes were plotted based on the estimated soil hydraulic parameters listed in Annex Tab. 1 and 2. As an example, the retention and hydraulic conductivity curves for the USDA textural class sand is plotted in Fig. 3 (see Annex Fig. 1 for all retention and hydraulic conductivity curves of all USDA textural classes). Figure 3 shows that the retention curves based on the 13 PTFs greatly differ along the entire pressure head range. For these curves, the saturated water content, θs , varies between 0.472 for Rawls MvG and 0.375 cm3 cm-3 for Cosby SSC with a mean of 0.417 cm3 cm-3 over all PTFs. The corresponding coefficient of variation (CV) is 7.5 % for the sandy soil. The residual water content, θr , varies between 0 for those PTFs setting θr to 0 such as Weynants, Clapp&Hornberger, Cosby SC, and Cosby SSC, to 0.061 cm3 cm-3 for the ‘Toth class’ PTF with a mean of 0.029 cm3 cm-3 and a CV of 80.1%, indicating a much higher variability in terms of CV inθr compared to θs . An even larger variability can be found in the saturated hydraulic conductivity,Ks , for which we found a maximum value of 1520.6 cm d-1 for Clapp&Hornberger and a minimum of 8.3 cm d-1 for the ‘Toth class’ with a mean of 315.8 cm d-1 (CV = 129.2 %). In general, the smallest CV values (data not shown) for all USDA textural classes were found forθs , with lowest ranging from 2.4% for the silty clay loam class to CV = 7.5 % for the sand class. (Larger variability was observed in θr with the lowest coefficient of variation in the loamy sand class (CV = 76.3 %) and highest in the silty clay loam class (CV = 106.6 %). As expected,Ks showed the largest variability, with lowest CV in the loam class (91.7 %) and largest in the clay loam class (215.2 %). The larger CV values for the Ks estimation is not surprising as a large uncertainty in predictedKs has been already widely reported (e.g., Jaynes & Tyler, 1984; Ahuja et al., 1985; Tietje & Hennings, 1996; Schaap et al., 1998). Furthermore, it has to be noted that the PTF developed by Weynants et al. (2009) predictsKs* instead ofKs , whereKs* is a hydraulic conductivity acting as a matching point at suction head h = -6. Therefore, some slightly lower Ks (hereKs* ) value will be predicted by Weynants’ PTF. On the other hand, there seems to be a clear grouping among the class PTFs, with regards to the estimation ofKs . Clapp&Hornberger predicted the highest values for six classes (sand, loamy sand, sandy loam, loam, silt loam, and silty clay loam), followed by the PTF of Woesten for five soil classes (sandy clay loam, clay loam, sandy loam, silt loam, and clay). Even more pronounced is the picture for the prediction of lowestKs , whereby the PTF of Rawls MvG estimated the smallest Ks values for 11 soil textural classes, except for sand, whereas the ‘Toth class’ PTF showed the lowestKs values. Unfortunately, the α and n (or 1/b ) values cannot be directly compared between the BC and MvG approaches, as both parameters have a slightly different physical meaning.
Numerical model performance
As numerical stability of the simulation is one of the crucial aspects in the choice and application of the PTF, especially for large scale modelling, we analysed each PTF with respect to numerical convergence, when using HYDRUS. For each PTF and the seven model scenarios (see Fig. 2), 44 individual model runs were performed: for each PTF, the three homogeneous soil layer model scenarios were modelled for each soil textural class (these are 36 model runs). In addition, the four layered configurations are run for a coarse and a fine soil layering, resulting in eight model runs; hence, a total of 44 model runs per PTF were obtained. Note that for the Clapp and Hornberger (1978) PTF, only 41 model runs were performed as no parameters were reported for the silt class. For 486 model runs, out of the total 569 model runs (i.e., 85 %), convergence was achieved. A total of 184 out of 217 (85%) of the model runs for the BC and 279 out of 352 (79 %) for the MvG parameterization converged, even though it has been reported that the BC type function sometimes prevents rapid convergence and might therefore cause numerical problems. This was deemed to be caused by the discontinuity present in the slope of both the soil water retention and unsaturated hydraulic conductivity curves (van Genuchten, 1980).
For those cases where the simulation did not converge for the MvG parameters, the air entrance value (the inverse of α ) was set to -2 cm, and the model was rerun. This procedure increased the total number of converged MvG simulation runs to 302 (86 %), which is a similar percentage to that obtained for BC. Looking at individual PTFs (see Tab. 3), we can see that the Rosetta SSC+BD and Cosby SSC seem to be numerically very stable with 100 % converged runs. The Woesten and ‘Toth continuous’ function converged for > 95 % of the runs, after setting the 1/α to -2 cm. On the other hand, the lowest convergence was found for the Rawls MvG and Rawls BC with 43 and 39 %. Unfortunately, using 1/α = -2 cm did not improve convergence for Rawls MvG.
The reason why some PTFs prohibited the HYDRUS model from converging is quite apparent for some cases. For example, Rawls MvG and Rawls BC yielded very low Ks values of 0.8 cm day-1 for the loam and sandy clay class and ≤ 0.3 cm day-1 for clay loam, silt, and silt loam class. Extremely low values were obtained for silty clay loam (0.04 cm day-1) and clay (0.004 cm day-1), almost allowing no infiltration at all. Another extremely lowKs value was predicted by the ‘Toth class’ PTF for silty clay, with 0.01 cm day-1; these unrealistically low values again led to numerical instabilities.
It has to be noted that the reported convergence here is only valid for the numerical model (HYDRUS-1D) used in this exercise with the given numerical (convergence) default criteria, vertical discretization and temporal resolution, and atmospheric boundary conditions. The performance of these PTFs may change if a different numerical scheme, e.g. solving the Richards equation in the mixed or diffusivity form were used, or a different spatial discretization and/or temporal resolution. Furthermore, the lack of certain processes in our simulations (e.g. coupled heat and water transport or evaporation from the wet canopy), or the nature of the atmospheric forcings (e.g. a difference in rainfall frequency and amount) will affect the likelihood of convergence.
Fluxes and outliers
Simulated fluxes over time
Firstly, the simulated cumulative fluxes were analysed.ETa (vegetated surface) orEa (bare soil) is a key flux as it indirectly contains information of the net infiltration into the soil profile (net daily infiltration = daily sum of precipitation – daily sum ofETa or Ea ), deep drainage (over long-run), and plant available water in the root zone. Furthermore, ETa or Eadetermines the return of water from the soil profile to the atmosphere, and as such affects the land surface energy budget. CumulativeETa or Ea data for each scenario/soil class combination was plotted and the arithmetic mean of all data (model ensemble mean, MEM) for each combination, as well as the spread of the data, was calculated by the 70 and 90 percent tolerance interval according to Eq. [8] and [9].
As an example of the high variability in simulated fluxes, the simulated cumulative Ea , for the homogeneous bare soil scenario of loamy sand texture, over the entire simulation period of 30 years, that ends on day 10988 (tend ), is plotted in Fig. 4a. There is a large variability between the various simulations based on the 13 PTFs. MEM at tend is 1692 cm (564 mm year-1). The smallest simulated cumulative Ea was 1273 cm (424 mm year-1) for the Carsel&Parrish PTF, and largest, with 2043 cm (681 mm year-1), for Weynants. The difference of 257 mm year-1 between the largest and smallest simulated Ea , and their deviation of 140 and 117 mm year-1 from MEM clearly indicates that the choice of PTF substantially affects the estimation of theEa for this soil class. In contrast, low variability was found for cumulative Ea for the bare homogeneous clay loam (see Fig. 4b). Notably, two out of the 13 simulations did not converge (Rawls MvG and Rawls BC), which potentially also impacts the variability in simulated fluxes. Nevertheless, for the remaining 11 simulations the lowest simulated flux was 1744 cm (581 mm year-1) for Rawls class PTF and largest for Weynants PTF with 2041 cm (680 mm year-1) (for this soil, MEM = 1893 cm or 631 mm year-1). Overall, the difference between the largest and smallest flux is only 99 mm year-1, i.e., 2.5 times smaller than the difference found for the loamy sand. As Ea will be also be influenced by the precipitation entering the soil (total precipitation over 30 years = 2479.7 cm (827 mm yr-1)), we also looked at the cumulative runoff. For most soil textural class/PTF combinations, still for the homogeneous bare soil scenarios, runoff is low or negligible, with zero runoff, or values < 1 cm over 30 years, for 121 model runs, which is equivalent to 88 % of runs. Nine simulations (7 %) returned a runoff >1 cm but <10 cm, and eight exceeded 10 cm over 30 years (7 %). The highest cumulative runoff was generated for the Rawls MvG/silt combination with 675 cm (27 % of total precipitation) followed by Rawls silt loam with 664 cm and Rawls MvG loam with 389 cm. These three combinations have also been flagged as outliers of the lower 0.15 percentile for total cumulative evaporation attend , which can be explained by the fact that less water enters the soil, and therefore less will be evaporated. The same holds for the Carsel&Parrish PTF and silty clay loam (runoff = 135.1 cm) and sandy clay (runoff = 70 cm) combinations. On the other hand, the combination ‘Toth class’ PTF/silt loam generated 122.9 cm runoff but was classified as an upper 0.95 percentile outlier, generating more Ea . Finally, Rawls BC/silt and ‘Toth continuous’/silty clay combinations generated runoff of 156.1 and 43.7 cm, respectively, yet are not classified as outliers. In general, runoff generation is linked to low Ks values (see Annex Tab. 1 and 2). An overview of all cumulativeEa fluxes at tend for the bare soil scenarios is plotted in Fig. Annex 2.
Our findings with regards to ETa for the vegetated scenarios (grass and wheat), still with a homogeneous soil profile (Annex Fig 3 to 4), were comparable. For some soil classes, such as clay, clay loam, and silty clay loam, variability between PTFs was low, for both grass and wheat, whereas the ETafor the sandy and sandy loam soils showed consistently high variability. In contrast, ETa for the loamy sand class exhibited relatively high variability for grass (as was also the case for bare soil) and a slightly smaller one for the wheat scenario configuration. For the other soil textural classes, the picture is less clear. Again, as was the case for the bare soil, there is a substantial number of soil class/PTF combinations that result in runoff. A slightly larger, compared to the bare soil scenario, percentage of simulations with runoff > 1 cm was found for the grass (18 %) and the wheat (20 %) scenarios. Moreover, maximum runoff att end value increased from bare (675 cm for Rawls MvG silt) via grass (859.2 cm for Rawls MvG silty clay) to the wheat scenario (999.2 cm for Rawls BC silty clay). Surprisingly, eight out of twelve soil class/PTF combinations yielding runoff > 100 cm were not flagged as outliers for the 90 % tolerance interval for the grass and four out of 10 for the wheat. There are some unexpected findings, namely that the simulation for the Toth continuous’ PTF yielded 46.1 and 11.7 cm runoff, respectively, for the silty clay and silty clay loam under wheat vegetation, despite the fact that theETa flux at tend was flagged as an outlier of the upper 0.95 percentile, indicating relatively high evaporation with respect to the model ensemble.
Finally, the simulation for the scenario of sandy loam overlying silt loam and loamy sand plotted in Annex Fig 5 showed much lower variability for Ea compared to Ea of the homogenous profile with the texture of the uppermost layer (silt loam in Annex Fig. 2). This indicates that soil layering will reduce the effect of the choice of PTF on the cumulative evaporation. This holds true even more for the layered bare soil scenario where silt loam overlies silty clay loam that is overlying silty clay. Again, the variability in Ea for the layered system is much lower than that of the homogeneous silt loam, that forms the first layer in the vertically heterogeneous soil profile. Besides, it can clearly be seen that when vegetation is introduced, variability increases slightly, which is reflected in the coefficient of variation (CV) of the flux attend , where for the layered profile topped by sandy loam the CV increased from 3.9, via 5.1 to 4.9 % for the bare, grass, and wheat vegetation scenario. For the profile with the first layer consisting of silt loam, CV values were 0.5, 3.7, and 3.7 % for the bare, grass, and wheat vegetation, respectively. Introducing a fluctuating ground water table increased the CV substantially to 13.6 % for both vegetated layered systems. Variability in simulatedEa or ETa for the layered scenarios can partly be explained by a large reduction in runoff. In total only two simulations (2 %) for the Carsel&Parrish (sandy loam topped layered profile under grass vegetation (290.2 cm) and sandy loam topped layered profile for the wheat vegetation and ground water fluctuation (302.6 cm)) exceeded runoff of 100 cm. A further four exceeded the runoff threshold of 1 cm (3 combinations for Carsel&Parrish and one for Rawls BC).
Overall, the choice of PTF substantially affects the simulated values ofEa or ETa for most soil classes, irrespective of the fact whether the soil was bare, where the water (vapour) can only leave the soil column via the pore-space at the soil surface, or vegetated, where a considerable proportion of the water being returned to the atmosphere consist of water taken up from the deeper rooted parts of the soil profile.
Outliers per scenario
As shown above, substantial variability in simulatedEa or ETa fluxes occurred for different PTFs and model scenarios. The fluxes exceeding the 70 or 90 % tolerance intervals, respectively, were marked as outliers and calculated for each scenario and soil class according to Eq. [8] and [9]. The number of outliers were counted for each scenario individually in a first step.
The number of outliers varies greatly between PTFs with regards toEa fluxes for the homogeneous bare soil (see Fig. 5a); these fluxes were shown in Figs. 4a and b. Naturally, more outliers are detected for the 70 than for the 90 % tolerance interval. For this scenario, Rosetta SSC, Weynants, and ‘Toth class’ exceed the upper 0.95 percentile, whereby Weynants exceeded this percentile for all soil classes where the model had converged, except for silt and silt loam. Rosetta SSC exceeded the upper 0.95 percentile for clay, whereas the ‘Toth class’ PTF exceeded it for silt and silt loam, respectively. Looking at the lower 5 percentile, Carsel&Parrish PTF exceeded this threshold for eight soil classes, and further outliers were found for Rawls MvG (N = 6) and Rawls BC class (N = 4). Two outliers were calculated for Cosby SSC, Rawls BC, and one for Rosetta SSC and Woesten PTFs.
Finally, only Rosetta SSC+BD, ‘Toth continuous’, and Cosby SC indicate no outliers for the upper and lower 70 % tolerance interval.
As some simulation runs did not converge (see discussion above), the comparison in terms of total number of outliers is limited. Therefore, the total number of outliers was normalized to the number of converged simulations for each scenario and PTF combination. Again, we present the relative number of outliers for the homogeneous bare soil profile simulations in Fig. 5b as an example (all others are shown in Annex Fig. 6 to 8). Here, the Weynants PTF shows the largest percentage of outliers for the upper 0.95 and 0.85 percentile with 82 and 100 % outliers, respectively. For the ‘Toth class’ PTF, we found 20 and 50 % outliers for the upper 0.95 and 0.85 percentile, respectively, indicating that also this PTF simulated larger fluxes with respect to the ensemble. On the other hand, Rawls MvG shows the largest percentage of outliers at the lower end (86 % for the 0.15 and 57 % for the 0.05 percentile) followed by Carsel&Parrish PTF with 73 % for the 0.15 and 45 % for the 0.05 percentile. However, Rawls BC and Rawls class also show substantial percentages of outliers for the 0.15 percentile. By comparing the relative (converged only) and absolute (all runs) number of outliers, it can be seen that despite equal or even lower or higher absolute number of outliers for different PTFs, the relative numbers differ due to non-converged simulation runs for some PTFs. For instance, ‘Toth class’ for the 0.95 percentile showed 2 outliers yielding 20 % relative outliers as two simulations (silty clay and silty clay loam) did not converge, whereby 1 outlier for Rosetta SSC yielded only 8 % relative outliers as all simulations converged.
As there is no clear trend in the analysis of the absolute or relative outliers for the individual scenarios (see Fig 5b and Annex Figs. 6b to 8b) which PTF generates most outliers, from here on the outliers over all scenarios for all soil textural classes were calculated for converged simulation runs only and expressed in relative terms. Figure 6a shows the outliers of the 90 % tolerance interval (sum of upper and lower outliers), combined for all textural classes, for the seven scenarios for the 13 PTFs for Ea andETa at tend . In this figure, the PTFs of the two main hydraulic formulations are clustered: those based on the Mualem van Genuchten (MvG) on the left and those based on Brooks Corey (BC) formulation on the right. Furthermore, two lines are added, dividing the results into three groups: i) those PTFs with relative number of outliers < 10 %, classified as ‘robust’, ii) those PTFs with 10 % ≥outliers ≤ 20 %, classified as ‘intermediate robust’, and iii) the PTFs with relative number of outliers >20 %, classified as ‘non-robust’. It has to be noted that these thresholds (10 and 20 %) were chosen arbitrarily, but may help to formulate the final recommendations for the choice of preferred PTF, to be used in land surface models, for example.
This classification shows that the Rosetta SSC, Rosetta SSC+BD, Woesten, ‘Toth continuous’, Rawls BC, Rawls class BC, Clapp&Hornberger, Cosby SC, and Cosby SSC PTFs are located below the 10 % threshold for the 90 % tolerance interval, and can be therefore classified as ‘robust’ with respect to the ensemble behaviour (spread). Interestingly, all PTFs using BC formulation show low relative numbers of outliers below 10%. Woesten PTF did not show any outliers at all, indicating that this PTF is very robust with respect to the PTF ensemble used. On the other hand, the ‘Toth class’ PTF was classified as intermediate robust, and three PTFs (Carsel&Parrish, Rawls MvG, and Weynants) were classified as non-robust, whereby Rawls MvG produced most outliers (32 %).
The results for the 70 % tolerance interval are shown in Fig. 6b and followed the same approach as for the 90 % tolerance interval discussed above. Four PTFs are characterised as robust (Rawls BC, Clapp&Hornberger, Cosby SC, and Cosby SSC). Again, all these four PTFs serve to produce parameters for the Brooks Corey hydraulic formulation. The intermediate robust grouping includes Rosetta SSC, Rosetta SSC+BC, Woesten, and ‘Toth continuous’, that provide parameters for the Mualem van Genuchten formulation. Finally, Carsel&Parrish, Rawls MvG, Weynants, ‘Toth class’, and Rawls BC class are those PTFs classified as non-robust. There are two class, rather than continuous, PTFs here, indicating that continuous PTFs are more likely to be robust. Also, the Weynants PTF was based on a relative small number of samples, for Belgium only.
Based on the results presented above, it can be concluded that the use of different PTFs results in different hydraulic properties that predict considerably different Ea orETa fluxes leading to different soil water contents in the root zone but also to differences in deep percolation (or ground water recharge). Furthermore, PTFs such as Carsel&Parrish, Rawls MvG and Weynants can be identified as systematically less robust. In contrast, others, such as Woesten or all PTFs using the Brooks Corey formulation (except Rawls BC class) seem to be robust with respect to the ensemble of PTFs used in this study.
To facilitate the identification of outliers, all outliers per PTF, scenario, and textural class combination were colour coded and plotted in Tab. 4. Again, Weynants overestimates Ea orETa fluxes (brown colour for dryer soil conditions) for nearly all textural soil classes except for clay, and silt. On the other hand, Rawls MvG shows underestimation (blue colour for wetter soil conditions) for loam and silt loam over all three homogeneous soil scenarios and for silt and sandy loam for two out of the three homogeneous soil scenarios. The Carsel&Parrish’ PTF, on the other hand, results in over- and underestimation, depending on soil class.
In the study of Zheng et al. (2020), who evaluated also a set of 13 PTFs using independent retention data (data not used in the development of the PTFs), the Carsel&Parrish PTF showed largest RMSE in predicted volumetric water content of the retention data points, which is in agreement to our findings that those PTF is characterized as less robust with respect to simulated fluxes. On the other hand, Weynants PTF was the one with lowest RMSE, whereas Weynants was characterized as less robust in our case. The reason why Weynants behaves differently between the study presented by Zheng et al. (2020) and our study, might be that only the retention characteristics were analysed by Zheng et al. (2020), whereas in the function sensitivity performed here, also the hydraulic conductivity function plays a virtual role. The reasons why Weynants is characterized as less robust is further discussed in the following sections.
3.3.3. Simulated spread with respect to scenario
We raised the hypothesis that differences (variability) in simulated fluxes from using different PTFs will be reduced with increasing model complexity. Increasing complexity was generated by introducing vegetation (grass or wheat), soil layering, or the assumption of a fluctuating ground water table, for the layered vegetated soil scenario only. As only the homogeneous scenarios (bare, grass, and wheat) used all soil classes, we restrict the analysis on these three scenarios.
For the analysis, again the simulated cumulative actualEa or ETa data attend was taken and the model ensemble mean (MEM) for Ea to ETa attend over all PTFs was calculated for each individual soil class and scenario. Based on the MEM value forEa or ETa , as well as the individual Ea or ETa value at tend for each model run, the % difference from the MEM (100/MEM*Ea @tend orETa @tend ) was calculated and visualized using boxplots in Fig. 7, where the red line indicates the median, the box indicates the 0.25 and 0.75 percentiles, the whiskers represent the most extreme data points not considered as outliers, and crosses represent the outliers (value is more than 1.5 times the interquartile range). From the boxplots, two types of information can be deduced: i), the variability of predictedEa or ETa over all PTFs for one soil class / scenario and ii), the change in variability (spread) resulting from a change in scenario complexity (bare, grass, or wheat vegetation).
In general, the largest variability in predictedEa or ETa was found for the bare soil conditions, which is most pronounced for the loam, loamy sand, sand, sandy clay, clay loam, and sandy loam class. Minor differences were found between bare and vegetated scenarios for the other soil classes. The silty clay soil class for the grass scenario showed the smallest overall spread between minimum and maximum predictedEa (or ETa ) with a value of 7 % (min = 97 % and max 104 %). On the other hand, the largest variability was found for the combination sandy soil/bare soil scenario with 53 % (min = 72 % and max 125 %). All spreads, throughout the 13 PTFs, for different soil classes and scenarios are provided in the final column of Tab. 4. Overall, bare scenarios show a mean spread of 30 %, whereby the grass and wheat vegetated scenarios have only 23 % spread over all soil classes. A possible explanation for the reduced spread in simulated Ea or ETa with increasing model complexity (in this case vegetation) is that for the vegetated profiles water is extracted from the rooted portion of the soil profile, whereas under bare soil the water can only leave the soil profile at the soil surface. In the latter case, differences in the soil hydraulic properties, especially in unsaturated hydraulic conductivity, which is highly variable (on the order of magnitudes) between PTFs, close to the surface will impact the Ea flux more substantially. As shown earlier, runoff will occur in both scenarios (bare and vegetated) and is even slightly larger for the vegetated scenario, and therefore, cannot explain the reduced variability.
Next, for the layered soil scenarios (bare, grass, and wheat, without fluctuating groundwater table) a clear reduction in the variability was observed, for both profiles (by sandy loam or silt loam). The bare and the vegetated scenarios showed nearly the same spread (mean 13.4 % for bare, 18.5 % for grass, and 15.5 % for the wheat). In general, the sandy loam overlaying silt loam and loamy sand showed always higher variability compared to the silt loam overlaying silty clay loam and silty clay, which is consistent to the finding that the sandy loam of the homogeneous soil profile also showed higher variability compared to the homogeneous silt loam scenarios.
Overall, the results indicate that adding vegetation reduces the variability in the simulated Ea orETa flux, even if runoff occurs more frequently. This conclusion also holds for adding more complexity in terms of soil layering, although the latter has to be regarded with some caution due to the low number of soil combinations selected for these model runs. However, taking into account that large portions of our global land surface is covered by vegetation, differences in predicted fluxes, as a result of differences in PTFs used to generate the hydraulic parameters, will most likely be smaller compared to an ‘unvegetated world’.
In contrast, adding a fluctuating ground water table to the layered wheat scenario greatly increased variability inETa flux, for both soil layering to 48 and 35 % for the sandy loam overlaying silt loam and loamy sand, and silt loam overlaying silty clay loam and silty clay, respectively.
3.3.4. Differences in instantaneous fluxes
Cumulative fluxes at tend will only provide long-term systematic under or overestimation, but will not provide information on how the instantaneous fluxes fluctuate compared to the MEM. Therefore, the instantaneous fluxes were also analysed. The same analyses as conducted for the cumulative fluxes were performed, i.e., calculation of the MEM and the 0.95 and 0.05 percentiles for time stepi , whereby i runs from day 1 to 10988. Secondly, the total number as well as the upper and lower percentile outliers were counted. As an example, the outliers of Ea for the sandy loam of the homogeneous bare soil scenario were plotted in Fig. 8, for the different PTFs. Carsel&Parrish PTF shows a substantial number of outliers for the lower 0.05 percentile (N = 2020 or 18 % of all days), indicating that for these days less water will evaporate and return to the atmosphere, which would have implications for the cloud forming processes of a numerical weather prediction or climate model if a LSM using this PTF were to be embedded within it. On the other hand, Weynants has 3053 outliers for the upper 0.85 percentile (28 %) but also a smaller number of outliers for the lower 0.05 percentile (N = 309 or 3 %), leading to larger Eaflux. A large number of 0.05 percentile outliers were also found for Rawls MvG (N = 2992 or 27 %), again combined with a lower number of upper 0.95 percentile outliers (N = 391 or 4 %). Cosby SC, Cosby SSC, and Rawls BC showed only low number of outliers (N< 10) for the upper and lower percentiles. Even though the model runs for which the hydraulic parameters were derived from Carsel&Parrish and Rawls MvG PTFs exhibit large numbers of outliers, both are not flagged as 90 % tolerance interval outliers when the cumulative flux at tend was analysed. This means that the non-flagged instantaneous Ea fluxes compensate for the lower fluxes determined as outliers in Fig 8, or that the outliers are close to the 0.15 percentile, which is reflected by the fact that the total sum of underestimated flux (outlier flux – flux for the lower 0.05 percentile for each outlier day) is low, amounting to 5.7 and 2.4 cm over the 30-
year period, respectively. Moreover, both PTFs show runoff exceeding a total of 1 cm in 60% (Carsel&Parrish) and 36 % (Rawls MvG) of all converged simulations, respectively. For the Rawls MvG the nine simulations with runoff even exceed the 100 cm threshold, with runoff ranging between 388.6 to 859.2 cm. Looking at all textural classes (data not shown) for the homogeneous bare soil scenario, 29 soil class / PTF combinations out of the total 151 do not exhibit any outliers at all for the instantaneous Ea flux. These outliers are clustered in three soil classes only (clay, silty clay, and silty clay loam). Interestingly, out of these 29 with zero outliers in instantaneous flux, five are flagged as outliers for the cumulative flux at tend (Rosetta SSC clay, Cosby SSC clay, Weynants silty clay and silty clay loam, as well as Carsel&Parrish silty clay loam), meaning that these PTFs over- or underestimate instantaneous Ea only very modestly, yet consistently throughout the simulation period.
The percentage of all 90 percent tolerance outliers (sum of upper and lower outliers) summed over all days for all three homogeneous soil scenarios (bare, grass, and wheat) for all soil classes and PTFs are provided in Tab. 4. Over all soil classes and PTFs, the bare soil scenario has the lowest total number of outliers (N = 119930 days or 6.5 % over all days and scenarios) followed by the homogeneous wheat configuration (N = 173961 days or 10.1 %) and the homogeneous grass scenario (N = 178249 days or 10.4 %). This finding is perhaps in contradiction to the finding that the percental spread in cumulative Ea or ETa attend was larger for the bare soil scenario, compared to the vegetated ones. Furthermore, for some texture classes the total number of outliers increased remarkably when vegetation was implemented, such as for the clay class, where the bare soil scenario has no outliers (0%), while the percentage of outliers increased to 14 % for the homogeneous grass and wheat scenario, respectively. This indicates that the differences in available root zone water, affecting actual transpiration, are the main driver for differences between PTFs, compared to fluxes over the soil surface Ea . On the other hand, only the silty clay and the silty clay loam showed no outliers at all for the instantaneous flux for all scenarios. Looking at all soil class/PTF/scenarios combinations, no clear trend in the total number of outliers in instantaneous evapo(transpi)ration flux, and flagged outliers for the cumulative Ea orETa flux at tend can be observed. This leads to the conclusion that the outliers in instantaneous flux alone do not necessarily sum up to a cumulative flux flagged as an outlier.
Explaining variability and outliers by soil physical properties
As has been shown, substantial variability exists in cumulative and instantaneous fluxes, and some PTFs are found to be more robust than others. In this section, we discuss in more detail the reasons for the differences between the predicted soil water fluxes, resulting from the use of different PTFs, by analysing the estimated hydraulic parametersKs , λ (MvG tortuosity parameter) and the soil physical characteristics. In general, variability between estimatedKs for the different PTFs is quite low (Fig. 9a), and values for Rawls MvG and BC only are significantly lower than all other PTFs. These lower values may explain the poor numerical convergence for these simulations, and the prevalence of lowerEa fluxes as well as a high number of lower 0.05 percentile outliers at tend as depicted in Tab. 4, especially for Rawls MvG. Clapp&Hornberger Ksvalues are significantly higher than those estimated by the Weynants PTF, ‘Toth class’, and Cosby SC and SSC, yet did not show any high outliers for Ea fluxes. Interestingly, Cosby SC and SSC were developed based on the same water retention andKs data as Clapp&Hornberger, as both used data from Holtan et al. (1968), nevertheless estimatedKs values are quite different. One reason might be that Clapp&Hornberger only used textural classes, and averagedKs for those classes, whereas Cosby SC and SSC is a continuous PTF. Coming back to the outliers listed in Tab. 4, those runs based on Weynants PTF indicate larger Eafluxes at tend and a large number of upper 0.95 percentile outliers, whereas their estimated Ksis not significantly different from most other PTFs. Here, it has to be noted that Weynants did not estimate Ks but rather estimated a near saturation hydraulic conductivityKs * that is mainly controlled by textural properties and which is lower that Ks . The results suggest that variability in Ks alone cannot explain the flux differences simulated.
Looking at the λ value used in MvG formulation, two different classes of PTFs can be distinguished, those setting λ to 0.5 as originally proposed by van Genuchten (1980) (Carsel&Parrish, Rawls MvG, and ‘Toth continuous’) and those who fitted λ as an additional free parameter (Rosetta SC and SSC, Woesten, Weynants, and ‘Toth class’). The variability in λ is plotted in Fig. 9b. It shows that the λ estimates of Weynants’ PTFs are significantly lower than those from the other four PTFs estimating λ, except for ‘Toth class’. ‘Toth class’ λ values are significantly lower than those calculated by Rosetta SC and SSC, and than those setting λ to 0.5, whereas Woesten is significantly lower than Rosetta SC and SSC, and < 0.5. The more negative λ values for Weynants appear strongly related to the larger number of upper 0.95 percentile outliers listed in Tab. 4, whereas the intermediately low λ values for ‘Toth class’ and Woesten PTF do not explain the number of flagged outliers. In general, λ is significantly correlated to the MvG parameter n for those PTFs setting λ ≠ 0.5 (R 2 =0.40, p = 0.05, data not shown) indicating a nonlinear behaviour which can be described as\(n=1.58\ e^{0.064\ \lambda}\) with an R 2 of 0.51. Looking at the ranges of λ, there is a systematic difference between PTFs, with largest λ values for Rosetta (-3.1>λ<0.62), followed by Woesten ((-4.46>λ<0.60), ‘Toth class’ (-5.5>λ<0.73), and Weynants (-7.87>λ<1.92). Rosetta and Woesten are characterized by low numbers of tolerance interval outliers, whereas ‘Toth class’ and Weynants are characterized by large number of tolerance outliers, both in the upper end (upper 0.95 percentile outliers). As λ is correlated to the n parameter, and n directly impacts the hydraulic properties and hence LG ,LC , τFC , andtFC , and to a less extend S , the correlation between λ and these soil characteristics was calculated. The results indicated (data not shown) that λ is not significantly correlated to LC , tgrav ,τFC , and tFC but moderately correlated to LG(R 2 =0.31, p = 0.05) and S(R 2 =0.30, p = 0.05), whereas λ is not correlated to the flux Ea attend .
For the calculated soil characteristics LG , Weynants shows large variability and high median and significantly differs from Woesten, Rawls MvG, Rawls BC class, Rawls BC and Cosby SC and SSC. In contrast, Rawls BC and BC class show lowLG , and Rawls BC class is significant different from Rosetta SSC and Clapp&Hornberger (see Fig. 10a). Here, it has to be kept in mind that LG solely depends on the water retention characteristics and hence the n and αvalues play a crucial role in the calculation. As n is positively correlated with λ, and Weynants shows the smallest λ values, the significant difference, with regards to LG , between Weynants and most other PTFs seems logical. LargeLG values occur for very fine textures, which are classically associated to low K values that limit water supply to the evaporating surface, which is reflected by the higher number of upper 0.95 percentile outliers for Weynants, leading to a drier soil profile. Lower Ea fluxes attend , and therefore, a wetter profile occurred frequently for Carsel&Parrish and Rawls (MVG and BC), whereby all these PTFs are also located in the low LG range.
The calculation of LC is based on knowledge ofLG and the actual hydraulic conductivity distribution above the evaporation front. Therefore,Ks plays also an important role in the calculation of LC . The impact ofKs on LC is clearly reflected in the high LC values for Clapp & Hornberger, which exhibit high Ks values across all soil classes compared to all other PTFs (see Fig. 10b). At the other end of the spectrum, the impact of Ks onLC is also apparent for Rawls MvG and Rawls BC which do not indicate much spread and are characterized by lowKs and hence low LC . Surprisingly, Clapp&Hornberger are not classified as outliers when looking at cumulative fluxes (see Tab. 4), whereas the lowLC for Rawls MvG corresponds to the number of outliers detected. On the other hand, Weynants, which was characterized as the PTF with most outliers at the upper 0.95 percentile, lies in the middle of the range of LC values depicted in Fig. 10b, indicating that LC might not be a good indicator for flagged outliers. As stated in Lehmann et al. (2008),LC longer than 1 m are considered as unrealistic (evaporative extraction of water by capillary flow across several meters is unlikely). Interestingly, only the Clapp&Hornberger PTF showLC > 1 m, while all other PTFs give realistic values.
The analysis of MFP shows a quite different picture (Fig. 10c). Here, the PTFs based on Brooks Corey group together and exhibit a higherMFP compared to the MvG based PTFs. Testing on significance showed that Rawls BC class, Clapp&Hornberger, and both Cosby PTFs are significantly different from all others and that only Rawls BC is not significantly different from those using MvG formulation, except for Rawls MvG. This is of interest, as Rawls MvG is only a ‘translation’ of the Brooks Corey to van Genuchten parameters from Rawls BC according to Morel-Seytoux (1986), while keeping Ks . As the Weynants PTF showed substantial outliers, as listed in Tab. 4, one would also expect Weynants to be different with regards to MFP as the λ value is much smaller compared to all other PTFs, whileKs does not differ (see Fig. 9a and 9b). One reason for the fact that MFP for Weynants does not differ from the other PTFs might be its relatively low n value, as λ andn are positively correlated. The impact of λ as opposed to the effect of MFP becomes clearer when we compare Weynants and Woesten, which show no significant difference in MFP , yet largerKs values for Woesten and lower λ for Weynants. Overall, the MFP cannot explain the outliers detected and depicted in Tab. 4 as only Rawls MvG is systematically different and exhibits large number of outliers, whereas Weynants MFP are in the centre of the range of values found for the different PTFs. On the other hand, MFP values for Clapp&Hornberger, as well as for both PTFs from Cosby, are significantly higher, yet do not stand out in Tab. 4.
With regards to the sorptivity S (Fig. 11a), there is a large variability in S for Clapp&Hornberger, which is significantly different from all other PTFs. Small variabilities in S , however, are found for Woesten, Rawls MvG and BC, Weynants, ‘Toth class’ and both Cosby PTF. In general, S is moderately correlated toLC (R2 = 0.40).
Rawls BC shows a high tgrav , which is significantly different from all other PTFs, except for Rawls MvG. Both Cosby PTFs and both Rawls continuous functions (Rawls MvG and BC) show relatively large variability (Fig. 11b). The highertgrav for Rawls MvG fits with the larger number of outliers listed in Tab. 4, whereas for BC this pattern is not clear, maybe due to the lack of numerical convergence. In general, largertgrav values are associated with more fine-grained soils such as loam and clays (Alastal, 2012), whereas the low tgrav of Woesten characterizes more coarse soils such as sands.
High τFC were calculated for Rawls MvG (Fig. 11c), whereby the large τFC is associated with extremely low predicted Ks values. Extremely high values were found for Rawls MvG with τFCexceeding 3 million days, whereby Rawls MvG hasKs values of 0.01 and 0.004 cm d-1 for the silty clay and clay class, respectively, and also did not converge. For the two soil classes, silt and silt loam, where the model run did converge τFC is also extremely large (>44,000 days) and for these soils again low Ks values of 0.2 and 0.3 cm d-1, respectively, were estimated. Additionally, these two model runs are also outliers at the lower 0.05 percentile. Clapp&Hornberger PTF resulted in the smallest τFC , whereas the Kspredictions are in general higher as for the other soils (see Fig. 9a) and none of the simulations were flagged as outliers. On the other hand, all other PTFs have comparable τFC values, and the outliers detected in Tab. 4 seem not to be linked withτFC .
Finally, tFC was analysed, which shows the same pattern as τFC , which is to be expected astFC and τFC are linearly correlated as also shown by Assouline and Or (2014).
As these soil physical characteristics were calculated to help explain differences in simulated Ea attend , all characteristics were correlated againstEa at tend (see Fig. 12). Only log10(LG ) shows a moderate correlation to Ea at tend(R 2 = 0.52, p =0.05) and a weak correlation was found for log10(LC ), withR 2 = 0.29 (p =0.05). AsEa , and also drainage D attend , will be biased if runoff is generated (because less water will infiltrate into the soil profile and be available for evaporation and drainage), Ea and drainage at tend were normalized (Ea_norm , Dnorm ) by dividing Ea or drainage attend by the difference of precipitation attend (2479.72 cm) and runoff attend . By doing so, the correlation between λ andEa_norm increased to R 2=0.31 (p = 0.05). For the derived soil characteristics the correlation also increased (to R 2 = 0.57;p =0.05) for log10(LG ) but decreased for log10(LC ), toR 2 = 0.10. On the other hand, the correlation slightly increased for tFC , fromR 2 = 0.09 to 0.22.
In a next step, a principal component analysis (PCA) using all converged model runs and soil hydraulic parameters available for MvG and BC (θr , θs ,Ks ) as well as all soil characteristics (Lc , LG , MFP ,S , and tgrav , tFC , and τFC ) and fluxes (Ea_tend , Ea_norm ,Dtend , and Dnorm ) was performed on log transformed data (except θr , θs , MFP ,Ea _norm, andDnorm ) and the results are plotted in Fig. 13. The first three components explain 76 % of the variability in the data and the important loadings on PC 1 (42.9 % of variability) aretFC (0.38), Ks (-0.35), and LG (0.33). PC 2 (24.5 % of variability) includes the important loadings LC (0.47),Ea at tend (0.40), andS (0.31). PC 3 explains only 8.6 % of the variability andtgrav (0.48) and θs (0.47) are the important loadings. The PCA triplot shows scatter of the individual PTFs around the origin of the triplot but also distinct PTF clusters, whereby Weynants (black circle) is oriented along the PC 1 in a fairly small volume and is positively correlated totFC and τFC and negatively to D at tend (as drainage D attend is negative per definition). Rawls (MvG and BC) is oriented in the same direction as Weynants but it exhibits larger scatter, whereas Clapp&Hornberger (red solid markers) is oriented along PC 2 and correlates positively with Ks .Ks values reported by Clapp&Hornberger are amongst the highest compared to all other PTFs as already discussed in relation to Fig. 9a.
Out of these 13 PTFs, three (Clapp&Hornberger, Weynants, and Rawls) can be identified as being distinctive from all others in the triplot as they do not cluster around the origin. Furthermore, they do not only differ considerably in their estimated soil hydraulic parameters (e.g., λ and n value for Weynants, and Ks for Rawls and Clapp&Hornberger) but also in the soil characteristics derived from these parameters, whereby in all soil characteristics either the n value (remember that n is correlated to λ) as well as Ks are directly or indirectly integrated. For example, the low LC values for Rawls PTFs indicate that the maximum extent of the flow region sustaining evaporation is much smaller than for all other PTFs. This results in lowEa at tend compared to other PFTs and larger number of outliers as depicted in Tab. 4.
Finally, a multiple regression was performed to test whetherEa at tend can be predicted by the soil hydraulic parameters and/or characteristics, whereby only one of those parameters or characteristics were used in turn, i.e. those that were available for MvG and BC. As per Fig. 13, all entries were log transformed except for θr , θs , and MFP , and the best regression was selected using bootstrapping. The best predictive model was found byEa @ tend = 1252.13 + 183.30 log10(LG ) + 367.88 log10(LC ) - 405.22 log10(S ) with an R 2 of 0.88 (see Fig. 14) pointing to the fact that the soil characteristicsLG , LC , and Sdescribe well the physical behaviour of soils with regards to actual evaporation. Using Ea_norm instead ofEa decreased the predictive power of the multiple regression (R2 = 0.75).
Summary and Conclusion
In this study 13 pedotransfer functions (PTF) were used to populate the hydraulic parameters required in the HYDRUS model that was then used to simulate the water fluxes for 12 USDA soil classes, for different model scenarios that varied in complexity (homogeneous or layered soil profile, with and without vegetation) over a period of 30 years. Plotting the hydraulic functions (water retention and hydraulic conductivity curves) for all PTFs revealed large differences, especially for the hydraulic conductivity curve, leading to the hypothesis that the different PTFs will also show substantial differences in simulated fluxes.
It turned out that some PTFs generated parameters that rendered the HYDRUS model numerically unstable, so that it failed to converge for certain soil class/configuration combinations, especially those reported by Rawls and Brakensiek (1985) (Rawls MvG) and by Rawls et al. (1982) (Rawls BC), which converged only in less than 44 % off all simulation runs. Surprisingly, PTFs using the Brooks Corey (BC) formulation resulted in higher convergence rates, compared to those based on Mualem van Genuchten, even though BC is in general perceived to be less stable.
In a next step, differences in simulated actual evaporationEa or evapotranspirationETa between the model runs were analysed, asEa and ETa indirectly contain information on the net infiltration, deep drainage (over long-term) and water stored in the root zone. Therefore, the cumulativeEa or ETa at the end of the simulation period (tend = 10988 days) was selected and the 90 and 70 % tolerance interval as well as the model ensemble mean were calculated. Fluxes exceeding the tolerance limits were flagged and counted. The results indicate that some PTFs (Rawls MvG, Weynants, and Carsel&Parrish) were classified as non-robust, as the fluxes generated by the parameters derived from these PTFs exceeded a defined threshold of 20 % of the 90 % tolerance interval outliers over all scenarios and soil classes. On the other hand, all PTFs using the Brooks Corey formulation (Rawls BC, Rawls BC class, Clapp&Hornberger, Cosby SC, and Cosby SSC) are classified as robust, as they generally result in a low percentage of 90 % tolerance outliers. The PTF of Woesten performed best, and it showed no outliers at all for the 90 % tolerance interval. A hypothesis raised at the beginning of the study was that increasing model complexity will reduce the variability in predicted fluxes. Therefore, the individual simulatedEa and ETa fluxes attend were compared to the model ensemble mean (MEM), and the relative spread of the individual simulations was calculated. The results show that the bare soil scenarios exhibit the highest mean percentage spread (30 %), whereas the grass and wheat vegetated scenarios had a reduced spread (23%), averaged over all soil classes. The reduction in relative spread with the inclusion of vegetation can be explained by the fact that for these runs the water leaving the soils can be extracted from the entire rooted soil profile (after which it gets transpired via the vegetation), whereas under bare soil conditions it can only leave the soil profile at the soil surface. In the latter case, differences in the soil hydraulic properties close to the surface, especially in unsaturated hydraulic conductivity, which is highly variable (in order of magnitudes) between PTFs, will impact the Ea or ETa flux more substantially.
The instantaneous Ea orETa fluxes over time were also analysed, whereby again the 90 % tolerance outliers were calculated and counted. The results indicate that some PTF/soil class/model scenario combinations showed substantial outliers in the instantaneous fluxes, yet were not flagged as outliers for the cumulative flux attend , indicating that the non-flagged instantaneous fluxes compensate these outliers. On the other hand, other PTF/soil class/scenario combinations showed no outliers for the instantaneous fluxes, but were flagged as outliers for the cumulative case, indicating that even small over- or underestimations in instantaneous flux can sum up to large errors in the long-run.
To explain differences in simulated Ea for the homogeneous bare soil scenario, different soil characteristics were calculated, and a PCA was conducted using all simulated fluxes, soil hydraulic parameters and soil characteristics available for both MvG and BC. The PCA revealed three distinct PTFs clusters, namely Weynants, Rawls, and Clapp&Hornberger, whereby Weynants and Rawls were also characterized by a large number of tolerance outliers. Weynants correlates positively to gravity time of infiltrationtFC and τ FC and negatively to drainage D at tend , whereby Clapp&Hornberger is oriented in the opposite direction and correlated with the saturated conductivity Ks . For Rawls a reasonable correlation with tFC andτFC is found, but due to the large scatter for this PTF the interpretation is less clear.
Finally, a multiple regression was performed, showing, that the gravitational length LG , characteristic length of evaporation LC and sorptivity S together explain almost 90% of the variability in simulatedEa at tend .
Overall, our results provide insights in the functional behaviour of the PTFs as a bases for the selection of PTFs in land surface modelling, but also for large scale hydrological or crop models, where considerations regarding the numerical stability, model behaviour and performance over the long run and instantaneously should be balanced against each other. Based on this, Rosetta SSC+BD, Woesten, and ‘Toth continuous’ seem to be the most robust PTFs for the Mualem van Genuchten function and Cosby SC for Brooks Corey. Note, however, that our study is in essence a sensitivity analysis; it does not include model verification using measured fluxes, and it employs one model only.
In any case, the results clearly demonstrate that the choice of PTF can substantially affect the simulated fluxes, and as a consequence, the water content stored in the soil profile with part of that available for root water uptake and crop growth. Therefore, we strongly recommend to harmonize the PTFs used in land surface, large scale hydrological, or crop model inter-comparison studies to avoid artefacts originating from the choice of PTF rather than from model structures. Additionally, our study should motivate future studies, where measured verification fluxes are available from lysimeters and or eddy covariance stations.
Acknowledgements
The authors would like to acknowledge Yakov Pachepsky and Attila Nemes for the help with Rawls PTF.