Influence of ecohydrologic feedbacks from simulated crop growth on integrated regional hydrologic simulations under climate scenarios

Hydrologic climate change modelling is hampered by climate-dependent model parameterizations. To reduce this dependency, we extended the regional hydrologic modelling framework SIMGRO to host a two-way coupling between the soil moisture model MetaSWAP and the crop growth simulation model WOFOST, accounting for ecohydrologic feedbacks in terms of radiation fraction that reaches the soil, crop coefficient, interception fraction of rainfall, interception storage capacity, and root zone depth. Except for the last, these feedbacks are dependent on the leaf area index (LAI). The influence of regional groundwater on crop growth is included via a coupling to MODFLOW. Two versions of the MetaSWAP-WOFOST coupling were set up: one with exogenous vegetation parameters, the “static” model, and one with endogenous crop growth simulation, the “dynamic” model. Parameterization of the static and dynamic models ensured that for the current climate the simulated long-term averages of actual evapotranspiration are the same for both models. Simulations were made for two climate scenarios and two crops: grass and potato. In the dynamic model, higher temperatures in a warm year under the current climate resulted in accelerated crop development, and in the case of potato a shorter growing season, thus partly avoiding the late summer heat. The static model has a higher potential transpiration; depending on the available soil moisture, this translates to a higher actual transpiration. This difference between static and dynamic models is enlarged by climate change in combination with higher CO 2 concentrations. Including the dynamic crop simulation gives for potato (and other annual arable land crops) systematically higher effects on the predicted recharge change due to climate change. Crop yields from soils with poor water retention capacities strongly depend on capillary rise if moisture supply from other sources is limited. Thus, including a crop simulation model in an integrated hydrologic simulation provides a valuable addition for hydrologic modelling as well as for crop modelling.


Introduction
In hydrologic models, vegetation characteristics are usually defined by "exogenous" parameters; a fixed dependency on the days of a year is assumed.This means that ecohydrologic vegetation feedbacks to the hydrologic system are neglected.The resulting limitations on the model validity become more poignant with the advent of climate change impact modelling using scenarios.These scenarios usually differ widely from current climate, which increases the necessity for endogenously simulating the weather-and climate-dependent vegetation feedback to the hydrologic system.
SWAT (Arnold and Fohrer, 2005;Neitsch et al., 2011) is an example of a regional hydrologic model that includes dynamic crop growth simulation and feedback to vegetationrelated parameters.However, the soil water submodel is a very simple two-layer type which cannot simulate capillary rise.This severely limits its applicability for simulations that require a feedback loop via groundwater.Another example of a regional integrated model is PROMET (Mauser and Bach, 2009).This approach to soil water modelling also lacks sophistication.It involves the repeated use of an analytical model with four sublayers.Compared to SWAT, the improvement is that capillary rise is modelled.However, the soil Published by Copernicus Publications on behalf of the European Geosciences Union.
water model is of the simple "bucket" type that neglects the soil water influence on the groundwater storage coefficient.
Other models that simulate the interactions between a soil column and crop development are the SWAP (Van Dam et al., 2008) and the Theseus (Wegehenkel, 2009) models.Both models are coupled to the crop growth simulation model WOFOST (Van Ittersum et al., 2003;Supit et al., 1994).However, in both models not all vegetation feedbacks have been included, and neither has a feedback loop via the regional groundwater system.
In this study, we use the regional hydrologic modelling framework SIMGRO (Van Walsum and Veldhuizen, 2011;Van Walsum et al., 2012).This framework hosts the coupling of the groundwater model MODFLOW (McDonald and Harbaugh, 1988) to the SVAT model MetaSWAP (Van Walsum and Groenendijk, 2008).The unsaturated flow simulation in MetaSWAP is based on a quasi-steady state schematization of flow processes in combination with water balances of control volumes.It is suitable for regions with moderate slopes with shallow or deep groundwater levels.It has been calibrated and successfully validated using results of the Richards-type SWAP model.Recently, the crop growth simulation model WOFOST has been coupled to MetaSWAP, including a two-way feedback.
To account for the increasing atmospheric CO 2 concentration, the maximum and initial angle of the CO 2 assimilationlight response in WOFOST are adapted according to Wolf et al. (2010).Since the increasing atmospheric CO 2 concentration and rising temperature directly affect growth and leaf development -and consequently the ecohydrologic feedbacks -we use leaf area index (LAI) values simulated with WOFOST to establish the fraction of radiation that reaches the soil, the crop coefficients, and the interception storage capacity.The root zone depth simulated by WOFOST is also included as feedback to MetaSWAP.
In this paper, we first present an outline of the models used.Special attention is given to the MetaSWAP-WOFOST coupling.Next, the test region is briefly described.Subsequently, the simulated effects of climate change scenarios on evaporation, transpiration and canopy interception evaporation are presented.Results obtained with endogenous vegetation parameters (i.e. with feedback from WOFOST simulation results) are compared to those obtained with exogenous parameters (i.e.static vegetation parameters).The model that uses crop simulation results is referred to as the "dynamic" model, whereas the model that uses static vegetation parameters is called the "static" model.
The static model parameters are derived from daily average crop variables that were simulated with the dynamic model for a 30-yr period.The advantage is that the simulated long-term averages of the water balance terms for current climate are then nearly the same for both the static and dynamic crop models.This facilitates the comparison between the two models when the climate scenarios are used as input.
2 Methods and materials

