Use of Field and Laboratory Methods for Estimating Unsaturated Hydraulic Properties under Different Land Uses

Adequate water management is required to improve the efficiency and sustainability of agricultural systems when water is scarce or over-abundant, especially in the case of land use changes. In order to quantify, to predict and eventually to control water and solute transport into soil, soil hydraulic properties need to be determined precisely. As their determination is often tedious, expensive and time-consuming, many alternative field and laboratory techniques are now available. The aim of this study was to determine unsaturated soil hydraulic properties under different land uses and to compare the results obtained with different measurement methods (Beerkan, disc infiltrometer, evaporation , pedotransfer function). The study has been realized on a tropical sandy soil in a mini-watershed in northeastern Thailand. The experimental plots were positioned in a rubber tree plantation in different positions along a slope, in ruzi grass pasture and in an original forest site. Non-parametric statistics demonstrated that van Genuchten unsaturated soil parameters (K s , α and n) were significantly different according to the measurement methods employed, whereas the land use was not a significant discriminating factor when all methods were considered together. However, within each method, parameters n and α were statistically different according to the sites. These parameters were used with Hy-drus1D for a 1-year simulation and computed pressure head did not show noticeable differences for the various sets of parameters, highlighting the fact that for modeling, any of these measurement methods could be employed. The choice of the measurement method would therefore be motivated by the simplicity, robustness and its low cost.


Introduction
The rubber tree (Hevea brasiliensis) has become a crop of high economic interest in the northeast of Thailand since the rise in price of natural rubber on the international market and the policy of the Thai government to extend rubber tree plantation.Despite climatic and edaphic conditions being very different to those of the original growing region (south of Thailand), rubber tree plantation extended by about 170 000 ha between 2004 and 2006 (Rantala, 2006).The introduction of the rubber tree in the northeast of Thailand may also contribute to important land use changes affecting soil and water resources.In this area, the average annual rainfall (1.1-1.2 m) does not completely meet the minimal requirements for rubber trees (1.3 m), and therefore it is necessary to design a wise water management to be able to ensure sustainable rubber tree farming.In order to achieve this practical goal, numerical modeling of water flow in soils is a valu-Published by Copernicus Publications on behalf of the European Geosciences Union.
able tool for quantifying precisely the water balance and the different mechanisms involved (Soares and Almeida, 2001;Antonino et al., 2004;Anuraga et al., 2006) and therefore for trying to forecast the evolution of the system, as the new land use may be leading to changes in the general water balance.In order to achieve accurate quantification of these changes with modeling, it is therefore necessary to estimate precisely the hydrodynamical properties of the vadose zone.
Several laboratory methods are commonly used to determine the hydraulic properties of soil, like the pressure plate method (Klute and Dirksen, 1986;Bittelli and Flury, 2009) to determine the retention curve, and hydraulic conductivity measurements like hot air, crust, one-step outflow, and sorptivity-based methods (Stolte et al., 1994).Some in situ methods have mainly been developed to measure hydraulic conductivity in the field, like the tension disc infiltrometer (Perroux and White, 1988;Akeny et al., 1991;Angulo-Jaramillo et al., 2000).In any of these cases, the determination of the hydraulic parameters of unsaturated soil is always tedious and time-consuming.Therefore, the quest for alternative methods of determination is a research issue.For example, the development of pedotransfer functions (PTFs), which provide relationships between soil hydraulic parameters and more easily measurable properties such as particle size distribution (PSD), has become successful amongst soil scientists and hydrologists (Bouma, 1989;Pachepsky et al., 2006).They are mainly based on the textural properties of soil, but show some inadequacy in describing structural aspects of soil.Consequently, specific methods taking into account structural properties have been developed to tackle this problem.For example, the Beerkan method (Haverkamp et al., 1996;Braud et al., 2005) is a field and laboratory measurement method that aims to determine retention curves and hydraulic conductivity curves by particle size distribution and single ring infiltration data, in order to describe the contribution of the texture and the structure to hydraulic properties.But, the main interest of this method is to provide an efficient and low cost estimation of soil hydraulic properties.These basic data are processed with the Beerkan Estimation of Soil Transfer parameters (BEST) algorithm according to this theory (Lassabatere et al., 2006) to provide acceptable estimations of hydraulic parameters to a complete characterization of hydraulic characteristic curves.We also compared it to a laboratory evaporation method inspired by the method proposed by Wind (1968), which is based on an evaporation measurement to estimate soil water retention and hydraulic conductivity.Wind's method has become very popular and is constantly improved by the introduction of inverse modeling (Tamari et al., 1993;Simunek et al., 1998a;Richard et al., 2001), and has been adapted for infiltration procedures (Bruckler et al., 2002).Though hydraulic conductivity is well predicted for the dry part of the curve, Wind's method is generally inadequate near saturation for hydraulic gradients that are too low (Tamari et al., 1993;Wendroth et al., 1993;Richard et al., 2001).In this study, we considered only part of the evaporation technique, in order to determine experimentally the retention curve.
Land use changes are known to modify soil properties and especially the unsaturated soil hydraulic characteristics (Zimmermann et al., 2006;Price et al., 2010), as the soil structure might be affected by different tillage methods, different root densities and sizes, and different biological activities of macro-fauna.Therefore, the hydraulic conductivity of the top soil will be affected and the infiltration capacity will be modified.In the context of tropical rain patterns with heavy rainfall followed by long dry periods, with these shallow soils, the infiltration capacity is an important factor in the soil water balance.The efficiency of the rainfall highly depends on the infiltration capacity.
Different land uses such as rubber tree plantation, pasture or natural forest in the same soil series are expected to show very different soil properties.The purpose of this study is to verify this hypothesis and to compare the performance of different measurement techniques for estimating the soilhydraulic properties of top soil in a typical mini-watershed of northeastern Thailand, in order to model soil water flow accurately.

