Evaluation of Hydraulic Parameters Obtained by Different Measurement Methods for Heterogeneous Gravel Soil

Knowledge of soil hydraulic parameters for the van Genuchten function is important to characterize soil water movement for watershed management. Accurate and rapid prediction of soil water flow in heterogeneous gravel soil has become a hot topic in recent years. However, it is difficult to precisely estimate hydraulic parameters in a heterogeneous soil with rock fragments. In this study, the HYDRUS-2D numerical model was used to evaluate hydraulic parameters for heterogeneous gravel soil that was irregularly embedded with rock fragments in a grape production base. The centrifugal method (CM), tensiometer method (TM) and inverse solution method (ISM) were compared for various parameters in the van Genuchten function. The soil core method (SCM), disc infiltration method (DIM) and inverse solution method (ISM) were also investigated for measuring saturated hydraulic conductivity. Simulation with the DIM approach revealed a problem of overestimating soil water infiltration whereas simulation with the SCM approach revealed a problem of underestimating water movement as compared to actual field observation. The ISM approach produced the best simulation result even though this approach slightly overestimated soil moisture by ignoring the impact of rock fragments. This study provides useful information on the overall evaluation of soil hydraulic parameters attained with different measurement methods for simulating soil water movement and distribution in heterogeneous gravel soil.


INTRODUCTION
The Richards equation, a foundation of the physicsbased models, was derived from the principle of the conservation of mass and Darcy's law (Lei et al. 1988) which can be extended to more complex conditions (Brunone et al. 2003;Pachepsky et al. 2003;Elmaloglou and Diamantopoulos 2008).Numerical models of the Richards soilwater kinetic equation are an efficient tool that can be used to describe variable saturated soil water flow such as rainfall infiltration and runoff, subsurface recharge, migration of nutrients, solution transport, and design and monitoring of irrigation and drainage systems (Bagarello et al. 2005).Unsaturated hydraulic conductivity K(θ), soil water diffusivity D(θ) and specific water capacity C(θ) as functions of water content θ are necessarily basal parameters which help understand soil water movement laws and quantitatively analyze soil water movement with mathematical modeling method.These three parameters are almost impossible to measure directly.Hence, soil water characteristic curves are proposed to indirectly solve this problem based on K(θ) = D(θ) × C(θ).Among the models fitting characteristic curves, van Genuchten model, describing the moisture absorption processes, was considered to be the most popular to fit soil water characteristics.Therefore, parameters of van Genuchten function, α, n, θ r and θ s , as well as saturated hydraulic conductivity (Ks) needed to calculate K(θ), become a key for solving the Richards equation.
Soil water characteristic curves can be produced using different methods (Lei et al. 1988).Each has its pros and cons.(1) The tensiometer method can obtain successive measurement of soil water potential and be easily applied in the field; however, this method has a limited measurement range and the results can be affected by soil spatial variability.(2) The pressure membrane method has a wide measurement range of soil potential, but is time consuming.(3) The tension table method is easy to operate and is low-cost, but is limited to a narrow measurement range.(4) The centrifugal method can save time, but error is easily incorporated during the centrifugal process through heating and compacting soil sample.The compaction effect is more profound for a loamy soil than a sandy soil (Reeve et al. 1980).
Different techniques have been applied for measuring K.They include the constant head-cutting ring method (Klute and DirKsen 1986), the Guelph permeameter method (Reynolds and Elrick 1985;Elrick and Reynolds 1992;Xu and Mermoud 2003), the single-ring infiltrometer method (van Es et al. 1999;Bagarello and Sgroi 2004), the doublering infiltrometer method (Starr 1990), and the inversed-auger-hole method (Messing and Jarvis 1990).In recent years, tension infiltromenters have been widely used for measuring the near-saturated K (Perroux and White 1988) in the field to investigate the spatial and temporal variability of soil hydraulic properties (Das Gupta et al. 2006;Hu et al. 2009).This method was considered to be water saving, repeatable and stable in space and time (Ventrella et al. 2005;Bagarello and Sgroi 2007).
Owing to the soil forming process and human activities, certain rock fragments or gravel layers are usually embedded in soil.The rock fragments can be ignored when their content is not high and has little influence on soil water movement.However, a large portion of rock fragments should be taken into account (EriKsson 1996).In general, the content of rock fragments can be investigated with four methods, including the cutting ring sampling method (Flint and Childs 1984), the digging method (Muller and Hamilton 1992), the probe method (Viro 1952), and the ray/radio wave method (Fleming et al. 1993;Rey et al. 2006).All of these methods are time and energy consuming and difficult to obtain satisfactory results in a heterogeneous medium.In addition, although the impact of size, content and distribution of rock fragments on soil water movement have been studied by many scholars, conflicting trends often exist in various reports revealing positive correlations, negative correlations, no correlations or a "turning point" between the size/content and infiltration ability (conduction) (Mehuys et al. 1975;Bouwer and Rice 1984;Ravina and Magier 1984;Brakensiek et al. 1986;Sauer and Logsdon 2002;Khaleel and Heller 2003).Consequently, it has been a challenge to understand the impact of rock fragments on soil hydraulic character, especially at a heterogeneous research site.Predicting soil water flow and distribution accurately and rapidly in heterogeneous gravel soil becomes a hotspot in recent years (Ma et al. 2009(Ma et al. , 2010)).
Numerical simulations are efficient tools to study soil water movement in various complex conditions (Lubana and Narda 2001;Gärdenäs et al. 2005).Researchers have been attempting to study the effect of rock fragments via the "black-box" model in complex conditions.Many simulation software packages have been developed in recent years, including HYDRUS-2D (Šimůnek et al. 1999), a well-known Windows-based computer software package used for simulation on water, heat and solute movement in a 2D variably saturated porous media (Cote et al. 2003;Skaggs et al. 2004).Based on a finite element method, the HYDRUS-2D was developed to solve the Richards's equation, and was proved to be in good agreement for observed and simulated data (Ben-Gal et al. 2004;Skaggs et al. 2004;Provenzano 2007).
The objective of this study was to evaluate soil hydraulic parameters by using the HYDRUS-2D to compare various methods.The parameters α, n, θ r and θ s were measured using the centrifugal method and the tensiometer method, and the parameter Ks was measured using the soil core method and the disc infiltration method, respectively.

