Spatio-temporal impact of climate change on the groundwater system

Given the importance of groundwater for food production and drinking water supply, but also for the survival of groundwater dependent terrestrial ecosystems (GWDTEs) it is essential to assess the impact of climate change on this freshwater resource. In this paper we study with high temporal and spatial resolution the impact of 28 climate change scenarios on the groundwater system of a lowland catchment in Belgium. Our results show for the scenario period 2070–2101 compared with the reference period 1960– 1991, a change in annual groundwater recharge between −20 % and +7 %. On average annual groundwater recharge decreases 7 %. In most scenarios the recharge increases during winter but decreases during summer. The altered recharge patterns cause the groundwater level to decrease significantly from September to January. On average the groundwater level decreases about 7 cm with a standard deviation between the scenarios of 5 cm. Groundwater levels in interfluves and upstream areas are more sensitive to climate change than groundwater levels in the river valley. Groundwater discharge to GWDTEs is expected to decrease during late summer and autumn as much as 10 %, though the discharge remains at reference-period level during winter and early spring. As GWDTEs are strongly influenced by temporal dynamics of the groundwater system, close monitoring of groundwater and implementation of adaptive management measures are required to prevent ecological loss.


Introduction
There is little doubt that the ongoing climate change will significantly influence the hydrological cycle worldwide (Kundzewicz et al., 2008;Maxwell and Kollet, 2008).Current observations show that already at this moment climate changes are influencing hydrological processes in certain areas (Rosenzweig et al., 2007;Kundzewicz and Döll, 2009).As the IPCC predicts that global atmospheric concentrations of greenhouse gases will continue to rise, it is expected that climate change will continue in the future (Solomon et al., 2007).Freshwater resources are among those systems that are particularly vulnerable to changes in climate (Solomon et al., 2007).
Recent research (Feyen and Dankers, 2009) showed that global warming is likely to amplify drought events over Europe.Especially during drought events groundwater is of vital importance for availability of water for food production and drinking water.Groundwater plays a vital role in maintaining the ecological value of many areas (Solomon et al., 2007;UN WWAP, 2009).Because groundwater is less visible and has a more complex relationship with the climate than surface water bodies it has been studied less than surface water bodies up till now (Kundzewicz and Döll, 2009;Scibek et al., 2007).However, there is an increasing awareness to protect the groundwater resources and to assess the impact of future land-use and climate changes (Solomon et al., 2007;Green et al., 2011).
In order to assess the impact of climate change on the groundwater system there is a need for reliable climate change scenarios and consistent methods to simulate water fluxes recharging and discharging the groundwater system.The uncertainty on climate change forecasts is still very high due to uncertainties in the future world visions, influencing for example the emissions of greenhouse gas, land-use changes, etc. and uncertainties caused by the General and Regional Circulation Models (GCMs and RCMs) (Murphy et al., 2004).In order to optimally incorporate the current knowledge on climate change, Kundzewicz et al. (2008) and Allen et al. (2010) suggest a joint analysis of ensembles of climate models driven by multiple emission scenarios.Hendricks Franssen (2009) emphasizes the importance of downscaling of future precipitation from GCMs for impact assessments on hydrology.Stoll et al. (2011) stress the large uncertainty introduced by those downscaling methods.
Previous studies show a large variety in complexity of approaches to simulate groundwater recharge.For example Chen et al. (2002), Hsu et al. (2006) and Serrat-Capdevila et al. (2007) apply a simple linear function including precipitation and temperature to simulate groundwater recharge, while Woldeamlak et al. (2007), Jyrkama and Sykes (2007), van Roosmalen et al. (2009), McCallum et al. (2010) amongst others apply a more complex approach.Holman (2006), Jyrkama and Sykes (2007), Hendricks Franssen (2009), Ferguson and Maxwell (2010) and Holman et al. (2011) advise a physically based approach that accounts for spatial and temporal variation of surface and subsurface properties of the study basin when simulating the impact of climate change on groundwater recharge.A majority of the current studies assessing the impact of climate change on the groundwater system estimate the impact on the annual or seasonal average spatially distributed recharge, e.g.: Dickinson et al. (2004), Scibek and Allen (2006), Scibek et al. (2007), Serrat-Capdevila et al. (2007) and Woldeamlak et al. (2007).Woldeamlak et al. (2007) stated that climate change impact studies based on steady-state groundwater simulation have limitations in representing boundary conditions and can only be used for assessing sensitivities.A few recent studies have applied transient methods to estimate the impact of climate changes on the groundwater system (Goderniaux et al., 2009; van Roosmalen et al., 2009;Goderniaux et al., 2011;Jackson et al., 2011;Stoll et al., 2011).However, van Roosmalen et al. (2009), Goderniaux et al. (2009), Goderniaux et al. (2011), Jackson et al. (2011) and Stoll et al. (2011) limit the analysis of the transient results to the predicted head change in some observation wells or to the average change in groundwater head of the basin without analyses of the spatial variation.Nevertheless, the groundwater dynamics within a year is of major importance for groundwater dependent terrestrial ecosystems (GWDTEs) (Naumburg et al., 2005).Groundwater dependent vegetations along with riverine landscapes have an important ecological function (Naumburg et al., 2005;Kløve et al., 2011) and should therefore be protected.Applying highly dynamic models also allows including more accurately changes in precipitation intensity and number of dry and wet days projections due to climate change.Precipitation intensity and number of wet and dry days have an important impact on the soil moisture content and consequently influence strongly the groundwater recharge.
This is one of the first studies analyzing the intra-annual response of a groundwater system to climate changes.These intra-annual changes determine the status of the groundwater resources as well as site conditions of GWDTEs (Naumburg et al., 2005).The climate for the reference period, 1960-1991, is compared with climate scenarios, predicted for 2071-2100.Due to the high variability of climate change predictions between different climate change models, an ensemble of 28 climate change scenarios is chosen from the European project PRUDENCE (Christensen and Christensen, 2007).By applying this ensemble of climate change models we obtain uncertainty bounds on the impacts of the climate change on the groundwater system.We limit the study to climate change impacts, disregarding other expected changes such as land-use change (Dams et al., 2008).
The Kleine Nete basin, situated in Belgium, was chosen as a study area.Due to the sandy soils and low slopes a large fraction of the effective rainfall in the basin percolates to the groundwater.The groundwater in the basin is extensively used for drinking water supply, and hosts important groundwater dependent wetlands.An impact assessment is therefore required to assess whether adaptive measures are essential to protect the groundwater system and related groundwater dependent natural vegetations from expected climate changes.

