Abstract
Soil-water characteristic curves (SWCCs) are important in terms of groundwater recharge, agriculture, and soil chemistry. These relationships are also of considerable value in geotechnical and geoenvironmental engineering. Their measurement, however, is difficult, expensive, and time-consuming. Many empirical models have been developed to describe the SWCC. Statistical assessment of soil-water characteristic curve models found that exponential-based model equations were the most difficult to fit and generally provided the poorest fit to the soil-water characteristic data. In this paper, an exponential-based model is devised to describe the SWCC. The modified equation is similar to those previously reported by Gardner (1956) but includes exponential variable. Verification was performed with 24 independent data sets for a wide range of soil textures. Prediction results were compared with the most widely used models to assess the model’s performance. It was proven that the exponential-based equation of the modified model provided greater flexibility and a better fit to data on various types of soil.
1. Introduction
The soil-water characteristic curve (SWCC) can be defined as the relationship between water content and suction in unsaturated soils. It can be viewed as a continuous sigmoid function describing the water storage capacity of a soil as it is subjected to various soil suctions. The soil-water characteristic curve contains important information regarding the amount of water contained in the pores in any soil suction situation and the pore size distribution related to the stress state in the soil-water [1].
Soil-water characteristic curves are important for groundwater recharge, agriculture, and soil chemistry. They are used to predict the soil-water storage, water supply to plants (field capacity), and soil aggregate stability [2]. These relationships are also of considerable value in geotechnical and geoenvironmental engineering. Unsaturated soil mechanics has primarily utilized SWCC for the estimation of unsaturated soil property functions which are subsequently used in numerical modeling solutions to geotechnical and geoenvironmental engineering problems [3]. SWCC field or laboratory measurements are difficult, expensive, and time-consuming; however, because of its importance and usefulness, many empirical models were developed to describe the SWCC.
A number of mathematical equations have been proposed in the literature to represent SWCCs. Fredlund and Xing [4] and Leong and Rahardjo [5] reviewed a range of proposed soil-water characteristic curve models along with parametric studies. Most of the equations described earlier are empirical in nature and based on the shape of the soil-water characteristic curve. Sillers et al. [6] also summarized numerous mathematical models of the SWCC. The mathematical models presented in their paper can be categorized in a number of ways to illustrate the characteristic equations such as parametric studies, as well as their advantages and disadvantages. Models such as those of Gardner [7], Brooks and Corey [8], Brutsaert [9], Tani [10], McKee and Bumb [11, 12], Van Genuchten [13], Burdine [14], Mualem [15], Kosugi [16], and Fredlund and Xing [4] were analyzed.
The most commonly employed classical retention models are the sigmoid function models by Gardner [7] and Van Genucthen [13] and the power function model by Brooks and Corey [8]. Recently, the lognormal distribution model by Kosugi [17] has gained popularity due to its great flexibility in terms of representing the water content in wet and dry ranges for all soil types. In contrast, the exponential-based models by McKee and Bumb [11, 12] are less commonly used. According to Sillers et al. [6] the problem with these models is that they are relatively less flexible. After conducting statistical assessment of a number of soil-water characteristic curve models, Sillers and Fredlund [18] concluded that the exponential-based models are the most difficult to fit. One reason for the difficulty in performing the fit is the overlap in the domain of each parameter. This is because the parameters of the exponential-based models affect both the shape and the position of the curve. The interdependence of these parameters causes uniqueness difficulties for the fitting routine which increase the number of iterations and the number of trials before convergence to the best-fit parameters. Tani [10], Russo [19], and Ross and Smettem [20] also proposed soil-water retention models. Their model equations are reproduced from Gardner’s [21] exponential function for unsaturated hydraulic conductivity when it is incorporated in Mualem’s [15] conductivity model. The equations provide reasonable plots for water retention curves; however, the models are less accurate than the other two- or three-parameter models as described by Kosugi et al. [22].
With this in view, therefore, an exponential-based model is proposed to describe the SWCC. The modified equation is similar to those previously reported by Gardner [7] but includes exponential variable. The Gardner model was chosen as a starting-point because of its simplicity and easiness to fit. Another advantage is that the model has parameters that have an independent effect on the soil-water characteristic curve. The performance of the modified Gardner model is then demonstrated on a variety of soil types and the fitting performance in terms of accuracy (RMSE) and linearity () is compared with the existing models suggested by Brooks and Corey [8], Van Genuchten [13], and Kosugi [17]. These three existing models were selected in this study because they are widely adopted and cited and also because of their relative simplicity. They also contain four parameters in their model equations.
2. Soil-Water Characteristic Curve (SWCC) Models
Models by Gardner [7], Brooks and Corey [8], and Van Genuchten [13], and the lognormal distribution model by Kosugi [17] are some of the notable models found in the literature. All these models are parametric models based upon a pore size distribution function and the capillary theory. They also contain four parameters in their model equations. The equations representing each model along with definitions of the variables used are given as follows.
In general, the normalized water content or a dimensionless water content term, , which is also called effective saturation, , will be used to represent the equations associated with the soil-water characteristic curve models. Consider the following: where is the volumetric water content and and are the saturated and residual volumetric water contents, respectively.
2.1. Gardner’s Model
The Gardner [7] equation was one of the first equations used to model the soil-water characteristic curve. It is a continuous function originally intended as a means of modeling the unsaturated coefficient of the permeability of soil. The equation has now been adapted, however, to model the soil-water characteristic curve. The equation uses two fitting parameters, namely, and . Parameter is related to the inverse of air entry value and the parameter is related to the pore size distribution [6]. Consider the following: where is the soil suction with a unit of kilopascals (kPa).
By substituting (1) into (2), we obtain the volumetric water content form of the Gardner model as follows: The Gardner model has a particularly simple form with few parameters; it is thus convenient to use and has a wide range of applications. However, it cannot accurately describe the soil-water characteristic curve for saturated and near-saturated soils [23].
2.2. Brooks and Corey’s Model
Brooks and Corey’s [8] model is among the earliest equations proposed for the soil-water characteristic curve and remains a popular model where it is in the form of a power-law relationship. The model is given by the following equation: The equation uses two fitting parameters, namely, and . Parameter is related to the air entry value of the soil. The parameter is termed the pore size index and is related to the pore size distribution of the soil. The model is assumed to be constant for suctions less than the air entry value. The soil-water characteristic curve is assumed to be an exponential decreasing function at soil suctions greater than the air entry value [6].
By substituting (1) into (4), the volumetric water content form of the Brooks and Corey model can be written as follows: The Brooks and Corey model is relatively simple and thus widely used [23]; but the model does not provide a continuous mathematical function for the entire soil-water characteristic curve [6].
2.3. Van Genuchten Model
The most widely adopted alternative to the Brooks and Corey model is that proposed by Van Genuchten [13]. The model uses three fitting parameters, namely, , and . The Van Genuchten model can mathematically be described as follows: To simplify and derive closed form equation for unsaturated conductivity based on Mualem [15], the and parameters in the SWCC equation can have a fixed relationship with . This suggestion therefore reduces the three-parameter equation of Van Genuchten to a two-parameter SWCC equation: Parameter is related to the inverse of air entry value, the parameter is related to the pore size distribution of the soil, and the parameter is related to the asymmetry of the model. By substituting (1) into (7), we can write the volumetric water content form of the Van Genuchten model as The Van Genuchten model has a complex form and relies on more fitting parameters than the models discussed above. However, it produces a continuous output in the unsaturated zone and provides a good description of the soil-water characteristic curve under most circumstances [23].
2.4. Lognormal Distribution Model
The last SWCC model considered is based on the model suggested by Kosugi [17]. This model was developed by applying a lognormal distribution law and its parameters are directly related to the soil pore radius distribution. The lognormal distribution model by Kosugi is described as follows: where is related to the complementary error function, erfc, and defined as The model uses two fitting parameters, namely, and . Parameter is a capillary pressure head related to the median pore radius and is a dimensionless parameter related to the width of the pore radius distribution. By substituting (1) into (9), we can write the volumetric water content form of the lognormal distribution by Kosugi [17]: The lognormal distribution model has a more complex form because of the complementary error function present and thus it is difficult to use. The model does, however, have greater flexibility in terms of representing the soil-water characteristic curve in the wet and dry regions for all soil types [6].
2.5. Modification of Gardner’s Model
In general, a good mathematical expression should have only a few parameters with clear physical meaning and easiness to use. Selection of the soil-water characteristic curve model and parameter determination was done by identifying a model that could accurately describe the soil-water characteristic curve under a broad range of conditions while being as simple as possible in order to facilitate its application on a region-wide scale. The Gardner model was judged to satisfy these criteria best and so was chosen as a starting-point.
Two basic parameters incorporated in Gardner’s model are the air entry value (AEV) and pore size distribution, which are denoted as and , respectively. Here, the parameter related to the air entry value of the soil is assumed to be where is a fitting parameter which is related to the air entry value. Substituting (12) into (2) yields To facilitate the application of the model, it was rewritten by taking the exponential function of the right sides of the model. Rearranging (13) as a function of the exponential gives The modified Gardner equation has two fitting parameters, namely, and . A parametric study was used to describe the fitting properties of the modified Gardner model. Figure 1 shows a plot of the modified Gardner model when one parameter is changed and the others remain fixed. Only that feature of the curve related to the parameter being varied is affected. As a result, parameter is related to the pore size distribution index of the soil (Figure 1(a)). The larger the value of , the more uniform the pore size in the soil, and the steeper the slope of the curve. Parameter locates the curve toward the higher or lower suction regions and has a unit of suction (Figure 1(b)). The parameter does not affect the shape of the curve but shifts the curve towards the higher soil suction region as increases.
| (a) Effect of parameter with kPa |
| (b) Effect of parameter (kPa) with |
The relationship between the volumetric water and suction can be obtained by substituting (1) into (14), yielding an expression for the soil-water characteristic curve: The capability of the modified Gardner (MG) model was demonstrated for 12 soil types ranging from sand to clay. Prediction results were compared with the most widely used models by Brooks and Corey (BC), Van Genuchten (VG), and Kosugi (LN) to determine the MG model performance. These models are expressed as (5), (8), and (11), respectively. Other than and , all models contain two fitting parameters. Therefore, for ease of scheduling, the fitting parameters are marked as P1 and P2. Table 1 summarizes the parameters for the SWCC models evaluated.
3. Materials and Methods
3.1. Sources of Soil-Water Characteristic Curve Data
Twenty-four soil samples with soil-water characteristic data selected from the unsaturated soil hydraulic database (UNSODA 2.0) [24] are used to demonstrate the performance of the modified Gardner model. The selected data sets represent 12 soil textural classes in which each soil class is represented by two samples of data to be tested. The soils comprise sand, loamy sand, sandy loam, silty loam, silt, loam, sandy clay loam, silty clay loam, clay loam, sandy clay, silty clay, and clay. All these soils are identified in Table 2. These enable the modified Gardner model to validate and identify its parameters and compare them with the three most widely used models: BC, VG, and LN. The SWCCs in this study are presented in terms of volumetric water content, plotted on an arithmetical scale, and soil suction plotted on a logarithmic scale.
3.2. Model Analysis
Optimization techniques are used to obtain the best-fit parameters for soil-water characteristic curve data sets. The curve fitting routine determines model parameters such that the mathematical function passes as close as possible to the experimental data points without necessarily going through any of the points. The fitting procedures for all 24 soil data sets were performed by nonlinear least-square analysis based on a trust-region algorithm method which employed the curve fitting tool in the MATLAB program. It is an iterative method starting with some initial values of the parameters. In the curve fitting tool program, the initial value of the parameters for the modified Gardner model is set to start with zero or 0.1 for all soil types. For the others (BC, VG, and LN model), the initial values for the model parameters in the iterative procedure were obtained by using reported literature values for the different soils.
There are various statistical measures which can be used to compare the fitting accuracy of the SWCC models. In this study, the root-mean-square error (RMSE) and the coefficient of determination () are used to help determine the best fit. The RMSE (m3m−3) statistic is an indicator of the overall error of the evaluated model function, with a value closer to zero indicating a better fit. The statistic is generally the best indicator of the fit quality. It is a measure of the linearity between observed and fitted data. The with a value approaching unity indicating that the observed and fitted data sets are linearly located around the line of perfect agreement or the fitted curve is of comparable shape as the observed curve.
The RMSE was expressed as where the SSE statistic is the least-squares error of the fit and defined as where is a parameter vector containing the unknowns that need to be estimated, is the number of soil-water characteristic data for each soil sample, is the soil water content corresponding to the data pair for each soil, and obs and fit denote observed and fitted values, respectively.
The value of reflects the proportion of the total sum of squares (SST) that is partitioned into the model sum of squares (SSM) since SST is equal to SSM plus SSE. Consider the following:
4. Results and Discussion
Table 3 shows the models’ fitted parameters for various soil textural classes. In general, the values of for BC are lower than for the LN, VG, and MG models. The values of are very close to each other for the LN and MG models and for the VG and BC models, respectively. The parameter values shown in Table 3 can serve as useful initial values for researchers attempting to use one of the models.
The statistical measure values are shown in Table 4. To assist and facilitate the evaluation in determining the best model, a weightage is given to each fit. In this case, a weightage of one is given to the model that gave the smallest RMSE value (best fit) and a weightage of four is given to the model that gave the largest RMSE value (worst fit). These weightages were chosen as there are four equations ((5), (8), (11), and (15)) for comparison. The weightage is given in parentheses in Table 4 and a summary of the weightage evaluation results is shown in Table 5. From the total weightage in Table 5, it can be observed that the MG model performs much better than the three other models, as it has a weightage of 26 out of 96 total weightages representing all soil samples. The statistical analysis in terms of RMSE shows that the MG model led to the best fit of all the soil samples except for the loamy sand in soil sample 1 and the silt in soil sample 2 in which the VG model was better (Table 4). The RMSE values associated with the MG model for all soil samples ranged from 0.0005 to 0.0263, whereas for the LN model they varied between 0.0010 and 0.0273. The results clearly show that the RMSE values associated with the MG model are well below those of the LN model which earned 51 of the total weightages. It was closely followed by the VG model with 67 weightages, where the RMSE values associated with the VG model for all soil samples ranged from 0.0015 to 0.0257 and the BC model gained 96 of the total weightages, in which the RMSE varied between 0.0039 and 0.0294. Note that the fitting errors for the MG model are smaller than those for the other models. This is because of the better representation of the MG model for all soil samples. Results suggest that the MG model is most suitable for describing the observed data.
As regards the values of , a similar trend could be observed. The MG model was the best model in terms of linearity, closely followed by the LN model, the VG model, and the BC model. In Table 4, the results obtained show that the values associated with the MG model are higher than those of the three other models for all soil samples except for the loamy sand in soil sample 1 and the silt in soil sample 2 for which the VG model result was better than that of the MG model. The values associated with the MG model for all soil samples ranged from 0.9585 to 0.9999, whereas for the LN model they varied between 0.9564 and 0.9997. The VG model’s values for all soil samples ranged from 0.9382 to 0.9994, and the BC model’s values varied between 0.8920 and 0.9961.
To illustrate further the behavior of the four models compared in this study when fitted to the SWCC data for various soil types, the observed data and fitted curves are compared in Figures 2(a) to 2(x). With regard to the sand samples in Figures 2(a) and 2(b), the correspondence between observed and fitted SWCCs exhibited deviation for all models where they did not accurately match several points near saturation. All models still gave relatively good and realistic fits; however, the discontinuous character of the BC model did not seem to be problematic for the sand soils. With respect to the loamy sand in Figures 2(c) and 2(d), the MG, LN, and VG models showed very good fits to the observed data and no significant differences were observed between the three models in terms of their curves, although the VG model performed slightly better than the MG and LN models in the loamy sand of sample 1 (Figure 2(c)). The BC model missed the shape of the data near saturation because of its discontinuous character. Almost the same results were observed in sandy loam (Figures 2(e) and 2(f)) and silty loam (Figures 2(g) and 2(h)) soils. MG, LN, and VG models are very well matched to the observed data, though the MG performs better than the LN and VG models for both soils. The BC model typically had problems matching data in the transition point and did not accurately match some points near saturation, which decreased the flexibility of the curve. Although the BC model did not achieve a rate of success as high as the other models, it was able to match the observed data adequately. A similar trend could also be observed in loam soil (Figure 2(k)), sandy clay loam (Figure 2(n)), silty clay loam (Figures 2(o) and 2(p)), clay loam (Figures 2(q) and 2(r)), sandy clay (Figures 2(s) and 2(t)), silty clay (Figure 2(v)), and clay soils (Figures 2(w) and 2(x)). Slightly different results were obtained for the silt soil (Figure 2(i)), loam (Figure 2(l)), sandy clay loam (Figure 2(m)), and silty clay soil (Figure 2(u)); all models exhibit deviations in matching data in the transition point. All models still performed very well for all the samples in which they produced comparatively realistic fits. The same trend could also be observed in the silt soil in sample 2 (Figure 2(j)), where the correspondence between observed and fitted SWCCs exhibited deviation for all models. As can be seen in Figure 2(j), the only model that showed a better match to the data in the dry region was the BC model. Although all models did not accurately match several points of the available observed data in their data sets, they were still able to match the observed data in the extrapolation region adequately with a reasonable fit.
(a) Sand 1
(b) Sand 2
(c) Loamy sand 1
(d) Loamy sand 2
(e) Sandy loam 1
(f) Sandy loam 2
(g) Silty loam 1
(h) Silty loam 2
(i) Silt 1
(j) Silt 2
(k) Loam 1
(l) Loam 2
(m) Sandy clay loam 1
(n) Sandy clay loam 2
(o) Silty clay loam 1
(p) Silty clay loam 2
(q) Clay loam 1
(r) Clay loam 2
(s) Sandy clay 1
(t) Sandy clay 2
(u) Silty clay 1
(v) Silty clay 2
(w) Clay 1
(x) Clay 2
Overall, Figures 2(a) to 2(x) indicate that all models performed very well for all data sets except for the BC model because of its discontinuous character. If we compare the fitted curve of the MG and BC models, we can see that the MG model showed a better fit over the entire range of available observed data. The BC model typically had problems matching data in the transition point near saturation and mostly missed the shape of the data near saturation in almost all the soil samples, which decreased the flexibility of the curve. Therefore, this suggests that the MG model is more consistent when applied to different soils than is the BC model. As can be seen in Figures 2(a) to 2(x), the difference in the quality of fitting data between the MG and VG models is only marginal. As a result, out of 24 soil samples, the MG model fitted the observed data marginally better than the VG model in 22 soil samples. The VG model was slightly better than the MG model at fitting the other two soil samples (Figure 2(c) loamy sand 1 and Figure 2(j) silt 2). When we look at the fitted curve of the MG and LN models, both models predict almost identical curves and show much closeness to each other. However, if we compare the MG and LN models in terms of accuracy, we can see that the MG model showed a better fit for all 24 soil samples. These results indicate that the MG model can perform better than the LN model in terms of the goodness of fit for all soil types. Based on the data sets of soil in this study, it is therefore suggested that the exponential-based MG model performs much better than the BC, VG, and LN models and it is suitable for describing the SWCC for the 12 soil types.
5. Conclusion
In this paper, a simple modification of the Gardner model to form an exponential-based model has been performed to describe the soil-water characteristic curve. The exponential form of the modified Gardner model was then tested to show its performance and compared with the three most widely used models (BC, VG, and LN). Using a limited group of 24 soil samples from the UNSODA database which represents 12 different soil types, we analyzed and compared the four SWCC models with four fitting parameters. The models were evaluated in terms of their goodness of fit by means of statistical indexes. Results obtained show that out of 24 soil samples, the MG model fitted the observed data better than the LN, VG, and BC models for 22 samples. Therefore, the modified Gardner model was found to be the most consistent model fitting the observed data for the 12 soil types. The major achievement of the MG model is its ability to predict the entire range of available SWCC data with a wider range of soils and texture data. These results indicate that the exponential-based MG model performed better than the other models in terms of goodness of fit and matching the observed data. Nevertheless, additional research on testing the model performance under a complete range of water contents, from saturation to oven dryness, will be needed to evaluate the predictive capability of the model for the entire regions of the SWCC.
In summary, the advantages of the modified Gardner model are the following: the model parameters have a physical meaning which is related to the shape of the SWCC; the effect of one parameter can be distinguished from the effect of the other parameter; the model provides a wide range of flexibility to better fit data from a variety of soil types; and the model contains only four parameters. Given the results of the study, it is therefore suggested that the exponential-based MG model should be used as an alternative to the soil-water characteristic curve.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to acknowledge Almushfi Saputra for his assistance in the preparation of this paper. The present work has been supported by the Universiti Kebangsaan Malaysia under Grant GUP 2013-051.