Study Site Description
The field research was conducted in a grape cultivated base, located in the eastern part of Shanshan county (42°54'N, 90°30'E, 416 m elevation), Xinjiang Uygur Autonomous Region, China.The climate was characterized as arid, with mean annual precipitation of 25.3 mm and potential evaporation of 2751 mm.The mean annual temperature was 14°C.Ground water table was 70 m depth below the surface.The soil had a deeper layer of Gobi gravel soil originating from Gobi gravel, and was covered by a layer of exotic sandy loam soil heterogeneously mixed with rock fragments in the top 50 cm.Soil particle sizes collected at the depths of 0 -5 and 60 -65 cm were determined using a MasterSizer 2000 laser particle size analyzer (Malvern, UK).These two soil layers had the texture of sandy loam soil (SI standard) (Table 1).The proportion of rock fragment in each soil layer was measured using the gravimetric method.Because the rock fragments were heterogeneously distributed in profiles, the rock fragment contents shown in Table 2 were the average values of the 0 -50 cm upper layer and the 50 -120 cm layer (Table 2).Thompsons Seedless grapes were planted in 1980.The cultivation inter-row space was 3.5 m, which included a 0.9 -1.2 m wide furrow part and a 2.3 -2.6 m wide ridge part within each row.The surface of cultivated furrow was treated as a reference level of 0 m.The ridge was about 40 -50 cm higher than the furrow (Fig. 1).The study site was drip irrigated during the study period.Three drip irrigation tapes were laid along the cultivated furrow: (1) at the bottom of furrow, (2) on the ridge and near the grapevine root, and (3) on the ridge, about 50 cm apart from the second as shown in Fig. 1.Water movement below each individual dripper is usually considered as a point source infiltration process (Chu 1994).Hence, water transfer under drip irrigation can be assumed to be a 2D infiltration process with multi-point source provided that soil was horizontally uniform along the furrow direction.The drip irrigation flow for each emitter was 3.3 L h -1 with 30 cm spacing, and the irrigation quota through all three tapes was 525 m 3 ha -1 in each irrigation event.