Study area
The study area is the Kleine Nete basin, which is a subbasin of the Scheldt basin (Fig. 1).The Kleine Nete basin has an area of 581 km 2 .The elevation ranges from 3 to 48 mTAW (Belgian reference height above sea level), the average slope is about 0.4 %.Interfluves are slightly elevated, the valleys broad and swampy.The dominant soil texture in the basin is sand, though in the valleys some loamy sand, sandy loam and sandy clay is present.The region has a temperate climate characterized by a warm summer and a cool winter with little snowfall.The average annual precipitation during the period 1960-1991 was 828 mm with a standard deviation of 136 mm.Precipitation is nearly equally distributed throughout the year and the different raingauges, indicated in Fig. 1, show similar annual precipitation amounts.Over the same period 1960-1991 the estimated average annual potential evapotranspiration (PET) is 664 mm with a standard deviation of 47 mm.hydrostratigraphy of the study area comprises the Miocene aquifer, the Pliocene clay layer, the Pleistocene and Pliocene aquifer, the Campine clay-sand-complex and the Holocene eolian sand aquifer.An overview of the formations is given in Table 1.Only the Miocene aquifer and the Holocene eolian sand aquifer are found throughout the basin, other hydrostratigraphic units are discontinuous as shown in Fig. 2. Figure 3 shows a 3-D view of the geological layers along a cross-section over the area.The Miocene aquifer has an average thickness of about 187 m and in the eastern part of the basin this aquifer reaches a maximum thickness of 410 m (Wouters and Vandenberghe, 1994).
The land-cover in the study area consists mainly of agricultural fields including meadows (60 %), coniferous and mixed forest (20 %) and urban areas (10 %).Groundwater is extensively used in the basin, in total there are 565 wells which extract a total of 54 291 m 3 day −1 of which about 30 200 m 3 day −1 is extracted by a single water production Fig. 2. Occurrence of the Campine clay-sand complex (0220) and Pliocene clay layer (0240), as described in Table 1.Other layers occur throughout the basin.The profile B-B'-B" is presented in Fig. 3. company for drinking water supply.Most important pumping wells are indicated in Fig. 1.
Within the Kleine Nete catchment several ecologically important areas are protected by the European Natura2000 network, set up for the protection of Europe's most vulnerable habitats.Several of these habitats depend largely on oligotrophic and mesotrophic site conditions, influenced by groundwater flow conditions.Typical habitats are Northern Fig. 3. Cross-section along profile B-B'-B" presented in Fig. 2 showing the different Tertiary formations.The HCOV aquifer codes are given in Table 1.

Data and method of analysis
This study compares the groundwater characteristics of a lowland watershed in Belgium for the reference period 1960-1991 with those subject to climate change conditions for the period 2070-2101.Figure 4 shows a conceptual overview of the applied spatial-temporal methodology.An ensemble of 28 climate change scenarios derived from multiple GCMs and RCMs and driven by multiple greenhouse emission scenarios is applied.

