Sensitivity of water stress in a two-layered sandy grassland soil to variations in groundwater depth and soil hydraulic parameters

M. Rezaei, P. Seuntjens, I. Joris, W. Boënne, S. Van Hoey, P. Campling, and W. M. Cornelis Department of Soil Management, Ghent University, Coupure Links 653, 9000 Ghent, Belgium Unit Environmental Modeling, Flemish Institute for Technological Research (VITO NV), Boeretang 200, 2400 Mol, Belgium Department of Bioscience Engineering, University of Antwerp, Groenenborgerlaan 171, 2020 Antwerp, Belgium Department of Mathematical Modelling, Statistics and Bioinformatics, Ghent University, Coupure Links 653, 9000 Ghent, Belgium


Introduction
Efficient water use and optimal water supply to increase food and fodder productivity are of great importance when confronted with worldwide water scarcity, climate change, growing populations and increasing water demands (FAO, 2011).In this respect, irrigation efficiency which is influenced by the type of irrigation and irrigation scheduling is essential for achieving higher water productivity.In particular, precision irrigation is adopting new methods of accurate irrigation scheduling (Jones, 2004).Various irrigation scheduling approaches such as soil-based, weather-based, crop-based, and canopy temperature-based methods have been presented (Jones, 2004;Mohanty et al., 2013;Pardossi et al., 2009;Evett et al., 2008;Nosetto et al., 2012;Huo et al., 2012).
Numerical models are increasingly adopted in water resource planning and management.They contain numerical solutions of the Richards' equation (Richards, 1931) for water flow and root water uptake (Fernández-Gálvez et al., 2006;Vrugt et al., 2001;Skaggs et al., 2006) or contain reservoir cascade schemes (Gandolfi et al., 2006).Hydrological models require determination of hydraulic properties (Šimůnek and Hopmans, 2002), upper boundary conditions related to atmospheric forcing (evapotranspiration and pre-Published by Copernicus Publications on behalf of the European Geosciences Union.cipitation) (Brutsaert, 2005;Nosetto et al., 2012) and groundwater dynamics at the lower boundary of the soil profile (Gandolfi et al., 2006).Numerical models such as Hydrus 1D (Šimůnek et al., 2013) have been used in a wide range of irrigation management applications -for example, by Sadeghi and Jones (2012), Tafteh and Sepaskhah (2012), Akhtar et al. (2013) and Satchithanantham et al. (2014).The tool has been combined with crop-based models for accurate irrigation purposes and for predicting the crop productivity for cotton (Akhtar et al., 2013), vegetables and winter wheat (Awan et al., 2012).The degree of soil water stress was used for irrigation management by coupling a hydrological model (Hydrus 1D) with a crop growth model (WOFOST) for maize (Li et al., 2012) and wheat (Zhou et al., 2012).The importance of correct average representation of the soilplant-atmosphere interaction in numerical models has been stressed by Wollschlager et al. (2009).A combination of crop growth model and the hydrological model makes it possible to calculate crop yield reduction based on soil water stress derived by the hydrological model.
Direct measurement of hydraulic parameters may be inaccurate for predictions at the field scale (Verbist et al., 2012;Wöhling et al., 2008).As an alternative, parameters can be determined by inverse modelling.A single-objective inverse parameter estimation using the Levenberg-Marquardt optimization procedures has been used in different studies (Abbasi et al., 2004;Jacques et al., 2012;Šimůnek et al., 2013).A typical challenge in parameter optimization is the nonuniqueness of the parameters, leading to parameter identifiability problems (Hopmans et al., 2002).Non-uniqueness can be reduced by decreasing the number of parameters to be estimated based on a sensitivity analysis.Sensitivity analysis has been used to optimize parameter estimation, to reduce parameter uncertainty (Rocha et al., 2006), and to investigate the effects of various parameters or processes on water flow and transport (van Genuchten et al., 2012).
In this study, we used a combination of soil moisture monitoring and modelling to estimate hydraulic properties and to predict soil water content in a two-layered sandy soil for precision irrigation management purposes.The objective of this paper is to investigate the impact of parameter estimation and boundary conditions on the irrigation requirements, calculated using a soil hydrological model in combination with a crop growth model.The effect of changing bottom boundary conditions on model performance was evaluated in a first step (see the Supplement).A systematic local sensitivity analysis was then used to identify dominant hydraulic model parameters.This was followed by a model calibration using inverse modelling with field data to estimate the hydraulic properties.Finally, the degree of soil water stress was calculated with different parametrization scenarios to show to what extent hydrological model parameter choice and boundary conditions affect estimations of irrigation requirement and crop yield.It is acknowledged that there is no stress in soil water, whereas the water stress is in the plant, indeed.But similar to a large bulk of papers and reports, we used the soil water stress term in the present paper instead of water stress in the plants.

