Pedotransfer Function for the Brunswick Soil Hydraulic Property Model and Comparison to the van Genuchten‐Mualem Model

Modeling soil hydraulic properties requires an effective representation of capillary and noncapillary storage and conductivity. This is made possible by using physically comprehensive yet flexible soil hydraulic property models. Such a model (Brunswick [BW] model) was introduced by Weber et al. (2019, https://doi.org/10.1029/2018WR024584), and it overcomes some core deficiencies present in the widely used van Genuchten‐Mualem (VGM) model. We first compared the performance of the BW model to that of the VGM model in its ability to describe water retention and hydraulic conductivity data on a set of measurements of 402 soil samples with textures covering the entire range of classes. Second, we developed a simple transfer function to predict BW parameters based on VGM parameters. Combined with our new function, any existing pedotransfer function for the prediction of the VGM parameters can be extended to predict BW model parameters. Based on information criteria, the smaller variance of the residuals, and a 40% reduction in mean absolute error in the hydraulic conductivity over all samples, the BW model clearly outperforms VGM. This is possible as the BW model explicitly accounts for hydraulic properties of dry soils. With the new pedotransfer function developed in this study, better descriptions of water retention and hydraulic conductivities are possible. We are convinced that this will strengthen the utility of the new model and enable improved field‐scale simulations, climate change impact assessments on water, energy and nutrient fluxes, as well as crop productivity in agroecosystems by soil‐crop and land‐surface modeling. The models and the pedotransfer function are included in an R package spsh (https://cran.r‐project.org/package=spsh).


Introduction
The accurate representation of unsaturated water fluxes in soils is important for reliable descriptions of water and energy fluxes ( van Genuchten & Pachepsky, 2011) at the soil-atmosphere boundary (Vereecken et al., 2019), for the prediction of plant growth and yields (Angulo et al., 2014;Gayler et al., 2009;Hoffmann et al., 2016;Lawless et al., 2008), including irrigation optimization (Elmaloglou et al., 2010), and for the quantification of the environmental fate of agrochemicals (Diamantopoulos et al., 2017;Vereecken et al., 2016). For this, the Richards equation, which is now known as the Richardson equation (Raats & Knight, 2018), is considered as de facto standard (Diamantopoulos & Durner, 2012). As necessary input, it requires functions to describe the soil hydraulic properties (SHPs), of which the most widely used is arguably the van Genuchten-Mualem (VGM; van Genuchten, 1980) model. However, a number of well-known shortcomings are summarized by Weber et al. (2019), motivating the authors to present a model framework. In the Weber et al. (2019) model, termed Brunswick (BW) model, the soil water retention curve (WRC) and the hydraulic conductivity curve (HCC) are ramified into a capillary and noncapillary pore space. This is achieved by introducing only one additional model parameter. The noncapillary part accounts for some underrepresented effects such as corner and film flow (Diamantopoulos & Durner, 2013;Tuller & Or, 2001). With the BW model, a water content of 0 at oven dryness is ensured. Further, the inclusion of the noncapillary part enables a linear reduction of the water content with increasing pF, as has been often observed in experimental data (Schelle et al., 2013). Moreover, the often observed change in slope in the mid-pF range in the HCC  can also be modeled with great accuracy. Finally, it should be noted that the BW model can be used with any capillary saturation function, which can be integrated in pF space (Streck & Weber, 2020). We define pF by pF ¼ log 10 |h|, that is, as the common logarithm of the absolute value of pressure head h (cm).
There are three principle methods of gaining knowledge about the SHP: first, by obtaining knowledge of the WRC and HCC from laboratory experiments (van Genuchten et al., 1997) and subsequent parameter estimation; second, through inverse modeling of time series of the relevant state variables (Romano et al., 2002); and third, by pedotransfer functions (PTFs) . In general, PTFs predict any difficult to be measured soil property from easily measurable soil properties. When considering soil hydrology, we should, in fact, explicitly refer to hydro-pedotransfer functions (hyPTFs), since other types of PTF exist, too (Goncalves et al., 2001;Martin et al., 2009;Selim, 1998). hyPTFs relate soil properties or soil texture class (and, if available, bulk density, soil organic carbon or matter content, cation exchange capacity, etc.) to SHP and SHP model parameters (e.g., Meskini-Vishkaee et al., 2014;Schaap & Leij, 2000;Tóth et al., 2015). While the first two methods may present more accurate knowledge, they are typically labor intensive, as well as being location and scale dependent. An advantage of hyPTF is that for field-scale application soil texture information may be derived directly from soil maps Tóth et al., 2017). Laboratory measurements can be effective in estimating large-scale effective SHP , but, in spite of the large number of data sets of laboratory measured SHPs included in the construction of hyPTF, concerns linked to scale and effectiveness may certainly prevail (Philip, 1980;Vogel, 2019). We can, nevertheless, expect that if the uncertainties during PTF development are contained (input uncertainty) and reported (output uncertainty), the application of hyPTF to the field scale can capture some problems of scale if adequate approaches (e.g., stochastic modeling; Beyer et al., 2009;Famiglietti & Wood, 1994) are used to represent the soil spatial variability .
The abundant literature on hyPTF contains a multitude of functional relationships between soil texture (properties) and SHP or SHP model parameters (Vereecken et al., 2010). Due to their widespread use in earth system models such as Hydrus 1D (Šimůnek et al., 2016), land atmosphere models (Vereecken et al., 2019), or soil-vegetation-atmosphere models such as the Daisy model (Abrahamsen & Hansen, 2000) and the agroecosystem multimodel library Expert-N (Priesack, 2006), there has been a strong emphasis on VGM model parameters. Some selected studies exist for other models, too . Again, it has been demonstrated numerous times that VGM and VGM-type parameterizations of the SHP suffer from well-known problems (Peters, 2013;Weber et al., 2019). However, no readily available hyPTF predicting parameters of more physically comprehensive models such as the BW model exist. In part, this might be due to the fact that the necessary large data sets with measurements of the WRC and HCC are typically not available to the wider scientific community, in spite of the proviso of replicable science.
Nonetheless, we are determined to make use of the fact that many numerical codes with implemented and tested hyPTF for the VGM model already exist. We can exploit this fact with the following approach: based on a data set of WRC and HCC data, compiled from literature sources, we established a constitutive 10.1029/2019WR026820 Water Resources Research relationship between the model parameters of the VGM model and the BW model in the VGM variant (BW-VGM). With this strategy, the direct use of existing codes is possible, while enabling a more physically comprehensive representation of the SHP. This leads us to the first hypothesis in which we state that there is a quantifiable relationship between the VGM parameters and the BW-VGM model parameters. Second, we hypothesize that the BW-VGM model, also in a data-based evaluation, is better than the VGM model. Therefore, our aims were to (i) compare the model performance based on a large data set, (ii) establish a parameter equivalence (Assouline & Or, 2013) between the VGM and the BW-VGM model parameters, and (iii) increase the utility of existing hyPTF by extending them to the BW-VGM.

