Modeling the Soil Water Retention Curves of Soil-Gravel Mixtures with Regression Method on the Loess Plateau of China

Soil water retention parameters are critical to quantify flow and solute transport in vadose zone, while the presence of rock fragments remarkably increases their variability. Therefore a novel method for determining water retention parameters of soil-gravel mixtures is required. The procedure to generate such a model is based firstly on the determination of the quantitative relationship between the content of rock fragments and the effective saturation of soil-gravel mixtures, and then on the integration of this relationship with former analytical equations of water retention curves (WRCs). In order to find such relationships, laboratory experiments were conducted to determine WRCs of soil-gravel mixtures obtained with a clay loam soil mixed with shale clasts or pebbles in three size groups with various gravel contents. Data showed that the effective saturation of the soil-gravel mixtures with the same kind of gravels within one size group had a linear relation with gravel contents, and had a power relation with the bulk density of samples at any pressure head. Revised formulas for water retention properties of the soil-gravel mixtures are proposed to establish the water retention curved surface models of the power-linear functions and power functions. The analysis of the parameters obtained by regression and validation of the empirical models showed that they were acceptable by using either the measured data of separate gravel size group or those of all the three gravel size groups having a large size range. Furthermore, the regression parameters of the curved surfaces for the soil-gravel mixtures with a large range of gravel content could be determined from the water retention data of the soil-gravel mixtures with two representative gravel contents or bulk densities. Such revised water retention models are potentially applicable in regional or large scale field investigations of significantly heterogeneous media, where various gravel sizes and different gravel contents are present.


Introduction
A large proportion of soils containing rock fragments are present in the world due to soil evolution and erosion [1,2]. The highly variable gravel content or size in soilscape greatly increases the variability of the soil properties [3][4][5]. Knowledge of soil water retention curves (WRCs) is a prerequisite for modeling the fluxes of water and solutes in the vadose zone and consequently it is necessary to determine their spatial variability [6][7][8]. Since direct field measurements of WRCs are time-consuming and expensive, laboratory measurements continue to be the most frequent means of characterizing the vadose zone [9,10]. Soil water retention data are typically obtained in laboratory for fine soils (,2 mm) using pressure cells, pressure-plate extractors, and centrifuge methods [11][12][13]. Several reports noted that water held by gravel cannot be neglected in the determination of water retention properties due to the significant porosity of gravel and the changed pore-size distribution [14,15]. For example, the ironstone gravel contained a large amount of available water ranging from 0.03 cm 3 cm 23 to 0.15 cm 3 cm 23 [16], while the sandstone fragments and shale fragments held 0.11 cm 3 cm 23 and 0.23 cm 3 cm 23 available water, respectively [17]. For soil-glass mixtures, the volume of coarse lacunar pore increased with glass content when glass content was less than 50% [18]. As a large number of water retention curves for fine soils have been obtained in laboratory, gravel corrections for moisture retention in soil-gravel mixtures have been developed on the basis of WRCs of fine soils. For example, Gardner [19] used mass-based approach while Bouwer and Rice [20] used volume-based approach to make gravel corrections for water retention of soil-gravel mixtures [20][21][22]. Correction procedures are available to determine water retention properties for soil-gravel mixtures, but they have limited utility for soil-gravel mixtures with weathered gravels especially on highsuction range [23]. Pedotransfer functions for predicting WRCs have been developed from more easily measurable and more readily available soil properties [24][25][26][27]. Scheinost et al. [6] predicted WRCs for soils with a wide range of particles which included 2-67 mm gravel with a new pedotransfer function including more textural fractions, but this did not greatly improve the prediction precision [25,28].
In order to quantify the effects of content and size of rock fragments on soil water retention and finally develop a new method for determining the WRCs of soil-gravel mixtures, especially for mixtures with weathered gravels and for WRCs at high-suction range, we investigated WRCs of a loess soil and gravels mixtures with various gravel contents and gravel sizes. The following constraints were soon evident in developing this new method: (1) To extend the practical use of the revised model, the water retention data of soil-gravel mixtures with weathered gravel should be obtained and the range of water potential should be extended to include very low values (high-suction). (2) To develop revised water retention models, the effects of gravel contents and gravel size on the effective degree of saturation (Se) of soil-gravel mixtures should be analyzed based on the measured WRCs data. The relationships between Se and gravel contents or bulk density of soil-gravel mixtures should be combined with a closed-form analytical equation such as the model of Brooks and Corey (BCfunction) [29] or the van Genuchten equation (VG-function) [30].
(3) To validate the effectiveness of revised water retention models and simplify the procedure of parameter-obtaining, the way that the gravel size affects the shape parameters of revised models should be analyzed, and the shape parameters of the revised models should be obtained from the representative WRCs data of soil-gravel mixtures, such as that with two extremes values of gravel content range. Those parameters and the method of parameter-obtaining may give references for determining WRCs of other soil-gravel mixtures.