Modelling framework and attached models
An overview of the SIMGRO modelling framework and its attached models (insofar relevant for this study) is given in Fig. 1: -MODFLOW: calculates water flows in the saturated zone on a regional scale; -MetaSWAP: calculates water flows in a Soil Vegetation Atmosphere Transfer (SVAT) column; and -WOFOST: calculates crop growth.

Calculation method for evapotranspiration terms
An evapotranspiration simulation method should be accompanied by a method that provides acceptable parameters for the intended user community.In our case we used the crop coefficients given by Feddes (1987).Those coefficients, however, are intended for an "all in one" evapotranspiration simulation method that does not account for variable feedbacks from crop growth.Those feedbacks affect each of the evapotranspiration components in a specific manner.We therefore model evapotranspiration as three separate and partly interdependent terms: crop transpiration, canopy interception evaporation, and soil evaporation.The parameters are calibrated on the coefficients of Feddes (1987).
As a starting point, we used an adapted form of the dual crop coefficient approach given by Wright (1982) and Allen et al. (2005): where T p is the potential crop transpiration (m d −1 ), E s,p is the potential soil evaporation (m d −1 ), K cb is the basal crop coefficient, K ew is the evaporation coefficient of a wet bare soil that is partly shielded by vegetation, and ET 0 is the reference crop evapotranspiration (m d −1 ).
The reference crop evapotranspiration is calculated with the simplified Makkink equation, presented by De Bruin (1981,1987): where s is the slope of the vapour pressure curve, γ is the psychrometric constant (kPa K −1 ), K ↓ is the incoming global radiation (W m −2 ), L v is the latent heat of vaporization of water (J kg −1 ), and c is a conversion factor (m d −1 / kg m −2 s −1 ).De Bruin (1987) gives the following considerations for choosing the Makkink equation:  -its behaviour is very similar to that of the Penman formula; -it is remarkably simple: it requires only air temperature and global radiation as input; -under dry conditions Makkink's formula performs better.
Soil evaporation also includes the evaporation from soil that is partly shielded by a canopy.The "bare" soil is assumed to be "diffuse", meaning that we do not distinguish between soil that is covered by vegetation and soil that is not.We assume that the net radiation inside the canopy decreases according to an exponential extinction function of the LAI and that the soil heat flux can be neglected (Goudriaan, 1977;Belmans, 1983): where κ gr is the extinction coefficient for solar radiation (-), LAI is the leaf area index, and K ew100 is the evaporation coefficient of a wet soil receiving full radiation.Ritchie (1972) and Feddes (1978) used κ gr = 0.39 for common crops.More recent approaches estimate κ gr as the product of the dimensionless extinction coefficient for diffuse visible light, κ df , which varies with crop type from 0.4 to 1.1, and the extinction coefficient for direct visible light, κ dr : We use the method from Boesten and Stroosnijder (1986) to calculate actual soil evaporation.This method only requires the time sequence of rainfall events and the reference crop evapotranspiration multiplied by the K ew coefficient; it has a proven track record in the Netherlands.Soil moisture conditions are assumed to limit transpiration according to Feddes et al. (1978): where T a is the actual transpiration rate (m d −1 ) and α rw is a dimensionless stress coefficient [0-1] that gives the influence of available soil water on transpiration (Fig. 2).A wet canopy leads to increased evapotranspiration due to interception evaporation, with a lower resistance for the moisture flow from the canopy to the atmosphere.Wright (1982) accounts for this lower resistance by temporarily adjusting the crop coefficient.We model this effect separately by explicitly modelling canopy interception evaporation.Interception evaporation from a vegetation canopy is not sensitive to soil moisture conditions as described by the α rw coefficient.This makes it relevant to model canopy intercerption evaporation separately from the transpiration.
To model canopy interception evaporation, Eq. ( 1) is expanded to: where E i,p is the potential canopy evaporation rate, and K iw is the evaporation coefficient of a wet canopy, with K iw > K cb .The actual canopy evaporation and dripping rate to the soil are determined with the method given by Rutter et al. (1971).However, this method includes a linear relationship between the canopy storage and the relative interception evaporation.As Rutter remarks, such a relationship entails that the complete drying out of a canopy takes an infinite time, which is not realistic.We have remedied this problem by assuming a discontinuous relationship with a "start-up" value instead of a linear one starting from zero, as described in Appendix A. Evaporation from a wet canopy is assumed to be a dominant process: as long as it is active, transpiration is assumed to be zero.This assumption is based on the notions that evaporation "has first choice" for using the energy flux and that due to a higher air moisture content transpiration is suppressed.To account for this dominance, the active time fraction W frac is calculated: where E i,a is the actual canopy interception evaporation rate (m d −1 ).This time fraction is used in the calculation of actual transpiration.However, to correctly calculate relative transpiration for the crop growth model (Sect.2.1.2),it is essential to account for interception dominance in the potential transpiration calculation step; otherwise, interception evaporation would lead to a reduced T a /T p , which is conceptually wrong.Therefore, Eq. 6 is finally modified to:

Coupling to WOFOST
The WOFOST principles have been discussed by Van Keulen and Wolf (1986) and Van Diepen et al. (1989).WOFOST's structure is described by Supit et al. (1994); a standalone version of WOFOST is described by Boogaard et al. (1998).It has been modified and applied for many purposes (e.g.De Koning and Van Diepen, 1992;Rötter, 1993;Supit, 1997).The crop growth model calculates daily crop photosynthesis based on absorbed radiation and water and/or salinity stress.After subtraction of the maintenance respiration, the remainder of the photosynthesis products are partitioned over leaves, stems, roots and reproductive organs.The most important internal driver is the leaf area index, which is the result of the leaf area dynamics controlled by photosynthesis, allocation of biomass to leaves, leaf age and development stage.In turn, LAI controls the daily rates of photosynthesis.
The temperature influence on crop development is simulated via the temperature sum, which is a summation of daily temperatures above a certain threshold value.This sum determines the germination process and the "phenological" development stage, which in turn affects CO 2 -assimilation.Apart from the indirect effect via temperature accumulation, there is also the direct effect of sub-optimal day temperature on crop growth and development.
WOFOST and MetaSWAP (Fig. 3) communicate on a daily basis.Crop water stress causes the leaf pores to close partly or completely, thus minimizing moisture loss.This also reduces CO 2 assimilation and crop growth.To model this reduction, WOFOST requires "relative transpiration" as input, i.e.T a /T p .During periods of interception evaporation, transpiration is assumed to be zero; then relative transpiration is set equal to the α rw -coefficient for a low value of T p (line for T low in Fig. 2).
WOFOST returns the root zone depth and LAI; the latter is subsequently used to calculate: -the fraction of radiation reaching the soil surface (Eq.3); -the transpiration crop coefficient; and -the rainfall interception fraction, the canopy interception storage capacity, and the interception evaporation crop coefficient.
Feddes (1987) gives transpiration crop coefficients for each ten-day interval during the growing season; these coefficients are based on field experiments.They are intended for a model that does not separately simulate soil evaporation and transpiration, let alone separate modelling of canopy interception evaporation.However, all these were of course active during the experiments, but were not to be measured directly.A method has been devised for back-calculating the separate coefficients for transpiration, and linking them to the LAI in the form of a table function (Appendix B).This method uses the LAI-simulations obtained from an unstressed MetaSWAP-WOFOST run, with a relative transpiration of 1.0 at all times.These results are then entered as data in the following adaptation of Eq. ( 1):  where K c,tot is the "all in one" crop coefficient given by Feddes (1987) for each ten-day interval of the growing season, K cb (LAI) is the crop transpiration coefficient as a function of LAI, E s,a is the simulated soil evaporation of the unstressed MetaSWAP-WOFOST run, and e is the residual error.The K cb (LAI) function has been implemented as a table in steps of 0.001 for the LAI, and with the constraint that it should be non-decreasing, and with a non-increasing first derivative (Appendix B).For grassland and arable land (potato), the derived relationships are given in Fig. 4. The simulation results for the unstressed run have also been used for deriving tables of daily crop parameter values that are applied in the "static" model.
The LAI affects the interception in various ways: -the fraction of rainfall that directly reaches the soil surface without passing through the interception reservoir is set equal to the fraction of radiation reaching the soil surface, as given in Eq. (3) for soil evaporation; -the interception storage capacity of a canopy is assumed to be linearly dependent on LAI according to S i,cap = s i,cap • LAI, where s i,cap is the capacity per unit of LAI (mm LAI −1 ); and -the interception evaporation crop coefficient K iw is set to 1.2 times K cb ; this factor has been set slightly lower than the ratio between the coefficients for open water and short grass (1.25).
Information about realistic canopy interception totals are hard to come by.Calder (1990) reported that grassland interception evaporation is about 15 % of the rainfall for sites in upland areas of Great Britain, as compared to 30 % for forests.Such figures should be treated with great care when using them for other locations.Dolman et al. (2000) reported roughly the same figures for forests in the Netherlands.This provides some degree of confidence that the grassland values of Calder are valid for this study.Applied to the annual mean rainfall of 800 mm (De Bilt, NL), the interception of grassland is ∼120 mm yr −1 .For this an interception capacity of 0.25 mm LAI −1 has been calibrated with the model given in Appendix A, using rainfall with a sampling interval of 1 h.For lack of information about interception capacity of arable land crops, this value has also been used for potato (the influence of this parameter uncertainty is reported in Sect.3).
In the climate change scenario W + , the CO 2 concentration is assumed to increase to 567 ppm.In WOFOST, the effects on crop growth are included by increasing the initial angle and the maximum of the CO 2 -assimilation light curve by respectively 6 % and 33 % (Wolf et al., 2010).There is also an effect on the transpiration crop coefficient.For both grassland and potato, the reduction is estimated at 6 %.