Description of the study site
The study site is located in a sandy agricultural area at the border between Belgium and the Netherlands (with central coordinates 51 • 19 05 • N, 05 • 10 40 • E), characterized by a temperate maritime climate with mild winters and cool summers.During the study period 2011-2013, the farmer cultivated grass.The farm is almost flat (less than 1 % sloping up from NW to SE) and runoff is not considered to be important.The measured depth of the groundwater table was between 80 and 155 cm and the Ap horizon thickness was between 30 and 50 cm below the soil surface at various locations across the field depending on the topography.The field is partly drained by parallel drainage pipes which are placed at 10 to 20 m intervals and at around 90 cm below the soil surface (as measured in the ditch).Drainage pipes are connected to a ditch in the northwest border of the field.Figure 1 shows the location and layout of the field.The apparent soil electrical conductivity, EC a , was measured at 5 m intervals between the measurement lines with a DUALEM-21S sensor (DUALEM, Milton, ON, Canada) corresponding to 0-100 cm depth of exploration.Then, EC a data were interpolated using ordinary point kriging (OK) to a 0.5 by 0.5 m grid to produce the field EC a map.More details about this methodology and its procedure can be found in Rezaei et al. (2016).Reel Sprinkler Gun irrigation (type Bauer rainstar E55, Röhren-und Pumpenwerk BAUER Ges.m.b.H., Austria) was used on a 290 m by 400 m field to improve crop growth in the sandy soil during dry periods in summer.The field was irrigated three times throughout each growing season (2012: 64.5 mm and 2013: 85.4

mm).
Figure 2 shows the soil profile at a sensors location, indicated by the star on the map in Fig. 1 (see also next section), a typical Podzol (Zcg-Zbg type according to the Belgian soil classification or cambisol according to WRB; FAO, 1998) consisting of a uniform dark brown layer of sandy soil (Ap horizon, 0 to 33 cm) with elevated organic matter content, followed by a yellowish to white sandy soil, including stones and gravels, (C1 horizon, 33 to 70 cm).A deeper horizon is light grey sandy soil (C2 horizon, 70 to 135 cm), including more stones and gravels (max 20 %), but having similar hydraulic properties as the C1 horizon.Maximum grass root density was found at about 6 cm and decreased from 6 to 33 cm (based on field observation during profile excavation).The properties of the two layers are summarized in Table 1.Table 1.Average of soil properties of soil profile: θ r , θ s are residual and saturated water content, respectively; α and n are shape parameters for the van Genuchten-Mualem equation.K s denotes the saturated hydraulic conductivity.

Field monitoring system
The site was equipped with two weather stations (type CM10, Campbell Scientific Inc., UT, USA), one in the study field and another 100 m away from the field.Soil water content was recorded (from 1 March until 25 November in both 2012 and 2013) using a water content profile probe (type EasyAG50, Sentek Technologies Ltd., Stepney, Australia), placed vertically, that measures soil water content at 10, 20, 30, 40 and 50 cm depths.The weather stations were connected to a CR800 data logger (Campbell Scientific Inc., UT, USA) and the water content profile probe provided the soil water content wirelessly.All measurements were taken on an hourly basis and an hourly reference evapotranspiration was calculated based on the Penman-Monteith equation (Allen et al., 1998) using weather station data.The amount of irrigation was derived by subtracting measurements of rain gauges of the field's weather station (i.e.rainfall and irrigation) and the local meteorological station (i.e.only rainfall) outside the study field.Grass yield was measured at each harvesting time (four times in each growing season) across the field (Fig. 3).At the sensor location (indicated by the star on the map in Fig. 1), duplicate undisturbed (100 cm 3 Kopecky rings, Eijkelkamp Agrisearch Equipment, Giesbeek, The Netherlands) soil samples were taken to determine the soil saturated hydraulic conductivity and water retention curve, and one disturbed sample to measure soil properties such as texture, dry bulk density and organic matter, from the Ap (topsoil) and C (subsoil) horizons in June 2013.Groundwater depth at the sensor location was measured four times on 4 June and 5 October 2012 (140 and 136 cm, respectively), and 24 June and 25 October 2013 (135 and 133 cm, respectively) using augering.
The saturated hydraulic conductivity (K s ) was determined using a constant head laboratory permeameter (M1-0902e, Eijkelkamp Agrisearch Equipment, Giesbeek, The Netherlands).The soil water retention curve (SWRC, θ (h)), was determined using the sandbox method (Eijkelkamp Agrisearch Equipment, Giesbeek, the Netherlands) up to a matric head of −100 cm and the standard pressure plate apparatus (Soilmoisture Equipment, Santa Barbara, CA, USA) for matric heads equal to or below −200 cm, following the procedure outlined in Cornelis et al. (2005).Bulk density was obtained by drying volumetric soil samples (100 cm 3 ) at 105 • C. Particle size distribution of the mineral component was obtained using the pipette method for clay and silt fractions and the sieving method for sand particles (Gee and Bauder, 1986).The organic matter content was determined by the method of Walkley and Black (1934).
Soil hydraulic properties were determined according to the van Genuchten (1980) and Mualem (1976) conductivity model (MVG model).The parameters of the water retention equation were fitted to the observed data set using the RETC, version 6.02 (van Genuchten et al., 1991).The MVG model (Mualem, 1976;van Genuchten, 1980) is given by where θ s , θ r and θ are the saturated, residual and actual volumetric water content respectively (cm 3 cm −3 ), α is the inverse of air entry value (cm −1 ), n is a pore size distribution index > 1, m = 1 − 1 / n (dimensionless) S e is the effective saturation (dimensionless), and l is a pore connectivity and tortuosity parameter in the hydraulic conductivity function, which is assumed to be 0.5 as an average for many soils (Mualem, 1976).