Climate change
Climate change scenarios are obtained from the PRUDENCE database and combine several GCMs: ECHAM4/OPYC, HadAM3H, HadAM3P, ARPEGE and HadCM3 that drive a range of RCMs: RCAO, RACMO, HIRAM, CHRM, HadRM3P, REMO, ARPEGE, CLM and PROMES (Christensen and Christensen, 2007).Baseline assessments for all scenarios used in this study have been performed by Baguis et al. (2010) and Ntegeka et al. (2008).All these scenarios are based on the A2 and B2 SRES greenhouse gas emission scenarios of the Intergovernmental Panel for Climate Change (IPCC) (Nakicenovic and Swart, 2000).In total, projections from 28 climate models runs were statistically analyzed, comparing the daily simulation results for precipitation and PET between the control period 1960-1991 and the scenario period 2070-2101.The precipitation results were obtained directly from the RCM outputs, and PET was calculated by Baguis et al. (2010) using the Penman equation based on the RCM outputs of mean sea level pressure, net terrestrial radiation, total solar radiation, cloud cover, temperature at 2 m, wind at 10 m and humidity.A change factor method (Deque et al., 2007;Diaz-Nieto and Wilby, 2005; The figure shows a watershed discretized using a rectilinear grid, surface water bodies are represented in blue.For every cell all water balance components are simulated daily and the runoff, interflow and groundwater flow are routed to the outlet of the catchment.Recharge and discharge are aggregated to half-monthly time steps.Two cells in this figure are highlighted.Cell A represents a typical groundwater discharge area: during most time steps the groundwater in this cell flows from the groundwater system towards the land surface where the groundwater can discharge to the surface water bodies or be used for evapotranspiration.Cell B represents a typical recharge area where the water table is recharged by water infiltrating from the land surface.The graphs on the right show how the groundwater discharge or recharge flux typically evolve over time.In this study the groundwater system is simulated for the reference condition (grey line) and several climate change scenarios (e.g.orange line).Prudhomme et al., 2002) was applied to downscale the RCM outputs.The change factor method involves the calculation of climate change signals from the RCM outputs, which are applied to an observed series to generate a climate change series.The main advantage of this method is that the climate change signals are transmitted robustly, irrespectively of absolute biases of individual models (Buonomo et al., 2007;Christensen and Christensen, 2007).
An adapted version of the traditional change factor method that takes into account changes in number of wet days and rainfall variability is applied in this study for the perturbation of rainfall series.To incorporate the change in variability, a quantile perturbation factor that accounts for the differences in rainfall intensities is used (Olsson et al., 2009).The climate change signal is calculated based on the difference between the cumulative distribution of the control and future RCM scenario.The difference in rainfall variability between the control and climate change scenario is applied on the observed historical series by multiplying the rainfall in each wet day by a unique change factor that is based on the exceedance probability of the rainfall intensity during that day.This method is applied separately for each month (all days in a month for all 32 yr), which results in an intensity perturbated observed series.Secondly, the observed series are adapted for the number of wet days.Again the control and climate change scenario are compared, this time for number of wet days (>0.1 mm).According to the result of the comparison, wet days are added or removed from the observed series based on a two-state second order Markov chain process.If the number of wet days increases, the four transition probabilities for wetness (dry-wet-wet, wet-dry-wet, wet-wet-wet and dry-dry-wet) are calculated on a monthly basis from the control and scenario series.The changes in transition probabilities are then calculated as a ratio of the control transition probabilities to the scenario transition probabilities.This factor is multiplied with the transition probabilities of the intensity perturbated observed series.The obtained transition probabilities are used to weight the probability of a dry day to be converted to a wet day during the random selection of new wet days.The intensity of the new wet days is estimated from a mixed exponential distribution as suggested by Wilks (1999) and Willems (2000).The above mentioned process is repeated one hundred times, from the generated series the series with the closest match of change in coefficient of variability with the considered climate model is selected.More details on the frequency-and quantile-perturbation procedure can be found in Ntegeka and Willems (2008) and Willems and Vrac (2011).
Because for PET no significant variations were found between the change factors for the different return periods, a simpler approach was followed.Monthly change factors for PET are calculated from a comparison of the control and scenario series.The observed PET series are multiplied by these change factors to obtain the perturbated PET series.

