Predicting the soil moisture retention curve , from soil particle size distribution and bulk density data using a packing density scaling factor

A substantial number of models predicting the soil moisture characteristic curve (SMC) from particle size distribution (PSD) data underestimate the dry range of the SMC especially in soils with high clay and organic matter contents. In this study, we applied a continuous form of the PSD model to predict the SMC, and subsequently we developed a physically based scaling approach to reduce the model’s bias at the dry range of the SMC. The soil particle packing density was considered as a metric of soil structure and used to define a soil particle packing scaling factor. This factor was subsequently integrated in the conceptual SMC prediction model. The model was tested on 82 soils, selected from the UNSODA database. The results show that the scaling approach properly estimates the SMC for all soil samples. In comparison to the original conceptual SMC model without scaling, the scaling approach improves the model estimations on average by 30 %. Improvements were particularly significant for the fineand medium-textured soils. Since the scaling approach is parsimonious and does not rely on additional empirical parameters, we conclude that this approach may be used for estimating SMC at the larger field scale from basic soil data.


Introduction
Increasing contamination of the groundwater resources has profoundly accentuated the need for accurate predictions of subsurface flow and chemical transport.Water flow and sub-sequent chemical transport are largely determined by the soil hydraulic properties, such as the soil moisture characteristic curve (SMC) (Wang et al., 2002;Mohammadi et al., 2009).Measuring the soil hydraulic properties is still difficult, laborintensive, and expensive.Therefore, many researchers have made an attempt to develop an indirect method as an alternative to the direct measurement of hydraulic properties.For the SMC, indirect methods are classified into conceptual methods (Nimmo et al., 2007;Mohammadi and Vanclooster, 2011), semi-physical methods (e.g., Arya and Paris, 1981;Haverkamp and Parlange, 1982;Wu et al., 1990;Smetten and Gregory, 1996) and empirical methods (e.g., Saxton et al., 1986;Schaap et al., 1998).
The semi-physical methods are mainly based on shape similarity between the SMC and the particle size distribution (PSD) curve (Zhung et al., 2001;Schaap, 2005;Haverkamp et al., 2005;Hwang and Choi, 2006), implying that the poresize distribution (PoSD) is closely related to the PSD (Arya et al., 2008).Arya and Paris (1981) did a pioneering work (AP model) for the development of semi-physical models.They showed that the pore size, which is associated with a pore volume, is determined by scaling the pore length, using a scaling factor, α.They demonstrated that an average value of 1.38 for α scales the pore lengths based on spherical particles to natural pore lengths properly.However, later investigations by Arya et al. (1982), Tyler and Wheatcraft (1989), Basile and D'Urso (1997) and Vaz et al. (2005) revealed that α value varies between 1.02 and 2.97 for fine-and coarsetextured soils, respectively.A slight error in the estimation Published by Copernicus Publications on behalf of the European Geosciences Union.
of α may result in considerable error in predicting the SMC (Schuh et al., 1988).Schuh et al. (1988) found that the value of α varies with soil texture and suction head, especially in the wet range of sandy soils.Using three formulations of α, Arya et al. (1999) modeled the parameter α as a function of particle sizes and showed that α was not constant.It decreased with increasing particle size, especially for the coarse fractions.Tyler and Wheatcraft (1989) showed that the parameter α is equivalent to the fractal dimension of a tortuous fractal pore.
Although the empirical methods have been developed extensively (e.g., Puhlmann and von Wilpert, 2012), the performance of an empirical method will depend on the databases being used for the model calibration and testing (Tietje and Tapkenhinrichs, 1993;Kern, 1995;Schaap and Leij, 1998;Schaap et al., 2004;Haverkamp et al., 2005;Hwang and Choi, 2006;Weynants et al., 2009).Moreover, direct measurements of SMC are often integrated as predictor variables of the continuous SMC function.Many attempts have been made to reduce the sensitivity of the indirect methods to empirical and database-dependent parameters.For instance, Mohammadi and Vanclooster (2011) proposed a conceptual robust model (MV) that does not include an empirical parameter and is independent of the databases that are being used.The disadvantages of semi-physical or conceptual models such as the AP and MV models are the use of "bundle of cylindrical capillaries" (BCC) concept to represent the pore space geometry and the lack of consideration of surface forces (Or and Tuller, 1999;Tuller et al., 1999;Mohammadi and Meskini-Vishkaee, 2012).These conceptual problems often lead to the underestimation of the dry range of the SMC (Arya et al., 1999;Hwang and Choi, 2006;Mohammadi and Vanclooster, 2011).Such underestimations would result in large modeling errors of hydraulic dependent soil functions such as mechanical resistance functions (Gras et al., 2010), plant water uptake functions (Ryel et al., 2002), and microbial activity functions (Jamieson et al., 2002;Santamarí a and Toranzos, 2003), in particular in arid environments.
To predict continuous SMC, Naveed et al. (2012) parameterized the van Genuchten model based on the SMC data points predicted from organic matter, clay, silt, fine sand and coarse sand content.Mohammadi and Meskini-Vishkaee ( 2013) integrated the MV model with the van Genuchten (VG) model (van Genuchten, 1980) to predict the continuous SMC curve (MV-VG model) from PSD data.They found that ignoring the residual moisture content (θ r ) is the main source of systematic error in the MV model.They further tested and compared four approaches to predict the θ r , and showed that the incorporation of predicted θ r will improve the MV-VG prediction results considerably.However, the estimation of θ r has some limitations, due to the lack of a conceptual underpinning and the poor predictability of θ r (Leij et al., 2002).Tuller and Or (2005) suggested that the introduction of θ r as a fitting parameter in most SMC models often makes the physical representation of key pro-cesses in the dry soils vague.Moreover, they pointed out that the dry range of the SMC shows remarkable scaling behavior.Arya et al. (2008) developed a procedure to scale natural pore lengths, directly from straight pore lengths.They showed that the scaling approach is less sensitive to uncertainties in model parameters and provides better predictions of the SMC, compared with the original AP model.Kosugi (1996) showed that the SMC can be expressed by a lognormal pore-size distribution function, while Kosugi and Hopmans (1998) found that the set of scaling factors is lognormally distributed when PoSD curve is lognormal.Havayashi et al. (2007) used the Kosugi model (Kosugi, 1996) to evaluate the effectiveness of three kinds of scaling factors obtained by the microscopic characteristic length, standard deviation of pore-size distribution and the porosity.They indicated that, in the natural forested hillslope soils, the variability in the SMC is characterized by variability in the effective soil pore volume.Nasta et al. (2009) concluded that the scaling of the PSD curves provides for adequate characterization of the mean and variance of SMCs, which allows for characterization of the soil spatial variability.
Many researchers developed empirical models for expressing the SMC since the parameters of these models do not address the physical significance of the medium.Hence the spatial variability in the pore structure of soils is not fully understood (Havayashi et al., 2007).Likewise, the conventional scaling approaches are based on empirical curve fitting, without considering the physical meaning of the scaling factor (Perfect, 2005;Millán and González-Posada, 2005).To apply these models, one needs to determine the scaling factor, where the complexity of measurements of the poresize and pore-volume distributions easily nullifies the estimation of the scaling factors.Nevertheless, some efforts have been made to relate the scaling factor to the soil texture (Tuli et al., 2001;Millán et al., 2003).
From this brief review we conclude that the scaling approaches improve the modeling and prediction of the SMC.Yet, most scaling approaches imply empirical parameters, and a robust fully conceptual approach for the estimation of the SMC from easily measurable properties is still lacking.
The MV model underestimates the moisture content in the dry range of the SMC because of the simplified conceptualization of the pore geometry.In particular the packing parameter does not effectively reflect the pore geometry.The general aim of this study is to improve the accuracy of the model proposed by Mohammadi and Meskini-Vishkaee (2013) using a scaling approach.
Therefore the objectives of this study are (i) to formulate a robust and physically based model to scale the SMC from the PSD and porosity, and (ii) to compare the model performance with the results from the existing MV-VG model, using soils documented in the UNSODA database (Nemes et al., 2000).We also evaluate the overall model performance with the results from the full empirically SMC prediction software ROSETTA (Schaap et al., 2001).