Simulation of leaf area index and grass yield
The simple generic crop growth model, LINGRA-N (Wolf, 2012), which can calculate grass growth and yields under potential (i.e.optimal), water-limited (i.e.rain fed) and nitrogen-limited growing conditions, was used to calculate the leaf area index (LAI) and grass yield.This tool was calibrated and tested for perennial rye grass and natural annual grass over Europe (Barrett et al., 2004;Schapendonk et al., 1998).LINGRA-N simulates the growth of a grass crop as a function of intercepted radiation, temperature, light use effi-ciency and available water (Wolf, 2012).The LAI and crop growth simulations were carried out from 1 January 2012 to 31 December 2013.The model calculated LAI and yield on a daily time intervals using daily weather data, solar radiation (kJ m −2 d −1 ), minimum temperature ( • C), maximum temperature ( • C), vapour pressure (kPa), wind speed (m s −1 ) and precipitation (mm d −1 ).A grass crop data file is available mainly derived from WOFOST.Soil data for our soil were produced using measured values of soil moisture content at air dry (pF = 6), wilting point (pF = 4.2), field capacity (pF = 2.3) and at saturation and also percolation to deeper soil layers (cm day −1 ) in the laboratory.The maximum rooting depth was adjusted to 40 cm.Irrigation supply was imposed at the specific applied times with optimal nitrate application.The simulated LAI was scaled to an hourly basis using linear interpolation between two adjacent simulated daily values of LAI.The model was run for optimal (no water limitation) and realistic conditions (actual water inlet i.e. irrigation and rainfall) for each growing season.Figure 3 represents predicted LAI and grass yield of 2012 and 2013.

Simulation of water flow
The simulated soil profile in the hydrological model extends to 150 cm depth and is divided into two layers: Layer 1 (0 to 33 cm) and Layer 2 (33 to 150 cm).Simulation of root water uptake and water flow, which is assumed to be in the vertical direction in the vadose zone, was carried out for two growing seasons (from 1 March until 25 November in 2012 and 2013) using Hydrus 1D version 4.16 which solves the 1-D Richards' equation: where θ is the volumetric water content (cm 3 cm −3 ), t is time (h), z is the radial and vertical space coordinate taken positive downward (cm), K(h) is the unsaturated hydraulic conductivity function (cm h −1 ), h is the pressure head (cm) and S(h) represents a sink term (cm 3 cm −3 h −1 ), defined as the volume of water removed from a unit volume of soil per unit time due to plant water uptake.
To solve Eq. ( 5), the MVG soil hydraulic model (Eqs.1-4) without air entry value and without hysteresis was used.The initial pressure head distribution was calculated using the inverse of Eq. (3), h(S e ),from the measured initial water content of each observation node.These point values were then interpolated linearly from the deepest observation node to the groundwater level (h = 0, GWL).The pore connectivity parameter of the MVG model was fixed at l = 0.5.The upper condition for water flow was an atmospheric boundary condition -based on rainfall and irrigation water supply, LAI calculated by LINGRA-N (see Sect. 2.3.1) and reference evapotranspiration (ET o ) -with surface runoff.The model performance was assessed using various implemented bottom boundary conditions, i.e. free drainage and incremental constant head conditions, as a manual sensitivity analysis (see the Supplement).The Feddes model (Feddes et al., 1978) without solute stress was used for root water uptake.The default grass parameters values provided by Hydrus 1D were used (Taylor and Ashcroft, 1972).