Site and soil description
The study area is set up in a small watershed located near Ban Noon Tun, Phra Yuen District, Khon Kaen Province, Thailand (16 • 19 43.90 N, 102 • 45 07.91 E) at about 19 km southwest of the city of Khon Kaen.This area is subject to a tropical savanna climate, with annual rainfall of 1.3 and 1.96 m in 2007 and 2008, respectively, and an average annual temperature of 29 • C. The site is located in the typically undulating landscape of northeastern Thailand, with an elevation ranging from 165 to 181 m above medium sea level and an average slope of 3.5 % (Fig. 1).The surface area of the mini-watershed is estimated to be 200 m × 104 m, where most of the surface is covered with rubber trees (Hevea brasiliensis) and ruzi grass (Brachiaria ruziziensis) cultivated for seeds.The area dedicated to rubber tree plantation is approximately 120 m × 104 m, and previously the land was planted with jute (Corchorus capsularis), cassava (Manihot esculenta) and ruzi grass.Following the governmental policy and the bright economical perspectives for the farmers, the soil use changed to a rubber tree plantation in 2004.
The soil is generally shallow; in the upper part of the toposequence, the soil thickness is about 0.70 m, 0.90 m midslope, and 2.5 m downslope.The soil texture is classified as sand to sandy clay loam, with a bulk density ranging from 1.38 to 1.54 × 10 3 kg m −3 .According to the USDA classification system, the soil fit into three subgroups, i.  layer with low clay and organic matter contents, overlying a 0.20 to 0.50 m thick grey clayey layer colored locally with red weathered material from the bedrock (Wiriyakitnateekul et al., 2009).The bedrock, which is also probably the parental material, is a fine red sandstone or siltstone, containing clay minerals and feldspars (orthoclase), thoroughly weathered in its upper part, and densely fractured.The physical and chemical properties of the soil measured in the laboratory (Table 1) show that the texture of the top soil is sandy for rubber trees (upslope and midslope) and ruzi grass and loamy sand for rubber tree (downslope) and forest.Moreover, they show slightly acidic soils with very low organic matter content, though slightly higher in the forest (Table 1).The first situation considered for measuring soil hydraulic properties was a toposequence of 300 m along a gentle slope (3 %) in the rubber tree plantation (i.e., RT site, Fig. 1) where three positions were selected: upslope, midslope, and downslope.Next to the rubber tree plantations, important extensions of ruzi grass pasture are still present in this small wa-tershed, which are ploughed once a year, during the rainy season, just before sowing.The second measurement site was located in a ruzi grass pasture plot close to the rubber tree plantation (i.e., RG site, Fig. 1).The third measurement site was located in the upper part of the watershed, where an original Dipterocarpus forest has been conserved (i.e., F site, Fig. 1).Experimental soil water flow measurement devices, namely, tensiometers with pressure sensors (SKT 850T, SDEC, Reignac-sur-Indre, France) and soil moisture sensors (Enviroscan probe, Sentek, Stepney SE, Australia), were installed at the different sites (RT upslope, midslope, downslope, RG and F) and monitored every 2 h continuously from 2007 to 2009.Tensiometric data were recorded at 0.10, 0.25, 0.40, 0.60 and 1.10 m and soil moisture was measured at 0.10, 0.30, 0.50, 0.70, 0.90, 1.10 and 1.40 m.Meteorological data (rainfall, temperature, air humidity, wind speed, solar radiation) were recorded continuously at a single site in the small watershed (Fig. 1).At each experimental site (RT upslope, midslope, downslope, RG and F), Beerkan infiltration experiments were performed in the field during the dry season, from November 2007 to February 2008, and soil cylinders were collected for laboratory evaporation method measurements.The soil volumetric water content at 0-0.10 m depth, measured when the Beerkan infiltration experiments were performed, was 0.01 (m 3 m −3 ) for RT upslope and F and 0.02 for the other sites.

Estimation of soil unsaturated hydraulic properties
Modeling and quantifying water flow in the vadose zone is mostly described by Richards' equation (Richards, 1931), and in order to solve this equation, one is required to deter-mine (i) the soil water retention function, relating the pressure head to the soil water content h(θ ), and (ii) the hydraulic conductivity based on the Brooks and Corey (1964) model (Eq. 3), or the van Genuchten relationship: These equations are defined by several parameters where h is the pressure head (L), S e is the effective saturation (-), θ (L 3 L −3 ) is the volumetric water content, θ r and θ s are the residual and saturated water content, m (-) and n (-) are shape parameters, and k is an integer usually chosen to be 1 (Mualem, 1976) or 2 (Burdine, 1953).α (L −1 ) is the scalefitting parameter, K s (L T −1 ) is the saturated hydraulic conductivity, and l (-) and η (-) are other shape parameters, with l usually considered to be 1/2.