Ethics statement
No specific permissions were required for these sampling activities because the location is not privately-owned or protected in any way and the field activities did not involve endangered or protected species.
The samples were excavated from the soil profiles in 30 cm depth of Yaoxianliang in Tongchuan (108.93 E, 35.28 N, at 1570 m altitude) and Weihe river bank in Yangling (108.10 E, 34.25 N, at 437 m altitude), Shaanxi province, China. The soil in Yaoxianliang is recognized as Aric Regosol (FAO). The rock fragments in these soils are schists and shales (S) with brownish green in color and sheet-like shape. The other type of rock fragments sampled in Weihe river bank was pebble (P) from alluvial sediment. The rock fragments in sizes ranging from 2 to 10 mm were sieved and washed to remove the soil particles from their surface, and then sieved into three different diameter classes: 2-3 mm, 3-5 mm, and 5-10 mm. The weathering degree of the rock fragments was defined by observing their weathering characteristics in the field and observing the surface fissures and coarseness under a magnifier in the laboratory according to the Geotechnical Engineering Handbook [31]. The gravel weathering degree is described in Table 1.
The texture of the air-dried fine soil (,2 mm), from which the rock fragments (. 2 mm) and plant residues were removed, was measured by the Laser Diffraction Particle Size Analyzer (Mastersizer 2000, Malvern Instruments Ltd. in England). According to the International Soil Classification System, the disturbed soil sample was a clay loam soil (labeled with CL). The mean density of the rock fragments (r r ) in 2-10 mm, determined   (Table 1). It should be noted that the rock fragments were saturated in water before determining the water volume substituted by the rock fragments in a container with scale when measuring r r . The saturated water contents (h s ) of the gravels soaked into water for three days were measured by oven drying at 105uC for 24 h and the measurements of each kind of gravels were replicated five times. The results showed that h s of shale clasts was 0.0960.02 g g 21 , while that of pebbles was almost zero. In the sample saturation process, the air bubbles enclosed in the gravels possibly caused the gravel in incomplete saturated and subsequently the underestimation of h s in this study. The clay loam soil was mixed with air-dried shale clasts or pebbles at the ratios of 10%, 15%, 25%, 35%, 45%, 55%, and 65% on a total mass basis in each of the three size classifications (2-3 mm, 3-5 mm, and 5-10 mm). We totally prepared 45 samples including one soil sample, two gravel samples, and 42 soilgravel mixture samples. The soil-gravel mixtures were packed uniformly by hand into soil containers (98.2 cm 3 with 5 cm high and 5 cm inner diameter) and the bulk density of clay loam soil without gravel predetermined 1.28 g cm 23 for all the soil and soilgravel mixture samples. Then the soil containers were posited inside the centrifuge rotor chamber (Hitachi-CR21G, Hitachi Ltd. in Japan) under dry conditions (i.e., water content was less than 0.5%). The bulk density of soil-gravel mixtures is presented in Table 2; it was calculated from the measured dry mass and the volume of the samples after finishing centrifuge rolling.
The WRCs of soil, gravels, and soil-gravel mixtures were measured separately by the centrifuge method [32], and each treatment had three replications. The tested suction head in the WRCs measurement was in the range of 102 cm to 10200 cm. The Hitachi-CR21G centrifuge was equipped to maintain air temperature constantly at 20uC. For the centrifugation method, soil water desorption was accomplished by applying a high gravity field (centrifugal force) to saturated soil samples. The suction heads were obtained by sequentially increasing the angular velocity or rotation speed of the centrifuge [33]. The rotation speeds and equilibrium times corresponding to the tested pressure heads are shown in Table 3. Water retention properties of the samples near saturation were not considered because the water retention of the artificially packed samples was conditioned by the samples packing and the geometry of soil samples, especially in the low matric suction region [34].
For those h-h data, BC-function and VG-function parameters describing WRCs of soil, gravels, and soil-gravel mixtures were determined by the computer program RETC.FOR [35] using the Marquardt nonlinear least-squares optimization algorithm [36]. BC-function and VG-function are presented as Eq. (1) and Eq. (2):  Table 4. BC-function parameters describing the water retention characteristics of soils, soil-gravel mixtures, and gravels.