Material and Methods
To achieve these aims, the following research approach was adopted, and a flowchart of the research approach is given in Figure 1.
1. Measured WRC and HCC from the literature were compiled in a data set, publicly available in the context of this work at the research data portal of the University of Tübingen, FDAT (https://fdat.escience.unituebingen.de/portal/). 2. Knowledge about the VGM and VGM-BW model parameters and their uncertainties was inferred by inverse modeling assuming identically and independently distributed residuals, the statistical properties of which were approximated by two separate Gaussian distributions: one for the WRC and one for the HCC. 3. The variance of the residual is unknown a priori. This is a direct consequence of the lack in information provided in the measurement error of the soil hydrological experiments and the unknown nature of the model error. Therefore, the variances of the combined error (residual, i.e., model + measured error) of the HCC and WRC were estimated as extra nuisance parameters. The estimate of the variance of the error can then be directly used to compare the models with each other. 4. The estimation was carried out in a Bayesian framework so that the uncertainty in the resulting model parameters is proportional to the variance of the likelihoods. The posteriors of the model parameters were sampled using a Markov Chain-Monte Carlo (MCMC) algorithm. 5. The functional relationship between corresponding VGM and BW-VGM model parameters was established for all parameters but one using weighted linear regression with individual weights in both the dependent and the independent variables. These weights were determined directly as the variance of the respective sample of the posterior.

SHP Models
Here, the two SHP models for the pressure head dependent WRC θ(h) and HCC K(h) are described. Some of the parameters exist in both models; we have chosen to use the traditional symbols but later in the study differentiate them by adding appropriate suffixes. Details of the robust inverse modeling scheme are given, and summary information describing the compiled data set is presented.

The VGM Model Soil Water Retention Function
The soil water retention θ (L 3 L −3 ) is a function of pressure head h (L) and given by the van Genuchten equation 10.1029/2019WR026820

Water Resources Research
where θ r (L 3 L −3 ) and θ s (L 3 L −3 ) are the residual and the saturated water contents, respectively, and the effective saturation function Γ(h) is given as (van Genuchten, 1980) where α (L −1 ), n (−), and m (−) are shape parameters, with the constraint of m ¼ 1 − 1/n.

Hydraulic conductivity curve
The closed-form expression for the HCC K(h) (L T −1 ) is given by where K s (L T −1 ) is the saturated hydraulic conductivity and τ (−) is a shape parameter. This leads to six adjustable model parameters which we will denote here as x VGM ¼ {θ r , θ s , α, n, τ, K s }.

The BW Model in the VGM Variant
In comparison with the VGM model, the BW framework model conceptually adopts the method of partitioning the SHP, introduced by Peters (2013), into a capillary and a noncapillary part (Weber et al., 2019). The obtained model for the WRC is continuously differentiable in h. In contrast to most other retention functions, a water content of zero at very high-pressure heads |h| is ensured. Finally, the desirable feature of a near linear decrease in water content on a log-linear scale at high-pressure heads |h| is obtained. The capillary and noncapillary conductivities are computed from the respective saturation functions, thereby accounting for more flow phenomena in the drier part of the HCC. All of this is achieved by one additional model parameter.
Here, the BW model framework is briefly explained. Full details can be found in Weber et al. (2019).

Soil water retention function
Again, the WRC is a function of h but given as a weighted sum of the capillary S c (h), and noncapillary S nc (h) saturation function, respectively, where θ cs (L 3 L −3 ) and θ ncs (L 3 L −3 ) are, respectively, the saturated water content of the capillary and the noncapillary part. In equivalence to Equation 1, the saturated water content, θ s (L 3 L −3 ) of the VGM model corresponds to the sum of θ cs and θ ncs .
The effective saturation function S c (h) in Equation 4 is given by rescaling Γ(h) (Equation 2) as (Iden & Durner, 2014) where the pressure head value at which oven dryness is attained (θ ¼ 0) is set to h 0 ¼ −10 6.8 cm, following the suggestion by Schneider and Goss (2012) and in agreement with literature values (Weber et al., 2017b). Other authors suggest to use h 0 ¼ −10 6.9 cm as anchoring point (Groenevelt & Grant, 2004). Soil water retention values in the very dry end may be determined experimentally with the dew-point method (Gee et al., 1992).
The effective saturation function of the noncapillary part S nc (h) is expressed directly as where h′ denotes the dummy variable of integration and the upper boundary of the integral −10 ϵ is a pressure head value very close to 0. From rescaling the Equation 6, an expression for S nc (h) is given as which ensures that S nc (h) scales between 0 and 1 when h varies between h 0 and ε cm. The expression for 10.1029/2019WR026820

Water Resources Research
Γ(h) can be replaced by any given capillary saturation function, such as those proposed by van Genuchten (van Genuchten, 1980), Brooks-Corey (Brooks & Corey, 1964), Kosugi (Kosugi, 1996), or Fredlund-Xing (Fredlund & Xing, 1994), demonstrating the flexibility of this approach for modeling the noncapillary retention function. For this study, we set ϵ ¼ −3. The integral in Equation 6 was solved numerically in h space as implemented in the R package spsh . The solution was derived by numerically integrating using 302 support points on the pressure head interval [10 −3 ; 10 6.8 ] and by adding measured pF values during parameter estimation for each sample individually. An analytical solution derived by Streck and Weber (2020) was derived after completion of this study and used as a benchmark for numerical solution. The results showed that no relevant differences exist between the two. Hydraulic conductivity curve Similar to the WRC, the liquid HCC is given as a sum of its capillary K c (h) (L T −1 ) and noncapillary K nc (h) (L T −1 ) conductivities. Under consideration of the isothermal vapor conductivity K ivc (h) (L T −1 ), the HCC K (h) is given similarly to Peters (2013) by where the HCC for the capillary part K c (S c (h)) is given by (Peters, 2014) where the functional argument h has been dropped from the notation and K sc is the saturated capillary conductivity (L T −1 ).
The HCC for the noncapillary part K nc (S nc (h)) is given by Weber et al. (2019) as where h r ¼ 1 cm ensuring matching dimensions, a (−) is a parameter, which governs the slope of the HCC in the part of the function, where noncapillary flow dominates and is set to 1.5 which is in agreement with direct measurements (Tokunaga, 2009) and as used in the literature (Weber et al., 2017a). The tortuosity model in K ivc we used was the Millington and Quirk model (Millington & Quirk, 1961). We give K ivc here for completeness but refer to Saito et al. (2006) or Peters (2013) for a full description.
This leads to seven adjustable model parameters denoted by x BW−VGM ¼ {θ snc , θ sc , α, n, τ, K cs , K ncs }. The shapes of the WRC and HCC of both models are shown in Figure 2, along with the specific soil water capacity function.

Data Set
A data set of 1,729 samples with WRC and HCC data was comprised from the following data sets: The Portuguese and German data sets from Weynants et al. (2013) (apart from the the BGR data, which we did not use as it was not freely available), the data in the supporting information of Rudiyanto et al. (2015), the UNSODA database from Nemes et al. (2001), and the "Vereecken database" from Vereecken et al. (1989) (Vereecken et al., 2017). The textural data are summarized in Figure 3 showing a coverage not untypical for studies on PTFs and the common underrepresentation of clay soils. The data were compiled and can be accessed publicly . Figure 4 summarizes the data set by reporting the fraction of samples for which measurements above and below certain pressure heads were reported. For example, 90% of the samples contained WRC data with pF ≤ 0.8 and 50% with HCC data (blue circles). In general, the wet end has more samples with measured WRC than HCC. In the dry end, this is reversed at pF ≥ 4.0. This is related to the adopted measurement techniques.

Parameter Estimation
In the following, two different approaches for the inference about the parameter are described. First, the method to estimate the SHP model parameters from Bayesian inference is described, in which we employ a converged MCMC sampler, to sample from the posterior distribution of the model parameters. Later, the functional relationships relating the parameters x VGM to selected parameters in x BW−VGM are described. Some parameters were transformed for the entire analyses as follows: log 10 (n − 1), log 10 (α), log 10 (K s ), log 10 (K sc ), and log 10 (K snc ). This means the prior and parameter bounds were specified accordingly.

Soil Hydraulic Properties
According to the continuous case of Bayes' theorem, the posterior probability density function (pdf) of the model parameters is (Box & Tiao, 1992) Figure 3. Texture information samples in the database (left) German soil classification system (right) USDA classification. The limits of the texture ranges can be found in the respective subfigures; in dark blue, the entire data set is shown, and in light blue, the reduced data set, DS1 (explanation in section 3.1).

Figure 2. The SHP models (a) WRC and (b) HCC with the BW-VGM parameters
, and τ ¼ 0.5 (−). Explanation of the parameters is given in the text. The continuous salmon-colored lines are identical to the VGM model, but the retention curve runs between 0 and θ sc (i.e., θ s − θ r ) and not θ s and θ r .
where f post (x k | y) is the posterior distribution of the model parameter set x k (i.e., x VGM or x BW−VGM ) conditional on the measured WRC and HCC data set, denoted by y. The term in the denominator of Equation 11, f(y), is independent of the model parameters and serves as a normalizing constant so that fpost (x k | y) integrates to unity. We explicitly note that, in formulating Equation 11, we do not explicitly formulate prior model weights (Schöniger et al., 2014). However, we do note that through the different number of parameters, different f prior (x VGM ) values are calculated, which are not equal for VGM and BW-VGM. This is an indirect result of the latter model carrying an additional parameter (K snc ), while the parameter bounds which correspond in both models were not altered. Scharnagl et al. (2011) demonstrated how the choice of prior may considerably influence the shape and the mode of f post (x k | y) for SHP parameters, but in the interest of a fair comparison, the prior was bounded, but no further shape assumptions beyond bounded uniform were included. Therefore, we are assuming here that f prior (x k ) is best described by a multivariate distribution with uniform distributed bounded marginals; the bounds are presented in Table 1. From this, we can compute ln f prior (x VGM ) ¼ −3.51 and ln f prior (x BW−VGM ) ¼ −5.3, implying prior weights of 6:1 in favor of the VGM model. The assumptions on the prior, here, render the posterior proportional to the likelihood function. The posterior sampling of f post (x k | y) by a MCMC is only known to a constant of proportionality (Gelman et al., 2014) so that MCMC sampling is possible without specifying f(y). The joint likelihood to observe the data is then given by where N θ is the number of data points in the observed WRC,N log 10 K is the number of the observed HCC, and y i,θ and y i; log 10 K are the ith measurements in the WRC and HCC data set, that is, the different data groups.
We point out that we have simplified the notation and subsequently dropped an explicit reference to L i,θ and L i; log 10 K . In the following we assume statistically independent and heteroscedastic errors (homoscedastic within θ and K) which are calculated by Note. Explanation of the parameters is given in the text. The parameters α, n − 1, K, and σ were estimated in the log 10 -transformed space. 10.1029/2019WR026820

Water Resources Research
where f(h i , x k ) denotes the model-predicted WRC or HCC value at pressure head h i and model for k, with N ¼ N θ þ N log 10 K . The errors in WRC and HCC are assumed to have a zero mean (E[e ik (h i , x k )] ¼ 0) and are normalized by the respective standard deviations σ lk , where l denotes the data group of WRC and HCC. The standardized residual e ik is obtained by The fact that σ lk can be different for WRC and HCC leads to the choice of estimating an individual σ lk ,. Thus, we explicitly write σ θ and σ log 10 K which were estimated as additional nuisance parameters in the Bayesian inference for each of the two models with the same bounds ( Table 2). As illustrated in Weber et al. (2018), the choice of the likelihood model L(y| x k ) influences f post (x k | y). With reference to the above, this can lead to better approximations of f post (x k | y i ) but requires close, that is, supervised examination of the inference result for each individual sample and model. Since this is not within the scope of this study, we assume that e ik is approximately normal and the likelihood model for HCC and WRC, under the consideration of Equation 14, is then formulated as For the estimation of the best parameter set, we used a parallelized global optimization algorithm (differential evolution; Storn & Price, 1997) as implemented in the R package DEoptim (Mullen et al., 2011) using a maximum of 500 iterations, default settings, and setting the number of population members as 10 times the number of parameters. The obtained best parameter estimates were used in subsequent analyses. To sample from the posterior, we used the delayed rejection adaptive Metropolis algorithm (Haario et al., 2006) as implemented in the FME package in R (Soetaert & Petzoldt, 2010). The sampling was performed using a single chain with 50,000 iterations. Convergence was attested by ensuring the Gelman and Rubin's convergence diagnostic b R (Brooks & Gelman, 1998;Gelman & Rubin, 1992) fulfilled the condition b R<1:2, for each element in x, calculated based on two additional chains. Of all metrics, the b R statistic provides the best guidance on exactly when convergence has been achieved (Vrugt, 2016). The sample was obtained using a maximum number of three tries for the delayed rejection procedure and an adaptation of the covariance matrix of the jump distribution every 100 iterations. The first 90% of the samples were treated as burn-in and discarded, and subsequent analyses were done on the last 10%. As starting point for the chain, we used the best parameter estimate for each given sample and model and added some small amount of noise to avoid the sampler getting stuck at the mode right from the outset.

Construction of the hyPTF
The construction of the hyPTF was based on establishing functional relationships between the estimated parameters (see section 2.3.1) for the two models (x VGM → x BW−VGM , except for K snc ). We used weighted least squares to account for errors in both the predicted and the explanatory variable (Glaister, 2001) for θ snc (θ r ), θ sc (θ s − θ snc ), α BW (α VGM ), n BW (n VGM ), K sc (K s ), and τ BW (τ VGM ), while for K snc , we simply used the median of a set of estimated b K snc . Next to the adopted weighted linear regression for the remaining parameters, we tested several other regression and machine learning models to predict b K snc based on different parameter combinations. None of the results were considered robust; thus, we suggest using median b K snc . The consideration of errors was achieved by using the calculated variance from the posterior's samples f post (x k | y) marginal for each parameter. The variance of parameter j in model k for the observed data set d was calculated from the MCMC sample of the posterior and used as a weight w jkd ¼ [var(x jkd )] −1/2 for the

Water Resources Research
York regression (York, 1966). We follow this procedure so that those parameters estimated in section 2.3.1 with a large variance, automatically receive a smaller weight in the regression, hereby reducing the biases in the hyPTF. In this way, less well-constrained parameters are assigned lower weights in the regression in contrast to well-constrained parameters.
For estimating unbiased confidence intervals for the regression, we used nonparametric bootstrapping (Davison & Hinkley, 1997)

Model Comparison and Evaluation
For the model comparison, we use several metrics to compare the model performance: (i) the mean error (ME) describing the average bias, (ii) the mean absolute error (MAE) describing the magnitude of the bias, (iii) the root mean squared error (RMSE) describing the spread of the error around the mean, (iv) the standard deviation of the residual for each sample and data group, which can be seen as a contributor to the variance of the model parameters, in the presence of a uniform f prior (x k ), and (v) the Akaike information criterion (AIC; Akaike, 1974) and the Kashyap information criterion (KIC; Kashyap, 1982), metrics rooted in information theory, used as model ranking metrics (Ye et al., 2008). ME and MAE are used as summary statistics over all results.
First, we directly compare σ θ and σ log 10 K (as introduced in 2.3.1) as a measure for the scatter of the observed data around the calibrated model. The assumption is that the total variance of the residual σ 2 lk; tot in data group (WRC, HCC) and of model k (VGM, BW-VGM) is described by the where σ 2 l; OBS is the measurement error and σ 2 lk; MOD is the model error, and σ 2 lk; tot is the estimated parameter from maximizing Equation 12. In a sample wise comparison of model performance, the σ 2 OBS can be seen as a constant for each individual sample, which is not known. Since it is constant, we can directly compare σ l,VGM,tot with σ l,BW,tot as a qualitative measure for the spread of the residual as a result of the model error of the corresponding model. The smaller value is to be favored.
The RMSE for data group (WRC and HCC) and model k (VGM, BW-VGM) is given by where N l is the number of observations in data group l. We acknowledge that there is an equivalence between the comparison of RMSE lk and σ lk .
Lastly, we used two different information criteria to evaluate the model performance by using the concepts behind the AIC and the KIC. According to Höge et al. (2018), AIC is the criterion for identifying the model with the largest predictive capability, and KIC is the criteria to identifying which model is the most likely data generating process.
The AIC for model k is given by 10.1029/2019WR026820

Water Resources Research
where N b x is the number of estimated parameters and L y ð jb x k Þ is the maximum likelihood estimate in Equation 12. The KIC is given by where ln f prior (x k ) is the prior pdf value, and |C b xb x | is the covariance matrix of the posterior, calculated from the MCMC sample. By comparison, models with lower values are to be favored. Interpreted in model weights, a difference in 10 signifies that the model with the lower value receives 100% of the weight, and the other model may be rejected outright (Burnham & Anderson, 2010).

Data Selection
From the parameter estimation of the SHP models, some initial analyses and subsequent sample removal were required, since it appears obvious to construct a hyPTF based only on data which contains sufficient information to constrain model parameters and where the sampler is likely not to have been influenced by the prior bounds. Samples with instances of estimated parameter close to the limits of the already wide prior bounds (Table 1) were removed. This was achieved by removing all samples for which the estimated parameters values did not fulfill the condition θ s < 0.69, τ VGM < −1.5, and θ snc < 0.29. After applying these criteria, the initial data set was reduced to 402 of which 10 samples were chosen at random and set aside for use in the hyPTF model performance (DS4), by way of example. The remaining 392 soil samples compose data set 1 (DS1; Figure 3, light blue crosses). Parameter τ is frequently poorly defined in the VGM model and lacks in physical meaning (Peters et al., 2011), as evident from literature frequently reported as having negative values (e.g., McCarter & Price, 2014; Schelle et al., 2011). In the presence of data from the dryer ranges of the HCC, this parameter can help achieve an improved match of model and data for the frequently observed change in slope in the HCC within the mid-pF range. In other words, in the VGM model, very small τ, that is, unphysical τ < 0, mostly enables a better match of the HCC where noncapillary conductivity becomes visible in the data. Honoring this, the functional parameter relation for τ was based on a reduced data set (DS2), filtering out the samples in DS1 which had τ VGM < 0. For the hyPTF of K snc , we used DS1 with additional filter criteria discarding K snc < 10 −5 cm d −1 resulting in DS3 with 359 samples. This precaution was taken to ensure that the MCMC sampler was hardly affected by the lower boundary of the prior. Moreover, the values <10 −5 showed a large variance in the posterior, indicating that the information in the data was, at best, limited to constrain the parameter.

Model Comparison
In Table 2, the median MEs and MAEs of the models are presented for DS1. The ME in all cases is essentially 0. The only significant difference in model performance could be attested for the MAE in the HCC, where the BW-VGM model was found to be significantly better than the VGM model (Wilcoxon rank sum test, p ¼ 1.4 E-8). The ME and MAE of the WRC for both VGM and BW-VGM are very similar. The MAE for the WRC is so small that in standard laboratory settings, these values would not be differentiable. Caution is called for when interpreting this result. As a reminder, the VGM comes with greater flexibility for the WRC, because it is not forced to match a water content of 0 at h ¼ 10 6.8 , considering that most of the samples contain data only up to a pF of about 3 (Figure 4). In other words, particularly for data sets containing observations only from limited pressure head ranges (Scharnagl et al., 2011), the VGM can match the data well. So from a utilitarian point of view, it is treated as an accepted model (Fatichi et al., 2020), proven by the myriad of applications it has found. This is a direct consequence of the VGM being more flexible at the cost of not maintaining physical comprehensiveness and is therefore not as rigid as the BW-VGM, which, in turn, is forced always to meet the physically necessary condition of a water content of zero at a fixed pressure around pF ¼ 6.8. The VGM can perform well and might even be superior when matching observed data. However, it does this by releasing the model from physical comprehensiveness, so it cannot claim to describe the properties of a semidry or dry soil, as it is the case for the BW-VGM. This points to the fact that the data are a 10.1029/2019WR026820

Water Resources Research
limiting factor, since it does not cover the entire pressure head range, highlighting the need for routinely obtaining SHP good data in the dry (and very wet) range of both HCC and WRC. Summarizing our findings, we can state that based on the comparison of the ME and the MAE for the WRC data, the model-data match for BW-VGM is, on average, equal to that of the VGM. On the other hand, for the HCC, the MAE is significantly larger for the VGM, while no significant difference could be distinguished for the ME. This means that on average, both models perform equally well, but the VGM model does not represent hydraulic properties of soils in the dry end of the WRC and HCC. With regard to the overall performance, the magnitude of the bias is significantly reduced by introduction of the physically more comprehensive model, which is mostly an attributed to the better performance in HCC.
This fact is further supported by Figure 5 where the RMSE values and estimated variances σ θ and σ log 10 K are shown. In the RMSE of the WRC ( Figure 5, top left), a slight scatter around the 1:1 line is observable. This is similar for the estimated variances ( Figure 5, bottom left), where, in the majority of cases, the variances are near equal ( Figure 5, bottom left). This is different for the HCC (Figure 5, top right panels). Here, almost all estimated RMSE and variances are plotted below the 1:1 line, and the BW-VGM model performs significantly better than VGM (Wilcoxon-Rank-Sum test p ¼ 1.68 E-10). This is linked to the VGM model's ability to describe data, which is in a pressure head range attributed to

10.1029/2019WR026820
Water Resources Research the noncapillary part, which can only be achieved by foregoing descriptive capability of data in the capillary part, a direct consequence of the lack in physical comprehensiveness. This is in contrast to the BW-VGM, which can describe this change in slope in the mid-pF range of the HCC (Fig. 2b), expressed in a smaller residual variance.
With regard to the evidence from the Bayesian selection criteria, the above-mentioned better performance is further corroborated by the histograms of the differences between the AIC and KIC (Figure 6 van Genuchten variant). The ΔAIC exhibits a slight shift to the left of the ΔKIC indicating that the AIC makes a clearer statement than the KIC. While both criteria are based on the maximum likelihood estimate −2ln L yjb x k ð Þ, the penalty terms vary. Both contain a term for the number of parameters, which is multiplied by a factor. This factor is 0.16 larger for the AIC (Equations 18 and 19). Following this, the AIC initially has a larger penalty than the KIC. However, the KIC contains extra penalty for extra parameters when they are not well defined through the prior, in particular if it is noninformative as in our case (Table 2). This difference, expressed in comparative weights, means a 6:1 favoring of the VGM model. This should not be confused with prior model weights; it is the calculated value of the prior at the best parameter estimate. This could have been slightly altered, by reducing the bounds of the θ cs to smaller values, but we chose to use the same bounds as for θ s , as we do not a priori know the value of θ ncs . At that point, we had no knowledge of realistic values for θ ncs ; it could also be 0, implying the saturated water content could equal the capillary water content θ s ≜ θ cs . This meant that the bounds were identical, and no further constraint could be legitimized. The introduction of badly constrained parameters is another penalty term, reflected in the term ln|C b xb x | . These two factors are larger than the 0:16N b x; k but not large enough to induce a ΔKIC > 0, and a systematic favoring of the BW-VGM is the result. This is noteworthy, since the KIC, initially, favors the VGM through the use of the prior, before either the data or parameter dependent terms are added. But from Figure 6, we can read that even if the prior favored the VGM over the BW-VGM 6:1, this is overcome "after the data has been seen by the models" in Equation 12. We note that the AIC does not account for the prior information and thus is also suitable to non-Bayesian methods but can demonstrate the similarity in results, if this study had been conducted using simple weighted least squares.
Information on the robustness of parameter estimates is presented in Figure 7. The lower triangle shows histograms describing the distribution of correlation coefficients between parameters. These were derived by calculating the parameter correlation matrix from the posterior sample of the parameter estimation for each individual soil sample in DS1. Note that for parameter τ, we reduced the data set a little further which does not, however, significantly change the result presented here. Correspondingly, the upper triangle gives the mean correlation coefficient between parameters. The diagonal shows the densities of the standard deviations of all marginals. This was done for both models individually. The density plots indicate small standard deviations and therefore, in almost all cases, clearly identifiable parameters. The mean correlation coefficients are small to moderate. While θ r and θ s hardly correlate (−0.21), their counterparts in the BW-VGM correlate more (−0.82). However, the correlations of all other parameters do not appear to be very strong (either close to 1 or close to −1). Overall, no clear picture can be derived as to whether or not correlation increased by adopting the BW-VGM model. This is significant because it means that adopting the physically more comprehensive model together with an additional model parameter does not come at the expense of losing robustness and identifiability. 3.3.1. θ sc , θ snc , α, n, K sc , and τ The predictions of all BW-VGM model parameters, except for K snc , were based on York regression. The corresponding regression parameters (Table 3) can be directly used to derive BW-VGM parameters from VGM parameters. We remark that the intercepts acquire values which are essentially 0. The regression results are depicted in Figure 8, and the corresponding sum of the variance in x and y is given alongside them. First of all, there is a perceivable difference between the weighted line regression and the ordinary least squares (OLS). This can be attributed to the fact that, in York regression, weights are considered leading to robust and unbiased estimates of hyPTF parameters. This is not the case for the saturated water content where OLS and the York estimate are nearly identical as well essentially identical to the 1:1 line. With regard to the physical interpretation of the model parameters, this appears to be consistent, as the modulation in the BW-VGM model affects the dry end and not the total saturated water content of a soil. For the θ r and θ ncs this is a little different, where for the hyPTF slope for θ ncs (θ r ) is 1.28, this means the underestimation of a residual water content, in other words the water content which is not part of the capillary water, is underestimated by about 30% in the VGM model. This may also be an artifact when considering that for many of the samples the data were obtained only from evaporation experiments. In evaporation experiments, the WRC data are usually available until about pF~3. This coincides with the pressure head region in which soils over a large range of textures only then show a decrease in the slope of the WRC on the semi-log scale as the soil dries. Exceptions are coarse-grained soils (e.g., sands) with low θ r , where the capillary water content, in the BW model sense, may have been nearly fully drained by pF~3. This means that in nonsandy soils, θ r often is used to enable a better model-data match, in particular for data in the pressure head range with the steep slope. In other words, θ r is much more of a fitting parameter than a physically consistent terminology and only delivers limited information about a "residual water content." Moreover, if WRC data do not exist for pressure heads up to pF > 6, like in Schelle et al. (2013), then changes in θ r lead to small changes in the objective function values, meaning that θ r is less well constrained. Therefore, the estimation of θ r often leads to larger uncertainties than for the other WRC parameters (e.g., Dettmann & Bechtold, 2016).

hyPTF: Functional Relationship Between Model Parameters
In BW-VGM, θ snc is larger (Figure 8) since the matching of the WRC data in the pressure head range of pF > 3 is achieved by the near-linear slope of the noncapillary retention curve, which appears only to be achievable when θ snc ≥ θ r . It is noteworthy that log 10 α and log 10 (n − 1) are nearly identical. We speculate that this stems from an interplay with τ, K sc , and K snc , and the information content in most observed data sets provides sufficient information content for the parameters of the capillary model to account for the air 10.1029/2019WR026820 Water Resources Research entry pressure head and the width in the unimodal pore-size distribution. Generally, the larger τ, the greater the steepness of the capillary bundle model describing the HCC. In Figure 6 (left), it can be seen that the slope of the York fit for τ is so large that hardly any negative τ BW − VGM remains. This is physically consistent (Peters, 2013) and an indicator that through inclusion of the noncapillary conductivity model part, τ is no longer used to account for HCC in the dry range, which can be achieved Figure 8. Relationship between the soil hydraulic property model parameters from top to bottom θ r : θ snc (here denoted as θ r ), θ s : (θ sc +θ snc ), α VGM : α BW−VGM , n VGM : n BW−VGM , K s : K sc , τ VGM : τ BW−VGM , τ VGM : τ BW−VGM , constrained (α i , n i − 1, and K i are shown as log 10 -transformed values on both axes), OLS: ordinary least squares, CI: confidence interval.

10.1029/2019WR026820
Water Resources Research by τ VGM ≪ 0. Thus, we decided to take a functional approach, and use τ BW ¼ min(τ VGM , 0); otherwise, the same approach is used as in Table 3. 3.3.2. K snc K snc has no corresponding parameter in the VGM model. Initially, we tested the York regression (and several other methods, e.g., artificial neural networks and multiple linear regression, results not shown) to find a method to predict K snc from VGM parameters. However, it was not possible to achieve a better predictive performance than with the median of b K snc , which we determined as 10 −1.72 cm d −1 . Therefore, using the median appears to be the parsimonious choice. Diamantopoulos and Durner (2013) presented a model for corner flow based on theoretical considerations at the pore scale, whose maximum unsaturated conductivity value was determined as being very close to the value calculated here. Based on Ahe works by Tuller and Or (2001) and  and the similarity in K snc to the weighted saturated conductivity for the film flow component in Peters (2013), it appears not farfetched to construct informative priors in the near future. In particular, by considering the existing knowledge of the relationship between soil properties and saturated noncapillary flow, we are confident physical constraints can be incorporated. At this stage, we encourage the use of the published estimated hydraulic model parameters and their associated uncertainties in the accompanying data publication (cf. acknowledgement) in order to explore this further. Given current knowledge, we deem it premature to consider a prior in the Bayesian inference scheme for the parameter K snc but encourage to use the presented data set to construct priors for future inference of parameters. It means that the value of K snc ¼ 10 −1.72 cm d −1 might change in the future and still needs further research.

Performance of the hyPTF
By way of example, the performance of the hyPTF is illustrated in Figure 9 for one of the samples (UNSODA 2573) of the randomly chosen set-a-side data set DS4. The results of the remaining nine samples from the validation data set can be inspected in Figure A1. Three different results of models are compared against the measured WRC and HCC data: (i) the VGM model as parameterized by inverse modeling, (ii) the same for the BW-VGM model, and (iii) the BW-VGM model parameterized by predicting the model parameters using the hyPTF established in the paper. Evidently, the inclusion of the noncapillary part in the WRC

Water Resources Research
and HCC not only leads to the mentioned physically comprehensive model but also leads to a formidable description for pF > 2. This is also true for the BW model whose parameters were predicted applying knowledge of the estimated VGM parameters and the hyPTF.
To test the performance of the intended use of the hyPTF developed in this paper, we gathered five hyPTF from the literature. The five models were CP (Carsel & Parrish, 1988), RM7 (model H2m in Schaap et al., 2004), RM8 (H3m in Schaap et al., 2004), WPTF &Weihermüller et al., 2017, and euptf (Tóth et al., 2015) and made a double parameter prediction. In the initial step, we predicted the VGM model parameters for the six soils in DS4 which had soil property information. Based on those predicted parameters, and using the hyPTF from this paper, we predicted the BW-VGM parameters. The summary statistics are shown in Figure 10 (top row). The calculated RMSE results for both the WRC and the HCC between the modeled and measured data are displayed. In each, the RMSE value for the VGM model is given on the x-axis, and the RMSE of the corresponding BW-VGM model is given on the y axis. Visual inspection of Figure 10 clearly shows that the plotted RMSE fall below the 1:1 line for almost all cases in both data groups and, again, with a more pronounced effect for the HCC. This demonstrates that, in spite of a double prediction step, the BW-VGM model is still better than the VGM. For one sample, UNSODA 2573, we explicitly show the modeled WRC and HCC for (i) five double predicted parameter sets for the BW-VGM, (ii) the maximum likelihood estimate of the VGM-BW model  Figure A2).