Groundwater system modeling
The impact of climate change on the groundwater system is simulated by applying a coupled WetSpa -MODFLOW approach.WetSpa (Liu et al., 2003), a physically based distributed hydrological model, simulates with a daily time step the river discharge at the outlet of the basin and the groundwater recharge for each 50 by 50 m model cell in the watershed.WetSpa updates the root zone water balance for all model cells during each timestep (Safari et al., 2011): where The evapotranspiration flux includes, evaporation, transpiration from the root zone and direct uptake of groundwater by plants.The routing scheme of the WetSpa model accounts for the spatial differences by integrating the slope, soil and landuse characteristics of the grid cells along the flow path.Based on the spatial characteristics of grid cells along the flow path the runoff, interflow and groundwater flow are routed directly from the grid cell to the watershed outlet.The rate of percolation (R rate ) or groundwater recharge in the WetSpa model is derived through the Brooks and Corey relationship (Brooks and Corey, 1964): where water content at saturation, θ r [L 3 L −3 ] residual soil moisture content, and B [-] is the soil pore size distribution index.The soil pore size distribution index B is obtained from an empirically derived univariate regression, based on the percentage of clay content (Cosby et al., 1984).The Brooks and Corey relationship assumes the pressure potential gradient in the soil to be zero, causing the percolation to be controlled by gravity only.Furthermore, the applied methodology assumes the percolation out of the root zone to pass directly to the groundwater reservoir without considering any delay.Nyenje and Batelaan (2009) applied WetSpa to assess the impact of climate change on groundwater recharge and baseflow.Daily spatially distributed recharge results are aggregated over half monthly periods to be compatible with the MOD-FLOW time step.Additionally, the results of a hydraulic model for the main rivers in the basin are used to obtain half monthly average river heads for every 50 m transect of those rivers, based on WetSpa simulated river discharge at the basin outlet.
The groundwater flow model MODFLOW (Harbaugh and McDonald, 2000) simulates the effect of the climate induced changes in river head and groundwater recharge on the groundwater level and flux.Previous groundwater flow models that included the Kleine Nete basin used a variety of concepts to represent the geological units playing a role in the local hydrology of the study basin: Batelaan and De Smedt (2004) and Gedeon et al. (2007) conceptualized all geological units above the Boom clay as one layer while Verbeiren et al. (2006) used five layers according to the geological units listed in Table 1.In both cases satisfying calibration results were obtained.Woldeamlak et al. (2007) applied a two layer conceptualization in which the top layer represents the Pleistocene and Pliocene aquifer, Holocene eolian sand deposits, Pliocene clay layer and Campine clay-sand complex and the bottom layer represents the Miocene aquifer.To incorporate the spatial variability of hydraulic conductivity and specific yield within the upper layer, the layer is sub-divided into seven zones.The initial horizontal and vertical conductivity in these zones are calculated using the weighted arithmetic and harmonic mean, respectively.Independent calibration multipliers were applied for each of the zones in the layer.The two layer model (Woldeamlak et al., 2007) gave similar calibration results than the five layer model of Verbeiren et al. (2006).Furthermore, the numerical stability of the two layer model is significantly better then the five layer model, which sometimes fails to converge during dry periods.Therefore, the conceptualization of the geological units for the MODFLOW model in this study is chosen as suggested by Woldeamlak et al. (2007) and represents a compromise between geological detail and optimal hydrogeological model conceptualization.The watershed boundaries of the applied model are assumed to be no-flow boundaries.All major rivers, canals and lakes are simulated as internal boundaries and parameterized with the RIVER package.The RIVER package controls the flux exchanged between the groundwater system and the river, based on the river stage, the elevation of the bottom of the riverbed, the riverbed hydraulic conductance and the hydraulic head calculated for the particular model cell containing the surface-water feature.The river stage of the main rivers, in the MODFLOW model, is based on an upstream river profile simulated by a hydraulic model according to the river discharge at the outlet calculated by WetSpa.The groundwater drainage from ditches, small streams and wetlands is simulated using the DRAIN package (Batelaan and De Smedt, 2004).Drain cells are defined for the whole model area except for the cells where river conditions are defined.Drain simulated by MODFLOW represents flow from the saturated zone towards the land surface.The drain flow simulated by the DRAIN package depends on the drainage level and conductance.The drainage level is set to the deepest location in the soil profiles where oxidation appears as suggested by Stuurman et al. (2002).
WetSpa is calibrated manually on seven global parameters: interflow scaling factor (i), groundwater recession coefficient (ii), initial soil moisture content (iii), initial active groundwater storage (iv), maximum active groundwater storage (v), moisture or surface runoff exponent (vi) and maximum rainfall intensity (vii).The calibration period for WetSpa is from 1 January 1992-31 December 2001 and val-idation period from 1 January 2002-31 December 2004.The calibration parameters for the MODFLOW model were selected based on a sensitivity analyses of a steady-state version of the model.The most sensitive parameters were used for the calibration of the transient model and include: the horizontal hydraulic conductivity of different zones in the two layers that are based on the occurrence of Tertiary formations as represented in Fig. 2 (i-vi), the river conductance of northern and southern part of the river (vii-viii) and the drain conductance (iv).The MODFLOW model is calibrated using 10 226 head observations measured between 1 January 1992 and 31 December 2001 from 113 observation wells (Fig. 1) more or less equally distributed over the basin.To make optimal use of the available groundwater head observations (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001) all observations are used to calibrate the MOD-FLOW model.Initial calibration parameters were obtained from an automatic calibration of a steady state version of the MODFLOW model using PEST (Doherty and Johnston, 2003).The transient MODFLOW model was further manually calibrated.Because at the moment there is no automatic coupling between WetSpa and MODFLOW, the models are calibrated separately while being loosely coupled.Because the ability of the models to simulate groundwater recharge and discharge is important in this paper, the baseflow is integrated in the calibration procedure.A measured baseflow time series is extracted using the baseflow filter developed by Arnold and Allen (1999).

Groundwater level and flux analyses
The mean highest groundwater level (MHGL), mean lowest groundwater level (MLGL) and mean spring groundwater level (MSGL) are calculated respectively as the three highest, the three lowest and the three groundwater level measurements around the 1st April per year, based on two weekly measurements, and averaging these values over at least eight years (Van der Sluijs and De Gruijter, 1985).In this study the MHGL, MLGL and MSGL for each model cell are estimated based on the half monthly groundwater level simulated by MODFLOW.In order to reduce the effect of the initial conditions, results of the first six half-monthly time steps are not used.The groundwater discharge frequency of each 50 by 50 m MODFLOW cell is calculated as the percentage of time steps in which a groundwater discharge to the soil surface is simulated by the DRAIN or RIVER package.

Model calibration and validation
The WetSpa model was calibrated using measured river discharges and estimated baseflow at the catchment outlet.During calibration a Nash-Sutcliff efficiency of 73 % was obtained for the river discharge and 87 % for the baseflow.For the validation period the Nash-Sutcliff efficiency is 74% for the total discharge and 57 % for the baseflow.Figure 5 compares the baseflow extracted by the baseflow filter with the simulated baseflow of WetSpa and MODFLOW obtained during the calibration period.It is shown that the baseflow simulated with the WetSpa model is very similar to the baseflow derived from the baseflow filter.The MODFLOW model, while using the WetSpa simulated recharge, tends to underestimate high baseflows.After calibration the MOD-FLOW model has an average bias between observed and simulated hydraulic head of −0.03 m, a mean absolute error of 0.59 m and a root mean square error of 0.81 m. Figure 6 shows a scatter-plot of all measured versus simulated groundwater heads.The Nash-Sutcliff efficiency of the baseflow calculated by MODFLOW compared to the filtered baseflow is 0.76 [-].