Study region
The Kromme Rijn is a fork of one of the main Rhine branches, the Lek (Fig. 5).The region surrounding it (33 610 ha) is part of the waterboard De Stichtse Rijnlanden (www.hdsr.nl).Along the Northeastern side, the region is bordered by what is left over of an end moraine, the Utrechtse Heuvelrug, rising to an elevation of 65 m above sea level.The Southeastern border of the region is formed by the Amsterdam-Rhine Canal with a water level of 0.40 m below sea level, which is several meters below the soil surface, thus causing a substantial regional drainage.This drainage is mainly balanced by infiltration from the Lek and the Kromme Rijn.Seepage comes from the Utrechtse Heuvelrug, but it is much reduced in comparison to the past.This reduction is caused by heavy groundwater pumping for drinking water   supply.An overview of the land use in the region is given in Table 1.
The SIMGRO framework and the coupled models (Fig. 1) have been implemented for the whole waterboard and beyond, using a grid of 100 × 100 m.There can be several SVAT columns coupled to a single grid cell to represent the areal fractions of vegetated soil, surface water, and impermeable surface ("tiles").The groundwater model has a schematization of 8 aquifers.
Model calibration on phreatic level time series is not yet possible due to the scarcity of usable data.Most of the available phreatic measuring points are on the field boundaries and/or near water courses, thus making them nonrepresentative.The comparison between measured and simulated values for one of the points is given in Fig. 6.This comparison is of course no more than a cursory "plausibility check".However, for the purposes of the current study it is seen as sufficient: the focus is here on the influence of the feedbacks from the vegetation simulation.In such a conceptual exploration, it is the model sensitivity that is of interest, and not so much the absolute predictions for the investigated scenarios and variants.
The available data deck of the model covers the 17-yr period 1989-2005.The dynamic vegetation simulation can be expected to differ most from the static simulation in relatively warm and dry years.An example of such a year is 2003; for this year the maximum and minimum groundwater level transects are shown in Fig. 7 for the cross-section AB indicated in Fig. 5.The net saturated flux (as simulated with the "BCF"package of MODFLOW) is shown in Fig. 8 as a time-average for the year 2003.Locations with positive values correspond to locations with active drainage media.

Scenarios and investigated variants
Climate scenarios for the Netherlands (Table 2) have been taken from Van den Hurk et al. (2006).The per cent changes given in Table 2 relate to climatic means involving a 30-yr period.The scenario weather series have been derived with deterministic transformation rules; they are not randomized as would have been the case if a stochastic weather generator had been used.Therefore, comparing results for a specific year in the original weather series to results for its pendant year in the climate scenario series is possible.
The four scenarios form a 2 × 2 matrix, with respectively two possible developments for the global temperature increase and two possible developments for the atmospheric circulation.In this study we use the W + -scenario with +2 • C increase in 2050, because it has the biggest impact on vegetation development and on feedbacks affecting the hydrologic system.As can be seen from Table 2, substantial potential evapotranspiration changes are expected in this scenario owing to temperature and wind changes.Little is known about the possible changes of the radiation due to changes of cloud cover patterns in the climate scenarios; for this reason the radiation has been assumed unchanged.In the current land-use situation, the agricultural land use includes both grassland and arable land.As can be expected, both land-use types have a different sensitivity to climate change.For each type separate runs have been made, one with all agricultural land as grassland and the other with all of it as arable.The investigated variants of climate and land use scenarios without/with endogenously simulated crop growth are listed in Table 3.
This study focuses on agricultural land use types.For grassland, the WOFOST parameters have been taken from Kroes and Supit (2011).Other arable land crops are modelled with "potato" using crop parameters provided by Wolf et al. (2010).In order to "fairly" compare the dynamic and static simulation results, a fixed sowing date for both models is assumed (day 111 for potato).

Introduction
The vegetation parameters have been determined in such a manner that the average evapotranspiration over 30 yr for the static and dynamic vegetation simulations are equal.In more extreme meteorological years, conceptual differences between the static and dynamic crop model may have larger impacts on the simulation results.To study the effects of warm and dry conditions, simulation results for 2003 have been analysed; in the climate scenario data series for "2050", the year 2003 has its pendant in 2063.based on the growth in a median year (Appendix B).In a warm and dry year like the shown 2003, the LAI development simulated by WOFOST can partly fail from mid-July until the end of September, as can be seen in the results for scenario C Gr Dyn (Fig. 9) with a retarded hay cutting.The dynamic simulation for the climate scenario shows a completely missed hay cutting.However, since the grass never completely dies off, the change of the grassland evapotranspiration year total with respect to the current climate (Table 4, column R cl ) predicted by the dynamic vegetation simulation does not significantly differ from the static one.

Arable land
In a warm year, potato develops faster in the dynamic simulations than in the static one, especially in the W + climate scenario for 2003/2063 (Fig. 10).The higher temperatures shorten the germination period in the dynamic model, in comparison to the long-term average for the current climate that is applied in the static model; furthermore, the dynamic model simulates a shorter life span due to the increased temperature.The increased atmospheric CO 2 concentration leads to a higher LAI peak.The peak marks the end of the vegetative phase, which is followed by the generative phase.In the latter phase the tubers are formed.The faster LAI development leads to a higher potential transpiration in the vegetative phase.In contrast, in the generative phase the dynamically simulated LAI quickly drops below that of the static model, which translates for the year 2003/2063 to a lower potential transpiration than the transpiration obtained with the static model.This happens to such a degree in the W + scenario that the year total potential transpiration of the dynamic model is 10 % less than that of the static model (Table 4).The recharge (P -ET a ) is slightly higher for the dynamic than for the static crop simulation ( R cl ) for the W + scenario.This is because actual transpiration at the evaluation point EP1 is slightly more sensitive to drier conditions if the vegetation feedback is included.
As can be seen from Fig. 11, recharge from arable land vegetation strongly varies along the cross-section AB (Fig. 5).To a large extent this is caused by soil type variations such as the peak at x AB = 4200 m, which has a podzolic coarse-textured sandy soil with a poor water retention capacity, leading to a low actual evapotranspiration and high recharge.This high recharge is not sensitive for the dynamic/static simulation, because the soil moisture supply is in this case the limiting factor for the actual evapotranspiration.
To investigate the groundwater depth influence on the recharge sensitivity for static/dynamic vegetation simulation ( R veg ), locations with the same podzolic medium-textured sandy soil were selected.For this selection, R veg was plotted against the mean of the groundwater level on day 4 and day 264 of 2063 in the W + scenario (Fig. 12).The apparent strong correlation is explained by the difference between the potential and actual transpiration of static/dynamic models, with the static model simulating a 10 % higher value of the potential transpiration (Table 4).This higher value can be better approached by the actual value at locations with a shallower groundwater level, which have an increased potential for moisture supply via capillary rise.The increased actual transpiration has a lowering effect on the recharge of the static model.Thus the gap between recharge simulated by dynamic and static models widens for conditions with shallower groundwater levels that are supported by seepage from deeper layers.At locations where the actual transpiration is not at all limited by moisture stress, the recharge difference is equal to the above mentioned difference of 10 % in the potential transpiration.In terms of effects on the groundwater level itself (cf.Fig. 7) and the saturated flux (Fig. 8), the difference between the static and dynamic model is too small to be shown in a figure; this is partly due to the fact that agriculture is not the only land use.
To show the dependency of crop yield on groundwater conditions, the relative yield has been plotted in Fig. 13 for the same selection of points as in Fig. 12, involving the same soil type.As can be seen, the relative yield varies from 0.5 to 1.0.This strong variation is caused by the limited water retention capacity of the podzolic soil; the moisture supply to the crop depends heavily on the capillary rise from groundwater.Thus, at locations with shallow groundwater the productivity is the highest.This corresponds to situations with high evapotranspiration and low recharge; the plot in Fig. 13 is therefore the inverse of that in Fig. 11.