Soil water stress and yield reduction
In the Feddes model (Feddes et al., 1978) the sink term of Richards' Eq. ( 5), S(h), is specified in terms of quantify potential root water uptake and water stress as where R(x) is the root distribution function (cm), T p is potential transpiration (cm h −1 ) and w (h) is the water stress response function (0 ≤ w (h) ≤ 1) which prescribes the reduction in uptake that occurs due to drought stress.Crop-specific values of this reduction function are chosen from the default Hydrus data set.The actual plant transpiration is calculated numerically, as where L r is the rooting depth (cm).By assuming root water uptake is equal to actual transpiration, the ratio of actual to potential transpiration by the root uptake was introduced as a degree of water stress, DWS (Jarvis, 1989), as The effect of the boundary conditions and parameter uncertainty on soil water stress was evaluated using the ratio between the calculated actual water uptake/actual transpiration and the potential transpiration provided by the model (Li et al., 2012;Zhou et al., 2012).In optimal and stress-free conditions, this ratio should be (close to) unity (> 0.90 of maximum reference evapotranspiration).
The ratio between actual crop evapotranspiration and potential evapotranspiration was introduced as a water stress factor equal to the crop yield reduction due to water shortage (Doorenbos and Kassam, 1979), given as where Y a is actual crop yield, Y m is the maximum crop yield in optimal condition, K y is the crop yield factor (for grass K y = 1), ET a is actual crop evapotranspiration estimated by the model.The Y m value was simulated using LINGRA-N in optimal condition (no water stress) for 2012 and 2013 growing seasons.ET p is potential evapotranspiration and can be calculated from the reference evapotranspiration by where K c is the crop coefficient and equal to 1, assuming that grass at our site did not differ much from the reference crop.Accordingly, crop yield reduction of each scenario was calculated using Eq. ( 9) for both periods to show to what extent different scenarios affect soil water stress and crop yield.