The Beerkan method
The experimental procedure of this method consisted of two distinct processes, i.e., the soil particle size distribution analysis and a single ring infiltration test.According to the theory developed in Braud et al. (2005) and Lassabatere et al. (2006), 3-D axisymmetric infiltration was realized with a zero pressure head at the soil surface.It was achieved with a pulse flux of small volumes of water, and the infiltration tests were performed in situ with a 0.102 m diameter cylinder.The choice of the cylinder diameter for the Beerkan method was motivated by the experience from the literature (Braud et al., 2005) and the availability of the material.As the soil in this area is sandy and has little or no structure, it was assumed that the diameter of the infiltration cylinder would not impact the representativity of measurement.
The infiltration cylinder was driven firmly into the soil for about 0.01 to 0.02 m in order to prevent lateral loss of water from the infiltration cylinder (Fig. 2b).A series of constant volume of water (120 × 10 −6 m 3 ) was poured into the cylinder and the time was recorded after each volume of water had been completely infiltrated.Each volume represented a maximum of 0.015 m of pounding pressure at the soil surface when the cylinder was freshly refilled.Haverkamp et al. (1998) showed that small variations in pounding pressure did not influence significantly the infiltration rate, and the surface pressure head in the present case could be considered to be nil.The soil surface micro-topographic irregularities in a cylinder of 0.102 m in diameter are therefore negligible as regards the surface pressure.The time reading for each added volume was recorded precisely when the last free water puddle vanished completely from the soil surface.A new volume was added immediately to ensure continuous water supply.This procedure was continued until the steady-state infiltration rate was reached.Undisturbed and disturbed soil samples were collected at the test spot to determine dry bulk density, PSD by the sedimentation method, and initial (θ 0 ) and final (θ s ) volumetric water content.At each site, when the results were homogeneous, a minimum of three infiltration tests have been conducted.When the variability of infiltration was higher, the number of infiltration tests was increased to seven in some cases.The results were analyzed with the BEST algorithm (Lassabatere et al., 2006) in order to obtain both, the water retention curve (Eq. 1) with the Burdine condition, k = 2 (Eq.2), and the hydraulic conductivity (Eq.3).The details of this calculation procedure can be found in Lassabatere et al. (2006).However, in order to compare the results with other methods and to be able to use these parameters in numerical models like Hydrus1D (Simunek et al., 2005), the curves obtained with BEST were then adjusted to the van Genuchten with Mualem conditions (k = 1 in Eq. 2).The retention curve was plotted according to the Burdine parameters and the equation with the Mualem condition (m = 1 − 1/n) was adjusted with a Marquardt procedure to fit the new parameters.

Disc infiltrometer
This common field technique to measure saturated hydraulic conductivity was used at the different sites with the SW080 B device (SDEC, Reignac-sur-Indre, France).The principle of the tension disc infiltrometer is based on maintaining the water in the apparatus under a controlled tension, so that only pores with lower matric potential can be filled.With this technique, the bio-macropores, cracks and other structures promoting preferential flow can be ignored, to measure hydraulic conductivity strictly in the soil matrix.A tension disc infiltrometer consists of a water reservoir, a Mariotte bubbling tower, and a contact disc of diameter 0.20 m covered with a microporous nylon membrane (with a pore diameter of 20 × 10 −6 m).The water reservoir outlet was connected to the center of the disc with a flexible plastic tube (Fig. 2a).The tension applied to water in the device is controlled in the Mariotte tube; the depth of the air inlet under the water surface controls the minimal tension necessary to draw water out of the infiltrometer throughout the disc.Infiltration is realized until it reaches the constant infiltration rate successively with decreasing tension.The soil surface has to be cleared of vegetation and leveled to ensure perfect contact with the infiltration disc.Usually, the soil surface was slightly covered with clean fine sand to get a smooth horizontal surface and to provide a good contact between the base of the disc and the soil below.The relative position of the infiltration disc with the water reservoir is not necessarily constant.Therefore, it is important to measure it in order to calculate the actual water potential head controlled with the immersed tube of the Mariotte device.In order to calculate saturated hydraulic conductivity with the multi-potential method (Perroux and White, 1988;Smettem and Clothier, 1989), the infiltration measurements have been realized for two different tension values and interpreted with Wooding's method (Wooding, 1968;Akeny et al., 1991;Angulo-Jaramillo et al., 2000).The tension values used for the experiments were not always exactly the same, as they were partly controlled by the soil microtopography.But, they were generally between −15 and −10 hPa for the higher tension and between −7 and −3 hPa for the lower tension.