Theory
Because of the close similarity between the shapes of the PSD and SMC curves, many researchers expressed an SMC model in terms of a PSD model (Haverkamp and Parlange, 1982;Fredlund et al., 2000;Zhuang et al., 2001).The SMC model developed by van Genuchten (1980) is very flexible, widely used and given by where θ (L 3 L −3 ) is the soil moisture content, S e (−) is effective saturation degree and θ s (L 3 L −3 ) and θ r (L 3 L −3 ) are saturated and residual soil moisture contents, respectively.The parameters n, m, θ r and α (L −1 ) are fitting coefficients, and h (L) is the suction head.The suction head, h i (L), corresponding to the particle radius of the ith fraction R i (L) is given by (Mohammadi and Vanclooster, 2011) where ζ (−) is a coefficient depending on the state of soil particles packing and is defined as where e (−) is the void ratio given by where the ρ s (ML −3 ) and ρ b (ML −3 ) are soil particle and bulk densities respectively.Arya and Paris (1981) suggested that the moisture content, θ i (L 3 L −3 ), can be obtained from PSD and θ s (L 3 L −3 ), as where w j is the mass fraction of particles (−) in the j th particle size fraction.Consider that would result in where S (−) is the saturation degree and P i (−) is the cumulative mass fraction of soil particles.It is obvious that if θ r = 0, then S e = S and subsequently S = P i .Arya and Paris (1981), however, ignored the residual moisture content, while it may be a considerable value for many types of soil and clayey soils in particular.Combining Eqs. ( 1) and (3) with Eq. ( 7) yields In Eq. ( 9), the cumulative mass fraction, P i , is substituted with the S e in Eq. ( 1).Hence, fitting Eq. ( 9) to the PSD data enables one to directly predict the SMC parameters (n, m and α).Moreover, these coefficients allow expression of the continuous form of predicted SMC.Since assuming that θ r = 0 would result in model underestimation in dry range of the SMC (Mohammadi and Meskini-Vishkaee, 2013), we developed a conceptual scaling approach to reduce the model bias.