Samples
Gravel content (%) Gravel size in 2-3 mm Gravel size in 3-5 mm Gravel size in 5-10 mm where S e is the effective degree of saturation, also called reduced water content; h s is the saturated water content (cm 3 cm 23 ); h r is the residual water content (cm 3 cm 23 ); h is the matric suction (cm); a is the empirical parameter whose inverse equals to h a and is referred as the air entry value or bubbling pressure (cm 21 ); and l is the pore-size distribution parameter; n, m, A = a 2l are empirical parameters affecting the shape of the retention curve, and m = 1-1/n.

Effects of rock fragments on WRCs
Soil matric suction (h) is a function of water content (h) in an unsaturated soil and WRCs express the relationship between them. The VG-function resolves the coherence problem of WRCs near soil saturation and therefore it has become one of best choices for the analytical model for h(h), especially for undisturbed field soil and many fine-texture soils [37,38]. While BC-function, leads to an air-entry value in the WRC above which soil is assumed to be saturated, has the ability of more accurate description for coarse texture soil with structural deterioration or compaction, especially in the dry end of WRCs [39]. However, in this study, they gave similar fine fitting performance ( Fig. 1 and Table 4).
The maximum relative standard error (RSEm) for the saturated water contents of each sample was lower than 2.5% (n = 3), and the RSEm for the unsaturated water contents of them at various pressure heads was lower than 2.2% (n = 3), which showed a good representation for each sample; the measured WRCs data by the centrifuge method had acceptable accuracy for soil-gravel mixtures. While, the measured data of water retention for the soil-gravel samples might be underestimated due to the experimental set up, in which the gravels was possibly not moistured completely under normal air pressure conditions. Parameter sensitivity at one pressure head for WRCs was evaluated by the average value of the total ratios of the relative changes of h to the relative changes of values of one parameter. The results of sensitivity analysis showed that the absolute value of sensitivity for parameter a in BC-function increased with soil suction till it reached 1/a, and then kept at that maximum value. Correspondingly, that value in VG-function increased with soil suction till it reached the dry end of curve (h = 10200 cm, except for the mixture samples containing 55% gravels, it reached at h = 407 cm). Parameter n became most sensitive at h = 10200 cm for both BC-and VG-function. The increasing sensitivity of a and n to BC- Table 6. Regression empirical parameters for curved surface functions with r obtained by fitting all the measured water retention curves.  and VG-function with the increased suction implied that the water content at h = 0 and h = 100 cm decided the trend and shape of WRCs in the wet end (h = 0-100 cm), which varied insignificantly when the parameters were changed. Those results indicated that the comparison between WRCs combining extrapolated and measured one for soil-gravel samples was acceptable. It seems that a reflected the order of air-entry values from those a values (lower, media, upper value at 95% confident limits) for various soil-gravel samples due to its high sensitivity to BC-function at h = 1/a. One thing should be noted that the measured h s and the fitted h r in BCand VG-function were as the fixed value in the analysis of parameter sensitivity. There were considerable differences between the water retention parameters of soil-gravel mixtures with different gravel contents ( Table 4). The air-entry values (1/a) of the soil-gravel mixtures decreased gradually with increasing gravel content when it was lower than 55%; the reason may be that the amount of coarse pores increased with gravel content (Although the gravels were packed by hand to increase uniformity, there was still some pores in the soil-gravel mixtures because these gravels possibly overlapped each other and functioned as skeleton). However, when the gravel content reached 65%, especially for the soil-gravel mixtures containing the larger gravels (5-10 mm), the air-entry values increased slightly. It was further noted that the stronger the weathering degree of rock fragments in the soil, the higher the airentry values, at least in the range of the higher coarse fragment contents for soil-gravel mixtures. Apparently, the shale clasts with smaller sizes had more void characteristics similar to fine soil medium than the pebbles. This observation was confirmed by the fact that the h s of clay loam-shale clast mixtures (CL+S) were larger at any gravel content than those of the clay loam-pebble mixtures (CL+P), while this effect was probably aggravated by the side-effect (side conditions in the container) in the experiment. It is reasonable to speculate that the measured h s of the soil-gravel mixtures decreased linearly with the increase of rock fragment contents because the majority of water was held by the fine soil of soil-gravel mixtures. The variation of h r in the soil-gravel mixtures in this research differed from that reported by Indrawan et al. [40], according to whom the residual water content decreased with increasingly coarse fragments.
The WRCs for fine soil, gravel, and soil-gravel mixtures with different gavel contents (15%, 35%, and 65%) are presented in Fig. 1. This figure only displays the results of soil-gravel mixtures with 3-5 mm gravel because of the similarity among the different size groups. It can be seen that the WRCs of the soil-gravel mixtures are located between the WRCs of the pure soil and those of the rock fragment media. According to the fitting curves of VGfunction, the slopes of the WRCs for the soil-gravel mixtures increased with the coarse grain contents at a low suction range (0-100 cm). That suggested an increase in the rate of the soil water volume changes with respect to the matric suction changes (hh/ hh). However, the trend to increase in slopes of the WRCs for the soil-gravel mixtures became less significant with increasing rock fragment contents at the high suction range (100-1000 cm), and became opposite at a higher suction range (.1000 cm), even lower than the slope of WRCs of the fine soils, especially for the soil-pebble mixtures (see CL+P in Fig. 1). This indicated that a soil containing a high amount of rock fragments at lower water contents releases less water than the others, which might impair plant growing.