Evaporation method
Two undisturbed soil samples were collected in the field at each location (RT upslope, midslope, downslope, RG, F) by driving a PVC cylinder (height 9 cm and diameter 15 cm) into the previously wetted soil surface.The soil cylinders were air-dried, ensuring that the soil structure was preserved, and slowly saturated with tap water to control the initial water content.Four micro-tensiometers (3 × 10 −2 m length × 1.5 × 10 −2 m diameter) filled with de-aired water were connected to pressure transducers (Honeywell 15 PSI model, Honeywell Aerospace Plymouth, Plymouth MN) and inserted into the soil sample at four depths: at 0.011, 0.034, 0.056 and 0.079 m from the top of the soil sample (Fig. 2c).The cylindrical soil sample was dried from the upper surface in a ventilated thermostatic oven at an air temperature of 40 • C, imposing very slow evaporation conditions with a constant evaporation rate of 2.3 × 10 −8 m s −1 .Pressure and the weight were recorded during evaporation at a regular time interval (30 min) on a data logger (CR10 model, Campbell Scientific, Logan UT).The measurements continued until the air entry limit of the ceramic cup of the tensiometer was reached.The soil was removed from the PVC cylinder and oven-dried at 105 • C during 24 h, and the final water content and dry bulk density were calculated.As the evaporation rate was set to be slow with a quasi steady-state flux, the pressure head gradient in the soil sample was close to zero, with a uniform pressure head profile.Therefore, it was legitimate to derive the retention curve from the average tension measured in the four tensiometers and from the water content variation during the drying of the soil cylinder.The relationship of van Genuchten (1980), Eq. ( 1) with the Mualem condition (k = 1), has been used to fit to the experimental water retention curves.

The inverse method
The original method of Wind is known for introducing some biases into the calculation of the hydraulic conductivity near saturation, especially as the hydraulic gradient is low (Tamari et al., 1993;Wendroth et al., 1993;Richard et al., 2001).Therefore, the same experimental data set as for the evaporation method was used to derive the soil hydraulic parameters by numerical inverse modeling with Hydrus1D (Simunek et al., 2005).The height of the simulated domain was set to 0.09 m, with four observation points corresponding to the positions of the tensiometers in the soil cylinder.The boundary conditions were set to Neumann conditions; namely, the lower boundary condition was set to zero flux and the upper boundary condition was set to the experimental evaporation flux.The inversion procedure is based on the minimization of an objective function describing the difference between observed and computed values.The estimated unsaturated soil characteristics were chosen to be described by the van Genuchten (1980) relationship (Eqs. 1 and 4) with the Mualem condition (k = 1, Eq. 2).In order to improve the fitting procedure, the values of θ r , θ s , n, and α obtained with the evaporation method and the values of K s obtained with the disc infiltrometer were used as initial guesses for the inversion procedure.The water content parameters θ r and θ s were kept fixed, whereas the parameters n, α, K s were fitted, and l was set to 0.5.

Pedotransfer function
Many pedotransfer functions relating simple PSD to the soil water characteristics (SWC) and especially the water retention curve have been developed in the last decades (Pachepsky and Rawls, 2003;Mohammadi and Vanclooster, 2011).Therefore, a well-established model was chosen for this study.
The model proposed by Arya et al. (1999) based on the original work of Arya and Paris (1981) is a partly physically based model, deriving from the pore size distribution of a soil from the PSD data, according to the following equation:

S. Siltecho et al.: Estimating unsaturated hydraulic properties
where r i is the mean pore radius and R i the mean particle radius for the ith particle size fraction, e is the void ratio of the natural structured soil sample, n i is the equivalent number of spherical particles in the ith particle size fraction, and α i is a scale factor.The calculation details for these parameters can be found in Arya et al. (1999).Finally, the pore radii, r i , are converted into equivalent pressure heads h i , using the capillary rise equation where η is the water surface tension at the air-water interface, is the contact angle, ρ w is the density of water, and g is the acceleration due to gravity.The volumetric water content, θ i (m 3 m −3 ), is obtained from successive summation of waterfilled pore volumes according to where ε is the total porosity (m 3 m −3 ), and S w is the saturation rate.

Evaluation of the methods
Experimental soil water flow measurement devices, namely, tensiometers and soil moisture sensors (Enviroscan probe, Sentek, Stepney SE, Australia) were installed in the different soil locations, and were continuously recorded every 2 h for several years.In order to evaluate the goodness of the soil parameter estimation methods, a reference was needed.The final use for these parameters was modeling of soil water flow, and no intrinsic correct values were known, so these parameters were evaluated by modeling with Hydrus1D, a robust and well-documented software (Simunek et al., 2005), and by comparing modeled results with experimental data.The calculation has been performed on a uniform 0.35 m deep domain with two tensiometric boundary conditions measured in the field at 0.10 m as the upper boundary condition and at 0.45 m for the bottom boundary condition.This choice has been motivated by the necessity to perform the simulation with well-constrained boundaries.The development of matric potential during an almost 1-year period was calculated for the intermediate tensiometer located at 0.25 m.These calculations were performed with the combination of the unsaturated soil parameters measured with the different methods (i.e., the Beerkan, disc infiltrometer, PTF, evaporation and inverse methods) and were compared to the experimental tension values measured in the five different locations (i.e., in the forest, ruzi grass, and the three positions in the rubber plantation).The four criteria generally proposed were used to evaluate the performance of modeling (Table 2), namely the root mean squared error (RMSE), the coefficient of determination (CD), the modeling efficiency index (EF) and the coefficient of residual mass (CRM) (Loague and Green, 1991;Kim et al., 1999).