10.1029/2019WR026820
Water Resources Research from maximizing Equation 12, and (iii) the hyPTF-predicted parameters of the BW-VGM model based on the maximum likelihood estimate of the VGM. We refrain from a discussion of the performance of the individual literature hyPTF at this stage. It follows that the RMSE shown in Figure 10 (top row) is within the range of other published PTFs (Nguyen et al., 2017;Román Dobarco et al., 2019;Zhang & Schaap, 2017).
We note that the presented double prediction hyPTF does not alleviate problems mentioned in the literature relating to unconsidered effects of stone content on SHPs (Nasri et al., 2015;Tetegan et al., 2015), on organic (Wallor et al., 2018), forest (Puhlmann & von Wilpert, 2012), or multimodal soils (Pollacco et al., 2017). In addition, the effects of dynamic nonequilibrium Šimůnek & van Genuchten, 2008) and macropore flow are not captured, since in the sense of the BW-VGM model, only an application to matrix flow is possible. Nevertheless, as the input hyPTF for the VGM model is region specific (e.g., Hollis et al., 2015;Román Dobarco et al., 2019), we are convinced our approach to work, here, too. Another advantage is that there are different VGM PTFs available , with different requirements in input, so that, depending on the case specific information available, the case-specific double prediction step can be made. However, it is possible that this might not replace a newly constructed hyPTF based on machine learning tools, deemed superior for these kind of predictions (Szabó et al., 2019). A difficulty for that is that existing large data sets such as the EU-HYDI (Weynants et al., 2013) and HYPRESS (Wösten et al., 1999) data sets are not publically available, which in itself hampers model improvement.