Measurement of α, n, θ r and θ s , Hydraulic Parameters of van Genuchten Function
The van Genuchten function was selected as the soil water characteristic curve (van Genuchten et al. 1980) as follows where θ s is the saturated volumetric water content (cm 3 cm -3 ); θ r is the residual volumetric water content (cm 3 cm -3 ), h is the soil water pressure head (cm), α is a parameter related to air-entry suction (cm -1 ), m and n are coefficients related to shape coefficient with a relation of m = 1 -1/n and n > 1, K(h) is hydraulic conductivity (cm min -1 ), S e is the effective water content, l is a pore-connectivity parameter which is often assumed to be about 0.5, an average value of many soils (Mualem 1976).
In this study, the hydraulic parameters of van Genuchten function, α, n, θ r and θ s , were obtained using two methods: the soil water characteristic function method (SW-CFM) and the inverse solution method (ISM).
For SWCFM, soil water characteristic function was obtained by a centrifugal method (CM) and tensiometer method (TM).For CM, 9 undisturbed soil cores were taken from the depths of 0 -5 and 60 -65 cm using a cutting ring (5 cm in height 200 cm 3 in volume) in 2008.Soil water potentials were measured using a thermostatic high-speed centrifuge (HITACHI, Japan).For TM, soil water potentials were manually measured using 9 tensiometers (2500S, France).For the top soil layer, 9 tensiometers were buried at a shallow depth of 20 cm in 2009.For the deep soil layer, 9 tensiometers were buried at a depth of 60 -65 cm in 2008.Soil moisture around the tensiometers was synchronously measured using Time Domain Reflectometry (TDR) method (TRIME, German).The hydraulic parameters were obtained by fitting the soil water retention curves using RETC V6.02 software (van Genuchten et al. 1991).
For ISM, the initial soil water distribution in a vertical cross section along the whole 3.5 m inter-row horizontal width was measured using 5 neutron tubes at the depth of 20 -120 cm with 20 cm intervals (503BD, China) (see Fig. 1) prior to the irrigation event.Soil water in the cultivated furrow surface and cultivated ridge was measured using the gravimetric method.After an irrigation event, a 3.5 m wide and 1.5 m deep vertical cross-section along the whole inter-row horizontal width was collected to measure soil water content by using the gravimetric method.These experiments were conducted from 2009 -2010.The Table 1.Soil particle size distribution at 0 -5 and 60 -65 cm layers for the study site.
Table 2.The gravimetric percentage of gravel content in each soil layer.
a The size less than and equal to 2 mm is defined as soil particle and the size greater than 2 mm is defined as stone according to soil particle definition.hydraulic parameters were obtained by the inverse solution of soil water distribution using the HYDRUS-2D V2.1 model.

Measurement of Ks
In this study, Ks values were measured using three methods, the soil core method (SCM) (Klute and DirKsen 1986), disc infiltration method (DIM) (Logsdon and Jaynes 1993), and inverse solution method (ISM).
For SCM, Ks was obtained according to Darcy's law.The undisturbed soil core was sampled using a cutting ring as described in 2.2.1 and a Marriote bottle to keep a stable hydraulic head of 4.0 cm.The discharge water volume Q (cm 3 ) was recorded every 30 min.Ks was calculated using stabilized Q that was an average of the last three values until the differences were within the 5% range for individual 30 min intervals.Then, Ks in cm min -1 was calculated using Eq.(3): where A is the cross sectional area of the soil core (cm 2 ), T is the time (min), H is the hydraulic head (cm), and L is the length (cm) of soil core.
For DIM, the configuration of disc infiltrometer was similar to that described by Ankeny et al. (1988) including an infiltrometer base with a radius of 7.5 cm and a reservoir tube with a radius of 1.7 cm.The disc infiltrometer was used to determine infiltration under pressure heads of -15, -6, -3 and 0 cm in an ascending sequence.Before each measurement, the flat and cleared soil surface was prepared using a knife; and, a fine layer of sand (about 1 mm) was evenly placed on the surface to ensure good contact between the infiltrometer base and the soil.For each infiltration mea-surement, the cumulative infiltration was recorded at 5 min intervals until steady infiltration occurred.After the first infiltration experiment finished, the following measurements with different pressure heads were repeated at the same location to reduce the spatial variability.The multi-tensions with a non-linear regression method as reported by Logsdon and Jaynes was believed to produce more stable results than other methods (Logsdon and Jaynes 1993) and hence was chosen in our study.The fitting equation is expressed as ^h is the steady infiltration rate (cm 3 min -1 ) under pressure head of h f (cm), R is the radius of the disc infiltrometer (cm), α 0 is the Gardner constant that characterizes soil pore size distribution (cm -1 ).
For the disc infiltrometer where r is the radius of the reservoir tube (cm), H 0 (cm) is the height of water drop in reservoir tube at time interval of t (min).Consequently, Ks can be calculated by combining Eqs.(4) and (5).
For ISM, Ks was obtained the same way as other hydraulic parameters using the HYDRUS-2D model as described in 2.2.1.

