Performance of the Groenevelt and Grant Model for Fitting Soil Water Retention Data from Brazilian Soils

The soil water retention curve (SWRC) is essential for vadose zone hydrological modeling and related applications. In 2004, Groenevelt and Grant (GRT) presented a mathematical model for describing the SWRC and reported its mathematical versatility and good fit to soils from a Dutch database. In order to evaluate the application of GRT to SWRCs of Brazilian soils, we aimed to analyze the performance of GRT for 72 soils from Brazil. Besides that, the obtained results with GRT for these soils were compared to the fitting performance of the most frequently used models: Brooks and Corey (BC) and van Genuchten-Mualem (VGM). The three models were fitted to available soil water retention data by minimizing the sum of square errors. The Pearson correlation coefficient (r) and the Root Mean Square Error (RMSE) were used to assess the goodness-of-fit. Results showed high correlation coefficients (r≥0.95) and small values of RMSE (RMSE ≤0.03 cm cm) for all fits. The goodness-of-fit was of similar performance for the three models with a positively correlation between them. The major difference in shape among GRT, BC, and VGM occurred in the near saturated range, while they were almost identical for low matric potentials. The exponent of GRT showed to be highly correlated with exponents of BC and VGM, but the correlation between the other shape parameters is not well defined, making a direct conversion still difficult.


INTRODUCTION
The phenomenon of water retention in the soil is driven by the action of capillary and adsorptive forces, which together give rise to the soil water matric potential (Dane and Hopmans, 2002). The hydraulic function that relates the volumetric ratio of water retained in the soil to its matric potential is the soil water retention curve (SWRC). Over the last decades, several models have been developed to better describe the SWRC, such as those proposed by Gardner (1958), Brooks and Corey (1964), Campbell (1974), van Genuchten (1980), and Broadbridge and White (1988).
Fitting to a wide range of soils, the equations of Brooks and Corey (1964) (to be referred to as BC) and van Genuchten (1980) with the parametric restriction of Mualem (1976) (to be referred to as VGM), are among the most frequently used models in literature. On the other hand, the Campbell (1974) model, as well as the exponential model, are very useful in analytical solutions of complex problems regarding water retention due to their mathematical simplicity. Other less common formulations are polynomial and exponential equations (Too et al., 2014).
The SWRC is used in soil physics as well as in related areas like hydrology, soil conservation, irrigation and drainage, among others. The SWRC directly links to the soil pore size distribution function, and is used in hydrological studies (Silva et al., 2017), soil physical quality evaluation (Reynolds et al., 2009;Armindo and Wendroth, 2016) as well as in the prediction of field capacity (Turek et al., 2018) and crop water availability (Feddes and Raats, 2004). The understanding of soil water dynamics is important in applications involving infiltration, water redistribution, evaporation, and root water uptake, and helps to promote management that allows an increase in water use efficiency (Prevedello and Armindo, 2015). Groenevelt and Grant (2004) proposed a SWRC characterization model (to be referred to as GRT) showing its fitting performance to water retention data of soils from The Netherlands. Like VGM and BC models, GRT allows the prediction of the unsaturated soil hydraulic conductivity based on soil water retention data making use of Mualem (1976) or Burdine (1953) theories (Grant et al., 2010).
Since then, the GRT model has not been systematically tested for soil databases. Most Brazilian soils are the result of a lengthy pedogenesis under tropical climatic conditions with a precipitation surplus under well-drained conditions, leading to a specific clay mineralogy and distinct structure. It is therefore imperative to evaluate the performance of any SWRC model, including GRT, for these soils. This study aimed to assess the fits of the GRT model to a SWRC database of Brazilian soils. The goodness-of-fit was also compared to the two models most frequently used in literature, BC and VGM.

Database
A set of 72 soil water retention curves extracted from the Brazilian Soil Hydrophysical Database (HYBRAS) (Ottoni et al., 2018) was used. For each soil, between 6 and 13 data pairs of soil water content (θ) versus matric potential (h) were available, from soil saturation (h = 0) up to dry condition (h = -15300 cm). This database also provides information on some soil physical properties as sand, silt and clay contents, bulk density and, total porosity. Each soil was classified according to the texture classes defined in the Brazilian system of soil classification (Santos et al., 2013). In this classification, clay is defined as particles with an equivalent diameter smaller than 2 µm, silt is between 2 and 50 µm, and sand has an equivalent diameter larger than 50 µm. Very clayey texture is defined as a clay content larger than 600 g kg -1 , clayey texture corresponds to a clay content between 350 and 600 g kg -1 , silty textured soils contain less than Rev Bras Cienc Solo 2019;43:e0180217 350 g kg -1 of clay and less than 150 g kg -1 of sand, whereas a sandy texture is defined by a sand content exceeding the clay content in more than 700 g kg -1 . All other soils are of medium texture (Figure 1). This figure shows the selected data to cover all textural classes, with fewer representatives for the silty class, very uncommon in tropical soils.