Sensitivity analysis
The effect of each input factor or parameter on the model output is determined by a local sensitivity analysis (SA), using a one-at-a-time (OAT) approach.We used this approach because it allows a clear identification of single-parameter effects.Relevant parameters have major effects on output variables with only a small change in their value (Saltelli et al., 2008).SA is, among other purposes, used to find the most relevant parameters which enable a reduction of the number of parameters that need to be optimized.In a local SA, only the local properties of the parameter values are taken into account, in contrast to global SA which computes a number of local sensitivities.Since the interest in this study goes specifically to the measured (parameter) values in the field, a local SA is chosen.Furthermore, an OAT approach (local or global) does not provide direct information about higher-and total-order parameter interaction as is provided by variancebased SA (Saltelli et al., 2008).However, by evaluating the parameter sensitivities in time, insight is given about potential interaction when similar individual effects are observed.The latter can be quantified by a collinearity analysis (Brun et al., 2001), but will be done graphically in this contribution.
A dynamic sensitivity function can be written as follows: where SF(t), y(t) and x denote the sensitivity function, output variable and parameter respectively.If an output variable (y) significantly changes (evaluated by calculating the variance or coefficient of determination or by visualizing in a scatter plot) due to small changes of the parameter of interest x, it is called a sensitive parameter.This partial derivative can be calculated analytically or numerically with a finite difference approach by a local linearity assumption of the model on the parameters.Local sensitivity functions evaluate the partial derivative around the nominal parameter values.The central differences of the sensitivity function are used to rank the parameter sensitivities and can be expressed as follows: where p f is the perturbation factor, x j is the parameter value and x j is the perturbation, CAS is the Central Absolute Sensitivity, CTRS is the Central Total Relative Sensitivity analysis, and CPRS is a Central Parameter Relative Sensitivity.Since the parameters and variables have different orders of magnitude for which the sensitivity is calculated, direct comparison of the sensitivity indices with CAS is not possible.Hence, recalculation towards relative and comparable values is needed.In order to compare the sensitivity of the different parameters towards the different variables, CTRS is preferred.CPRS is sufficient when the sensitivity of different parameters is compared for a single variable, i.e. soil water content.Here, a dynamic (time-variable) local sensitivity analysis was conducted by linking Eqs. ( 11)-( 14), programmed in Python software (https://www.python.org/) to Hydrus 1D.
Given the output accuracy of Hydrus 1D (0.001), a perturbation factor of 0.1 was chosen.To carry out the SA, each hydraulic parameter (K s , θ r , θ s , α, and n) in each layer was varied (measured value ± perturbation factor) and its CTRS was calculated (Eqs.13-14), while the values of other parameters were fixed to the measured values.The model was run in forward mode 20 times, i.e. 10 runs for each layer and two runs for each parameter.A weak direct effect of a parameter in SA is denoted by low absolute values close to 0. A positive effect is expressed by a positive value and a negative effect by a negative value.

Model calibration
For accurate parameter estimation, a longer period such as a growing season (i.e.2012) with several drying and wetting events was selected.This was also suggested by Wöhling et al. (2008Wöhling et al. ( , 2009)).Therefore, the period between 1 March 2012 (00:00 CET) and 25 November 2012 (23:00 h) was used as the calibration period.We used a time interval of 2 h, resulting in 12 960 soil water content records for four depths (as data for inverse solution), based on hourly precipitation and evaporation input data.Based on our experience this number of data is sufficient for optimization purposes.The objective functions were soil water content and water retention data for both soil layers with unit weighting.In the calibration, we optimized only the values of the most sensitive parameters (K s , n and α) of the two layers, taking initial values of hydraulic parameters for each layer equal to the values estimated by the RETC program for the independent field samples, while keeping the insensitive hydraulic parameters (θ s , θ r ) fixed to the measured values.Thirty-seven parameter optimization scenarios were selected and analysed to identify correlations among optimized parameters and to identify the most influential parameter sets on soil water stress and water content in different lower boundary conditions.The 37 scenarios comprised optimizing all six parameters simultaneously (one scenario), four parameters (nine scenarios), three parameters (18 scenarios) and two parameters (nine scenarios).Finally, the best-performing parameter set -based on performance criteria, the correlation between optimized parameters (non-uniqueness of the parameter sets) and the visual inspection of simulated and observed soil water content -was selected for validation using independent data from 2013 (from 1 March until 12 September 2013).

Model evaluation and statistical analysis
The performance of models can be evaluated with a variety of statistics (Neuman and Wierenga, 2003).It is known that there is no efficiency criterion which performs ideally.Each of the criteria has specific pros and cons which have to be taken into account during model calibration and evaluation.It is suggested to use a combination of different efficiency criteria to assess of the absolute or relative volume error (Krause et al., 2005).The root-mean-square errors (RMSE), the coefficient of determination (r 2 ) and the Nash-Sutcliffe coefficient of model efficiency (C e ) (American Society of Civil Engineers, 1993) are popular and widely used performance criteria to evaluate the difference between observed and modelled data (Wöhling and Vrugt, 2011;Verbist et al., 2009Verbist et al., , 2012;;Gandolfi et al., 2006;Vrugt et al., 2004;Wollschlager et al., 2009;Nasta et al., 2013).They are calculated as follows: where O and S are observed and simulated values at time/place i, respectively.C e and r 2 are considered to be satisfactory when they are close to 1, while RSME should be close to 0. C e may result in negative values when the mean square error exceeds the variance (Hall, 2001).

Irrigation scheduling
The value of soil water stress, and the number and the duration of stress periods was calculated for two growing seasons (2012 and 2013), as an indicator for the performance of the irrigation scheduling (van Dam et al., 2008).To optimize the irrigation scheduling (timing of application), the actual water supply (all irrigation events) was deleted from the model input of the hydrological model.Secondly, the LAI simulated with the LINGRA-N for optimal conditions (no water stress) was used as a variable in the hydrological model.Then, the hydrological model with a constant bottom boundary condition was run with the new input variables to elucidate water stress without actual water supply (see the Supplement).Subsequently, the required irrigation was added to the precipitation at the beginning of each water stress period to exclude water stress from the simulations.To simulate crop yield at the optimized condition, the new precipitation variables (rainfall and required irrigation) were used in LINGRA-N model.The optimal yield obtained using the optimized irrigation scheduling was compared to the actual (simulated and measured) yield of current irrigation management practices.

Parameter sensitivity analysis
Due to the variable rainfall, irrigation, evapotranspiration and drainage, the soil water content changes in the soil profile, and, consequently, parameter sensitivities are time dependent.The soil water content has a low sensitivity to θ s and θ r , especially for the second layer.Low sensitivities to θ r have been reported by others (Kelleners et al., 2005;Mertens et al., 2006;Wöhling et al., 2008).
Figure 4 illustrates the results of the sensitivity analysis as a function of time for the most influential parameters α, n and K s , and for both soil layers as depicted by the suffix 1 for layer 1 and suffix 2 for layer 2. A weak direct effect of a parameter is reflected by low absolute values (close to 0).
The results show for all parameters a general change in sensitivity with time with the seasonal changes in irrigation application and rainfall.Generally, all soil hydraulic parameters showed higher sensitivity in dry periods as compared to wet periods.On the other hand, there is a clear effect of parameter variability in layer 1 on water content estimation at 10 cm, and the effect is slightly declining at 20 and 30 cm, which suggests the great importance and influence of upper boundary variables, especially evapotranspiration.Similar results were observed by Rocha et al. (2006).They found that soil water content and pressure heads were most sensitive to hydraulic parameters variation in the dry period near the soil surface using local sensitivity analysis of Hydrus.
Soil water content is sensitive to variations of α, n and K s in both layers.The sensitivity is the largest for n, α and less so for K s in the first layer.For the second layer, soil water content was most sensitive to α followed by n and K s .Abbasi et al. (2003) reported that n, θ s and K s were the most sensitive parameters in their study and that this sensitivity was more pronounced in deeper parts, however they also observed some sensitivity near the soil surface during the drier conditions.The most sensitive parameters were θ s , n and α and least sensitive parameter was K s in the study by Schneider et al. ( 2013) using Hydrus 1D.They found large interaction (correlation) among sensitive parameters.In contrast, Wegehenkel and Beyrich (2014) reported that soil water content predictions were most sensitive to θ r and θ s and least sensitive to α, n and K s input parameters using Hydrus 1D.Similarly, Caldwell et al. (2013) found that θ r , n and l were sensitive and θ s , α and K s were insensitive to water content simulation.In dry periods, there is a general negative correlation between n and α on the one hand and soil water content on the other hand, whereas a positive correlation exists between K s and soil water content (Fig. 4). Figure 4 shows that in the first layer, the soil water content is more influenced by rainfall at 10 cm than at 30 cm (higher and lower sensitivity for observation nodes 10 and 30 cm, respectively, within first layer).
The fact that the model predictions in the upper part of the soil profile are extremely sensitive to variations in hydraulic parameters in dry periods, is of great importance to irrigation management.To improve the timing of irrigation in these crucial periods, numerical soil models that are used to determine irrigation requirement, need to be well parametrized for α, n and K s .

Model calibration
Since soil water content prediction was insensitive to the parameters θ s and θ r , they were fixed to the measured (initial)  values (Table 1).Similar strategies were used by Verbist et al. (2012) and Schwartz and Evett (2002).
The model was run inversely using time series of soil water content with values for α, n and K s being optimized for the two layers (i.e.six-parameter optimization scenario).A significant correlation appears between optimized α and K s for both layers (layer 1: r = 0.85; layer 2: r = 0.95 constant head; and layer 1: r = 0.82; layer 2: r = 0.80 free drainage) and between optimized n and α (both layers: r = −0.99constant head; and layer 1: r = −0.83and layer 2: r = −0.84free drainage) within each layer, but not between layers.On the other hand, there is a significant correlation between n and K s in both layers (layer 1: r = −0.85;layer 2: r = −0.94constant head; and layer 1: r = −0.75;layer 2: r = −0.98 free drainage).This means that α, n and K s within one layer cannot be determined independently and different sets of correlated parameters lead to very similar predictions of soil water content.The high correlation between optimized parameters within a layer leads to a large uncertainty of the final parameter estimates (Hopmans et al., 2002).To avoid non-uniqueness of the inverse solution (Šimůnek and Hopmans, 2002), 36 additional systematic four-, three-and twoparameter optimizations were conducted.All optimizations resulting in correlations among the optimized parameters were removed and only the optimization scenarios with the uncorrelated parameters were kept.This resulted in parameter values as shown in Table 2 for a constant head corresponding to a groundwater depth of −140 cm and free drainage.For comparison purposes, the six-parameter scenario (all parameters optimized) and only the best-performing optimization with two parameters is presented for the other boundary condition (i.e.GWL = −120 cm).
The performance results of the parameter optimizations according to the performance criteria for all scenarios with uncorrelated parameters and different boundary conditions are presented in Table 3, together with the performance of the six-parameter scenario.The results show that a twoparameter optimization (optimizing only K s in both layers) performs equally well as compared to a six-, four-or threeparameter scenario for all performance criteria and observation depths.However, parameters in the six-parameter scenario are considered unidentifiable due to their correlations.In this case, the model was not able to find a global minimum but found a local minimum (Levenberg-Marquardt method) due to the high dimensionality of the problem (Ritter et al., 2003) and the large uncertainty of the optimized values.
Large differences in model performance were obtained when using free drainage or constant head conditions (Table 3).After optimization, the r 2 for different free drainage and constant head conditions and various optimization scenarios was similar, while C e and RSME were different.Overall, the performance of the model to predict soil water content at 40 cm was lowest.The model performs well for the 10, 20 and 30 cm depths where the plant roots are concentrated and which are consequently the most critical in terms of irrigation optimization.The model with a constant head (−140 cm) clearly performed better than the free drainage boundary condition.The smallest differences were detected at the top node (10 cm) compared to deeper nodes in constant head and free drainage conditions.The optimization approach showed that the free drainage condition was unsuccessful to predict soil water content sufficiently well, in agreement with observations, even using different parameter estimations.
The two-parameter scenario requires fewer parameters (one parameter for each layer) to be optimized, performs better as compared to the uncalibrated model (see Supplementary Material) and is therefore to be preferred.Large confidence limits indicate uncertain estimations of a particular parameter (Šimůnek and Hopmans, 2002).The optimized K s with 95 % confidence limits (CL) for the first and second layer were 1.20 (1.15-1.24)cm h −1 , and 2.17 (2.06-2.26)cm h −1 , respectively, in the two-parameter scenario with −140 cm GWL.Therefore, this optimization result was considered the best and was chosen for the evaluation run.

Model evaluation
The validation results (using the same hydraulic parameter values as in the calibration period) under different upper (rainfall and water supply, ET o , LAI) and lower (groundwater depth, i.e. −135 cm) boundary conditions, show that model performance during the calibration was superior to the validation period at all observation depths (Fig. 5, Table 3).The same result was reported by Wöhling et al. (2008Wöhling et al. ( , 2009)).Similar to the calibration period, soil water content was predicted better during the rain and irrigation period than in the dry period.Specifically, soil water content was overpredicted during summer months (June-August) and underpredicted during winter and spring.Wöhling et al. (2009) explained that the differences can be partly attributed to nonuniqueness of the optimization process, inadequacy of the model structure, the large number of optimized parameters, different information content in the calibration and evalua- tion data, and seasonal changes in soil hydraulic properties.The extent to which the soil water content prediction affects the calculated irrigation requirements is dealt with in the next section.4).Despite the similarity, the extent of soil water stress was larger in 2013 as compared to 2012.This can be attributed to the fact that the first water stress event in 2012 with about 328 h duration is not related to soil water availability but is also due to climate limitations (low temperature and light-radiation limitation).No significant reduction or increase in yield and LAI was achieved during this first water stress event in current and optimum conditions (Fig. 3).
There was a significant effect of the bottom boundary condition on the calculated water stress.A free drainage condition resulted in a larger number, longer duration of stress con-ditions (Fig. 6 and Table 4) and overestimated water stress due to excessive recharge to the groundwater (more than 148 mm).On the other hand, a shallower imposed groundwater level (−120 cm) creates less estimated water stress (Fig. 6 and Table 4), because this boundary condition allows inflow (upward flow) from the groundwater table.When the groundwater level was −140 cm the outflow of the bottom flux increased from the six-optimized parameter scenario (−4.6 mm) to the two-parameter scenario (−15.4 mm) in the calibration period, while upward flow increased with increasing number of optimized parameters in validation period (63.3 to 76.9 mm).But these inflows did not meet the crop water requirement (see next section).Huo et al. (2012) reported that the maximum contribution of groundwater level to crop water requirement occurred when the groundwater level was less than 100 cm.Overall, to overcome the water stress effects on crop yield, additional irrigation should be supplied for different optimization scenarios and boundary conditions.During water stress, yield reduction would be in range of 0 to 33 % for different optimization scenarios (Table 4).In addition, two-to six-parameter optimizations showed a similar value in yield reduction (16 % for two-and 13 % for three-to six-parameter in calibration; and 13 % for two-and 11 % for three-to six-parameters to be optimized in validation periods).The maximum yield reduction occurred in the free drainage condition among different boundary conditions and parameter optimization scenarios.Different parameter optimization strategies (two-, three-, fouror six-parameter optimizations) do not affect the calculated water stress as significantly as does the bottom boundary.Therefore, these results suggest that simultaneous optimization is needed for irrigation management purposes, i.e. optimize/choosing boundary conditions to accurately describe recharge to or from groundwater and, in second order, opti-  mize hydraulic parameters to accurately describe soil water content variation in the topsoil.