Model Description
The HYDRUS-2D model is based on the two-dimensional Richards equation.The program can be used to nu- merically solve the Richards equation for variably saturated water movement for a given soil, and has been used to simulate soil water movement by many researchers (Cote et al. 2003;Skaggs et al. 2004).The root uptake and evapotranspiration were not taken into account in this study.The water movement equation is described as where θ is the volumetric water content, h is the soil water pressure head, t is the time, x is the horizontal coordinate with the origin at ridge-furrow interface beyond grape vine (X 0 , see Fig. 1), z is the vertical coordinate with the origin at the soil surface (positive upward) (Z 0 , see Fig. 1), K( h ) is the unsaturated hydraulic conductivity.
The simulation zone (Fig. 2) with 3.5 m width, 2.0 m height for the furrow part and 2.5 m for the ridge part (a positive value represents above reference surface and a negative value stands for below the surface).The simulation zone was set for two soil types (above 50 cm and below 50 cm depth), then the finite element meshes of simulation zone were constructed by dividing the flow region into irregular triangular elements.
For the present study, both the initial condition and the upper boundary condition were , , , h x z h x z where , h x z i ^h is the initial soil water pressure head, and h 0 is the soil water potential at soil surface.These values were set according to observed results.The upper boundary condition of infiltration points under the dripper was set as constant flux.The free drainage was considered as a lower boundary condition.
In our study, the simulation conditions were the same for different hydraulic parameters.The simulation time was set to 900 min in accordance with field sampling experiment.The HYDRUS-2D program was performed with 5 hydraulic parameters of α, n, θ r , θ s and Ks by 5 treatments: (1): CM-SCM; (2): CM-DIM; (3): TM-SCM; (4) TM-DIM; (5): ISM.An additional filed irrigation experiment was conducted to evaluate the accuracy of simulation results with the observed data obtained with different methods.

