Interactive comment on “ Using an inverse modelling approach to evaluate the water retention in a simple water harvesting technique ”

The article is of great interest particularly, as the authors rightly say, few efforts have been made to quantify and publish the hydrological significance of improvements in small-scale agriculture. This is also partly evidenced by few very recent cited articles in the manuscript. The HYDRUS model is indeed an appropriate tool to explore this area of science. The article is well written which makes it very easy to follow. There are, however, some areas which, in my view, require attention for further improvement:


Introduction
Arid and semi-arid zones are characterized by an important deficiency in water availability for plant growth.On the other hand, precipitation often comes in the form of short bursts of high intensity rainfall, causing rapid saturation of the uncovered soil surface and promoting soil erosion, flash floods and mud flows in extreme cases.In Andean arid lands, a range of agricultural solutions for these conditions were implemented by a large number of indigenous communities, as extensively documented by Denevan (2001).These technologies were designed to improve the crop environment, increase labour efficiency, enhance sustainability, improve productivity, and to minimize risks from changing environments, especially unpredictable climatic conditions.The fact that many of these systems are still used under present day conditions is a strong indication that they are sustainable.This has led to the belief that ancient Andean populations managed natural resources in those areas better than we manage them today (Murra, 1983).Browman (1987) stated that the principal economic strategy of former Andean arid land producers was risk management, focusing their efforts on water management, frost reduction, erosion control and soil accumulation.Pandey et al. (2003) investigated the apparent link between climate changes resulting in droughts and the increase in the use of water harvesting techniques (WHT) throughout past civilizations, indicating that WHT can partially alleviate the negative climatic conditions.They also stated that traditional systems would become more efficient if scientific attempts be combined to enhance the productivity of local knowledge systems.
More recently, modern versions of these WHT have become an important tool in the efforts against the desertification in Latin America and, more specifically, in Chile.A common WHT used for this purpose is the infiltration trench, a depression in sloped soils to capture runoff water, comparable to ancient indigenous dug ditches that are still used in northern Peru (Denevan, 2001).Implementation of these WHT in semi-arid and arid zones of Chile has increased dramatically over the last years, from 52 hectares in 2001 to 2200 hectares in 2003 for reforestation purposes (Pizarro et al., 2004), due to strong incentives.Although various demonstration projects were realized (CONAF, 1988(CONAF, , 2000;;Pizarro et al., 2003), some even dating back to 1975, few efforts were observed to quantify the effect of runoff water harvesting techniques on water retention and soil conservation.
Recent developments in soil hydrology allow the use of complex distributed models to describe hydrological processes at the field scale.Due to increased computational power, for example, the numerical approximation of the Richards equation (Richards, 1931) has only recently been implemented, thereby obtaining a more accurate modelling of the infiltration-runoff process, compared to more empirical or analytical methods that were previously used (e.g. Green and Ampt, 1911).Examples of such models equipped with calculation procedures to solve the Richards equation numerically are TOUGH2 (Pruess et al., 1999), FEHM (Zyvoloski et al., 1997) and HYDRUS ( Šimunek et al., 2006) for 2-D and 3-D applications.Due to their dependency on physically meaningful parameters, these models are equipped with parameter optimization techniques that minimize a suitable objective function which expresses the discrepancy between the observed values and the predicted system response ( Šimunek et al., 2006).Originally, in soil science inverse modelling was applied only under laboratory conditions on soil columns using a one step (e.g.Kool et al., 1985) or multi-step drainage approach (e.g.Van Dam et al., 1994;Hopmans et al., 2002), and later used in field scale applications (e.g.Zijlstra and Dane, 1996;Mertens et al., 2006), with the early field application by Dane and Hruska (1983) as an exception.An important advantage of inverse procedures, if formulated within the context of a parameter optimization problem, is that a detailed error analysis of the estimated parameter can be obtained.However, while parameter optimization methods provide several advantages, a number of problems related to computational efficiency, convergence, and parameter uniqueness remain to be solved, especially when many hydraulic parameters should be estimated simultaneously (Van Genuchten and Leji, 1992;Zijlstra and Dane, 1996).
These advances in soil hydrology make it now feasible to use these complex models to evaluate and improve ancient irrigation structures and related soil-management and irriga- tion practices (Van Genuchten and Šimůnek, 2004).In this paper, we describe a method to evaluate the water balance of a simple WHT using a parameter estimation procedure which combines a Marquardt-Levenberg nonlinear parameter optimization with a numerical model solving the variably saturated flow equation.Objective functions are formulated in terms of measured cumulative infiltration, soil-water content and soil-water retention data or are based on a combination of these measurements.
In the first part of the paper hydraulic properties of the study area are discussed in view of its usage inside the HYDRUS-2D model, and different measurement strategies are compared.Based on this best estimate, soil physical parameters are allowed to optimize using an inverse modelling approach with three different and independent data sets.The final model with the calibrated and validated parameter set is then used to evaluate the water harvesting processes for the infiltration trench under study.

Field site
All the field measurements were performed on a hillslope near the town of Quebrada de Talca (30 • 00 45 S, 71 • 02 37 W), in the greater arid Elqui Valley, Chile.At the field site, soil and water harvesting techniques were installed in 1997 to diminish soil losses, to reduce flash flood hazards to the town centre, to improve local infiltration and to increase the soil-water storage and plant growth locally (Pizarro et al., 2003).Apart from the wooden dikes and stone dams installed in the existing preferential pathways, infiltration trenches, which are the subject of this study, were dug at the lower part of the slope, perpendicular to the flow lines and parallel to the contour lines (Fig. 1).A total of 184 infiltration trenches were constructed at the site, with an average volume of 0.23±0.07m 3 and a spacing along slope between trenches of 5.14±0.82m (Table 1).
The study area has an arid Mediterranean climate, and, based on data from 1980 to 2000, characterized by (1) an average annual precipitation of 99.2 mm, of which more than 70% is produced in the southern winter season; (2) moderate temperatures, with an absolute minimum of 2 • C (June) and an absolute maximum of 30 • C (March); (3) a high relative humidity (80%) with frequent cloudiness; (4) an average annual solar radiation of 4075 cal cm −2 d −1 , resulting in a water deficit of 800 to 1000 mm year −1 (Miller, 1976).The climate is categorized as "arid" using the Aridity Index proposed by UNEP (Middleton and Thomas, 1997), and scarce precipitation is often concentrated in short bursts of high intensity rainfall, leading to a high Modified Fournier Index (Fournier, 1960; also called the Climatic Aggressivity Index) in 10% of the years.Due to these climatic conditions, vegetation cover in the area is limited and mainly composed of shrubs, herbs and cacti (Miller, 1976;Olivares and Squeo, 1999), leading to a high exposure to runoff and erosion risks.Due to climatic limitations and steep slopes, dominant land use is pasture, although hi-tech irrigated agriculture, especially grapes, is found on similar slopes with access to water.

Field measurements
In 2005, eight years after its construction, the water balance of one infiltration trench was studied in more detail.Therefore, a field plot of 6×2 m was selected on the hillslope (indicated as a grey square in Fig. 1) and consisted of a trench with its impluvium (10 m 2 ), under a slope of 23%.Twenty two 30-cm long Time Domain Reflectometry (TDR) probes, connected to a TDR100 device and datalogger (Campbell Scien-tific, Loughborough, Logan, USA 1 ) were installed horizontally at a depth between 7 and 45 cm below the soil surface.The equation proposed by Topp et al. (1980) was used to convert measured relative dielectric constants into volumetric water contents, since it fitted well (R 2 =0.96) to a prior established calibration data set.A rainfall event with an intensity of 120 mm h −1 was simulated for 20 min -the total amount of rainfall effectively applied during the event was hence 40 mm -using a rainfall simulator similar to the one described by Verbist et al. (2003).In brief, it consisted of a sprinkler boom on which nozzles (type TeeJet TG SS 14w) were fixed at a distance of 1 m from each other and were positioned at a height of 1.8 m above the soil surface.The characteristics of the rainfall simulator, the drop size distribution, average drop velocity and the overall kinetic energy of the simulated rain are further described in Gabriels et al. (1997) and Erpul et al. (1998).During the simulated 20-min long rainfall event, the advancing of the wetting front was monitored with the TDR probes and soil-water content was further monitored after the rainfall event until no significant changes were observed for a 24 h period, which was 4 days from the start of the rainfall simulation.The water level in the infiltration trench was monitored, reaching its maximum volume at 17 min (overflow was observed), and was continued until all water had infiltrated 61 min after the start of the rainfall simulation.Immediately after the rainfall simulation, the field plot was covered with a plastic cover, to prevent evaporation.
After the rainfall simulation and moisture content monitoring, soil texture, bulk density and the soil-water retention characteristics were determined on twenty five undisturbed soil samples taken using standard sharpened steel 100 cm 3 Kopecky rings at various depths (0-5 cm, 15-20 cm, 25-30 cm and 35-40 cm) at 1 m interval distance along the slope.Soil texture was determined with the pipette method (Gee and Or, 2002), whereas organic matter measurements were based on the Walkley and Black (1934) method.The basic physical and chemical characteristics of the soil are summarized in Table 2.In the shallow soil profile (0-0.45 m) texture is homogeneous (sandy loam), with a slight increase in the sand fraction when approaching the parental material.Due to the coarse sandy loam texture, drainage is optimal and salt content is low (EC<2500 µS cm −1 ).Soil bulk density is rather high, as is typical for arid soils.The soil-water retention curve was established using the sand box apparatus (Eijkelkamp Agrisearch Equipment, Giesbeek, the Netherlands) for pressure potentials between −1 and −10 kPa, and pressure chambers (Soilmoisture Equipment, Santa Barbara, CA, USA) for pressure potentials between −20 and −1500 kPa.Infiltration measurements were performed at the selected field plot to quantify the field unsaturated and saturated hydraulic conductivity (K sat ) as an input parameter for the modelling phase.Kool et al. (1987) stated in their review of parameter estimation techniques, that poor identifiability of K sat could often be observed, resulting in large confidence intervals on the estimated value, related mainly to the lack of independent measurements of the parameter.Therefore, special attention was given to identify K sat at the field site, using four different measurement techniques.First, the field saturated hydraulic conductivity of the soil surface was determined using a single ring pressure infiltrometer (Soilmoisture Equipment, Santa Barbara, CA, USA) at a distance of 1 m, 2 m, 3 m, 4 m and 5 m from the infiltration trench, whereas a well permeameter (Soilmoisture Equipment, Santa Barbara, CA, USA) was used to determine K sat at the depths of 15 cm, 30 cm and 45 cm at the same locations on the selected field plot.For each experiment two positive pressure heads were applied, 5 and 10 cm, which allows for the use of both the simultaneous equation approach (Reynolds, 1993), as well as the single head analysis (Reynolds and Elrick, 1990;Elrick and Reynolds, 1992) for calculating K sat .At each location, readings were taken until infiltration rates were steady for at least three consecutive time intervals.K sat was calculated as the geometric mean of measured K sat values (Reynolds et al., 2000).
Additionally, a tension infiltrometer with a diameter of 20 cm (Soilmoisture Equipment, Santa Barbara, CA, USA) was used as a second means of K sat determination, following the parameter estimation method outlined in Šimůnek and van Genuchten (1997).They found unique solutions for K sat when using multiple tension infiltration measurements in combination with knowledge of the initial and final water content.At the start of each measurement, the tension infiltrometer's porous plate was placed on a previously prepared contact sand layer and a series of decreasing pressure heads was applied: −11.5 cm, −9 cm, −6 cm, −3 cm, −1 cm and −0.1 cm.Initially, to achieve a homogenous water content, the soil was saturated shortly before the start of the measurement by allowing the infiltration of 15 cm of water prior to the tension infiltration measurement, thus avoiding problems with the identification of the initial water content, which was reported to be of high importance to the identifiability of K sat ( Šimunek and van Genuchten, 1997).The water level in the water reservoir was monitored 15 min for each pressure head applied, resulting in a total duration of 90 min for each measurement.The Disc software ( Šimůnek and van Genuchten, 2000) was used to obtain estimates for K sat as well as the K(h) relationship as described by van Genuchten (1980; see below).Initial estimates of the van Genuchten coefficients θ s , θ r , n and α were obtained from an independent water retention analysis on undisturbed soil cores in the laboratory, using the RETC software (Van Genuchten et al., 1991).
Finally, two rainfall simulations were performed near to the selected field plot to determine the field saturated hydraulic conductivity and the cumulative infiltration characteristic under rainfall using the rainfall simulator as described above.During each 20-min simulation, a rainfall intensity of 120 mm h −1 was applied on two bare 6 m 2 plots simultaneously (plots A and B).The rainfall amount and intensity was measured by placing 20 cups with an inner diameter of 9.6 cm in and at the border of the plots.Runoff discharge was measured continuously using calibrated beakers, obtaining infiltration rates by subtracting the observed runoff data from the rainfall intensities.K sat was then computed by fitting the well-known Philip equation (1957) to the cumulative infiltration curve: where I (t) is the cumulative infiltration (m 3 ), S is the sorptivity (m s −0.5 ) and K a constant (m s −1 ) which approaches the hydraulic conductivity K sat at high values of time t (s).

Modelling
The HYDRUS-2D ( Šimunek et al., 2006) software package was used for simulating infiltration and water flow in the soil domain represented in Fig. 2. The program numerically solves the Richards equation for saturated or unsaturated water flow and implements the soil-hydraulic functions of van Genuchten (1980) who used the statistical pore-size distribution model of Mualem (1976) to obtain a predictive equation for the unsaturated hydraulic conductivity function in terms of soil water retention parameters: (2) where θ r and θ s denote residual and saturated water content, respectively (m 3 m −3 ), h the matric head (kPa), α is the inverse of the air-entry value (m −1 ), n is a pore size distribution index >1, λ is a pore-connectivity parameter (-), S e is the effective water content (-) and m=1-1/n.The parameters α, n and λ in HYDRUS-2D are considered to be merely empirical coefficients affecting the shape of the hydraulic functions.The pore-connectivity parameter λ was estimated to be about 0.5 as an average for many soils (Mualem, 1976) and is often fixed at this value (e.g.Schwartz and Evett, 2002;Šimunek et al., 1998), resulting in five independent parameters: θ r , θ s , α, n, and K sat .
The software program includes a Marquardt-Levenberg type parameter optimization algorithm for inverse estimation of soil hydraulic parameters from measured data for twodimensional problems.In order to solve the Richards equation in time and space, initial and boundary conditions must be specified by the user.Three different boundary conditions are considered for the modelling problem, as outlined in Fig. 2: an atmospheric boundary condition located at the soil surface mimicking rainfall, a variable pressure head boundary condition in the infiltration trench used for simulating the filling up of the trench during runoff accumulation and subsequent drainage, and thirdly, a free drainage boundary at the bottom of the soil domain.It should be noted that HYDRUS-2D does not model surface runoff processes directly, requiring accumulated runoff to be specified as a variable pressure head inside the infiltration trench, which was calculated from the water-depth measurements during and after the rainfall event.As an initial soil-water boundary condition, we used the average soil-water content measured for the TDR probes at the start of the experiment (0.11±0.02 m 3 m −3 ).Ten TDR probes were selected for calibration purposes and used in the inverse modelling optimization procedure.The selection assured us of a good distribution of soil-water data, both below the impluvium and around the infiltration trench, as presented in Fig. 2.
The unknown, or estimated, soil hydraulic parameters were determined by matching simulated and observed data from TDR water content measurements, (t), cumulative infiltration measurements from the rainfall simulation experiments, I (t), and soil water retention measurements, (h).values obtained from single head analysis and multiple head analysis for pressure head infiltrometer measurements at four depths along the field slope a .

Single head analysis
Multiple head analysis Depth The objective function, ( (t), I (t), (h)), describes the difference between observed and measured values and is minimized by the model to yield the best possible fit between them.In the optimization procedure, six different objective functions are considered: (I (t)), using only cumulative infiltration data; ( (t)), using only moisture content data; both of these data sets together, ( (t), I (t)), or a combination with the water retention data, ( (t), (h)), (I (t), (h)), and with ( (t), I (t), (h)) incorporating all data sets.The latter combination was written as: where β represents the optimized parameter set; u i , v i and w i , v j and u k are weighing factors for individual measurements and measurement types and were determined using the method described in Hopmans et al. (2002), m is the number of TDR probes in the calibration data set and n 1 to n 3 are the number of observations for each of the data sets.
Model performance was evaluated with two functions: the Root Mean Squared Error (RMSE) and the Index of Agreement (IA), according to Willmott (1982), which calculates the deviations between observed and estimated values and which indicates a perfect agreement when approaching 1.

RMSE
where O i and X i represent measured and simulated values of a certain variable, and O i and X i their mean value.

Initial parameter estimation
Table 3 shows the hydraulic conductivity K sat obtained with the pressure infiltrometer at the soil surface (depth of 0 m) and with the well permeameter (depths of 0.15, 0.30 and 0.45 m) both with the single and multiple head analysis.Note that in case of the multiple head analysis, only those measurements resulting in positive K sat values were used in the calculation of the geometric mean.Negative K sat values have been reported in a number of cases (e.g.Elrick and Reynolds, 1992;Logsdon and Jaynes, 1993) and are mostly related to non-homogeneous initial water contents and heterogeneous soils.In the multiple head analysis, 30% of the measurements had to be discarded due to negative K sat results, against only 2.5% in the single head analysis.As is typical for K sat measurements (see e.g.Reynolds et al., 2000), a large range of K sat values could be observed per depth and per method.When comparing both calculation methods, the single head analysis resulted in lower coefficients of variations within the same soil layer, compared to the multiple head analysis.Furthermore, the variation in mean K sat among soil layers was smaller for the single head analysis compared with the multiple head analysis, where large differences in mean K sat between layers were observed.Nonparametric statistical analysis showed that for both calculation methods, three of the four layers were not significantly different, although the layers differed for both approaches.In Hydrol.Earth Syst.Sci., 13,[1979][1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990][1991][1992]2009 www.hydrol-earth-syst-sci.net/13/1979/2009/ the case of the multiple head analysis, this lack of difference between layers was probably due to the large standard deviation observed in the calculated K sat for each soil layer.These results suggest that the single head analysis is a more appropriate technique to determine K sat from well permeameter and pressure infiltrometer measurements, since the obtained values have lower confidence intervals compared to the multiple head technique and each pressure head measurement results in an independent K sat value.
Table 4 shows the K sat values as computed from infiltration measurements with the tension infiltrometer and inverse modelling with the Disc software ( Šimůnek and van Genuchten, 2000) in combination with the saturated water content at the end of the experiment.The K sat values were calculated for each individual measurement.Also shown are the 95% confidence interval, the Root Mean Squared Error (RMSE), the coefficient of determination (R 2 ) and the mass balance error (MBE) observed during simulation.In all cases, the model achieved convergence on the solution, and resulted in very high correlation coefficients and low RMSE values, indicating a very good agreement between observed and simulated values.The confidence intervals on the K sat values had a rather narrow range, showing that the model was sensitive to the K sat value and thereby most likely resulted in a global minimum of the objective function and a unique value for K sat .The mass balance error gives an additional indication that no model errors were generated during the simulations, since reported values are well below 1% for all model runs.As an example, Fig. 3 shows the measured and optimized cumulative infiltration curve and the residuals for the soil surface at a 1.5 m distance upward from the infil- a S is the sorptivity, K sat is the saturated hydraulic conductivity, RC is the runoff coefficient, expressed as the percentage of precipitation that appears as runoff tration trench.A very good fit was found between observed and fitted cumulative infiltration data, showing that the model was perfectly capable of reproducing the measurement conditions.This is also reflected in the differences between observed and fitted values that oscillate around zero and never attain values higher than 1.5% of the total infiltrated volume.Similar results were obtained for the other infiltration measurements.A non-parametric Kruskall-Wallis test showed that K sat can be considered equal in all layers, with a geometric mean value of 7.72×10 −6 m s −1 .
Table 5 shows the regression coefficients of the Philip equation ( 1957) as applied to the rainfall simulation data.A very good consistency was found between the observations and the fitted curves (Fig. 4), resulting in high coefficients of determination.The geometric mean K sat value from this approach was 7.19×10 −6 m s −1 .To allow the comparison between the various approaches, a box plot summarizes all K sat results (Fig. 5).Large differences between the approaches are observed for the surface layer and the 45 cm depth, whereas similar values were found at 15 and 30 cm depth.The multiple head approach on the pressure infiltrometer and well permeameter data re- sulted in high fluctuations and less reliable K sat estimates, probably due to a lack of larger data sets.The single head approach gave better results, which were comparable to those measured with the tension infiltrometer, except for the surface measurement.This can be explained by the disturbance of the soil surface when inserting the ring which, according to Reynolds et al. (2000), influences the infiltration rates causing preferential flow on the sides and, thus, higher K sat values.Pressure infiltrometer and well permeameter measurements in general were more variable compared to the tension infiltrometer measurements, with outliers influencing the mean K sat value.As a result, field saturated hydraulic conductivities from the tension infiltrometer measurements were preferred as an initial estimator for the modelling phase.The rainfall simulation results are used as an independent data set for inverse modelling purposes and therefore, can not be used as initial estimates.However, they gave similar saturated hydraulic conductivities at the surface when compared to those obtained with the tension infiltrometer, confirming the usability of these TI measurements as a first estimate of the saturated hydraulic conductivity in the modelling process.
Figure 6 represents the soil-water retention data from 25 soil cores, five for each of the different soil layers.As shown in Table 6, the van Genuchten (1980) equation fitted the data rather well, resulting in low RMSE and high coefficients of determination (R 2 ).Table 6 also lists the values of the van Genuchten parameters as obtained by nonlinear least squares analysis using the Marquardt-Levenberg algorithm.The differences between layers proved to be rather small, resulting in similar van Genuchten parameters when using all water retention data together and an R 2 of 0.9.Therefore, it was preferred using a uniform water retention curve for the entire profile in the modelling phase, with its corresponding initial estimates for θ r , θ s , α and n.

Parameter optimization
Starting with the field measurements discussed in the previous section as initial estimates, the physical soil parameters were optimized by minimizing objective functions containing independent measurements of soil-water content, cumulative infiltration and the soil-water retention characteristic.In Table 7, optimization results are given for the six different objective functions.Two parameters, θ r and λ, were intentionally fixed at their estimated value to prevent nonuniqueness of the solution.
With the exception of the model run using all data sets, the objective functions incorporating cumulative infiltration data performed poorly, evidenced by large confidence intervals on all optimized parameters and large RMSE values compared to the other objective functions.The index of agreement clearly marks the lack of correlation between observed and simulated water content values and a bad fit between measured and simulated water retention data for those objective functions.Due to the relatively large differences in the cumulative infiltration curves between the four replicates (Fig. 4), no exact parameters could be obtained during the optimization process for these objective functions.In the case where water retention data was added to the infiltration a θ r , θ s , α and n are the van Genuchten parameters, K sat is the saturated hydraulic conductivity, θ is the water content, h is the matric head, I is the cumulative infiltration, t is time, IA is the index of agreement and RMSE is the root mean squared error b not optimized c values preceded by ± give the standard deviation on the optimized value d the model did not converge to a solution for this objective function data, ( (h), I (t)), the model did not converge into a solution, indicating that a (global) minimum could not be reached for this combination of data sets and, therefore, parameters could not be optimized.The model run using both measured moisture contents and infiltration data did not result in an acceptable match to the water retention curve either (Fig. 7), with large confidence intervals on all estimated parameters, indicating non-uniqueness of the solution.
The model runs optimized with the water content data alone in the objective function, (θ(t)), and also where water retention data was added as well, (θ(t), (h)), resulted in well-defined minima, with small confidence intervals on the parameters, and a good match with the measured water content and retention data.The objective function including all independently measured data set, (θ(t), (h), I (t)), performed equally well, with a simulated water retention curve closely following the measured water retention data, as shown in Fig. 7.
When comparing the optimized van Genuchten parameters with their initial values, it can be seen that θ s was reduced during the optimization procedure for the objective functions including water content data.It is a common observation in field soil-water monitoring that θ s becomes smaller in comparison with the laboratory value, because of incomplete saturation of the soil profile due to air entrapment during the wetting front advancement (e.g.Schwartz and Evett, 2003).The optimized K sat value was slightly higher than the initial value, but remained within the range of values measured by all instruments (Fig. 5), suggesting that the K sat measuring methods can be used to determine initial estimates for water balance studies, but that parameter optimization is necessary to obtain an effective parameter set that results in good agreement between measured and simulated values of infil-tration and water redistribution.Although this seems to be a deficiency of the model at first, differences between in-situ measurements and model parameters have been observed by many researchers (e.g.Mertens et al., 2005).This is most probably due to the different temporal and spatial resolutions at which measured values (K sat in this case) are obtained, compared to the scale at which the processes are modelled.
These results also indicate that water content readings contribute considerably to the identifiabilty of hydraulic properties under field conditions and should be included in obtaining a correct parameter optimization, when using simulated rainfall as an input to inverse modelling.Other researchers (e.g.Šimunek and van Genuchten, 1997;Schwartz and Evett, 2003) also reported that including water content measurements or pressure head measurements was essential to obtain correct parameter values with narrow confidence intervals for their applications.Since root mean squared errors were the least on the objective function including all prior knowledge, the optimized parameters from this model run were used in subsequent analysis.
As expected from the small RMSE values observed for the water content data used in the optimization process, the observed soil-water content values versus simulated water content values yielded a close fit and a coefficient of determination of 0.93 (Fig. 8).With the data points scattered around the line of perfect agreement, no bias in the simulations was observed.Since the initial water content of the soil was set to a fixed value of 0.11 m 3 m −3 , no values below this water content are simulated, even though the observed initial soil-water content ranged from 0.08 to 0.125 m 3 m −3 .For the twelve probes not included in the optimization process, a good agreement between observed and simulated values could still be found, with an R 2 of 0.78 (Fig. 9a).The overall coefficient of determination for all water content data was 0.85 (Fig. 9b).
When considering the TDR probes individually, very good matches were obtained.Examples are shown in Fig. 10a for the simulations with the highest coefficient of determination.All peaks were satisfactorily reproduced by the model, as well as the slower rising limbs, corresponding to probes at greater depths.Even for probes not present in the calibration phase, such as probe 1, good matches and high correlations were observed.Four probes were modelled with lower accuracy, as illustrated in Fig. 10b.For those probes that did not show any soil-water changes during the duration of the experiment, i.e. probe 3,4,9,10,20,21,22,23 and 24, the model showed no changes with respect to the initial soil-water content, indicating that the deepest penetration of the soil-water front was modelled as observed during the experiment.
Possible sources of variation influencing the overall correlation between observed and simulated soil-water contents in time are: (1) the assumption of a homogeneous soil profile, which is a simplification of the reality; (2) the exact physical location of the TDR probe, since relative small changes in the location of the probe in the model give relative large changes in soil water values in function of time; (3) the assumption of an equal initial soil water content for the whole profile that gives rise to a considerable amount of unexplained variability, such as observed for probe 1 in Fig. 10a, thereby reducing the correlation coefficient; (4) the measurement error associated with the TDR probes and the indirect method of measuring water content using the TDR principle, which is rather small, and identified to be less than 0.015 m 3 m −3 for the cable lengths used in this study (Bilsky, 1997).

Runoff water harvesting
The soil water content changes during the rainfall event and model simulations were used to visualize the infiltrationrunoff and water harvesting process and the effect of the infiltration trench on water availability for potential plant growth.Figure 11 shows the distribution of water content at different time steps: (a) at 20 min, when the rainfall simulation experiment had stopped; (b) at 61 min, when all water had infiltrated in the infiltration trench, (c-e) at the intermediate time steps of 150, 800 and 3000 min and (f) at 5800 min, when soil-water measurements had ceased.A clear increase in water storage can be observed due to the presence of the infiltration trench.
During the simulated rainfall, infiltration excess was observed, resulting in runoff that was partly captured by the infiltration trench.The water balance of the trench can be expressed in terms of its boundary conditions, with the source fluxes in the soil domain given by the rainfall amount and the variable pressure head, with values of 0.46 m 3 and 0.09 m 3 , respectively.The sinks are the total runoff over the whole slope (0.25 m 3 ), the accumulation of water content in the soil domain (0.28 m 3 ) and the drainage (0.02 m 3 ).The runoff coefficient calculated by the model indicates that 56% of the rainfall was lost to runoff, which corresponds to the field measurements (Table 5).The 10 m 2 impluvium area contributing to the infiltration trench generates an overland flow of 0.18 m 3 during the 20-min long simulated rainfall event, of which 0.09 m 3 (48%) was captured by the infiltration trench.The runoff water that could be collected by the infiltration trench, infiltrated slowly until 40 min after the precipitation had ceased and resulted in a significant increase in soil water storage in the direct neighbourhood of the infiltration trench.After 4 days, soil-water was redistributed and the highest water content was observed near the soil surface and within a radius of 25 cm surrounding the infiltration trench (Fig. 11).These results suggest that, at least under the experimental conditions, the best location for reforestation would be inside the infiltration trench, where plants can optimally benefit from the infiltrated water.

Summary and conclusions
In this study, a combination of rainfall simulations, soil water content monitoring and model simulations was used to evaluate the water balance of a simple water harvesting technique.Through inverse modelling, a calibrated model was constructed to allow visualizing the water harvesting and storage process of an infiltration trench often used in semiarid zones in Chile and other similar regions in the world.
In the first phase, hydrophysical parameters were obtained from independent field and laboratory measurements.Four techniques were used to determine K sat , which allowed the comparison of measuring methods as effective estimators of hydraulic properties.The tension infiltrometer measurements with multiple heads were selected as the best estimator, and these results were found comparable to those obtained from rainfall simulations.Pressure infiltrometer and well permeameter measurements showed more scatter in comparison with the other two methods, with the single head analysis of the data given preference over the multiple head analysis.Large variations in the K sat values obtained were observed for the multiple head approach and negative values were calculated regularly, probably due to soil heterogeneity and soil disturbance due to ring insertion or hole drilling.
By minimizing objective functions with different data sets, we learned that water content data was essential in providing effective model parameter sets and to attain narrow confidence intervals on their values.Infiltration data alone obtained from rainfall simulations were found to be too Hydrol.Earth Syst.Sci., 13,[1979][1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990][1991][1992]2009 www.hydrol-earth-syst-sci.net/13/1979/2009/ variable for these purposes, and should be combined with water content and water retention measurements to reach a unique optimized parameter set.After correct optimization, measured values of the soil's physical properties deviated little from the final model calibrated values, indicating the model was capable in representing the hillslope conditions adequately.
Model calibration, performed with 10 out of 22 soil water probe readings, proved successful, bringing simulated soilwater contents in close agreement with observed values.An overall coefficient of determination of 0.85 using all 22 TDR probes indicated that the model was capable of reproducing the observed water infiltration in the field experiment.The heterogeneity of the initial water content boundary condition in the soil profile prior to the measurements was found to be an important contributor to the unexplained variance.
The moisture content profile of the slope showed a marked increase in the water content just below and around the infiltration trench during and after the rainfall event, confirming the water harvesting potential of these techniques for reforestation purposes.Nevertheless, it should be stressed that the increase in soil moisture in the root zone is spatially very limited, a conclusion which is not widely known among users at present and which should eventually be taken into account to improve the effectiveness of existing Chilean financial incentives to farmers adopting the technique for reforestation.
However, to fully utilize the potential of these techniques, the determination of correct spacing of the trenches under varying soil and slope conditions remain a further objective for study, for which modelling tools such as the one described in this paper are especially appropriate.In this respect, additional efforts to upscale these results to the hillslope level would be a logical step forward.
Additionally, for a full evaluation of water retention efficiency of this water harvesting technique, natural rainfall events should be used to confirm actual results.Since HYDRUS-2D does not allow for simulating rainfall-runoff processes directly at present, but needs to be modelled using variable pressure heads inside the infiltration trench, this approach can only be used when detailed field measurements are performed during these rain storms.However, a different approach could be developed by using fully integrated surface-subsurface models, such as the one recently described by Sudicky et al. (2008), that reduces the necessity of monitoring water heights in the infiltration trenches, although they still would be useful for calibration-validation purposes.

Fig. 1 .
Fig. 1.Digital terrain model of the hillslope under study, with the installed infiltration trenches indicated in black and the selected field plot as a grey square.

Fig. 2 .
Fig. 2. Position of the 22 TDR probes with respect to the infiltration trench; the probes used in the calibration phase are shown as white circles; boundary conditions are indicated as arrows.

Fig. 3 .
Fig. 3. Observed and fitted cumulative infiltration values measured at the surface for different (a) pressure heads and (b) their residuals.

Fig. 4 .
Fig. 4. Cumulative infiltration curves from simulated rainfall measurements and the fitted Philip equation (1957); Different symbols indicate different replicates.

Fig. 5 .
Fig. 5. Boxplots showing the mean K sat values, their 25% and 75% percentiles and the outliers observed with four different measurement techniques: a ring pressure and well infiltrometer, using single head analysis (PI-SHA) or multiple head analysis (PI-MHA), a Tension Infiltrometer (TI) and Rainfall Simulations (RS).

Fig. 6 .
Fig. 6.Water retention data and fitted van Genuchten equations for four depths.

Fig. 7 .
Fig. 7. Measured water retention data for all layers and the fitted van Genuchten equation compared to the van Genuchten equations obtained by minimizing each of the selected objective functions.

Fig. 8 .
Fig.8.Observed and fitted water content data for 10 selected TDR probes used in the optimization process; the line of perfect agreement is also indicated.

Fig. 9 .
Fig. 9. Regression line between observed soil-water contents versus simulated values at the same time step for (A) the 12 TDR probes not used in the optimization process and (B) all 22 TDR probes; the line of perfect agreement is also indicated.

Fig. 10 .
Fig. 10.Examples of observed and simulated water content values in function of time, for the regressions with (A) the highest and (B) the lowest accuracy.

Table 1 .
Physical characteristics of the infiltration trenches at the Quebrada de Talca experimental site.
a values preceded by ± give the standard deviation

Table 2 .
Physical and chemical characteristics of the soil profile at the experimental site.

Table 3 .
Comparison of K sat N is the number of measurements, K satGM is the geometric mean of saturated hydraulic conductivity K sat , K satmin and K satmax , respectively, are minimum and maximum K sat value observed in the data set b Geometric Mean values within the same column followed by the same symbol are not significantly different at P <0.1 a c Coefficient of variation

Table 4 .
Optimized K sat for the 12 tension infiltrometer experiments at three locations and four depths a .
a K sat is saturated hydraulic conductivity, K satlower and K satupper refer to the boundaries on the 95% confidence interval, RMSE is root mean squared error, IA is index of agreement and MBE is mass balance error

Table 5 .
Parameters from the Philip equation fitted to the cumulative runoff data from four rainfall experiments a .

Table 7 .
Initial and optimized parameter sets for six different objective functions a .