Scaling approach
Following Havayashi et al. (2007), we suggest that the porosity is an appropriate property for inferring a characteristic scaling factor.Since the soil porosity is linked to the packing parameter, ζ , in the MV model (Eq.4), we hypothesize that ζ is the characteristic scale of the soil.We assume that the reference soil is the one that consists of uniform-size spherical particles that are arranged in random close packing state, leading to minimum porosity (known as the Kepler conjecture in literature of crystallography).Literature suggests that the porosity of this packing state is 0.259 (Hopkins and Stillinger, 2009).Subsequently, the maximum value of packing parameter, ζ max , would equal 1.41432 for reference soil.Hence the scaling factor, λ, for each soil sample can be suggested by In general, the values of the pore-size distribution index, n (Eqs. 1 and 9), and ζ are large for the coarse-textured soils and small for the fine-textured soils.We suggest that the λ can scale the parameter n, obtained from fitting Eq. ( 9) to the PSD data, to the n parameter in the SMC model (Eq. 1) (hereafter n * ) as follows: where n * is scaled to the PoSD index in VG model.Hence, the modified model is In summary, given a known θ s we can calculate ζ and subsequently λ using Eq. ( 10).The soil parameters α and m are obtained from fitting Eq. ( 9) to the PSD data, and n * is estimated by Eq. ( 11), and consequently the SMC is predicted directly by Eq. ( 12).

Material and methods
A total of 82 soil samples, with a wide range of physical properties that contained at least five PSD data, were selected from the UNSODA hydraulic properties database (Nemes et al., 2000).UNSODA is a database of basic soil and hydraulic properties from 790 samples, gathered from all over the world, and compiled by the US Department of Agriculture.All soils are summarized in Table 1.
In this procedure, volumetric moisture contents corresponding to the ith fraction were computed using Eq. ( 6), and suction heads were predicted using Eq.(3), in which the parameter ζ was calculated with Eq. ( 4).In this study, we assumed that the porosity is equivalent to θ s .For soils that neither provide a porosity nor a θ s , the first point of the SMC data that corresponds to the lowest suction head was used as θ s (Chan and Govindaraju, 2004).
We fitted Eq. ( 9) to the PSD data.We used nonlinear regression analysis to fit Eq. ( 9) to the PSD, using Matlab 7.1 software (Matlab 7.1, The Mathworks Inc., Natick, MA, USA) and the Marquardt-Levenberg algorithm (Marquardt, 1963).We calculated for each soil the scaling factor, using either the bulk density or the available saturated soil moisture content, and predicted the SMC.
For each prediction, the agreement between the predicted moisture content θ i(p) and measured moisture content θ i(m) was expressed in terms of the root mean square errors (RM-SEs), given by in which N is the number of observed data points.The relative improvement (RI) resulting from the scaling approach rather than MV-VG model was calculated as follows (Minasny and McBartney, 2002): where RI is the relative improvement, RMSE M and RMSE s are RMSE values of the MV-VG model and the current scaled model, respectively.Obviously, a negative RI value indicates that the scaling approach would diminish the accuracy level of the prediction of the SMC in comparison with the MV-VG model.We also fitted a cubic polynomial function to the overall predicted data and calculated the area between the fitted polynomial and the 1 : 1 line from the difference of the numerical integrals of these two functions (do Carmo, 1976).
Moreover, to consider and compare the reliability of the scaled MV-VG model with a fully empirical SMC prediction model, we compared the estimations of scaled MV-VG model with the estimations of the ROSETTA software.In this software, we used the SSCBD model option; that is, we used textural (sand, silt, clay) percentages and bulk density as model predictors (Schaap et al., 2001).

