Modeling the unsaturated flow associated with a border irrigation event on an alfalfa plot

HYDRUS-1D, a one dimensional water movement and solute transport model, was applied to analyze a border irrigation event in an alfalfa plot located in the irrigation district of "La Violada" (Spain). The research encompassed monitoring of the irrigation event, the setting-up of the numerical model, model calibration and evaluation, and sensitivity analysis of the van Genuchten-Mualem parameters. Both the model predictions and field measurements clearly indicated that the soil water content was high during and following the irrigation event. The model even predicted a small evapotranspiration decrease induced by water logging and a significant contribution of capillary rise subsequently the water table rise in the adjacent plot caused by an irrigation event prior to the monitored event in the study plot. The research revealed the need for (i) the improvement of the schematization of the soil profile and the better experimental definition of key model parameters; and (ii) the derivation of more accurate data series for respectively the upper and lower boundary conditions.


INTRODUCTION
Irrigation brings along many benefits, however, the need of producing enough food for the growing world population and industry-related agricultural products such as bio-fuels, leads to an intensification of land and water use and an increase in the application of agro-chemicals with the consequent imbalance between water demand and availability, as well as the contamination of the environment.Mitigation of the negative effects of irrigation requires, among others, the introduction of suitable management strategies that improve water use efficiency and minimize the negative environmental impacts (Ragab, 2002;Causapé et al., 2003).The latter requires understanding of the main processes governing soil-water-plant interactions as a function of the time and space variable boundary conditions (Behera et al., 2003).Basic understanding of the soil-water-plant interactions will help irrigators to efficiently manage their irrigation systems.In this respect, numerical models are very valuable tools for better interpreting field data and deriving knowledge about the main soil-plant-water interaction processes, and to better underpin and guide decision-making (Ragab, 2002;Vázquez, 2003;Šimůnek et al., 2005;Nathan et al., 2005).
In the semi-arid Mediterranean basin of the Ebro River traditional irrigation systems such as border irrigation are used.The concept of efficiency is not taken into account by the farmers who receive water on a fixed schedule regardless the actual soil water status and crop growth stage.This often results in over irrigation and poor crop response (Isidoro et al., 2004), rise of the water table and leaching of nutrients with the consequent contamination of soils, surface water and groundwater (Isidoro et al., 2006).To assess the efficiency of the current irrigation practice and to use the interpreted information as guide for improving the system operation an irrigation event was monitored and field data analyzed using a numerical model.Emphasis was given to analyze the soil water dynamics and evaluate the model sensitivity as a function of the main model parameters.
Given the particular properties of the lowlands on which the numerical modeling experiment took place, as a first step in the study, the ability of the numerical model to simulate the soil water dynamics during and after the irrigation as a function of the soil properties and the boundary conditions was examined.Additionally the sensitivity of the model performance to the individual parameters was explored using a combined deterministic-stochastic approach (Beven and Binley, 1992;Vázquez et al., 2009).
As the success of any modeling activity depends on the availability of sufficient and accurate data (Stephens, 1996;Šimůnek et al., 2005;Vázquez et al., 2008), a further objective of this research was establishing a series of simple and cheap approaches for the monitoring of the soil-plant-water interactions in the study zone as to gather information for the applied numerical model.One final objective was to adapt the collected information over the whole plot to the input requirements of a one-dimensional model.The applied monitoring-modeling experiment differs from similar ones (Nathan et al., 2005) in that the study site, far from being a research field, is more representative of local irrigation practices since it is a plot run by a local farmer.