Intra-annual impact of climate change on groundwater characteristics
Figure 7a shows the reference and forecasted average precipitation within the basin for each time step.The total annual precipitation decreases on average by 50 mm, from 821 to 771 mm yr −1 with a standard deviation of 35 mm between the scenarios.As observed from Fig. 7a the change in precipitation varies in time: from October to April the precipitation increases on average 50 mm but from May to September the precipitation decreases about 100 mm. Figure 7b illustrates the projected intra-annual change in PET, obtained from averaging the 32 yr simulation period.The average yearly PET of 664 mm yr −1 measured during 1960-1991 is predicted to increase almost 30 % with a standard deviation between the scenarios of 91 mm yr −1 .The increase in PET occurs almost completely between April and October.Due to restrictions on the water availability the actual evapotranspiration (AET) is only a fraction of the PET.The intra-annual variability of the AET is represented in Fig. 7c.During winter the scenarios predict no significant changes in AET.In spring we notice a small increase in AET, however the largest changes are observed for the summer period where a significant decrease in AET is predicted.The largest decrease in AET is predicted for August where the average AET of the scenarios is 30 % below the AET of the reference scenario.The average yearly AET is predicted to decrease by 3.5 % from 532 mm yr −1 for the reference scenario to 513 mm yr −1 for the average of the climate change scenarios with a standard deviation between the scenarios of 23 mm yr −1 .The difference between the precipitation (Fig. 7a) and the AET (Fig. 7c) is available for surface runoff or groundwater recharge.The mean average half-monthly runoff represented in Fig. 7d is the sum of the interflow and surface runoff simulated by WetSpa.The results show that the runoff increases during winter and decreases during summer.The average yearly runoff decreases with 6.5 % from 103 mm yr −1 to 96 mm yr −1 with a standard deviation between the scenarios of 5.7 mm yr −1 .Figure 7e represents the average intra-annual variability of soil moisture storage, one of the controlling factors of the percolation rate.During winter season the average soil moisture storage is around 90 mm for both reference and future scenarios.During summer however, the average soil moisture storage decreases to around 65 mm for the reference period, while for the future scenarios the soil moisture storage decreases to around 55 mm. Figure 7f illustrates the average annual groundwater recharge pattern.The already low groundwater recharge during summer decreases further due to higher evapotranspiration and lower precipitation, on the other hand additional precipitation during winter causes the groundwater recharge to increase.On average the groundwater recharge is predicted to decrease about 40 mm during summer and increase about 20 mm during winter, resulting in an annual decrease from 278 to 258 mm yr −1 or 7.2 %.The standard deviation of the change in yearly groundwater recharge calculated from the different climate scenarios is 20 mm.
Figure 8a shows the intra-annual variability of the groundwater head for the reference and climate change scenarios.In April the average simulated future groundwater head is close to the average reference head.During summer however, the climate scenarios predict a larger seasonal groundwater storage depletion.The maximum average groundwater depth simulated for the reference period is 2.2 m below the topography and is reached during the first half of September, while the maximum average future groundwater depth predicted by the climate scenarios is 2.3 m and occurs later at the end of September.The timing of the minimum average groundwater depth also shifts, from late December-January to early February-late March.The maximum average difference in simulated groundwater depth between the reference period and the future scenarios occurs in November, when the average simulated future groundwater depth increases with about 15 cm.Over the entire year the average reference simulated groundwater head declines about 7 cm.Figure 8b displays the average intra-annual groundwater discharge simulated for the reference period and future scenarios.Similar to the groundwater head we observe that the average future groundwater discharge decreases at a faster rate during summer compared to the reference scenario, however, also the increase of the average groundwater discharge during autumn is more profound for the future scenarios.The average future groundwater discharge from February until May fluctuates around the groundwater discharge simulated for the reference climate.On the other hand, from August until December the groundwater discharge is predicted to decrease by more than 10 %.

Spatial impact of climate change on average groundwater recharge and extreme groundwater heads
Figure 9 shows the spatial distribution of the temporally averaged change in groundwater recharge between the average of the climate change scenarios and the reference scenario.The figure clearly shows that in this basin the change in groundwater recharge is largely determined by the land-cover and soil texture.Urban land-use seems more sensitive to the climate changes in comparison to agricultural land and forest.Areas with a sandy loam soil texture show larger decreases in groundwater recharge than soils with a sand texture.Open waters in Fig. 9 are shown in blue because in these areas the groundwater-surface water interactions are simulated by the MODFLOW model.Figure 10 illustrates the spatial impact of the climate change scenarios on the groundwater head.Compared with the reference scenario the groundwater head decreases most on the interfluves and near the fringes of the watershed where the average groundwater level can be as much as 30 cm lower.In the valleys the average groundwater level decrease is generally less than 5 cm (Fig. 10a).For GWDTEs, especially the yearly extreme groundwater depths (MHGL,  Fig. 10b and MLGL, Fig. 10c), and the MSGL (Fig. 10d) influence the plant species distribution.Both the MLGL and the MHGL show a generally decreasing trend.Similar to the average groundwater levels, the interfluves are more sensitive and show the greatest decrease in yearly extreme groundwater levels.From Fig. 10b-d we notice that the largest decrease is obtained for the MLGL, for which an average decrease of 6 cm is simulated, with a standard deviation of 3 cm between the scenarios.The MHGL and MSGL decrease on average 3 and 1 cm, respectively.The standard deviation between the different scenarios is 5 cm for both MHGL and MSGL.