Results and discussion
Table 2 gives the comparison between the MV-VG model, the ROSETTA model and the scaled MV-VG model in terms of RMSE, R 2 and RI.Table 2 demonstrates the significantly improved accuracy of the scaled MV-VG approach as compared to the original MV-VG model and ROSETTA.The RMSEs of the predicted and measured moisture contents ranged from 0.0223 to 0.1502 for the original MV-VG model (average 0.086), from 0.0169 to 0.1122 (average 0.0601) for the scaled approach and from 0.0188 to 0.2453 (average 0.745) for the ROSETTA software.In terms of RM-SEs, the scaled approach performed better than Schaap et al. (1998) and the Schaap and Leij (1998) models with similar predictor variables.The results showed that there is a significant difference between performance of scaled MV-VG   approach and ROSETTA (p = 0.05).Despite the pure statistical and empirical nature of the ROSETTA approach, it provides worse results than the approach based on the current scaling technique.The improvement of the scaled approach is also reflected by RI in Table 2. Except for soils no.3033 (clay loam) and 3090 (silt loam), the scaling approach resulted in more accurate predictions for all soils.Table 2 also indicates that the scaling approach can improve the model estimations of the original MV-VG model by 30 %.
For the fine-and medium-textured soils, the values of RI are larger than for the coarse-textured soil.This result was expected, because the MV and MV-VG models underestimate the dry range moisture content for the fine-texture soils (Mohammadi and Vanclooster, 2011;Mohammadi and Meskini-Vishkaee, 2013), and subsequently the scaling approach was more effective for these soils.
We examined the possible relations between the RI and soil physical properties.Among all parameters, the saturated moisture content and scaling factor show strong relations with the RI. Figure 1a shows that the RI values increase significantly with the saturated moisture content of the soils; that is, the scaling approach would more effectively improve the model accuracy for the fine-texture soils with higher θ s .This result can be confirmed with Fig. 1b, which shows that the scaling factor is inversely correlated with the IR factor (Fig. 1b).Indeed, the soils with high porosity commonly have an abundant amount of clay materials and organic matter, characterized with high surface energy.These attributes are the main sources of errors of the MV and MV-VG models.
Typical examples of measured vs. predicted SMCs with the MV-VG model, the scaling approach and ROSETTA for clay, sandy loam, loam, and silt loam textures are presented in Fig. 2a-f.For the clay (codes: 2340 and 4681), sandy loam (code: 3180 and 3200), loam (code: 3191) and silt loam (code: 3090) soils, the scaling approach fits the data well and performs better than the MV-VG model in the entire range of the SMC.For the silt loam soil (code: 3090), the scaling approach slightly overestimates the moisture content through the entire range of suction heads and the MV-VG model underestimates the moisture content at low suction heads.Figure 2a-f show that the ROSETTA software performs worse in the wet part of SMC.Overall, the scaling approach performs better than MV-VG model and ROSETTA software for all soil samples (Table 2), but the performance of scaling approach did not suitably respond for two soil samples (codes: 3033, 3090).The residual model error may be related to the simplified representation of the total porosity which is considered equal to the saturated volumetric moisture content.The swelling properties and high organic carbon content of these soils (> 4 %, 3.85 %, respectively) may partially be a source of these errors.We further suspect that the complexity of the relationship between PSD, PoSD and pore connectivity can be effective in the model performance (Zhuang et al.,

2001).
The assumption of the similarity between PSD and PoSD does not perform equally well to all soils.
We tentatively conclude that the scaling of the PSD curves using the parameter ζ generally performs better in predicting the SMC as compared to the original MV-VG model.The unscaled MV-VG model underestimates the moisture content at high suction heads.
The most semi-physical based methods for predicting SMC rely on the use of empirical parameters to improve the SMC estimates from PSD (Lilly and Lin, 2004).Hydraulic properties are indeed affected by both the soil texture and the soil structure (Haverkamp et al., 2002).The MV-VG model uses the packing parameter, ζ , derived from soil bulk density as a metric of soil structure.Moreover, the scaling parameter that is inferred from the packing state is integrated in the scaled MV-VG model.Hence the soil structural features are integrated in the MV-VG model at two levels: first at the MV-VG model to convert moisture into pressure head and, second, to correct the SMC model prediction.The good performance of scaling approach in the wet range and dry range of SMC suggests a convenient of soil structural features in the SMC prediction.comparing the experimental data and the predicted soil moisture content on a 1 : 1 plot.Linear regression of the measured and estimated moisture contents, using the two methods for all the soil samples, showed that the slope values were 0.7675 and 0.8484, and the coefficients of determination (R 2 ) between the estimated results and measured data for all soils were 0.765 and 0.8565 for MV-VG model and scaling ap-proach, respectively.Hence, the MV-VG model and the scaling approach underestimate the moisture content by about 23 and 15 %, respectively, while the bias of the scaling approach is smaller than the MV-VG model.Since it has been reported that usual measurement method of SMC (pressure plate apparatus) is susceptible to some errors at high soil suction heads (Campbell, 1988;Gee et al., 2002;Cresswell Figure 3. Comparisons of the measured and estimated moisture contents for 82 selected soils using the MV-VG model (Eq.9) and scaling approach (Eq.12).Dashed lines: the 1 : 1 line.Solid lines: linear-regression line.Solid green line: nonlinear-regression line.et al., 2008).We suggest that a part of the underestimation in the dry range of SMC of our method is partially related to limitation of this method for measuring the SMC (Solone et al., 2012).Regarding the R 2 , the scaling method still remains the most preferable method.Comparing the overall predictability of the two methods, the correlation coefficient of linear regression should not replace the visual examination of the data.We, therefore, use the cubic polynomial function to adequately express the data variations.The fitted polynomial functions are drawn as shadowed green curves in Fig. 3.We further suggest that the area between the fitted curve and the 1 : 1 line (AE) is an expression of the systematic error.The values of AE were 0.0369 and 0.0250 for the MV-VG model and the scaling approach respectively.It confirms that the level of systematic errors of the scaling approach is about 33 % less than that of MV-VG model.This result can be confirmed by the comparison of the R 2 values obtained when predicting the SMC for each soil with two methods (Table 2 columns 5 and 6).We conclude that the scaled PSD curve will result in a more accurate prediction of the SMC as compared to the unscaled PSD data.
To scale the SMC, Tuller and Or (2005) used the soil specific surface area (SA) and the thickness of water film to express the moisture content in dry range of the SMC.Despite their reasonable model performance, the application of their procedure was limited due to difficulty in the measurement or the estimation of SA.Havayashi et al. (2007) found that, in natural forested hillslope soils, the variability in the SMC is scaled and characterized by the variability in effective porosity.Nevertheless, the determination of the effective porosity is also difficult.
The scaling factor proposed in the current study is defined by using the index of packing state which can be determined easily from the bulk density of soil and particle density.As compared to prediction models that rely on the measured attributes as suggested above, or prediction models that rely on measured SMC data, our approach is based on a robust metric of the soil structure: the packing density.It does not rely on any other additional empirical parameter.The scaled MV-VG model is therefore very parsimonious and robust.We therefore conclude that the scaled MV-VG model may be appropriate for predicting SMC from basic soil data.