The study site
The study site is an irrigation plot located within "La Violada" irrigation district, in the north-eastern province of Huesca (Aragón) of Spain.The irrigation district is situated in the "La Violada" gully watershed.The elevation in the district ranges from 345 to 414 meters above sea level (m.a.s.l).The drainage network of the irrigation district discharges directly to "La Violada" gully, in several locations upstream the main outlet.The district is limited by the "Monegros" canal in the northeast, "La Violada" canal in the west and the "Santa Quiteria" canal in the south.The three canals embrace a surface of 5282 ha, but only about 4000 ha are irrigated.
The zone is characterized by a dry sub-humid mesothermic climate, with precipitation concentrated in spring and autumn.The mean annual precipitation is 479,9 mm, with 76 days of rainfall per year and the mean annual temperature is 13,3 ºC (data for the period 1964-1997;Isidoro, 1999).Border irrigation is dominant in the district, with most plots being blocked-end but with some, as the study site, having free drainage at the end.The irrigation system was initially designed to provide water to winter grains, wheat and barley, but currently it irrigates crops with higher water demands such as maize and alfalfa, which is nowadays the main crop in the district.In the irrigation district, alfalfa can be sown in November with the winter cereal, in March/April in plots already sown with wheat or barley, or in September in fallow plots.The habitual sowing dose is 70 to 80 kg ha -1 .No fertilizer is applied if the alfalfa is sown along with cereal; whereas single sown plots are fertilized with doses around 300 kg ha -1 of low Nitrogen (N), high Phosphorous (P) and Potassium (K) fertilizers (i.e. 8-15-15 or 5-20-20 complexes).The average number of cuttings of the perennial crop varies from five to six per year; usually only 3 to 4 in the first year.Generally, the first cut is carried out by the end of April and the last in mid-September.The cutting interval is approximately one month.Dry matter yield (18% of humidity) is about 14000 kg ha -1 .
For this study, alfalfa cropped plots were identified around the town of "Almudévar", in the catchment lowlands along the main drains where the soils, known as "Vales", have developed from alluvial parent material.Soils are deep (> 120 cm), normally stone-free, silty texture, and within some locations low internal drainage and high salinity (Torres, 1983).Two neighboring plots, F1 and F2 (Fig. 1) were selected for monitoring.Plot F2, 8517 m 2 in size (about 170 m x 50 m), was selected for testing the numerical model because (1) no significant contributions could be observed from other overland sources than irrigation, and (2) the water inlet and outlet could easily be monitored.Figure 1 shows the main geometric characteristics of both fields, and the location of the main water inlet and drainage features.To the east of F2, from north to south runs the "Valsalada" gully, one of the two branches of the "La Violada" gully.The different plots in the zone are irrigated successively from west to east, beginning with the higher ones and ending in the last plot immediately to the east of F2.Irrigation inlets are located in the northern edge of the plots.The plots are considerably longer than wide, with a mean slope along the longitudinal direction of 0,1%.
The advancement of the water flow in the border is controlled by the inflow rate, the smoothness of the surface, and the soil infiltration characteristic.Fields are separated by raised earth check banks preventing water to flow to neighboring fields.The surface outlet of field F2, "D2" (Fig. 1) is connected to a drain pipe situated 150 cm below the surface, running west-to-east.A similar drain pipe is located along the northern edge of the plot intercepting the groundwater flow from the neighboring plots in the north.All drains, respectively at the northern and southern border of the fields, discharge in the "Valsalada" gully.Given the configuration of underground drain pipes, it is likely that the water table height in the study plot F2 is only affected by groundwater from the neighboring plot at the western side, plot F1.