Impact of climate change on groundwater discharge
Climatic conditions also influence the groundwater interaction with surface water and groundwater discharge towards the land surface.Because sufficient groundwater exfiltration is crucial for the presence of GWDTEs changes due to climate change are further analyzed in this section.Figure 11 plots the average change in groundwater discharge flux simulated for the different scenarios versus the reference groundwater discharge flux of the same grid cells.The scenarios predict for most cells a decrease in average groundwater discharge.Figure 11 also shows that the maximum decrease is about 50 %.The highest decrease in groundwater discharge occurs in cells with a reference groundwater discharge flux of less than 1 mm d −1 where the average decrease is about 15 %.Groundwater discharge cells with a reference flux between 1 and 10 mm d −1 seem to be buffered quite well to the predicted climate changes.The average total groundwater discharge from the basin decreases from 5.0 m 3 s −1 (reference scenario) to 4.8 m 3 s −1 (average of the future scenarios).
In addition to the magnitude of groundwater discharge, the temporal availability of groundwater is important for GWDTEs.Figure 12 plots the change in groundwater discharge frequency versus the reference groundwater discharge frequency.The groundwater discharge frequency of a cell is the temporal frequency that groundwater discharge occurs from this cell.Figure 12 shows that there is an average decrease in the frequency of groundwater discharge.For cells with a reference groundwater discharge frequency between 0 and 80 % the climate change scenarios predict an average reduction of the frequency of about 20 %.Cells with groundwater discharge frequencies above 80 % have a lower sensitivity that further reduces towards permanent discharge cells.

Discussion of results
Both the calibration and validation results of the WetSpa runoff and calibration results of MODFLOW groundwater heads are satisfying.A closer analysis of the MODFLOW results for the different observation wells revealed that in some cases there is a systematic shift between measured and simulated groundwater head, although the groundwater dynamics is modeled relatively well.This systematic shift could be caused by differences in reference head between the measurements and the groundwater model.The systematical underestimation of the winter baseflow by MODFLOW (Fig. 5) can be explained by the fact that the drain output simulated by MODFLOW is not added to the baseflow but is assumed to be lost by ET.The baseflow from MODFLOW was calculated as the net river leakage from the groundwater system towards the river.During summer it is very likely that a major part of the drain is indeed lost by ET, however, during winter when PET is low it is likely that a significant part of the output of the drains contribute to the baseflow of the river.
The climate change scenarios predict that more than in the reference scenario the availability of water will be the limiting factor for (eco-)hydrological processes in summer.Results shown in Fig. 7 show that although the PET is predicted to increase significantly in the future the AET is simulated to decrease due to the decrease in summer precipitation.The decrease in precipitation during summer, although slightly compensated by less ET, causes a stronger depletion of the soil moisture storage and a decrease of the summer runoff and groundwater recharge.
Figure 8 shows that the groundwater system is less variable than the highly dynamic precipitation and groundwater recharge, due to the moderating effect of the flow systems.The changes in groundwater recharge and river head cause the groundwater head and discharge to decrease especially towards the end of the summer.This predicted reduction in groundwater availability and groundwater discharge could cause an increased risk of water shortage towards the end of the summer in the future.
Figure 9 shows that groundwater recharge in urban areas is most vulnerable to climate change.A land-use map of the catchment can be found in Dams et al. (2008).The relatively large decrease in groundwater recharge in urban areas is mainly caused by a significantly higher increase in AET simulated by the WetSpa model due to the predicted climate change, in comparison to non urban areas.An explanation of the increase in AET is that because of the rise in PET predicted by the climate change scenarios a larger fraction of the water stored in the relatively large depression storage of the urban cell is actually evaporated, causing an increase in AET.However, because parameterization of urban areas in hydrological models is difficult due to large heterogeneity in for example impervious surfaces density, drainage density etc. further research is required to verify the impact of climate changes on the urban hydrology.
Comparing the spatial distribution of changes in groundwater recharge and head due to climate change, we can conclude that the groundwater level change is mainly topographically controlled but amplified by the position of the urban areas which are generally also located on the interfluves.As shown in Fig. 10 for both the temporally average groundwater as well as the MHGL, MLGL and MSGL the largest decrease in groundwater head occurs upstream and on the interfluves and less in valley's and wetland areas where GWDTE's are located.The interfluves and catchment fringes are most sensitive to changes in groundwater recharge because their groundwater level is primarily controlled by recharge.In valleys, on the other hand, also the groundwater flow from upstream areas and river head boundary conditions have an important impact on the groundwater level.
The buffering of groundwater discharge cells with a flux between 1-10 mm d −1 (Fig. 11) can be explained by the fact that cells with a small groundwater discharge flux (<1 mm d −1 ) are situated on the borders of the groundwater discharge areas.Due to climate change the discharge zones will become smaller in most scenarios and as such lead to drastic changes in the groundwater discharge at the fringe of these areas.Within the center of the wetlands where the groundwater discharge is higher, MODFLOW predicts a similar decrease of groundwater discharge for all cells.This is due to the fact that the changes in recharge of the upstream infiltration zones are spatially averaged.Closer to the river where the highest groundwater discharges are found the groundwater discharge is also influenced by the river heads which causes more diverse changes.
The impact of the climate change scenarios on the groundwater discharge frequency (Fig. 12) shows that the highest changes in frequency occur in cells that during the reference period change from groundwater recharge cells (generally during summer) to groundwater discharge cells (generally during winter).The majority of the cells that are infiltrating water during all time steps also remain groundwater recharge cells in the future.The same applies for cells that are groundwater discharge zones throughout the year.Cells with a reference groundwater discharge frequency above 80 % seem to have a lower sensitivity which could be explained by the relatively large amount of discharge in these cells which might be increased or decreased but will nevertheless remain discharging groundwater to the soil surface and as such not influence the frequency.For a large number of cells with a low reference groundwater discharge frequency, a small increase in discharge frequency is simulated under the climate change conditions (Fig. 12).This increase is caused by the fact that the presented future groundwater discharge frequency is the average of the frequencies simulated for the different climate scenarios.A few scenarios predict a significant increase in groundwater discharge area.Although after averaging the groundwater discharge frequency for the different scenarios the frequency becomes very small, the frequency is still larger than zero.