Conclusion and Outlook
The aim of our study was to rigorously compare the VGM model with the physically comprehensive SHP model proposed by Weber et al. (2019). In their model framework, some important deficiencies in the VGM model are overcome, which is a result of a larger physical comprehensiveness. By using a large data set of 402 soils samples, it could be shown that this framework, the BW model, outperforms the VGM model in describing measured hydraulic properties. This is particularly visible in the HCC with a significant decrease in average MAE of the model-data matches. We could also demonstrate that it would be helpful for model development and assessment if WRC and HCC data were determined over a much wider range of pressure heads than those commonly included in available data sets.
Many environmental modeling codes, which adopt the Richards/Richardson equation to simulate water fluxes, are often implemented with the VGM and, also hyPTF, for the VGM. With the presented work, it will now be easy to use the better SHP models and the here presented new hyPTF. With this hyPTF, readily available data from soil maps (Dai et al., 2019) can be used to parameterize the BW model. The approach of the double prediction can also be used to permute input and algorithmic uncertainties (e.g., Carsel & Parrish, 1988;Deng et al., 2009;McNeill et al., 2018) through the hyPTF. In addition, Europe wide maps of BW-VGM model parameters based on available VGM maps (Szabó et al., 2019;Tóth et al., 2017) can easily be created, facilitating spatially explicit soil moisture modeling. This means that, for the first time, the application of a physically more comprehensive models may be applied for landscape-scale research.
In the worst case, the implementation of the BW model will lead to no significant improvement in simulated water redistribution in soils. In the best case, it can lead to more realistic water uplift from the groundwater or deeper soil zones. From this, we expect to have an effect on simulated plant growth, yield development and evapotranspiration when considering the soil-vegetation-atmosphere nexus. However, this requires a separate investigation. Finally, we are convinced that the presented hyPTF will (a) strengthen the utility of the BW model, (b) enable field-scale simulations of water and nutrient flux estimations in agroecosystems also under dry conditions, and (c) enable better impact assessment of warmer and drier meteorological conditions caused by climate change (at least for Western and Central Europe where these conditions have been predicted). In other words, it is possible that in modeling global food production, considerable differences may result from adopting a new model for the soil water fluxes. One step toward achieving this is by mapping model parameters based on soil property

Water Resources Research Data Availability Statement
The data sets  used in this work can be freely accessed under http://hdl.handle. net/10900.1/c2f45822-86ba-4338-b3e0-2d068f024e07.