Calculation method for evapotranspiration terms
Each method to compute the reference crop evapotranspiration has its limitations.This also applies to the much used Penman-Monteith equation for short grass or alfalfa (Wright, 1982;Allen et al., 1998), which requires actual vapour pressure measurements that are sensitive to the circumstances in the direct vicinity of the gauging station.If this method is used under circumstances that differ from the defined reference situation and the measurements are not first interpreted using an advanced regional meteorological model, the results will contain an error of unknown value.Also, the "full" Penman-Monteith equation used with a crop resistance of the actual crop considered (instead of the crop coefficient) is not necessarily superior to the Makkink equation.Jacobs and De Bruin (1998) showed that Makkink's equation, when applied to unstressed maize in combination with a crop coefficient, explained experimental results better than the Penman-Monteith equation.An additional reason for using the Makkink equation is that the climate scenarios for the Netherlands ( Van den Hurk et al., 2006) have been specified in terms of modifications of the Makkink values.For the present study, the method used is not essential to estimate the ecohydrologic feedbacks from crop growth.In practice,  there is not much difference between the two methods: the long term average for 1971-2000 (De Bilt, NL) for Makkink is less than 2 % higher.

Water uptake function
Water uptake by plants is known to be influenced by temperature, as shown for cucumber plants by Yoshida and Eguchi (1989).The pressure-head dependent reduction function (Feddes et al., 1978) lacks this dependency.Bartholomeus et al. (2008) proposed a process-based simulation method that includes such a dependency, in combination with a model for the influence of too wet conditions on the oxygen stress of plants.This method has, however, not yet found its way to modelling practice, partly due to lack of both soil-specific data and experimental facilities to validate the approach: There is not a single lysimeter operational in the Netherlands.