Irrigation scheduling scheme
The simulated results further showed that, to avoid drought stress during summer, a more accurate irrigation schedule  would be needed in the drier part of the field.It would be better to supply water in June and July instead of a huge amount in late summer or at an inappropriate time (see Figs. 6 and 7).
Results revealed that the actual water supply exceeded crop demand but did not meet the crop requirement (Fig. 7 and Table 5).Irrigation volume affects soil water fluxes.In the "no irrigation" scenario for 2012 the upward/inflow fluxes from groundwater were larger than current and guided irrigation scenarios (Fig. 8).The upward flow of water was not sufficient to meet the crop requirement.For guided irrigation, recharge from groundwater was larger than current irrigation in 2012 and 2013 -which means some part of crop water de- Results show that, despite reducing water supply throughout the growth period by about 22.5 % in 2012 and 12 % in 2013, yield would have increased about 4.5 % in 2012 and 6.5 % in 2013 on average (Table 5, Fig. 3), by rescheduling irrigation at the precise time when the crop is exposed to water stress.The number of irrigation events would remain similar to realistic applications (three times in each growing season).At the field scale, non-uniform irrigation distribution (water supply in drier parts with groundwater level below 120 cm) would be necessary.

Conclusions
The results of this study have demonstrated clearly the profound effect of the position of the groundwater table on the estimated soil water content and associated water stress in a sandy two-layered soil under grass in a temperate maritime climate.Indeed, field-scale variations in soil water content can be very large, due to topography and variable depth of the groundwater.Furthermore, the model performance was affected by the spatial variability of hydraulic parameters such as K s .Results show that the uniform distribution of water using standard gun sprinkler irrigation may not be an efficient approach since at locations with shallow groundwater, the amount of water applied will be excessive as compared to the crop requirements, while in locations with a deeper groundwater table, the crop irrigation requirements will not be met during crop water stress.
The results show that the effect of groundwater level was dominant in soil water content prediction, at least under conditions similar to those in our study.This reflects the need for accurate determination of the bottom boundary condition, both in space and time.In a subsequent field experiment in an adjacent field, the temporal fluctuations of the groundwater table based on diver (Mini-Diver, Eijkelkamp Agrisearch Equipment, Giesbeek, the Netherlands) measurements in boreholes revealed changes in groundwater depth of about 10 cm.The temporal changes were smaller than the expected variation due to topography which may well range more than 100 cm even for relatively flat areas.This has important consequences for precision irrigation management and variable water applications at sub-field scale.The use of detailed (cm scale) digital elevation models, geophysical measurement techniques such as electromagnetic induction or ground-penetrating radar as proxies for hydraulic parameters will serve as valuable data sources for hydrological models to calculate variable irrigation requirements within agricultural fields.The parametrization scenarios in the calibration and validation stage of model development should be kept simple in view of the information they generate.We have shown that it is sufficient to estimate a limited amount of key parameters for which the temporal variant informa-tion of the sensitivity is crucial, and also that optimization strategies involving multiple parameters do not perform better in view of the optimization of irrigation management.We have shown that a combined modelling approach could increase water use efficiency (12-22.5 %) and yield (5-7 %) by changing the irrigation scheduling.However, these efficiencies can only be achieved if rainfall is known a priori-while the soil water status could indicate when to irrigate, it would be impossible to know how much to irrigate if the rainfall cannot be accurately predicted.Therefore, the results of the study call for taking into account accurate weather forecast and water content data in irrigation management and precision agriculture.The combination of accurate and spatially distributed field data with appropriate numerical models will make it possible to accurately determine the field-scale irrigation requirements, taking into account variations in boundary conditions across the field and the spatial variations of model parameters.The information gained in this study with respect to dominant parameters and the effect of boundary conditions at the plot scale (1-D) will be scaled up in a 2-D approach to the field scale using detailed spatial information on groundwater depth and hydraulic conductivity K s .
The Supplement related to this article is available online at doi:10.5194/hess-20-487-2016-supplement.