Direct calculation of the saturated water content of soilgravel mixtures
Knowing the saturated water content is a prerequisite for calculating the effective degree of saturation. At or near saturation, the moisture of the soil-gravel mixtures equals the sum of water amount hold by the fine soil and that by the rock fragments. As a result, the water content of the soil-gravel mixtures could be predicted according to the saturated water content of the fine soil and the rock fragments, respectively, as well as their content in the soil-gravel mixtures. The formula expressing the saturated water content of soil-gravel mixtures is given in Eq. (3):  whereh b s , h f s , and h r s are saturated volumetric water content of the soil-gravel mixture, fine soil, and gravel, respectively (cm 3 cm 23 ); M r is the gravimetric gravel content (g g 21 ); r b , r r , and r f is the bulk density of the soil-gravel mixture, fine soil, and gravel, respectively (g cm 23 ); and both superscripts and subscripts refer to bulk sample (soil-gravel mixture) (b), fine soil alone (f), rock fragments (r), and at saturated state (s), respectively.
Reasonably good results were obtained using Eq. (3) to predict the saturated water contents of all the samples, even those containing the weathered rock fragments. The values resulting from Eq. (3) and those actually measured were almost identical for the soil-gravel mixtures with all sizes of gravels (see Fig. 2); a small deviation between the measured water contents and the predicted values for CL+P with relative low gravel content due to measured error was recorded (see CL+P in Fig. 2). Of course, when the h r s value of the pebbles was assumed to be zero (i.e., Eq. (3) was the same as Bouwer-Rice equation) [20], the predicted results for CL+P with relative low gravel contents could not be improved.