Crop coefficient of transpiration
Allen et al. (1998) give the following K cb (LAI) relationship for non-pristine agricultural vegetation: where K cb,mid is the estimated mid-season basal K cb when plant density and/or leaf area are lower than for full cover conditions, K cb,full is the estimated mid-season basal K cb (at peak plant size or height) for vegetation having full ground cover or LAI > 3, and K cb,min is the minimum K cb for bare soil (K cb,min ≈ 0.15 − 0.20).To compare this method to the ) and the difference between the recharge that is simulate 3 scenario with dynamic and static vegetation models, with ∆R veg = R dyn -R stat ( 4shown data points are for the same podzolic soil type.5 6 Fig. 12. Relationship between the mean groundwater level on days 4 January 2063 (day 4) and 24 September 2063 (day 264) and the difference between the recharge simulated for the W + scenario with dynamic and static vegetation models, with R veg = R dyn -R stat (mm yr −1 ).All shown data points are for the same podzolic soil type.method described in this paper, we have set K cb,min to zero, which results in the function shown in Fig. 14.More recently, Jayanthi et al. (2007) developed a canopy reflectance-based crop coefficient for potato; they also give a relationship between the "soil adjusted vegetation index" (SAVI) and the leaf area index.Rearranging the given relationships and scaling for reaching K cb,full at an LAI of 3.5 (the value for which they find an alfalfa-based coefficient of 0.8, the maximum value in their field observations), gives the K cb (LAI) function shown in Fig. 14   The simulation results were compared (Table 5) for the 1971-2000 weather series (DeBilt, NL) for a crop with an optimal water supply (T a = T p ).The results show a significant difference with both the function of Allen et al. (1998) and with that of Jayanti et al. (2007), which both simulate a potential transpiration that is roughly 9 % higher and a total evapotranspiration about 6 % higher.We explain this difference by the fact that the other methods are actually based on soil cover instead of on LAI.A full cover is reached at a LAI value of 3-3.5.However, the crop coefficient sensitivity does not stop there: the conductance of stoma will increase further as the crop continues to develop, reaching a maximum LAI of about 5.

Interception capacity
As indicated in Sect.2.1.2,information about the interception capacity is indirect and thus uncertain.A value of 0.25 mm LAI −1 has been assumed; it yields for grassland interception evaporation a long-term average (1971-2000, De Bilt, NL) of ∼120 mm yr −1 .If the interception capacity is doubled, then this increases to ∼160 mm yr −1 .The results for evaluation point EP1 for the year 2003 (cf.Table 4) show a similar percentage increase.The simulated change of recharge for the W + scenario ( R cl ) is slightly less for the static model (14.7 vs. 15.2 mm yr −1 ), but the change for the dynamic model increases from 14.5 to 18.5 mm yr −1 .This means that our results are only moderately sensitive to the uncertainty in the assumed interception capacity.

Conclusions
Correct parameteriszation of ecohydrologic parameters is of key importance when using regional hydrologic models for climate change analysis.Increasing temperatures and rising CO 2 concentration not only directly affect evapotranspiration, but also the crop development changes.Consequently, hydrologic processes that depend on crop developmentsuch as soil evaporation, transpiration, and interceptionare also indirectly affected.In this study, we investigated whether the regional hydrologic modelling with MetaSWAP can be improved by using vegetation parameters coming from WOFOST.We used climate change scenarios from the Dutch Meteorological Institute (KNMI) as input and compared results of standard MetaSWAP with an extended version that includes dynamic vegetation modelling.As crops we used grass and potato.
The dynamic potato model simulates accelerated growth and a shorter growing cycle, owing to increasing temperatures in the climate change scenarios.In combination with the increasing CO 2 concentration, this causes a 10 % lower potential transpiration in the example year compared to the potential transpiration of the static model.Similar to the dynamic potato simulation, potential transpiration in the dynamic grassland simulations is lower than the static model.The effect on the water balance terms of grassland, however, is minor, because the soil cover is always >50 % in both model options, thus reducing the impact of any changes in the simulated growth.If the moisture supply is ample, the higher potential transpiration of the static model translates directly into a higher actual transpiration.The degree to which this happens determines how much lower the recharge is when the static model is used.
The productivity of soils that have relatively poor water retention capacities strongly depends on capillary rise in a warm dry year.This dependency can be simulated by coupling a groundwater model to the used modelling framework.
Climate change and extreme weather events lead to changing ecohydrologic feedbacks in the hydrologic cycle.These changes are difficult to capture in the static model parameters.This study demonstrates that these feedbacks can be endogenously accounted for by coupling a crop growth model to a hydrological modelling framework.This avoids the necessity of repeating the static model parameterization procedure for each climate change scenario.

A2 Case 2: Semi-continuous dependence of evaporation rate on canopy saturation
If f i,min in Eq. ( A4a) is less than unity, there is a semicontinuous dependence of the evaporation rate on the canopy saturation.After inserting the expression for the evaporation rate as given in Eq. (A4a) into Eq.( A2) and rearranging, we get for an unsaturated canopy (S i < S i,cap , D i =0): S i,cap and Assuming that β is non-zero (see the above Case 1 for β = 0) and that the interception reservoir does not become full (S i < S i,cap , D i = 0), the solution to the differential equation for S i (t) is given by: If the reservoir does not become full (S i < S i,cap ), the evaporation is given by Eq. ( A6).
If S i (Eq.A10a) yields a value > S i,cap , then we first determine the time at which this happens by solving for t cap ( t cap ≤ t): where t cap is the time to fill the reservoir during which there is no dripping from the canopy.The canopy evaporation can thus be found from a balance like the one given in Eq. (A6); for the remaining part of the interval ( t − t cap ) the evaporation rate is equal to the potential value.So the total canopy evaporation rate can be found from: The dripping rate then follows from Eq. (A8).

Appendix B Parameterization of the crop coefficient as a function of the leaf area index
The K cb (LAI) parameterization is based on Eq. ( 9) (main article text): where K c,tot is the "all in one" crop coefficient given by Feddes (1987) for each ten-day interval of the growing season, K cb (LAI) is the basal transpiration crop coefficient as a function of the LAI, ET 0 is the reference crop evapotranspiration, E s,a is the simulated soil evaporation of an unstressed MetaSWAP-WOFOST run, e p is the positive residual error (≥0) and e n is the negative residual error (≥0).The potential transpiration T p (=K cb (LAI) ET 0 ) plus the actual soil evaporation E s,a is here denoted as ET opt .
The K c,tot coefficients for grassland and potato are presented in Table B1.These coefficients apply to long-term averages of the ten-day intervals within a year, for a climate that is defined by a weather series from 1971-2000.For a wet bare soil, the K ew100 coefficient (Eq.3, main article text) is assumed to be 1.0.Coefficients less than unity in Table B1 are caused by the soil evaporation reduction in the field studies that they were based on.This reduction occurred even though the crops themselves were optimally supplied with water.
The first step in the parameterization procedure is to execute an unstressed MetaSWAP-WOFOST run for the period 1971-2000 using the best Dutch soil with a fixed groundwater level at 1 m b.s.s.The field measurements used for parameterization were obtained with sprinkler irrigation.However, irrigation only took place at the height of the growing season when the bare soil evaporation played a minor role.So the situation used here sufficiently emulates the field conditions.
The simulation yields for each day (dy) and year (yr) values of: -LAI pot (dy, yr), the leaf area index of an unstressed crop, and the value of the exponential extinction function in Eq. (3) (main article text); and -actual soil evaporation E s,a (dy, yr).-the table function K cb (LAI) as an array of 10 000 decision variables, for the LAI in steps of 0.001 from zero to 10; constraints are included for: -the starting point of the function (K cb (0) = 0); -enforcing non-decreasing values; -enforcing a non-increasing first derivative; Table B1.Crop coefficients for grassland and potato, per ten-day interval of the growing season as given by Feddes (1987).These K c,tot coefficients are for the "all-in-one" transpiration, canopy interception evaporation and soil evaporation models.
-calculation of ET opt,Feddes directly from K c,tot times the reference crop evapotranspiration ET 0 , per day of every year , and then averaged per decade of every year, and then averaged per decade over the years;  ; E s,a -actual soil evaporation, T p -potential crop transpiration, E i -canopy interception evaporation.
-calculation of total ET opt,LAI as the sum of K cb (LAIpotX1000) • ET 0 (dy, yr) and E s,a (dy, yr); the values of ET opt,LAI per day of every year are averaged in the same fashion as described for the calculation directly from K c,tot ; -the long-term average of the total ET opt,LAI computed with the K cb (LAI) method is set equal to that computed directly from Feddes-factors time ET 0 ; -the objective function to be minimized is formulated as the sum of the absolute deviations between the values per ten-day interval (values averaged over the years), as computed with the above two methods.
Figure B1 presents the fit for the two calculation methods for a potato crop.Figure B2 shows the plot on a cumulative time scale of the ten-day intervals of a year.The deviation between the two curves is never more than 6 mm d −1 .The long-term averages of the ET-terms are given in Table B2, run 1.
In run 2 (Table B2) the canopy interception evaporation has been added as a model feature, by setting the interception storage capacity to 0.25 mm LAI −1 and the crop factor for interception evaporation equal to that of transpiration.The effect of this addition is primarily that part of the transpiration simulated in run 1 is substituted by interception evaporation, as a consequence of Eq. ( 8) (main article text).A secondary Table B2.Evapotranspiration terms of the simulation runs used for K cb (LAI) parameterization.The meteorological data of De Bilt (NL) for the period 1971-2000 are used, consisting of hourly precipitation data with a mean yearly total of 794 mm yr −1 .The mean reference crop evapotranspiration (Makkink) is 543 mm yr −1 .Where E s,a is the actual soil evaporation plus ponding evaporation, T p is the potential transpiration, E i is the interception evaporation from canopy and ET opt is the total evapotranspiration of an unstressed crop.des (1987) for grassland and values of K cb (LAI) for the year with median precipitation surplus (year 2000, in 1971-2000).

Run
effect is that the soil evaporation is slightly reduced, owing to less rainfall reaching the soil surface.
In run 3 (see Table B2), the interception evaporation crop coefficient is set equal to 1.2 times that of transpiration to account for the reduced resistance of water returning to the atmosphere directly from the canopy surface, compared to water via root uptake.This increased coefficient has two effects: -the interception evaporation is increased, caused by the accelerated emptying of the interception reservoir, and thus increasing the fraction of the rainfall that can be stored on the canopy at some moment; -the transpiration is increased, caused by the accelerated interception evaporation process, leaving more time for the transpiration to be active.For this run, the simulated total evapotranspiration was found to be a few per cent higher than without the interception enabled, as can be seen in Table B2 (run 3).For run 4, the LP-code for calibrating the K cb (LAI) function was first rerun; in this run, the results of run 3 were used to anticipate that addition of the interception evaporation adds a few per cent to ET opt .As can be seen from the last column of Table B2, the ET opt of run 4 approaches that of run 1 within 0.1 %.The partitioning of ET opt into separate terms is illustrated in Fig. B3.
To establish the "static" model coefficients, the following procedures were followed: -for grassland, the year with the median precipitation surplus (2000, of 1971-2000) was used to obtain a representative LAI time series and dependent parameters; Fig. B4

Figure 1 .
Figure 1.Overview of SIMGRO modelling framework and attached models used in this study.

Fig. 1 .
Fig. 1.Overview of SIMGRO modelling framework and attached models used in this study.

Figure 2 . 4 Fig. 2 .
Figure 2. Reduction coefficient for root water uptake, α 2 head and potential transpiration rate T p (after Feddes et a 3 4 Fig. 2. Reduction coefficient for root water uptake, α rw , a function of soil water pressure head and potential transpiration rate T p (after Feddes et al., 1978).

Fig. 4 .
Fig. 4. Fitted relationships between the leaf area index (LAI) and the crop factor K cb (of Makkink reference crop evapotranspiration), for potato and grassland.The shown relationships are for Eq.(1).

Figure 6 .
Figure 6.Measured and simulated phreatic levels for the monitoring point 39BL0020 indicated in Fig. 5.

Fig. 6 .
Fig. 6.Measured and simulated phreatic levels for the monitoring point 39BL0020 indicated in Fig. 5.

Figure 9
Figure 9 presents the simulated LAI obtained with the static and dynamic models.Dips in the LAI values are caused by hay making.The LAI dynamics of the static model are

Figure 9 .
Figure 9. Simulated Leaf Area Index of grassland in evaluation point EP1 (Fig. 5), for three of the variants (Table 3): C_Gr_Statstatic vegetation model current climate, C_Gr_Dynsimulation with dynamic vegetation model, current climate; W + _Gr_Dyn -simulation with dynamic vegetation model, W + -scenario.The simulations are for the year 2003 (warm year, current climate) and 2063 (W + -scenario).

Fig. 9 .
Fig. 9. Simulated leaf area index of grassland at evaluation point EP1 (Fig. 5), for three of the variants (Table 3): C Gr Stat -static vegetation model current climate, C Gr Dyn -simulation with dynamic vegetation model, current climate; W + Gr Dyn -simulation with dynamic vegetation model, W + -scenario.The simulations are for the year 2003 (warm year, current climate) and 2063 (W + -scenario).

Figure 11 .
Figure11.Simulated recharge R = (P -ET a ) along the cross-section AB (Fig.5), for the part that is in use by agriculture.The plot is for the arable land scenarios, for the static/dynamic simulations, for the W + climate scenario.

Fig. 11 .
Fig.11.Simulated recharge R = (P -ET a ) along the cross-section AB (Fig.5) for the part that is in use by agriculture.This plot is for the arable land scenarios, for the static/dynamic simulations, and for the W + climate scenario.

Figure 13 . 5 Fig. 13 .
Figure 13.Relative yield of potatoes along the cross section 2 use by agriculture and has the same podzolic soil type.3 4 5 Fig.14

Fig. 14 .
Fig. 14.Comparison between the K cb (LAI) function from the present study with those from Allen et al. (1998) and Jayanthi et al. (2007).

Fig. B1 .
Fig. B1.Total evapotranspiration of unstressed potato (ET opt ) calculated directly from the "all-in-one" K c,tot factor multiplied by ET 0 and calculated with K cb (LAI) multiplied by ET 0 plus the actual soil evaporation E s,a , per ten-day interval, averaged over 30 yr. Fig.B2

Fig. B2 .
Fig. B2.Cumulative total evapotranspiration of unstressed potato (ET opt ) calculated with Feddes-factor multiplied by ET 0 , and calculated with K cb (LAI) multiplied by ET 0 plus the actual soil evaporation E s,a , per ten-day interval, averaged over 30 yr.

Fig. B4 .
Fig. B4.Comparison between K c,tot coefficients as given byFeddes (1987) for grassland and values of K cb (LAI) for the year with median precipitation surplus(year 2000, in 1971-2000).

Fig. B5 .
Fig. B5.Comparison between K c,tot coefficients as given byFeddes (1987) and values of K cb (LAI) averaged for ten-day intervals of a calendar year.

Table 1 .
Land use in the Kromme Rijn region.

Table 2 .
Climate scenarios for the Netherlands for two possible developments in global temperature change for 2050 ( T 2050 ) and two possible developments of the atmospheric circulation ( A).The given changes are for the (climatic) mean precipitation ( P ) and the Makkink reference crop evapotranspiration ( ET 0 ).For this study the scenario W + has been used.

Table 3 .
Investigated variants of climate and land-use scenarios without/with dynamically simulated crop growth.

van Walsum and I. Supit: Influence of ecohydrologic feedbacks on hydrologic simulationsTable 4 .
Water balance terms (mm yr −1 ) of monitoring point EP1 in Fig.5(+/− = in/out) for the simulation year of 2003/2063.The listed terms are: P -precipitation; ET 0 -Makkink reference crop evapotranspiration; T p -potential transpiration; T a -actual transpiration; E i -interception evaporation; E s -bare soil evaporation; T rel -mean relative transpiration; ET a -total evapotranspiration; R -recharge (P − ET a ); R veg -change of simulated recharge of dynamic versus static vegetation model; and R cl -change of simulated recharge with respect to current climate.

Table 3
): C Gr Stat -static vegetation model current climate, C Gr Dyn -simulation with dynamic vegetation model, current climate; W + Gr Dyn -simulation with dynamic vegetation model, W + -scenario.The simulations are for the year 2003 (warm year, current climate) and 2063 (W + -scenario).

Table 5 .
Jayanthi et al. (2007)potranspiration terms simulated with the K cb (LAI) function from the present study with those fromAllen et al. (1998)andJayanthi et al. (2007).The results are longterm averages for 1971-2000 (De Bilt, NL) for a crop that is optimally supplied with water (T a = T p ). Canopy interception evaporation has been disabled for the comparison.Used symbols: E s,a -actual soil evaporation; T a/p -actual/potential transpiration; and ET opt -total evapotranspiration of unstressed crop.
Method for K cb (LAI) E s,a (mm yr −1 ) T p (mm yr −1 ) ET opt (mm yr −1 ) These results are not significantly influenced by the initially assumed K cb (LAI) table function.LAI pot(dy, yr)values are multiplied by 1000 and rounded to the nearest integer, here denoted as LAIpotX1000.These integers are used as tableindices in realizations of the table function K cb (LAI) in the LP-code that is used to determine the optimal fit.The LPcode has been implemented with the Xpress Optimization Suite (www.fico.com), in the Mosel language 1 .It contains the following main elements: K cb (LAI) of Crop E s,aT p E i ET opt ET opt / LP-run nr (mm yr −1 ) (mm yr −1 ) (mm yr −1 ) (mm yr −1 ) ET opt,run1 (-) K c,tot (Feddes, 1987)K cb (LAI) median year

van Walsum and I. Supit: Influence of ecohydrologic feedbacks on hydrologic simulations -for
potato, run 4 of TableB2was used to obtain daily average values of coefficients; to compare the "static" model coefficients to the K c,tot -values given by Feddes (1987), the coefficients have been averaged for each ten-day interval of a year, as shown in Fig.B5.