Characterization of the study site and monitoring of the irrigation event
The field survey included a topographic mapping and a soil prospecting campaign.The characterization of soil profiles was based on the analysis of the vertical variation of soil color, structure and texture of the subsequent soil layers and results from laboratory tests.For the textural characterization, soil samples were collected in each of the observed soil layers: 0-30, 30-55, 55-65, 65-90, 90-115 and 115-125 cm.A less permeable soil layer was identified at 1,35 m below surface.
For the monitoring of the water table build-up during each irrigation event on top of the less permeable soil layer, 5 shallow observation wells ("w 1 " to "w 5 " in Fig. 1) with a depth of 1,2 m were installed at plots F1 and F2.The temporary water table depth was measured manually.To examine the fluctuation of the piezometric level below the limiting substratum a deeper piezometer ("p 1 " in Fig. 1) was installed.For measuring the runoff leaving the field, a square notch weir/orifice was installed just before the horizontal orifice connects to the drainage pipe in the lower end of the field (Fig. 1).The weir/orifice system consisted of a 5x5 cm 2 square orifice bored in a vertical metallic plate that works as a sharp crested weir when the water level is between the lower and higher edges of the orifice (i.e.throughout a vertical distance of 5 cm) and as a submerged orifice for water levels that are higher than the upper edge of the square orifice.The lower edge of the orifice was aligned with the neighboring ground level of the study plot F2 but the 1,5 x 1,5 m 2 upstream land swath was lowered a few centimeters so as to ensure that the approaching water fulfilled the hydraulic conditions for the discharge measuring structure.A water height recorder was placed in situ for the continuous monitoring of the water level on the discharge measuring system.
The studied irrigation event took place on field F2 on the 15 th of September 2007.The total volume of water infiltrated was calculated as the difference between the in-and out-flow volumes.The inflow rate, Q a , in the semi-circular irrigation channel was frequently measured using a current meter.Measurements revealed that the inflow rate was nearly constant throughout the monitored irrigation event.
The time advancement of the overland water sheet during the irrigation event was monitored with the aid of a precision GPS, and complemented using the Surface Irrigation Simulation Model (SIRMOD; Walker, 1993).SIRMOD is generally used to simulate, evaluate and design surface irrigation systems.Due to the lack of sufficient data on the measured mean water height during the irrigation event, water height during the event was reconstructed using SIRMOD.Inputs to SIRMOD included the surface slope, the inflow rate and duration of the application.The SIRMOD model was calibrated by varying the parameters of the empirical Kostiakov infiltration equation (Chow et al., 1988) until the predicted runoff volume and duration matched approximately the respective observations.
Five soil sampling campaigns were carried out, respectively one campaign before the irrigation event and the remaining 4 campaigns after the event, to determine gravimetrically the soil water content (θ wh , M M -1 ) at four depth intervals: 0-30, 30-60, 60-90 and 90-120 cm.The volumetric water content, θ (L 3 L -3 ), was derived multiplying θ wh with the mean soil bulk density (θ b, M L -3 ) of the corresponding soil layer.The time-variant mean water table levels before, during and after irrigation were derived from the water level records in the observation wells in the study fields F1 and F2 (Fig. 1).The water level hydrograph, WT tot , is the result of the combined contribution of both, the irrigation event that took place at the upstream neighboring plot F1 (i.e.WT F1 ), and the posterior irrigation event on F2 (i.e.WT F2 ).WT F1 was separated from WT tot to define the lower boundary condition (BC) of the numerical model of the study site, by applying the procedure that is depicted in Fig. 2a.Hereby, the recession limb of WT tot was approximated by an exponential curve, assuming that the groundwater storage in the area is linear (discharge proportional to the stored volume).This exponential curve was in turn adopted for approximating the corresponding recession limb of WT F1 (Fig. 2a).The rising limb of WT F1 is the observed rising limb until approximately the end of the irrigation event on F1 (Fig. 2a).
The lower BC was estimated then as the rise of the water table in the study site with regard to the mean water table level after irrigation, z m (i.e.WT F1 -z m ).
The reference evapotranspiration (ET o ) was calculated using the meteorological data from the Tardienta station, belonging to the SIAR network of the Spanish Ministry of Agriculture, situated at a distance of 6 km from the study site.ET o was calculated using the Penman-Monteith equation (Allen et al., 1998) using 30-minute interval measurements of air temperature, relative humidity, solar radiation and wind speed.Crop evapotranspiration (ET c ) was calculated as the product of K c and ET o , with K c being the daily derived crop coefficient according to Allen et al. (1998), taking into account the days past since the last cutting in the alfalfa crop (Sayah, 2008).