Revised formulas for water retention processes of soilgravel mixtures
BC-function with simple power equation, introduces a welldefined air-entry value, which is associated with a largest pore-size through the relation of Young Laplace, assuming complete wettability [41]. In addition, the parameter a in BC-function could reflect the variation of air-entry value of varied soil-gravel samples even lack data of wet end of WRC. Due to these reasons, BC-function was selected to be the revised model for soil-gravel mixtures in this study. To examine the moisture content of the unsaturated soil-gravel mixtures, we generated water retention curved surfaces (curve surfaces which reflected the water retention variation of soil-gravel mixtures) by adding one variable to BCfunction, and then fitted the empirical parameters of the curved surface based on the measured WRCs data of the soil-gravel mixtures. We found that the effective saturation (S e ) of the soilgravel mixtures at any soil potential had a reasonable power or linear correlation with the gravimetric gravel content (M r ) or bulk density of the soil-gravel mixtures (r b ). Thus M r and r b were added respectively as further variables of BC-function and the revised water retention functions were established. Eq. (4) and (5) express power relations:  Table 7. Regression empirical parameters obtained by fitting two measured water retention curves of soil mixtures with 10% and 65% gravel content. where A M , A r , b, b', and l are empirical parameters. Eq. (6) and (7) express power-linear relations.
When Eq. (6) and Eq. (7) were simplified to decrease the number of the parameters, the water retention curved surfaces were expressed as Eq. (8) and Eq. (9): whereA' M , A' r , B' M , and B' r are empirical parameters. The water content of the soil-gravel mixtures can be calculated according to the revised water retention functions and the definition of S e . By taking both Eq. (5) and Eq. (8) as examples, the water content of the soil-gravel mixtures can be calculated by the formulas listed as follows.
The saturated water content of soil-gravel mixtures (h b s ) in Eq. (10) and Eq. (11) could be obtained by direct measurement in the laboratory, while it can be calculated indirectly using Eq. (3) above if pertinent detailed information is known, such as gravel content, sample bulk density, gravel density, saturated water content of fine soil, and water content of gravels.
Parameter determination for the proposed revised formulas through nonlinear regression for the soil-gravel mixtures Two variables (Mr and rb) which presented a linear or power relation with Se were added to BC-function. In order to test if those revised water retention functions were practical and applicable for obtaining the WRCs parameters of the soil-gravel mixtures, we used the measured data to fit the revised formulas and obtain the parameters of curved surfaces by a nonlinear regression.
Parameter determination for separate classes of the gravel sizes. The fitting results are given in Table 5 and Table 6. It is shown that the values of coefficient of determination (R2), estimated standard error (SEE), and square root of the residual sum of squares (RSME), which all reflect the nonlinear regression matching level, indicate that the fitting results were in considerable agreement for both the power function and the linear-power function curved surfaces with either selected variable (M r and rb). Though the values of R2 and SEE for the different curved surfaces of the various soil-gravel mixtures were similar ( Table 5 and Table 6), the difference of Norm indicated that the power function with the variable of sample bulk density (Eq. (5)) and linear-power function with the variable of gravimetric gravel content gave a better fitting (Eq. (8)). The regressed curved surfaces of power function with the variable ''bulk density'' taken as an example here are drawn in Fig. 3 showing the regression results based on the measured data. The two different types of revised formulas, with the variables bulk density or gravel content, provided options for spatial heterogeneous research of regional soil hydrology according to the types of the practical measured data. The above formulas and the empirical parameters listed in Table 5 and Table 6 concur to simplify the determination of variability in water retention properties for soil-gravel mixtures with highly variable gravel content in soilscape.
Model effectiveness and parameter determination for soil-gravel mixtures with distribution of the three classes of gravel sizes. Table 5 and Table 6 show that the empirical parameters for various soil rock fragments mixtures vary with gravel type and size. The shape parameter l for soil mixtures containing the same type of gravels with three size groups showed little difference, while by comparing the l values between the two kinds of gravels it was observed that the stronger the weathering degree of the rock fragments in the soil-gravel mixtures, the smaller the l values. The shape parameter b or B in the power function or linear-power function, respectively, was decreasing with increasing gravel sizes, while the differences in b or B between different size groups were not significant. Therefore, it seemed reasonable to consider obtaining the ''trade-off'' and appropriate parameters for the curved surfaces to account for the water retention properties of soil-gravel mixture with a large range of gravel size (2-10 mm). The fitting parameters, obtained by fitting the measured data of soil-gravel mixtures within all size groups (2-3 mm, 3-5 mm, and 5-10 mm) for the same kind of gravels, are given in Table 5 and Table 6. The fitting results of the soil-gravel mixtures containing the large range of the gravel sizes from 2 to 10 mm were also rather satisfactory since the R2 values of the two typical soil-gravel mixtures were both 0.97. The match effect between measured and calculated data is presented respectively in Fig. 4 and Fig. 5 using the trade-off parameters according to Eq. (5) and Eq. (8). Both Fig. 4 and Fig. 5 show that there was a good match between the measured and simulated data. In addition, the two kinds of the revised WRCs formulas (Eq. (5) and Eq. (8)) had almost the same fitting level. The well fitted results for the soilgravel mixtures with the gravel sizes of 2-10 mm lead to conclude that the effect of gravel sizes was much lower than that of the gravel contents and types. This conclusion gives support to the effectiveness of revised models for large scale field investigation of the soil-gravel mixtures with highly variable gravel size.
Parameter determination with two data sets of representative gravel contents. As shown above, the water characteristic properties of the soil-gravel mixtures can be well determined based on the seven WRCs data sets corresponding to the seven different gravel mass contents (10%, 15%, 25%, 35%, 45%, 55%, and 65%) by fitting the regression surfaces, but the WRCs measurement work for seven gravel mass contents might be time-consuming. It was explored consequently if two measured water retention curves of soil mixtures with minimum and maximum gravel contents (10% and 65%) could represent well all the measured data within the whole range of gravel contents (10%-65%) in determining the parameter of revised water retention model. To examine that, the measured water retention curves of the soil-gravel mixtures with 10% and 65% gravel contents were used to fit Eq. (5) and Eq. (8), respectively, by the nonlinear regression. The regression results are reported in Table 7 showing that high R2 values were achieved. Importantly, the values of the empirical fitting parameters from the two measured curves were rather close to those from the whole measured curves, which indicated that the curved surface obtained based on the two measured curves were acceptable. The comparison between fitting parameters calculated from the two data sets and those from the seven data sets is presented in Fig. 6. The parameters of A and A' for Eq. (5) and Eq. (8) obtained from fitting the two curves were slightly larger than those from fitting the whole measured curves. However, the small deviation did not bring about a considerable difference between the two groups of the curve surface parameters. Therefore the method of determining the parameters of the revised water retention model based on the measured data from two representative gravel contents can be recommended, probably being more practical for the mixtures with embedding gravels in fine soil and for poorly sorted sediments with small size gravels in regional or large scale field investigation.