Discussion of methodology
Hydrological impact assessments are subject to large uncertainties.The uncertainty of the climate change scenarios is incorporated in this study by applying a range of scenarios from the PRUDENCE database.It should, however, be noted that all scenarios in the PRUDENCE database are based on the A2 and B2 emission scenarios and therefore not cover the whole range of scenarios as set out by the IPCC (Christensen and Christensen, 2007).The integrated scenarios also include only the direct consequences of changes in PET and precipitation, indirect consequences of the projected climate change such as increased future groundwater extraction during dry summers, changes in greywater recycling or rainfall water harvesting, etc. are not taken into account in this study.Also uncertainties of the hydrological models are not included in this study.Assessing both the individual and combined impact of the different sources of uncertainty is important for the selection of appropriate mitigation and adaption measures.
The one-way coupling between the WetSpa model that simulates the groundwater recharge and river head, and the MODFLOW model has both advantages and disadvantages.The most important reason for choosing the one-way coupling in this study is that the calculation time is significantly lower in comparison to a fully integrated surface-subsurface model.To investigate the spatial and temporal impact of climate change on the groundwater system both a high spatial and temporal resolution is required and simulations should be performed on basin scale.In order to incorporate climate variability time series of at least 30 yr should be applied, which leads to considerable calculation times.Furthermore, to incorporate the climate change uncertainty, simulations should be repeated for different scenarios, in our case 28 times.The disadvantage of one-way coupled surface and sub-surface hydrological models compared to fully coupled models is that the groundwater recharge calculation is not updated by the groundwater head and flux information available from the MODFLOW model.WetSpa simulates the groundwater recharge based on the soil moisture content of the soil layer independent from the groundwater depth.Although WetSpa allows direct AET from the groundwater system if the soil moisture is below the field capacity, the groundwater recharge output of WetSpa does not integrate this groundwater discharge towards the soil surface.The water flux from the groundwater system towards the soil surface is simulated as drain in the MODFLOW model.MODFLOW does not differentiate between drainage that is lost by AET or that contributes to the river baseflow.Fully integrated hydrological models are better suited to simulate these interactions between the saturated zone, unsaturated zone, land surface and vegetation in a physically based way (Hendricks Franssen, 2009;Holman et al., 2011).
One of the consequences of the applied one-way coupling on the results is that AET could be underestimated in zones with shallow groundwater depths (Maxwell and Kollet, 2008).This is most likely partly compensated during calibration of the drain conductance in the reference MOD-FLOW model, however winter-summer differences remain.Additionally, when due to future climatic changes grid cells change from recharging the groundwater system to receiving groundwater discharge or vice versa the drain flux might be under-or overestimated.Since the decrease in groundwater level obtained in this study, especially in the valleys where wetlands are concentrated, is relatively small, the impact of the climate change on the ET from the groundwater system is most likely also small.The small change in groundwater level in the valleys justifies the one way coupling applied in this paper.However, in studies predicting a larger change of the groundwater level in areas where surface water processes are influenced by the groundwater level, a fully integrated surface water-groundwater simulation of the hydrological processes is required.