Statistics
Firstly, the estimated variables, namely saturated hydraulic conductivity (K s ), parameter α and parameter n, were analyzed using standard statistics to calculate their mean and coefficient of variation values.The statistical R package (R Development core team, 2008) was used to test the normality of data frequency distribution with the Anova test.However, as the number of repetitions were uneven and sometimes very low (from 3 to 10 samples per method tested in each site), Kruskal-Wallis non-parametric tests were also used to check the hypothesis of the same continuous distribution for the different cases.This test determines equality between means or medians of different groups of data, without presuming any specific hypothesis on the distribution.

Beerkan method
With a simple 3-D infiltration experiment combined with PSD analysis, the BEST software (Lassabatere et al., 2006) derived the main shape parameter n and scale parameter K s and α describing the unsaturated hydraulic properties.The results displayed in Table 3 show relatively low variation between the different sites.Saturated hydraulic conductivity (K s ) ranged from 5.59 to 10.23 × 10 −6 m s −1 within the different sites and showed a clear gradient along the slope in the rubber tree plantation.The highest values of K s were found upslope and on the lowest downslope, probably due to the accumulation of finer particles in the lower part, as mentioned by other authors (Heddadj and Gascuel-Odoux, 1999;Jing et al., 2008) describing topographic variation of hydraulic conductivity.Average K s values were found to be very similar for the soil under ruzi grass and under forest (5.84 and 5.93 × 10 −6 m s −1 , respectively), though the standard deviation was much higher for ruzi grass than for forest (3 × 10 −6 and 9 × 10 −7 m s −1 , respectively).Along the slope, the relative standard error on K s was also quite high, with values ranging from 49 to 62 %.The retention curve parameters, α and n, also showed little variation between the different sites: for most of the cases, the average value of α was 1.4 or 1.5 × 10 −4 m −1 , except for the forest (F), where α displayed higher values (2.8 × 10 −4 m −1 ).The shape parameter n followed a similar trend as K s , namely, decreasing values along the slope (2.33, 2.26 and 2.06, respectively, in the upslope, midslope and downslope positions).Outside the plantation, the average value for n in ruzi grass (R) was within the range previously observed, namely 2.23, but it was significantly Indicates that the degree of deviation between the experimental determinations and calculated values tends to zero when the calculated and experimental values tend to be equal.
Describes the ratio between the dispersion of experimental determinations and the dispersion of the calculated values, tending towards unity when the experimental and calculated values are consistent.
Indicates whether the model provides a better estimate of experimental determinations than the mean value of these determinations.The expected value for EF tends towards 1.
Indicates whether the model tends to overestimate (CRM < 0) or underestimate (CRM > 0) compared to experimental values.The optimal value for CMR tends towards zero.lower in the forest situation, with a value of 2. Normalized or scaled retention curves (Fig. 3d) were established with dimensionless values of water content (S e ) and of the matric potential (h * = h α), in order to simplify the retention curve h(θ ), and to compare the shape parameters.They showed very similar shapes for all the locations, with a small shift for the F soil samples, indicating differences in shape parameters for the soil in forest resulting from a slightly higher clay content.As hydraulic conductivity is mainly ruled by shape parameters like texture that vary less in local than in scale parameters, scaled hydraulic conductivity curves (not shown here) were all grouped together and were therefore less informative.

Disc infiltrometer
Results obtained for K s with the disc infiltrometer showed a behavior similar to that with the Beerkan method, namely, a progressive decrease down the slope in the rubber tree plantation, though the actual values for K s were systematically lower, namely, 6.9, 5.1 and 3.5 × 10 −6 m s −1 from upslope to downslope positions, respectively (Table 3).The average values of K s in the forest and ruzi grass were larger than those for the Beerkan method (6.65 and 6.93 × 10 −6 m s −1 , respectively).The relative standard errors for this method were significantly larger than for the Beerkan method, with values ranging from 66 to 90 %.The higher dispersion for hydraulic conductivity values derived from the disc infiltrometer could be explained by the fact that measurements are very dependent on the quality of the contact between the disc and the soil surface.Despite all the efforts to meet this requirement, it can be quite difficult to fulfill, and therefore affects the kinetics of infiltration.The precise infiltration surface area was also subjected to uncertainty when sand was applied to improve the contact, as the sand could overlap the actual disc surface.Moreover, in sandy soils, local hydrophobicity can occur (especially in rubber tree plantations, where natural rubber can modify water repellent properties of soil), and therefore affects infiltration dynamics and more specifically the measurements with the disc infiltrometer, as soil suction was the driving force.The presence of macro-pores related to biological activity (earthworms, ants, termites) could also partly explain the differences between the Beerkan method and the disc infiltrometer, as with this latter method, water can only infiltrate into the soil matrix.

Evaporation method and the associated inverse method
Theoretically, the water retention curve obtained from the evaporation method could be considered a reference, for it corresponds to direct experimental data with no model as-sumed a priori.In order to compare the results more easily with the other methods, the van Genuchten model with Mualem conditions was used to fit to the experimental retention curves.Parameter α varied only in a very narrow interval between 1.8 × 10 −4 and 2.5 × 10 −4 m −1 , and was found to be equal for both methods.Nevertheless, the values were found to be slightly larger than for the Beerkan method.
Values for shape parameter n are significantly different between the two methods, with higher values for the evaporation method.By comparison, the n values obtained with the Beerkan method are much lower than those derived from the evaporation and inverse methods.The scaled retention curves depicted in Fig. 3c for the evaporation method clearly show slight variations between the different situations.A single soil sample from the RT in the midslope position stood out, with a noticeably larger value for n.Saturated hydraulic conductivity obtained with the inverse method was much more important than when determined with other methods, as it ranged from 10 to 22 × 10 −6 m s −1 .Regarding K s along the slope, the same pattern was found, namely, showing an increase along the slope from the bottom to the upper position.
Unlike the disc infiltrometer and Beerkan methods, K s for ruzi grass and for forest soils was found to be very different, with a very high K s value for ruzi grass: 22 × 10 −6 m s −1 instead of 6 × 10 −6 m s −1 with the other methods.