Conclusion
Based on this investigation it can be seen that the saturated water contents of the soil-gravel mixtures can be directly calculated from the amount and saturated water contents of both fine soil and rock fragments. The experimental results achieved from the different samples containing typical rock fragments with various sizes (2-10 mm) confirmed that the given equation was applicable for computing the water content of the soil-gravel mixtures. Revised formulas for water retention processes of the soil-gravel mixtures were proposed. Revised water retention models were developed by adding one variable (gravimetric gravel content or bulk density of soil-gravel mixtures that reflect the change in rock fragment contents) which had linear or power relation with the effective saturation of soil-gravel mixture to BCfunction forming curved surfaces. Furthermore, the laboratory results indicated that the revised water retention models calibrated either by using the experimental data from the soil-gravel mixtures with a small range of gravel sizes or with a large range of gravel sizes, would be acceptable. However, determination and application of the regression surface model based on the data sets from a few typical different gravel sizes with two representative gravel contents would be more practical. Finally, it seems that the revised water retention models, the procedure of calculation parameters and the specific parameter values of the curved surfaces obtained in this study may be applied to different kinds of soil-gravel mixtures, because the whole analyses in this study were based on the measured water retention in soil-gravel mixtures containing two different types of rock fragments. However, their applicability needs to be validated through the water retention data from soils containing different rock fragments.