Determination of soil hydraulic parameters
The HYDRUS-1D model used to describe the water dynamics in the unsaturated zone is based on the Richards equation (Stephens, 1996;Šimůnek et al., 2008), for which, owing to its non-linearity, only few simplified analytical solutions exist and most of its practical applications involve a numerical approximation.The numerical solution of the Richards equation requires the availability of the hydraulic functions describing the soil water retention and the soil hydraulic conductivity The analytical expression of the soil water retention function necessitates the knowledge of parameters such as the saturated soil water content (θ s ) [L 3 L -3 ] and preferably some measurements of the relationship between θ [L 3 L -3 ] and the pressure head (h) [L; or pF = log10(h), with h in cm].On the other hand, the analytical representation of the hydraulic conductivity curve demands the knowledge of at least θ s and the saturated hydraulic conductivity (k s ).
In this study only the pF-range 2,0-4,2 log10 (cm) (i.e.0,1-15,5 bar) was measured using the ceramic plate extractor.The soil water content was determined for the following pressure values: 0,1, 0,3, 1,0, 3,0, 5,0, 10,0, and 15,5 bar.Soil samples were collected at five locations within the study plot F2 and from the following depth ranges: 0-30, 30-60, 60-90 and 90-120 cm.The Inverse Auger Hole method (van Beers, 1976;Oosterbaan and Nijland, 1994;Vázquez, 2003) was used to estimate the average k s of the soil profile above the water table.Tests were carried out at 7 different locations within plot F2 and at two different depths so that average values were obtained for the soil columns in the depth intervals: 0-55 and 55-120 cm.
The parameters of the water retention curve (i.e. the parameters of the van Genuchten-Mualem model) θ r , θ s , k s , the tortuosity (l) and the curve shape parameters n [--] and α [cm -1 ], were estimated by means of the following protocol that was implemented as well to cross check on the congruency of the experimental data and the model predictions.A first step consisted on the application of the artificial neural network computer code ROSETTA (Schaap et al., 2001) using as input soil textural information (i.e.percent sand, silt, and clay), θ b and water content at field capacity (θ FC ) and at wilting point (θ WP ).
A second step was based on the use of the curve-fitting/parameter-optimization code RETC (van Genuchten et al., 1991) that fits the parametric models of Brooks and Corey (1964) and van Genuchten to actual soil water retention data (Hollenbeck et al., 1999).Because, only few experimental pF data were available, which affects the statistical degrees of freedom, RETC was used to refine the values of only the parameters θ s , n and θ, while still using the values of θ r and k s determined through the previous approach.Finally, a third step involved the calibration of the numerical model of the study site by manually varying the values of all of the van Genuchten-Mualem parameters.

Model calibration and evaluation
The HYDRUS-1D modeling system used to simulate the water movement in the study site is a wellknown one-dimensional deterministic model that has been widely described and applied (Inoue et al., 2000;Bumgarner and McCray, 2007;Heatwole and McCray, 2007;Leij et al., 2007;among others).It was developed to simulate the soil water and solute dynamics in the root zone, root growth, root water and solute uptake, carbon dioxide transport and heat transport (Šimůnek et al., 2005).HYDRUS-1D uses the Richards equation for saturated-unsaturated water flow and convection-dispersion type equations for heat and solute transport.The joint deterministic-stochastic modeling protocol with HYDRUS-1D included three steps: calibration, evaluation and a posterior sensitivity analysis (SA).Model calibration was based on a step by step manual trial and error variation of the van Genuchten-Mualem parameters by considering both visual inspection of the agreement among the observed and predicted variables and evaluation of a model performance index, the Coefficient of Efficiency (Nash and Sutcliffe, 1970;Legates and McCabe, 1999): where 2 obs = the observed variance and 2 i = the error variance for the model.Although this index is better tailored for model evaluation than other correlation-based indexes such as the Coefficient of Determination (usually reported in literature as either "r 2 " or "R 2 "; Vázquez et al., 2002), it is still affected by some oversensitivity to the simulation of peak values (Legates and McCabe, 1999;Vázquez et al., 2002).Nevertheless, it was selected for the current study as it measures both the average systematic and random discrepancies among simulations and observations (Vázquez et al., 2008).
As only one irrigation event was monitored, model calibration was carried out upon the time series of the total water content in the profile (θ T ), whilst model evaluation and the posterior SA were implemented using the depth distributed time series of θ (4 depths) as well as the time series of observed water table depth.θ T was calculated by integrating and averaging the vertical distribution of θ measured in 3 observation points in plot F2 in each of the 5 sampling campaigns.
The feasible, i.e. physically based, van Genuchten-Mualem parameter intervals considered during model calibration and the SA are depicted in Table 1 and were derived from literature sources (Yates et al., 1989;Bohne et al., 1993).The protocol used for the SA was based on Monte Carlo simulations (Beven and Binley, 1992).A total of 25000 sets of parameter values were sampled randomly from the multidimensional parameter space, considering uniform parameter distributions for all of the parameters (Vázquez, 2003;Vázquez et al., 2009).The same van Genuchten-Mualem parameters that were included in the trial and error model calibration were included in the SA.The sampling parameter intervals were the same as those used for the model calibration (Table 1).The model performance statistic that was used throughout the SA is the EF 2 , calculated on the predictions of θ, distributed vertically for every observation time.
Table 1.Feasible variation intervals for the van Genuchten-Mualem parameters.
The sensitivity of the model performance to the individual parameters was explicitly explored using a form of Global Sensitivity Analysis (GSA; Freer et al., 1996;Vázquez et al., 2009).In this method the EF 2 -based weights of the behavioral simulations are grouped into intervals.A cumulative distribution of parameter values is then calculated for each one of these intervals.The criteria used for assessing the relative sensitivity of the individual parameters are: (i) strong differences among the cumulative distributions of behavioral values indicate a sensitive parameter; whilst (ii) similarities among the cumulative distributions suggest an insensitive parameter (Spear and Hornberger, 1980;Freer et al., 1996).Given the extensive computational work required in the SA, HYDRUS-1D was run in batch mode.Furthermore, the subroutines needed to sample the 25000 parameter sets, to prepare the respective simulation files, to process the HYDRUS-1D results and to calculate the EF 2 values were written in PERL (Practical Extraction and Reporting Language), following a protocol outlined in Vázquez (2003).