Figure 1 .
Figure 1.Geographical location of the experimental field and the map of the apparent soil electrical conductivity (EC a ) of the study site corresponding to three different zones of groundwater levels (GWL).The black star on the EC a map indicates the sensor location.

Figure 3 .
Figure 3. Predicted leaf area index, LAI and grass yield using the LINGRA-N model for 2012 and 2013.

Figure 4 .
Figure 4. Parameter sensitivity as a function of time.The numbers 1 and 2 correspond to the first and second layer, respectively.

Figure 5 .
Figure 5. Observed and simulated time series of soil water content with calibration using the two-parameter K s scenario for 2012 and validation results of 2013.

Figure 6 .
Figure 6.Degree of water stress at potential reference evapotranspiration in 2012 and 2013 for various scenarios and bottom boundary conditions.

Figure 7 .
Figure 7.Comparison degree of water stress between farmer's conventional irrigation (current irrigation), without irrigation and optimized irrigation scheme for calibration and validation periods.

Figure 8 .
Figure 8. Actual flux of farmer's conventional irrigation (current irrigation), without irrigation and optimized irrigation scheme (guided irrigation) for 2012 and 2013.

Table 2 .
Optimized values of hydraulic parameters for the optimization scenarios yielding uncorrelated parameters (except for reference scenario with six optimized parameters).Values indicated in italic are values fixed to the measured values close to the sensor location.Numbers in parentheses represent the standard errors of optimized parameter.

Table 3 .
Calculated performance criteria describing the correspondence between measured and simulated soil water content for each scenario for various boundary conditions.RMSE, C e and r 2 are the root-mean-square deviation (cm 3 cm −3 ), the Nash-Sutcliffe coefficient of efficiency and the coefficient of determination. *

Table 4 .
Total duration, number and extent of water stress for different boundary conditions and scenarios (from 1 March to 12 September).Total rainfall and irrigation amount were 398.2 and 64.5 mm in 2012 and 343.3 and 85.4 mm in 2013 respectively.Number between parentheses represents the duration of first water stress event due to light-radiation and temperature limitations.

Table 5 .
Comparison of optimized irrigation schedule with farmer's conventional irrigation schedule.