Models
The measured data of matric potential (h) and soil water content (θ) were fitted to the nonlinear models proposed by Groenevelt and Grant (2004) (GRT), Brooks and Corey (1964) (BC), and van Genuchten (1980) with parametric restriction of Mualem (1976) (VGM). The GRT model was originally written as: in which θ is the volumetric soil water content (L 3 L -3 ), θ s is the saturated soil water content (L 3 L -3 ), h is the soil water matric potential (L) and p, k 1 , and k 0 are model fitting parameters. The parameter k 1 has same physical dimension as soil water content (L 3 L -3 ). The parameter k 0 has the same physical dimension as |h| (L) and corresponds to the value of at the inflection point of the SWRC, as confirmed by De Jong van Lier (2014) and Grant and Groenevelt (2015). However, its physical meaning is not clear, as occurs with parameter α of VGM (De Jong van Lier and Pinheiro, 2018). Equation 1 can be rewritten as: taking k 1 = (θ s -θ r ) and k 0 = k. The θ r is the residual soil water content (L 3 L -3 ) and Θ is the effective saturation (L 3 L -3 ), which is found by the expression Θ = (θ -θ r )/(θ s -θ r ).
Sand (g kg -1 ) C l a y ( g k g -1 ) S i l t ( g k g The BC model is defined by: in which h b is the absolute value of the air-entry pressure head (L) and λ is a fitting parameter.
The VGM model is given by: in which n and α (L -1 ) are model fitting parameters.

Calibration and validation
The parameters of the respective models were calibrated by fitting equations 2, 3, and 4 to the measured data of θ(h). The Sum of Square Errors (SSE) was minimized to obtain the best fit for each model. The fitted parameters were θ s , θ r , k, and p for GRT; θ s , θ r , h b , and λ for BC; and θ s , θ r , α, and n for VGM. During the fitting procedure, values of all fitted parameters were restricted to non-negative values, according to their physical or mathematical meaning. Then, for each fitted model to each SWRC, the goodness-of-fit was evaluated by metrics that quantify model precision and accuracy to estimate the function θ(h).
Model precision was assessed using the Pearson correlation coefficient (r), defined as: in which θ i-mea is each value of measured soil water content, θ i-est is each value of estimated soil water content, θ ̅ mea is the mean of measured values, and θ ̅ est is the mean of estimated values, all with dimension L 3 L -3 . The value of r represents a measure of the linear correlation between the measured and estimated values of θ. The closer to one, the greater is the model precision.
Model accuracy was analyzed by the Root Mean Square Error (RMSE), defined as: in which N is the number of data pairs. The RMSE expresses the difference between measured and estimated values and thus, the closer to zero, the greater is the model accuracy.

RESULTS
The Brazilian soil texture triangle with all 72 data points is presented in figure 1. This database is composed mostly of soils of medium and clayey texture and some few soils of silty, very clayey, and sandy classes. The data number identification (ID), the number of measured θ(h) pairs (N) of each data from HYBRAS database are shown in table 1 together with their respective textural classes information (T). Furthermore, fitted parameters for SWRC models of GRT (θ s , θ r , k, and p), VGM (θ s , θ r , α, and n), and BC (θ s , θ r , h b , and λ) are also exhibited with respective values of r and RMSE.  -Mualem, Brooks and Corey (1964), and Groenevelt and Grant (2004) together with Pearson's correlation coefficient (r) and Root Mean Square Error (RMSE). ( The GRT model presented the lowest value of θ s (0.272 cm 3 cm -3 ) for soil ID-995 and its highest value (0.862 cm 3 cm -3 ) for soil ID-1027. The maximum value for θ r (0.347 cm 3 cm -3 ) was found for soil ID-1000. The parameter k was lowest (1.65 × 10 3 cm) for soil ID-1027 (silty texture) and highest (2.49 × 10 4 cm) for soil ID-449 (clayey texture). The lowest value for parameter p (0.078) was found for soil ID-423 (clayey texture) and the highest value (2.963) for soil ID-153 (sandy texture).
For the BC model, like for VGM, the lowest value for θ s (0.260 cm 3 cm -3 ) was also found for soil ID-995, whereas its highest value (0.833 cm 3 cm -3 ) occurred in soil ID-401. The maximum value for θ r was found for soil ID-398 (0.303 cm 3 cm -3 ). The soil ID-157 (sandy texture) presented the highest value for parameter h b (35.24 cm), whereas the lowest value for parameter λ was 0.028 for soil ID-1000 (medium texture) and the highest value 2.434 for soil ID-157. ( Mean values for the Pearson's correlation coefficient r were highest for GRT (0.996), closely followed by VGM (0.995) and BC (0.993), showing a slight superiority of precision for the GRT model. On the other hand, mean values for RMSE were smallest for BC (0.028 cm 3 cm -3 ), closely followed by both GRT and VGM (0.030 cm 3 cm -3 ), showing a slightly higher accuracy for the BC model. The major difference in shape among the three models occurs in the near-saturated range, as shown for the cases with the best (soil ID-545; figure 2a) and the worst fit (soil ID-282; figure 2b) with GRT among the evaluated soils. In these examples, BC (black curve in figure 2) is the only one that remains constant for h in the near-saturated range (h between -10 and 0 cm). In case of figure 3, the values of Pearson's correlation r among the three assessed models for all 72 measured θ(h) data points are presented, in which GRT and VGM models exhibited larger values. Lastly, an important finding of linear correlation between exponents p (GRT) and n (VGM) shows up in figure 4.

DISCUSSION
Since all fits, regardless of the used model, resulted in very high precision (r≥0.948) and high accuracy (RMSE≤ 0.030 cm 3 cm -3 ), we conclude that the three studied models fit well to the 72 measured θ(h) data points. Based on all found measures of r and RMSE, the goodness-of-fit was slightly better (larger r and smaller RMSE) for the GRT model in the case of 35 of the evaluated soils (48.6 %), followed by VGM for 20 of the soils (27.8 %), and BC for 17 of the soils (23.6 %).
About the difference in curve shapes, the number of fitting parameters is the same for all three models and thus curve shapes are almost identical in the best fit for values of |h| larger than 40 cm (soil ID-545), showing almost equal goodness-of-fit among the three analyzed models. Possibly, one or two additional measured values between 0 and 40 cm of |h| might reduce the uncertainty of the non-linear fitting procedure near the saturation point. Even though more measured points were obtained for soil ID-282, a worse performance of the three models to fit these points together with a larger difference between their estimates was observed due to the incongruence between the measured values of this SWRC and the curve shapes. This is illustrated in another way in figure  HYBRAS database. The strong correlation between (1 -r) for the different models is clear (Figure 3), in other words, the goodness-of-fit among the three models correlates positively, and data that allow a better fit for one of the models tend to a better fit for the other models as well.  (1 -r) VGM The equation by Groenevelt and Grant (2004) may be considered mathematically more convenient than VGM, allowing straightforward integration of θ(h) to obtain the integral water capacity (Groenevelt et al., 2001;Grant and Groenevelt, 2015) and K(h) to obtain the matric flux potential (Raats, 1977;Pullan, 1990;Grant and Groenevelt, 2015). Furthermore, the exponent p is linearly correlated to the slope of the SWRC, with |h| on a log-scale, sometimes referred to as the S-index (De Jong van Lier, 2014). Nevertheless, most databases on soil hydraulic properties report the VGM parameters. A correlation between parameters  of both equations would allow to transform databases in VGM to GRT. We verified the correlation between exponents p (GRT) with n (VGM) and also with λ (BC), obtaining a strong linear correlation (r = 0.985) between p and n and moderate linear correlation (r = 0.781) between p and λ for the evaluated database ( Figure 4) according to: p = 0.969n -0.869 Eq. 7 p = 1.217λ + 0.162 Eq. 8 The same figure shows that correlations between parameters α and k as well as parameters k and h b are not well defined. Analyzing these correlations for each texture class separately did not generate promising results either. This is somehow unexpected, as 1/α and k apparently have a similar role in the equations. The correlation between parameters of GRT with VGM and BC models could support the exchange of information related to SWRC between these models providing several applications due to the higher mathematical versatility of the GRT model. Therefore, a further investigation of the correlations for other soils may be of interest.

CONCLUSION
An analysis of water retention data for 72 Brazilian soils allowed to conclude that soil water retention data can be fitted with equal quality to the equations by Groenevelt and Grant (2004) (GRT), van Genuchten (1980) with Mualem restriction (VGM), and Brooks and Corey (1964) (BC), suggesting the use of the GRT model for Brazilian soils to be of interest. The major difference in shape among the three models occurs in the near saturated range. Exponents from GRT are correlated with exponents from BC and VGM, but the other shape parameters (k for GRT, with h b for BC, and α for VGM) do not show clear correlation, making a direct conversion between the equations difficult.