Arya method
Van Genuchten water retention parameters were fitted on the experimental curves and showed little variation for parameter α between the different sites, reflecting a low variability of PSD.On the other hand, the values for α were found to be significantly lower than for the other estimation methods.The shape parameter n showed very high variability, with average values ranging from 1.82 for RG to 3.35 for RTmid and an important standard deviation running from 3 to 67 %.The scaled retention curves (Fig. 3b) for this method were very similar, though the forest soil samples showed a slightly particular behavior.Nevertheless, when all are represented on the same graph, the slight variations of scaled retention curves corresponding to the location of the soil samples seemed negligible compared to the differences due to the measurement methods (Fig. 3a).The value of the shape parameter n represented by the slope of the retention curve was lowest for the Beerkan method, intermediate for Arya method and highest for the evaporation method.Retention curves obtained with the Arya and Beerkan methods were, respectively, strictly based on PSD, and on PSD partly influenced by infiltration data.As for these methods, the values for n were lower, and the pore size distribution was not as sharp as for the evaporation method.One can therefore suppose that part of the porous volume was not drained with the evaporation method.In fact, by imposing a low evaporation rate to avoid overly important hydraulic gradients inside the soil sample, the energy necessary to draw out water from the smallest or less accessible pores is probably not sufficient.

Statistical analysis
The number of replicates for each method was not equal, because it depended on (i) the time necessary to perform the measurement, and on (ii) the quality of the measurement.For example, only a few evaporation measurements have been performed, as each measurement took up to 2 weeks to be completed.On the other hand, more Beerkan and disc infiltrometer measurements have been performed, as the infiltration experiments lasted only between 30 and 60 min, and between 2 and 6 h, respectively.The number of results was also controlled by the data processing of the infiltration experiment, as some experiments had finally to be discarded.
The statistical analysis of the van Genuchten unsaturated soil water parameters for the different locations with the different experimental techniques was performed with the Kruskal-Wallis non-parametric method.Both κ 2 and p values (Table 4) showed clearly that measured unsaturated soil water parameters were highly dependent on the measurement techniques.As already pointed out on the scaled retention curves, parameter n seemed to be the most dependent on the type of method, considering the high value of κ 2 and the extremely low p value.The shape parameter α showed the highest κ 2 value, but saturated hydraulic conductivity K s is also highly dependent on the measurement method.The influence of location was found to be secondary or even negligible, with a p value of 0.11, 0.94 and almost 1 for α, n and K s , respectively, except for α showing a relatively high value for κ 2 (7.15), but still under the acceptable limit (9.49) for the corresponding degrees of freedom (4).Saturated hydraulic conductivity K s was therefore not a discriminating parameter for the different situations, whereas parameter n tended to show an evolution with the different sites.When data were considered globally, regrouping all the measurement methods, the variability of the results was higher considering the measurement method rather than the measurement location.When studied separately for each measurement method (Table 4), some parameters seemed to differentiate clearly the different sites.For example, with the Beerkan method, the shape parameter n appeared to be significantly different for the different sites.With the pedotransfer function of Arya, scale parameter α was the most discriminating parameter, with an extremely low p value, though parameter n also showed a low p value (< 0.05).In both cases, it was shown that the particle size distribution and the derived parameters were significantly different at the various sites for these two measurement methods.With the others, no clear distinction between the different sites was observed, as the p values were very high, between 0.15 and 0.86.Table 4. Medians and parameters of Kruskal-Wallis, chi-squared (κ 2 ), degree of freedom df, and p value for the van Genuchten parameters for different methods and different sites, in general and in detail for each method, K s × 10 −6 (m s −1 ) and α (m −1 ).B: beerkan, I: inverse, A: Arya, E: evaporation, D: disc.* indicates a statistically significant difference between values.