Conclusions
Using a new scaling approach, the current study showed that the continuous form of SMC curve can be predicted from knowledge of PSD, as modeled by the van Genuchten (1980) model and particle packing state.In this approach it was assumed that the scaling factor can be defined as the ratio of packing state of a soil sample and the packing state of a reference soil.Results showed that the proposed approach can adequately predict the SMC of 82 soil samples selected from the UNSODA database.It was further found that the scaling approach provides better predictions of the SMC than MV-VG model and ROSETTA software, especially in the dry range of the SMC.For soils for which the error was important, we attributed the proposed scaling approach error to high organic carbon content and swelling properties of the soil.Indeed, in these soils the soil pore structure and porosity is changing in time, leading to uncertainty in the scaling factor based on the soil porosity.
In summary, we concluded that the main advantages of the proposed scaling approach as compared to many SMC prediction models are the following: (i) the applied scaling factor is determined easily from soil bulk and particle densities; (ii) the scaling factor has physical meaning, which does not depend on soil database and empirical parameters; (iii) the proposed approach predicts a continuous form of the SMC; and (iv) this approach estimates the SMC more appropriately in comparison with many other models.Considering that there is no further need for empirical parameters, we conclude that this approach may be useful in estimating the SMC for regional-scale soil hydrological studies.
), and relative improvement (RI) compared to the MV-VG model and hydraulic parameters for each soil textural group, with standard deviations in parentheses.The lowercase letters a and b indicate significant differences at P <

Figure 1 .
Figure 1.The efficiency of scaling approach, % RI, defined with RMSE (Eq.14) as function of (a) the saturated moisture content and (b) scaling factor for all soil samples.* * : significant at P = 0.01.

Table 1 .
Textural classes and UNSODA codes for soils used for testing and evaluating the approach.

Table 2 .
Average root mean square error (RMSE), coefficients of determination (R 2

Table A1 .
Symbols and abbreviations.