RESULTS AND DISCUSSION
The field survey revealed the presence of alfalfa roots to a depth of 90 cm.Soil color showed a prominent change at a depth of 55 cm.The oxidation of iron (Fe)/magnesium (Mn) was apparent below this depth (reddish-black dots), caused by the alternate saturated/unsaturated conditions.Figure 3a shows the vertical variation of the sand, clay and silt fractions, and Fig. 3b the variation of the soil organic matter (SOM).There is a significant vertical variation of clay (decrease with depth from about 42% to about 7%) and silt (increase with depth from about 4% to about 18%) fractions around 55 cm depth.The A horizon (0-55 cm) classifies as a silty clay, while the B 1 (55-105 cm) and B 2 (105-135 cm) horizons as silt loam.The vertical variation of the SOM content also marks the differentiation between the A and B horizons.
The abrupt textural change at a depth of 55 cm led to the conceptualization of the modeling domain composed by two soil layers, i.e. an upper 55 cm thick silt clay layer ("layer 1") and a deeper silt loam horizon reaching a depth of 135 cm ("layer 2").This schematization of the soil profile was confirmed by the fact that initial HYDRUS-1D simulations considering the soil profile of the study area consisting of two horizons (A and B) performed equally well as the case whereby the profile was mimicked by a sequence of three layers (A, B 1 and B 2 ).The average value and the range of the saturated hydraulic conductivity (k s ) and the bulk density (θ b ) are listed in Table 2.
Table 2. Experimental values for the hydraulic conductivity (k s ) and bulk density (θ b ).n t is the number of replicate measurements for k s ; SD stands for the standard deviation.θ b is the average of two measurements per soil horizon.