Modeling validation
Unsaturated soil water parameters are generally determined for use in mathematical models simulating soil water flow.
Their evaluation was therefore performed with Hydrus1D, and the computed matric potentials at 0.25 m were compared to the experimental data.The different locations could be divided into two distinctive groups where (i) infiltration was strictly 1-D with vertical water flow, corresponding to the forest and ruzi grass, and (ii) with a possible lateral flow component due to the slope in the rubber tree plantation.The results obtained in the forest and at the ruzi grass sites showed good conformity with the experimental data (Fig. 4) and with RMSE values less than 40 hPa and EF values close to unity, especially for the forest site (0.83-0.95).The numerical simulations obtained with unsaturated soil water parameters issued from the different methods show very few differences.
In both cases (F and RG), the best fit was obtained with parameters derived from the combination of the Arya method for retention curves and the disc infiltrometer for hydraulic conductivity.Simulations performed with parameters from the Beerkan method showed a slightly higher RMSE value, whereas the results obtained with the evaporation-methodbased parameters seemed to fit worse to experimental data.These trends were confirmed by the other modeling evalua-tion indexes (DR, EF).CRM showed a systematic small negative value for the forest site, regardless of the measurement method, indicating an underestimation in the computed matric potential values.The scatterplots of computed vs. experimental matric pressure heads (Fig. 5) showed good agreement for all models in the forest and a slight underestimation for the high values in the ruzi grass.For the rubber tree plantation, concordance between computed results and experimental data was much worse (Figs. 4 and 5).The general trend was preserved but often overestimated compared to experimental data, as depicted in Table 4 and quantified by CRM in Table 5.During the rainy season, the fit between computed and experimental matric potential was generally better than during the drier period (Fig. 5).This discrepancy can be explained mainly by the slope in the rubber tree plantation and the inevitable subsurface lateral flow.Indeed, in the soil profile, an impervious clayey layer lying over the bedrock promoted the generation of a perched water table during the rainy season.Though the slope of the terrain was only 3 %, the actual slope of the interface on which the water table built up was more important, as the soil profiles increased from upslope to downslope.It has been shown that 40 % of the annual rainfall actually contributed to lateral flow along the slope (Seltacho et al., 2013).Therefore, the water flow in the field was not 1-D, and could not be com-  puted adequately; in upslope and midslope positions, water was lost laterally and accumulated in the downslope position.Nevertheless, when lateral flow was taken into account in a 2-D simulation with Hydrus2D (Seltacho et al., 2013), computed and experimental values fitted well.In any case, when contemplated from a yearly timescale, the differences between the different methods for simulating the water flow in soil seemed negligible in terms of water stock.Regardless of whether the computed results fitted the experimental data, as for forest and ruzi grass, or did not, as with the rubber tree plantation, no method could clearly stand out as being more efficient than any other one in modeling water flow for these different situations.These results are informative about the non-uniqueness of the parameters for modeling water flow in soil.Despite not describing exactly the development of soil water potential, different combinations of soil parameters lead to very similar results when used in Hydrus1D on a yearly timescale.The efforts to determine the soil hydraulic parameters precisely can therefore be seriously questioned.In any case, to improve the modeling performance for a longer time series or to forecast different scenarios, the parameters need to be adjusted by inverse modeling (Seltacho et al., 2013).

Discussion
A major result of this study is the apparent high dependence of unsaturated hydraulic properties on the measurement method.Several factors can explain this discrepancy, especially the size of the soil sample on which the measurement has been performed.The measurements have been performed on soil samples of different sizes for the different methods, depending on the availability of equipment and materials.For example, the diameter of the disc infiltrometer was almost twice the diameter of the infiltration cylinder.However, according to Anderson and Bouma (1997) and Bouma (1980), the representative elementary volume of sandy soil for measuring hydraulic conductivity is usually considered to be around 1 × 10 −4 m 3 .Consequently, considering the texture (mainly sandy) and especially the lack of structure of the soil in the different locations (except in the forest), the volumes of the soil samples exceed the representative elementary volume.Therefore, despite not having been measured on strictly the same volumes or areas, but still on the order of magnitude (Beerkan 8 × 10 −4 m 3 , disc infiltrometer 2 × 10 −3 m 3 , evaporation 1.5 × 10 −3 m 3 ), the results for the different methods should not be affected by the scale.
On the other hand, the differences can be explained by the specific properties of each method with their limitations and inherent assumptions for deriving the parameters.
The Beerkan method is popular for the straightforwardness of the experimental setup and rapid infiltration process; a constant infiltration rate is generally reached in less than an hour.However, the derivation of the unsaturated parameters is based on a rigorous hypothesis about the unsaturated hydraulic properties; namely, they are supposed to follow strictly the retention model of van Genuchten with Burdine's condition and the expression of Brooks and Corey for hydraulic conductivity.Even though this assumption agrees well in most of the cases, situations like bimodal porous networks, for example, are not taken into account.
With the disc infiltrometer, the experiment is more difficult to set up, as it needs a perfectly flat contact between the soil surface and the disc and is prone to many technical failures (leaks, etc.).Moreover, each experiment usually takes a very long time to reach a constant infiltration rate (sometimes several hours for finely textured soils).Wooding's model used to derive saturated hydraulic conductivity from disc infiltration measurements assumes an exponential relationship between hydraulic conductivity and matric potential that is quite different from the van Genuchten function.
The evaporation method is the only one for which the retention curves were actually measured without any a priori model.The slow evaporation rate imposed at the surface generated a very slight tension gradient inside the soil, with a uniform water content distribution.Average water content and pressure head values could therefore be calculated.The drawbacks of this method are the length of time needed for a soil sample to dry completely (up to 2 weeks) and the costly equipment (oven, computer, balance, microtensiometers, pressure gauge, data logger).
The pedotransfer function is an extremely easy method to derive the retention curve based only on particle size distribution.The Arya relationship is physically based, deriving the size of the voids between the grains by assuming a packing model.Nevertheless, as this model, unlike the Beerkan method, is exclusively governed by the PSD and the bulk density of the soil, little information about the the soil structure is available in the computed retention curve.
These technical and theoretical differences between the two infiltration methods can explain the contrast between their results.The shape of the hydraulic conductivity equations, Brooks and Corey for the Beerkan method and an exponential relationship for the disc infiltrometer, are different, especially near saturation.Saturated hydraulic conductivity determined through an exponential function with the Wooding method takes into account pressure head values near saturation, whereas in the BEST procedure, the Brooks and Corey equation is derived over a wider range of pressure head values.Therefore, saturated hydraulic conductivity determined with the disc infiltrometer is higher than when derived with the Beerkan method (Simunek et al., 1998b).When K s is derived by inverse modeling from the evaporation experiment, it is systematically higher than with the infiltration methods.The use of the van Genuchten equation (Eq.4) with Mualem conditions in the fitting procedure seems to be responsible for the overestimation of K s .More generally, the use of different type of equations, valid in different domains of pressure head or soil water content to derive K s , invariably leads to a wide range of values.The same conclusions could be drawn for the scale and shape parameters (α and n), as equivalent procedures were used to derive them.
Nevertheless, for each measurement method, a decrease in hydraulic conductivity down the slope in the rubber tree plantation was systematically observed, and could probably be related to the translocation of finer particles downslope (Wiriyakitnateekul et al., 2009).Despite showing slightly higher levels of organic matter, and unlike what would commonly be expected, forest soil had the lowest hydraulic conductivity.The higher content in finer particles (clay and silt) was most probably responsible for lower K s values at the forest site and to some extent at the ruzi grass site.During the early stages of the rubber tree plantation, the soil surface was not covered and was therefore more vulnerable to erosion.The fine soil particles were translocated downslope and down the soil profile, whereas under forest and ruzzi grass clay translocation was not as active.
The final use of these van Genuchten (VG) parameters was to use them in Hydrus1D to compute the water balance at these different sites.This study showed clearly that the uniqueness of VG parameters did not apply in this case, as many different combinations provided similar or equivalent results when a yearly timescale was considered.The simulation results obtained with the parameters derived from different methods were generally equivalent when the different validation parameters were considered.In any case, to perform efficient and reliable modeling of water flow in soil, these parameters should be carefully adjusted by inverse modeling procedures on experimental data over a time series describing contrasted situations and then validated.The measured VG parameters should be used as appropriate first guess values for the inverse modeling procedure.No single method could be considered without doubt as producing better results than the others.It should therefore be recommended to use the easiest and cheapest methods for the experimental evaluation of the first set of VG parameters and, in this case, the Beerkan method would surely be the best option.