Conclusions
This paper discusses how climate changes alter the spatiotemporal dynamics of the groundwater system.Until recently hydrological impact assessment of climate change has been focused primarily on peak flows and flood events (Solomon et al., 2007;Maxwell and Kollet, 2008).However, most GCMs predict that global warming is likely to amplify drought events over Europe.Consequently, there is a growing concern on the future availability of water for drinking water supply, crop growth and natural vegetations throughout the year.Hence, there is an urgent need for more research on the impact of those drought events on low flows and on the groundwater system.
Our paper is one of the first that analyzes the impact of climate change on the groundwater system with a high spatiotemporal resolution at the watershed scale.Applying this high spatial and temporal resolution showed that the impact is highly variable both in space and time.We found that for our study area, situated in Western Europe, the ensemble average of 28 climate change scenarios predict a decrease in summer groundwater recharge causing reduced groundwater heads and lower groundwater discharge fluxes especially in late summer-early autumn.Because of the increasing precipitation during winter the groundwater head and flux during spring are expected to decrease only slightly.Groundwater level changes are shown to be more pronounced on the interfluves and upstream in the catchment.The MHGL, MLGL and groundwater discharge frequency are likely to decrease at most places.The results also indicate the importance of applying transient climate change impact assessments as seasonal variations of the changes are significant.
Additionally, our research shows the importance of applying an ensemble of climate change predictions.By applying 28 different climate scenarios obtained from different GCMs and RCMs we indicate the uncertainties associated with the results.As the uncertainties of the climate scenarios are large the additional uncertainties from the hydrological and groundwater flow models are not taken into account.Due to the large uncertainties in the predictions of climate variables, especially precipitation, the predicted impact on the groundwater system obtained in this research should be considered as trends and order of magnitudes rather than exact predictions.
To reduce model calculation time and increase the model stability a loose coupling is applied between the surface water model Wetspa and the groundwater flow model MOD-FLOW.Further research should examine how models could be improved for assessing the impact of climate changes on the groundwater system, for example by including vegetation growth, physically based ET calculation, hourly time discretization, further coupling of surface-subsurface processes without increasing the data requirements and computation time too excessively.
Although it is advisable to mitigate climate change as much as possible it has become clear over the past decade that we will also have to adapt to climate change.To prevent the loss of groundwater dependent vegetation and reduced crop growth due to drought problems, resource managers should consider adaptive measures as soon as possible.An important message from the results is that GWDTEs are especially vulnerable due to too low summer groundwater levels and reductions in the magnitude and frequency of groundwater discharge to the landscape.
Because climate models predictions are highly variable spatially (Solomon et al., 2007;Hendricks Franssen, 2009) similar research should be done for different hydro-climatological and hydrogeological type locations to gain insight into the meteorological and basin characteristics controlling the impacts of climate change on groundwater systems.

Fig. 1 .
Fig. 1.Location and topography of the study area including the geographical position of the observation and most important pumping wells and rain gauges.

Fig. 4 .
Fig. 4. Conceptual overview of the applied spatial-temporal methodology.The figure shows a watershed discretized using a rectilinear grid, surface water bodies are represented in blue.For every cell all water balance components are simulated daily and the runoff, interflow and groundwater flow are routed to the outlet of the catchment.Recharge and discharge are aggregated to half-monthly time steps.Two cells in this figure are highlighted.Cell A represents a typical groundwater discharge area: during most time steps the groundwater in this cell flows from the groundwater system towards the land surface where the groundwater can discharge to the surface water bodies or be used for evapotranspiration.Cell B represents a typical recharge area where the water table is recharged by water infiltrating from the land surface.The graphs on the right show how the groundwater discharge or recharge flux typically evolve over time.In this study the groundwater system is simulated for the reference condition (grey line) and several climate change scenarios (e.g.orange line).

Fig. 5 .
Fig. 5. Comparison of the filtered baseflow, the baseflow simulated by WetSpa and the baseflow simulated by MODFLOW.

Fig. 6 .
Fig. 6.Comparison of the measured groundwater heads and heads simulated by the MODFLOW model.

Fig. 7 .
Fig. 7. Average intra-annual variability of (a) precipitation, (b) potential ET, (c) actual ET, (d) runoff, (e) soil moisture storage and (f) groundwater recharge for the reference climate (1960-1991), 28 climate scenarios (2070-2101) and the average of the climate scenarios.One year is divided into 24 half monthly time steps, for every time step the average of 32 yr simulation is presented.Error bars represent one standard deviation between the climate scenarios.

Fig. 8 .
Fig. 8. Average intra-annual variability of (a) groundwater head and (b) groundwater discharge for the reference climate (1960-1991), 28 climate scenarios (2070-2101) and the average of the climate scenarios.One year is divided into 24 half monthly time steps, for every time step the average of 32 yr simulation is presented.Error bars represent one standard deviation between the climate scenarios.

Fig. 9 .
Fig. 9. Spatial distribution of the simulated change (average future minus reference scenario) in temporally averaged groundwater recharge.Negative changes indicate a decrease in groundwater recharge from the reference status to the average future state.

Fig. 10 .
Fig. 10.Spatial distribution of the simulated change (average future minus reference) in temporally averaged: (a) groundwater level; (b) mean highest groundwater level (MHGL); (c) mean lowest groundwater level (MLGL); and (d) mean spring groundwater level (MSGL).Positive changes indicate an increase in groundwater level, negative changes indicate a decrease in groundwater level from the reference status to the average future state.Rivers are shown in white.

Fig. 11 .
Fig. 11.Scatter-plot showing for all MODFLOW grid cells, which have an average groundwater discharge flux (drain and river leakage) between 0 and 20 mm day −1 (x-axis), the difference in groundwater discharge flux between the reference period (1960-1991) and the average of the 28 climate change scenarios (2071-2101) (yaxis).The moving average is calculated over a range of 500 successive values.The lines indicate the 50 % and 15 % decrease in groundwater discharge flux.

Fig. 12 .
Fig. 12. Scatter-plot of reference groundwater discharge frequency (x-axis) versus the change in groundwater discharge frequency from the reference period (1960-1991) to the average future period 2070-2101 (y-axis).The groundwater discharge frequency is the percentage of time that a certain model cell has a positive groundwater discharge flux, the quantity of this flux is not taken into account.Every point in the graph represents a model cell in the watershed where at least during one time step of the reference period groundwater discharge occurs.
The subsurface of the model area is limited to the Quaternary and Tertiary sediments which are deposited on the Boom clay aquitard during the Oligocene epoch.From depositionally oldest to youngest the Hydrol.Earth Syst.Sci., 16, 1517-1531, 2012 www.hydrol-earth-syst-sci.net/16/1517/2012/

Table 1 .
Overview of the hydrostratigraphy of the study area