Depth range
[cm] The mean values of both soil parameters, although only slightly different, reflect the composition difference between the two soil layers.Both the range and the coefficient of variation (CV) are higher for the lower layer ("layer 2").It is noticed that the measured values are different from the values typically suggested in literature for similar soil textures (Anderson and Woessner, 1992;Yates et al., 1989;Bohne et al., 1993).
Figure 2b shows the time series used for the upper and lower boundary conditions.The modeled irrigation event in F2 started on 15 of September of 2007 at 16:12 hr and lasted 4,6 hr.The total irrigation water volume applied was 1125 m 3 at a rate of 0,068 m 3 s -1 .The figure depicts the time series predicted with SIRMOD.The plot-end drainage started 17 minutes before the end of the irrigation event and lasted about 18,3 hr.Approximately 1,67 cm (12,3% of the applied irrigation volume, 13,67 cm) was lost as runoff, and 11,91 cm (87,7%) of the applied irrigation water infiltrated.Although, these results are important for the assessment of the irrigation application efficiency, they cannot be extrapolated directly to the rest of the irrigation district, since most of the plots located in the "La Violada" district do not present a surface drainage outlet at the plot-end.
Figure 2a shows, on the other hand, the average water table record at the study site F2 as well as the instants of the beginning and end of the irrigation events in both F1 and F2.The WT level in the absence of irrigation, z m , was measured at 133 cm below surface.Both, the estimated z m and the measured average WT tot evolution are confirmed by records obtained prior to the installation of the runoff gauging system in F2.The figure shows the clear rise of the WT tot hydrograph prior the beginning of the irrigation event in F2, due to the irrigation event in F1.The effect of the irrigation in F2 on the WT tot hydrograph was removed to obtain the transient lower BC applying the approach depicted in Fig. 2a.
The values of the van Genuchten-Mualem parameters adopted throughout this study along the parameterization process are shown in Table 3, along with the respective model performance assessment which was carried out by calculating the EF 2 -values upon the simulation of the average volumetric water content in the soil profile (θ T ) [L 3 L -3 ].The HYDRUS-1D parameterizations, considering the estimates generated with the codes ROSETTA and RETC, produced non acceptable model performances (i.e.EF 2 < 0,0).However, the calibrated parameter set produced a very good EF 2value.The calibrated k s was very close to the experimental k s for the lower layer, but the measured k s in upper layer was about five times higher than the calibrated value.This overestimation could be the result from the fact that the auger hole used for determining k s of the upper layer (layer 1) was too deep; stretched beyond the bottom boundary of the horizon into the underlying horizon (layer 2) with a higher value of the saturated hydraulic conductivity.As a consequence, it is very likely that during the test a considerable fraction of the normally horizontal flow took place through the bottom of the auger hole, leading to an overestimated k s value.Table 3. Experimental, calculated and calibrated values for the hydraulic conductivity (k s ) and the respective model parameters.
Figure 4 shows the vertical distribution of the observed and simulated (after calibration) volumetric soil water content at four subsequent observations during the recession phase following the irrigation.Although good simulations of the averaged θ T were obtained, the observed and simulated vertical distributions of θ started to deviate increasingly with time, particularly in the lower horizon.Figure 5 presents the result of the model evaluation analysis: the vertical distribution of the simulated soil water content (θ) is shown in conjunction with the observed average WT tot (horizontal dashed line) and the average θ s of the B horizon (vertical gray line) for the three available observation times.The depth to the WT, approximated as the intersection of the average θ s and the simulated vertical distribution of θ, was well predicted in the beginning and end phase of the simulation period that is, before and after the flooding of plot F2.Both, the field measurements and the model predictions revealed that θ was in general high throughout the monitoring period (17 days); in every sampled layer larger than 57% of θ s , suggesting that the alfalfa crop in this period was not affected by water stress.One day before the start of the irrigation event in F2 an increase in θ was observed in all layers, which given the absence of surface inputs (rainfall and irrigation) suggests a significant capillary rise contribution, in average terms estimated at 1,4 cm d -1 , and apparently due to the precedent irrigation event in the upper neighboring plot F1 (Fig. 1).
Figure 2a presents the cumulative time series of the potential (ET c ) and actual (ET a ) evapotranspiration rate, the potential (Q i_pot ) and actual (Q i_act ) infiltration rate, and the soil drainage rate (Q soil ), all expressed in cm hr -1 .Furthermore, the figure depicts the average measured (WT tot ) and estimated (WT F1 ) water table depth of plot F2 and F1, and the beginning of the irrigation event on plot F1 and F2, respectively.The figure shows a difference between Q i_pot and Q i_act during the advance and pounding phase of the irrigation event.As already stated, apparently the soil hydraulic properties adopted for the upper most soil layer, avoided the pounded water to infiltrate in that period (defined by the duration of the observed runoff drainage).HYDRUS-1D also predicts a bottom upwards flow into the root zone (modeling domain) before the irrigation in F2 and especially during the irrigation event in F1, up to a total of 31,8 mm (Fig. 2a).The net flux at the bottom of the root zone took place downwards during the irrigation of F2, remained almost nil during the period of inundation in F2 and decreased by the end of the infiltration period in F2 for a period of 7 days, to increase slightly thereafter (Fig. 2a).In the full modeled period, there was a net leaching of about 13 mm.The figure reveals a difference between the potential cumulative (ET c ) and the actual cumulative (ET a ) evapotranspiration with an average ratio of ET a /ET c = 0,92.The predicted reduction in ET occurred during the flooding phase of the irrigation event, resulting to anaerobic conditions in the root zone and a reduction in the root water uptake (Feddes et al., 1978).This high θ does not cause salt accumulation in the soil of the study plot F2 (as evidenced by the analysis of the EC of 1:5 soil-water solutions; data not shown herein), owing to the significant presence of drainage infrastructures in the plot (Fig. 1), but it might explain soil salinization problems taking place in other locations of the irrigation district that do not have comparable drainage facilities.With respect to the joint deterministic-stochastic sensitivity analysis, about 4600 parameterizations of the model structure (out of 25000) failed due to numeric instabilities.A brief visual examination of the results of the successful model runs indicated very few cases with systematic deviations of the model predictions relative to the observations.Figure 6 shows the scatter plots of the performance measure across the sampled parameter ranges.Every single dot in the scatter plots represents a sampled parameter set and as such a model run.The range of axes is the same as the sampled range.Since only positive values of the EF 2 index are meaningful, the parameter sets having an associated positive EF 2 value (i.e.EF 2 > 0,0) were retained for further analysis.The parameter set that was obtained through the manual trial and error calibration is highlighted in every plot of the figure by means of a square symbol.
Although there were not many behavioral parameter sets (i.e.only 2427 simulations with an associated EF 2 > 0,0), the scatter plots show that model performance was most sensitive to the saturated water content in the upper [(θ s ) upp ] and lower layers [(θ s ) low ], as their respective plots had a "peaked" cloud of dots approximately between 0,38 and 0,49 cm 3 cm -3 for (θ s ) upp and between 0,42 and 0,53 cm 3 cm -3 for (θ s ) low , both in the range of found experimental values (Table 1).The other scatter plots are flat topped across the whole parameter ranges, suggesting that model performance was not sensitive to their corresponding parameters.Figure 6 also shows that although the performance associated to the calibrated (manually) model is one of the highest; there are several parameter sets that produce model performances not only equivalent, but also superior to it throughout most of the extent of the sampled parameter intervals.
With regard to the GSA, in this research the behavioral simulations were grouped into 4 classes using a fixed width of the EF 2 -intervals (i.e.0,25).Only 4 classes were considered owing to the very limited number of behavioral simulations (i.e.only 2427 out of 25000 simulations).Class 1 includes the highest likelihood values, whilst, class 4 includes the lowest likelihood values (still higher or equal than the cut-off value 0,0).The analysis revealed that most of the parameters are insensitive because of the similarity between their associated cumulative distributions.However, the parameters (θ s ) upp and (θ s ) low seem to be sensitive.These suggest that the parameter response depicted in Fig. 6 indeed seems to be a consequence of the high insensitivity of the model performance to most of the inspected parameters, in addition to plausible parameter interactions, errors in boundary conditions and input data and model structure uncertainties.