Conclusions
In order to determine the unsaturated hydraulic properties of soil for different land use, several experimental methods have been used in a small watershed in northeastern Thailand.They included laboratory methods like the evaporation method, inverse methods and PTF (Arya et al., 1999), and also field evaluations like the Beerkan method and a disc infiltrometer.
Statistical analysis of the results obtained during this study showed clearly that significative differences in VG parame-ter measurements appeared, depending on the measurement method.Though the impact of land use and position was not completely negligible, the primary factor for measurement variability was found to be in the experimental methods.It was stated that the actual values of VG parameters depend on the method employed to determine them.However, when each measurement method was considered separately, the discrimination by site was significant for the Beerkan method and Arya's PTF method.As both methods are based on particle size distribution, parameters n and α, respectively, could distinguish the different sites.Unsaturated soil properties of this sandy soil seemed to be governed mainly by textural parameters related to the pedogenesis rather than by structural properties associated with land use and management.
Amongst the different methods tested in this study, none could clearly be considered superior to the others in terms of providing better parameters for modeling soil water flow.Therefore, the cheapest and easiest method to derive the VG parameters, like the Beerkan method, should be used to fulfill this task.The important land use changes taking place in northeastern Thailand can therefore be evaluated easily with numerical modeling and the consequences of rubber tree plantation for the water balance on different scales can therefore be predicted (Seltacho et al., 2013).
The Supplement related to this article is available online at doi:10.5194/hess-19-1193-2015-supplement.

Figure 1 .
Figure 1.Location of the experimental watershed with the different different land uses and measurement sites; F: forest, RG: ruzi grass, RT-up: rubber tree upslope, RT-mid: rubber tree midslope, RT-down: rubber tree downslope, and location of meteorological station MS.

Figure 2 .
Figure 2. Photographs of some of the experimental measurement methods: (a) disc infiltrometer; (b) Beerkan method; and (c) evaporation method.

Figure 3 .
Figure 3. Scaled retention curves for (a) the Beerkan method, (b) the evaporation method, (c) the method based on Arya's PTF, and (d) the different methods.

Figure 4 .
Figure 4. Development of the matric potential at 0.25 m depth in the different locations (black dots) and results of modeling with Hydrus1D (continuous lines).

Figure 5 .
Figure 5. Scatterplots of computed vs. experimental matric potential values in the different locations for different combinations of parameters.

Table 1 .
Physical and chemical properties of the topsoil at the different sites.

Table 2 .
Evaluation criteria for performance of modeling.

Table 3 .
Means and standard errors for the main van Genuchten parameters with Mualem conditions (n = 1/{1 − m}) determined for the different methods at the different sites.The number of samples is indicated in parentheses.

Table 5 .
Four modeling evaluation criteria for matric potential for the different situations with the different parameters used in Hy-drus1D.