Method of Analysis
Analysis of variance (ANOVA) and Post Hoc Tests (using a least significant difference test when equal variance occurred and Tamhane's T2 when equal variance did not occur) were performed to analyze the difference of hydraulic parameters obtained with different methods.Paired T test was employed to compare the hydraulic parameters of the layers of 0 -50 cm and below 50 cm.The coefficient of variation (CV) was used to evaluate the magnitude of soil water variability.Heterogeneity was considered weak if CV < 10%, moderate if 10% ≤ CV ≤ 100%, and strong if CV > 100%.These analyses were performed using SPSS V13.0 software.
The root mean square error (RMSE) was calculated to provide a quantitative comparison of the goodness-of-fit for the measured and simulated data, which can be expressed as where N is the total number of observations, O i and P i are the observed and predicted values, respectively.

Hydraulic Parameters of van Genuchten Function
Table 3 shows the typical values of the hydraulic parameters for various soils given by the manual for the HY-DRUS-2D program (Šimůnek et al. 1999), where α takes the mean values of 0.008, 0.016, 0.036 and 0.145 cm -1 , for clay, silt, loam and sand soil, respectively; n takes the mean values of 1.09, 1.37, 1.56 and 2.68, for clay, silt, loam and sand soil, respectively; θ r takes the mean values of 0.068, 0.034, 0.078 and 0.045 cm 3 cm -3 , for clay, silt, loam and sand soil, respectively; θ s takes the mean values of 0.38, 0.46, 0.43 and 0.43 cm 3 cm -3 for clay, silt, loam and sand soil, respectively.Table 4 summarizes the statistical parameters  In this study, Table 4 shows that, in the top 0 -50 cm, the average values of α for CM > ISM > TM.Moderate variation for CM (CV = 14.0%),TM (CV = 11.7%) and weak variation for ISM (CV = 3.6%) were observed.The average values of n for TM > ISM > CM.Moderate variation for CM (CV = 15.1%) and TM (CV = 11.5%) and weak variation for ISM (CV = 4.1%) were observed.The average values of θ r for TM > CM > ISM.Moderate variation for CM (CV = 12.4%) and weak variation for (CV = 9.9%) TM and ISM (CV = 3.6%) were observed.The average values of θ s for TM > ISM > CM.All the three measurement methods showed weak variation (with CV of 5.5%, 9.8% and 3.3% for CM, TM and ISM, respectively).
Below 50 cm, the average values of α for ISM > TM > CM.Moderate variation for CM (CV = 16.1%),TM (CV = 22.0%) and weak variation for ISM (CV = 3.6%) were observed.The average values of n for TM > ISM > CM.All the three measurement methods showed weak variation (with CV of 9.1%, 4.0% and 4.4% for CM, TM and ISM, respectively).The average values of θ r for TM > CM > ISM.Moderate variation for CM (CV = 11.5%),TM (CV = 22.9%) and weak variation for ISM (CV = 1.7%) were observed.The average values of θ s for CM > TM > ISM.All the three measurement methods showed weak variation (with CV of 5.9%, 5.6% and 2.9% for CM, TM and ISM, respectively).
One-way ANOVA was performed to explore possible impacts of measurement methods on the hydraulic parameters in van Genuchten function.Figures 3 and 4 show the averaged parameters α, n, θ r and θ s for various methods.Visually, some differences in α, n, θ r and θ s existed among various methods.ANOVA analysis showed statistically significant difference for α, n, θ r (P < 0.001) in 0 -50 cm depth and for n and θ r (P < 0.002) below the 50 cm depth.Results of Post Hoc Tests showed significant difference between TM and CM as well as between TM and ISM for α, n, θ r (P < 0.001) in the 0 -50 cm depth.Significant difference between TM and CM as well as between TM and ISM was only observed for θ r (P < 0.001) and between TM and CM for n (P < 0.001) below the depth of 50 cm.Moreover, differences might exist but not reach a significant level (P < 0.01) between TM and CM for α (P = 0.042) and between TM and ISM for n (P = 0.020) below the depth of 50 cm.
All three measurement sets for each soil layer were pooled together for each van Genuchten parameter, and then Paired T Test was performed to explore the impact of soils on van Genuchten parameters between the 0 -50 cm and below 50 cm layers.Table 4 also shows statistically results of van Genuchten function parameters of the two soil layers.Results showed significant differences for α, θ r , θ s (P < 0.001) between the two layers.ANOVA analysis showed that a significant difference in α values for TM, ISM (P < 0.001), and in θ r values for TM (P < 0.001) between the two soil layers.Differences in n for CM (P = 0.0305), θ r and θ s for ISM (P = 0.0317 and 0.0350) were also observed between the two soil layers.

Soil Saturated Hydraulic Conductivity
Values of soil saturated hydraulic conductivity (Ks) were obtained using the soil core method (SCM), disc infiltration method (DIM), and inverse solution method (ISM) in top 0 -50 cm and below 50 cm as shown in Table 4.
ANOVA analysis showed significant differences for Ks that were obtained using different methods (P < 0.001) in both the top 0 -50 cm and below 50 cm soil layers (Fig. 5).Results of Post Hoc Tests showed differences existed among all the three methods in the two layers, but the significant differences (P = 0.001) were observed between DIM and SCM in the top layer and between DIM and SCM as well as between DIM and ISM in the lower layer (Fig. 5).Comparing Ks values of soil layers between the 0 -50 cm and below 50 cm layers, Paired T Test indicated significant differences for Ks between the two layers (P < 0.001) (see Table 4).ANOVA analysis showed only significant difference for DIM and ISM (P < 0.001) between the two layers.

Analysis of Simulation Results
Numerical simulations of the water flow were performed for obtaining hydraulic parameters that were obtained using different methods.Simulation results of different methods were compared with observed data as shown in Fig. 6.Visually, CM-DIM and TM-DIM overestimated the soil water infiltration process, while CM-SCM and TM-SCM underestimated the infiltration process.Methods of DIM with higher Ks values have higher infiltration speed than the SCM methods, indicating that impact of Ks was more significant than other hydraulic parameters.Compared with other methods, ISM showed better simulating results, but it overestimates soil water content especially for the layer below the depth of 50 cm.This may be caused by the high content of rock fragment in the whole soil profile, especially for the deeper layer.Without considering the impact of rock fragments in the HYDRUS-2D simulation, the vertical soil cross section was more uniform and assumed to be "sandy soil," which could hold more water than rock fragments in the simulation process.Therefore, higher simulation values of soil water content were observed in the deeper layer.
Values of the simulation were compared to the observed data at the same location.Analytic results showed that (Table 5) the best simulation result with a higher correlation coefficient of 0.8369 for ISM, and 0.4068 for RESM.The relative error was less than 15%.However, the CM- SCM model produced the worst simulation result with a relative error of 50 % and RMSE of 1.1336.

DISCUSSION AND CONCLUSIONS
It is difficult to accurately and directly investigate soil hydraulic parameters in heterogeneous gravel soil.The HYDRUS-2D numerical simulation model was chosen in this study to evaluate the accuracy of hydraulic parameters that were indirectly measured using several methods.
The hydraulic parameters as obtained with the CM are significantly different from that obtained with the TM, so were from the CM and the ISM for α, n, and θ r .The Ks parameter was significantly changed from the DIM to the SCM and from the DIM to the ISM.The hydraulic parameter of α by TM, ISM obtained in 0 -50 cm soil layer are different from that obtained below 50 cm, and were similar from the 0 -50 cm and below 50 cm for θ r and Ks.The hydraulic parameters of n by CM obtained in 0 -50 cm soil layer are different from that obtained below 50 cm and the parameter of θ s by DIM and ISM obtained are also different between the two layers.The simulation results of hydraulic parameters obtained using various methods showed that the ISM had a relatively smaller error (15%) compared to other methods.The simulation results of parameters with the DIM were overestimated soil water infiltration process and results with the SCM were underestimated water mobility compared to the actual field results.
These analyses imply that the soil samples collected by the cutting ring sampler for the CM and SCM did not accurately reflect the impact of gravel on soil hydraulic parameters because the soil core cylinder was almost impossible to contain gravel due to the cutting ring size.Therefore, both CM and SCM yielded smaller Ks value than other methods and did not show difference between the two soil layers.In addition, the significant difference of Ks values as obtained with DIM, ISM and SCM was probably caused by two factors.First, the infiltration process measured by a disc infiltrometer could only reach a near-steady state under near-saturated condition rather than saturated condition.Second, the existence of rock fragments could easily lead to macro-pore flow and result in higher Ks values.The ISM yielded the best simulation results among all methods used in this study because the effect of gravel on the whole profile by a "black-box" scheme was taken into account.However, this method discounted the effect of rock fragments on soil water distribution especially in the deeper soil layer, hence the simulated soil water values were higher than the values actually observed in the deeper soil layer.Ankeny, M. D., T. C. Kaspar, and R. Horton, 1988

Fig. 1 .
Fig. 1.The schematic diagram of the study site.

Fig. 2 .
Fig. 2. The schematic diagram of the simulation zone.

Fig. 3 .
Fig. 3. Values of soil hydraulic parameters of the van Genuchten function above the depth of 50 cm by different methods.The characters above the columns indicate the significance.Columns with the same character do not show significant differences for P < 0.05.

Fig. 5 .
Fig. 5. Values of saturated soil hydraulic conductivities by different methods.The characters above the columns indicate the significance.Columns with the same character do not show significant differences for P < 0.05.

Fig. 6 .
Fig. 6.The isoline graph of filed observed data and simulated results of soil water distribution in the vertical cross section.

Table 3 .
The typical values of soil hydraulic parameters for the van Genuchten function.

Table 4 .
Statistical summary of hydraulic parameters of the van Genuchten function obtained by different methods.

Table 5 .
Statistical summary of soil water content for observed data verse different methods.