CONCLUSIONS, RECOMMENDATIONS AND FUTURE WORK
A basin irrigation event on a plot located in the central part of the Ebro River basin was modeled with the one dimensional soil water model HYDRUS-1D.With respect to the characterization of the soil properties, the study revealed that the Inverse Auger Hole method yielded adequate, congruent estimates of k s .However, a higher density of tests than the one considered in this study is highly desirable for future analyses.The experimental h-θ data collected to fit the water retention curve were also congruent, but the h-range that could be used because of the available equipment (2,5 < pF < 4,2) was most probably not wide enough.
The model was capable of simulating the phreatic level and the soil water content in the soil profile only when the van Genuchten-Mualem parameter values were assessed through model calibration.Some of these calibrated parameter values (for θ s and k s ) differed significantly from the experimental ones.This is likely the consequence of the inappropriate definition of initial/boundary conditions besides the mismatch between the scale at which the differential equations that are included in the structure of HYDRUS-1D were derived, the scale at which the experimental data were collected and the modeling scale (Vázquez, 2003).In this context, the process followed to approximate the one dimensional variation of the water sheet level on top of the soil surface using the predictions from another numerical model (i.e.SIRMOD), likely introduced additional model uncertainties, not explicitly addressed in this study.In addition, the simple averaging approaches used in this research to derive one dimensional (1D) data sets likely also introduced some numerical inaccuracies and as such might need to be improved in future.
Furthermore, the worse predictions obtained for the lower soil layer and the uncertainty attached to the approach, used to estimate the groundwater contribution from the previous irrigation event in plot F1, suggest considering the root zone as the modeling domain rather than the whole vadose zone.Indeed, the appropriate definition of the lower boundary conditions seems to be important.The fact that the WT level was not always correctly predicted might be as well due to air entrapment and preferential flow.Those phenomena can only be captured through more complex modeling, e.g.dual porosity models.Notwithstanding, model performance assessment suggests that the model acceptably captured the overall hydrological dynamics of water movement during the simulated irrigation event.Hence, both the observed records as well as the model predictions indicate that the basin irrigated plots in the lowland of La Violada district present high water content during the irrigation season, apparently resulting from lateral flows and capillary rise induced by prior irrigation events in uphill vicinity plots.The apparent contribution of capillary rise is certainly not to be neglected and could reduce the need for irrigation.In this context, the model even predicted a small evapotranspiration decrease induced by water logging.
The sensibility analysis revealed the insensitivity of the model performance to most of the van Genuchten-Mualem parameters, except for the saturated soil water content.The analysis showed that comparable or even better model performances were obtained for parameter values that are significantly different from the calibrated ones.That is, even when using the parameter values that were obtained experimentally for k s and θ s , it is likely that there exists at least one parameter set, i.e. a combination of certain parameter values, that would lead to model predictions that are comparable to the ones obtained through the implemented trial an error calibration.This fact could be useful for refining the calibration approach, by leaving the experimental parameter values out of the calibration process.The equifinality issue identified through the SA implies not only that there is not a clear single optimal parameter set but also that the prediction limits associated with the current modeling might be significantly wide.
Future activities could include not only the refinement of the model domain and initial/boundary conditions, but also the application of a two dimensional model, such as HYDRUS-2D (Šimůnek et al., 2008).Since the problem treated herein is typically bi-or three-dimensional.The modeling exercise also suggests the need of (i) generating better experimental parameter values (such as water retention curve data, k s data, etc.), and (ii) gathering better (i.e. less uncertain) time series records.Herein, the use of soil water content automatic measuring instrumentation (e.g.four electrode probes) could help establishing more correct θ distributions and consequently could help achieving a more appropriate evaluation of the HYDRUS-1D predictions.The final aim of the research, of which the current study is only a part, is to apply at catchment scale hydraulic/water quality modeling and to integrate in such a modeling protocol the knowledge acquired at plot scale, through studies such as the present one.

Figure 1 .
Figure 1.(a) Layout of the main irrigation and drainage structures and monitoring devices in the study plot (after Sayah, 2008); and (b) detail of the square notch weir/orifice that was installed on the field for measuring the runoff leaving the irrigation plot.

Figure 2 .
Figure 2. (a) Cumulative time series predicted by HYDRUS-1D for the study site: soil drainage rate (Q soil ), crop potential (ET c ) and crop actual, modeled, (ET a ) evapotranspiration; potential (Q i_pot ) and actual, modeled, (Q i_act ) infiltration into the soil.Downward flows are negative whilst upward flows are positive (HYDRUS-1D signs convention).The measured water table depth (WT tot ) and the estimated recession of the water table for F1 (WT F1 ) are also shown; and (b) time series used for the upper model boundary conditions (after Sayah, 2008).

Figure 3 .
Figure 3. Vertical distributions of (a) the soil textural properties; and (b) the soil organic matter (SOM).

Figure 4 .
Figure 4. Vertical distributions of the observed (θ obs ) and simulated (θ sim ) volumetric soil water content at increasing time periods following the irrigation event.

Figure 5 .
Figure 5. Vertical distribution of the volumetric soil water content predicted by HYDRUS-1D (black continuous line), the calibrated saturated volumetric soil water content (θ s = 43%) (gray vertical line), and the position of the observed phreatic water level (dashed horizontal line).

Figure 6 .
Figure 6.Scatter plots of the model performance statistic EF 2 as a function of the parameter ranges considered in the sensitivity analysis for the van Genuchten-Mualem parameters.EF 2 was calculated by considering the vertical distribution of the volumetric water content (θ) along the soil profile [L 3 L - 3 ] as a function of time.Manually calibrated parameter values are highlighted in the plots through a square symbol.