Hydrogeological evaluation of managed aquifer recharge in a glacial moraine complex using long-term groundwater data analysis

Managed aquifer recharge (MAR) is a strategy to address dwindling water availability that is used to recharge stressed groundwater systems for recovery or adaptation purposes. Glacial moraine complexes can host large groundwater volumes in thick coarse-grained sandy and gravelous aquifers, and therefore are often suitable for MAR. In this research, the impact of 20 years of MAR (infiltration) in the Veluwe glacial moraine complex in the Netherlands is evaluated through time-series analysis and water quality modeling. The research enhances the understanding of hydrological processes in a glacial moraine complex and thus supports effective MAR design. The results show that MAR in the Veluwe area has raised nearby groundwater heads, and that 20 years of infiltration did not significantly deteriorate the groundwater quality. In addition, the analysis revealed an unexpected hydraulic interaction between the glacial complex and the underlying aquifer. The results suggest that MAR in a glacial moraine complex can be an effective strategy for storage of surplus surface water, compensation of groundwater abstraction, or water quality improvement. Monitoring and evaluation of groundwater heads and quality will improve the understanding of the hydrology and hydrochemistry of the water system and the used water resource, which is essential to design effective MAR systems.


Introduction
Managed aquifer recharge (MAR) is used to recharge stressed groundwater systems (Dillon et al. 2018;Geelen et al. 2017), and can be used as a strategy to reverse a potentially dwindling water availability. Water obtained from aquifers that are recharged via MAR can be used for various purposes such as drinking water supply, industrial water supply and irrigation. Multiple examples of MAR systems have been developed to reduce the overexploitation of aquifers and restore natural hydrological conditions (Dillon et al. 2018). The potential of MAR is widely acknowledged by water resources specialists and is estimated to contribute to 10% of the global groundwater abstraction (Dillon et al. 2018).
Glacial moraine complexes often contain thick, coarse grained, sandy and gravelous aquifers, and are suitable for MAR because of their high hydraulic conductivity (Brun and Jensen 2001). However, these complexes are also known for their variation in hydraulic conductivity (Bense et al. 2013), which could limit the recovery efficiency of a MAR system (NRMMC et al. 2009). Good knowledge of sitespecific aspects and aquifer properties determine the MAR design (Ringleb et al. 2016). Hydrological models can be used to predict groundwater flow at MAR sites. In areas with geological heterogeneity such as glacial moraine complexes, Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10040-020-02145-7) contains supplementary material, which is available to authorized users. detailed hydrogeological mapping is essential to support the hydrological modeling (Brun and Jensen 2001).
Potential water sources for MAR are rainwater, surface water, groundwater from a different aquifer, reclaimed water, or runoff (Page et al. 2018). The quality of the water source and soil processes such as dispersion, cation exchange, or degradation of contaminants, determine the impact of the MAR to the groundwater quality (Laws et al. 2011). To ensure a hydrologically sustainable MAR design, it is therefore important to have a detailed understanding of the groundwater flow (Maliva et al. 2015) and of the contaminants' transport and fate in the groundwater and soil system (Gurjar et al. 2018).
The Veluwe area in the Netherlands (Fig. 1) is a glacial moraine complex with a large unconfined, coarse grained, sandy and gravelous aquifer, containing lateral barriers of glaciotectonic origin (Bense and Kurylyk 2017;van Engelenburg et al. 2017). For the last two decades, smallscale MAR has been applied in the north east of the Veluwe area. The research question is: What is the impact of 20 years of managed aquifer recharge (MAR) as compensation for the hydrological impact of the local drinking-water abstraction at Epe in the Veluwe glacial moraine complex in the Netherlands? To answer this research question, the available monitoring data over 20 years of infiltration at Epe have been analyzed through time-series analysis and one-dimensional (1D) modeling of groundwater quality. The results will enhance the understanding of hydrological processes in a glacial moraine complex to support effective MAR design.

Veluwe, The Netherlands
The Veluwe area is part of a glacial moraine complex originating from the Saalien glacial period in the Pleistocene (approximately 100,000 years BP; Rutten 1960;Gehrels 1999). The depth of the (freshwater) aquifer reaches up to 200 m, and the unsaturated zone has a maximum depth of approximately 70 m from ground surface. There is an absence of surface water in the elevated part of the system and the Veluwe area forms a large infiltration zone, where groundwater is recharged by infiltration of rainwater (Gehrels 1999;Kumar et al. 2016). Natural groundwater-level elevation across the area is characterized by fluctuation on the frequency time-scale of decades and with an amplitude of several meters, whereas around the edges of the elevated ridge the groundwater level is much shallower, and groundwater levels show annual fluctuations with an amplitude in the order of 1 m or less (van Engelenburg et al. 2017). Along these edges with shallow groundwater there are seepage zones and (manmade) brooks and springs that are vulnerable to groundwater level decline. The brooks and springs are mainly discharging groundwater from the Veluwe aquifer, supplemented with local rainwater.
The land use of the central part of the Veluwe has historically been limited to woodland, heath and sand drifts (Witte et al. 2019). The groundwater quality in the aquifer therefore is hardly affected by anthropogenic contamination. Because of the large availability and natural quality of groundwater in the Veluwe area, it has been used as a drinking water resource for centuries. Currently, Veluwe groundwater is abstracted for drinking water at 21 sites (Fig. 1).
The groundwater flow patterns in the Veluwe area can be erratic, due to the presence of glaciotectonic thrust zones, along which clay layers have been dragged into a subvertical position. These imbricated layers were formed when a glacier pushed up the sediment in the area during the Saalien glacial period (Bakker and Van Der Meer 2003;Gehrels 1999;Verhagen et al. 2014). The location of these tilted clay layers is barely visible in the landscape, but the impact of some layers on the groundwater flow is noticed through significant jumps in hydraulic heads (Van Engelenburg et al. 2012), similar to those found along tectonic fault zones (Bense et al. 2013). Based on these hydraulic jumps, the location of some apparent tilted clay layers is indicated (Fig. 1b), but the depth, length and slope of the tilted clay layers is still uncertain. Due to the origin of the layers, they are expected to descend from west to east, but the actual form has not been determined yet.

MAR in the Veluwe area
The drinking-water abstraction at Epe, near the eastern edge of the Veluwe area ( Fig. 1), has been abstracting a maximum of 6 million m 3 of groundwater per year since 1958 (Fig. 2) from a depth of 30-80 m below surface level. The raw water quality is good with low water hardness (mmol/L Ca 2+ +Mg 2+ ). The abstraction indirectly affects the wetland area Wisselse Veen and the Tongerense Beek brook system (Fig. 1b).
Due to the large unconfined coarse-sanded aquifer with a thick unsaturated zone, the Veluwe area is suitable for MAR, enabling spontaneous infiltration of water applied to the surface. Since 1998 untreated surface water has been infiltrated in ponds as part of a local-scale MAR project (Fig. 1). The total pond size is approximately 45,000 m 2 , and the maximum filling depth in the ponds is approximately 3 m. The intake of surface water takes place from a local brook. This brook mainly discharges groundwater from the Veluwe area towards the nearby River IJssel. By the infiltration, a part of this discharging Veluwe groundwater is infiltrated back into the Veluwe system instead of being discharged to the river. In the first stage of the infiltration scheme (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014), maximal 2 million m 3 /year surplus water was infiltrated (Fig. 2). In the second stage, since 2015, the infiltration system was extended to reach a goal of infiltration capacity 6 million m 3 /year, to fully mitigate the local drinking water abstraction at Epe. Due to limitations at the intake, this goal has not been met yet (Fig.  2). Through adjustments in the water management and technical optimization of the intake station in 2019 and 2020, it is expected that the goal of zero net abstraction will be met from 2021 onward.
The aim of the infiltration scheme is to restore groundwater levels that are influenced by the Epe drinking-water abstraction, and to restore ecological values and seepage flows towards Wisselse Veen (Fig. 1b). Because the land use and surface-water management were also adapted to restore ecology and seepage in Wisselse Veen, the actual impact from the infiltration on the ecology and seepage cannot be identified. A precondition for the infiltration is that the groundwater quality may not deteriorate; therefore, the impact of the infiltration is measured by analysis of groundwater heads and quality.

Methodology
In this research the impact of infiltration to the groundwater system is evaluated, to gain knowledge on the hydrogeology of the Veluwe glacial moraine area in the Netherlands. This will support the design of the MAR system for groundwater replenishment in glacial moraine complexes. To achieve the aims of the research, available groundwater head data and water quality data were analyzed, using time series analysis and 1D geochemical modeling. Data were collected on groundwater heads, groundwater quality, precipitation, evaporation, groundwater abstraction and surface-water infiltration (Table 1). To measure groundwater heads on different depths, two or more monitoring wells with various lengths with a screen at the end of the well are situated within one borehole. The measured groundwater heads in each monitoring well represent the groundwater head at the depth of the screen. For water quality analysis, groundwater is abstracted from the monitoring well for a period of time to refresh the water in the well before a sample is taken to ensure the water origins from the screen depth. Because the time series are long, the data collection methodology has developed and changed over time-for instance, the groundwater data used to be collected by hand measurements once every two weeks, and this has been replaced by automated diver measurements. The analytical methods used for water quality parameters have accordingly developed over the years. The groundwater data collection was performed and validated by Vitens (this study), and groundwater quality analyses were performed and validated Fig. 1 Overview: a Veluwe area, elevation, drinking water abstractions of Veluwe groundwater, meteorological stations, research area (purple rectangle) and cross section A-A′ (red line); b Research area (purple rectangle, see part a) with used monitoring wells, cross section A-A′ (red line) and well field/infiltration ponds (yellow rectangle); c Well field with drinking-water abstraction wells and infiltration ponds Epe (yellow rectangle, b) Fig. 2 Gross groundwater abstraction (blue), surface-water infiltration (orange) and net abstraction (gross groundwater abstraction minus surface-water infiltration volume) (green). Actual volumes at the drinking water abstraction and infiltration Epe, the Netherlands from 1958-2019 (data from Vitens, The Netherlands) Parameters analysed (this study) Step 1 Step 2 Step 3 Groundwater  Table 2 Water quality data: Average concentrations of macro parameters just before the infiltration started (1997)(1998)(1999) and recently (2014)(2015)(2016) at intake, in monitoring well 27B-356, and in three groups of abstraction wells: highly influenced (B5), mildly influenced (A9, B7, B9) and not influenced (A6, A7, A8, B8). For the location of the intake area and wells see Fig. 1c  TDS total dissolved solids (calculated); EC electrical conductivity; ++/--strong effect of infiltration; +/-moderate effect of infiltration; 0 no effect of infiltration; a = impact of acidification by the certified laboratory of Vitens, the drinking water company that operates the Epe site, or external laboratories of equivalent certification.

Groundwater heads
Time series analysis is a statistical modeling method, using time series of measured groundwater heads and different explanatory series to calibrate groundwater-head simulation. In time series analysis, the explanatory series are considered stress models, using a response-function to transform the stress into the contribution to the simulated groundwater heads. The parameters of each explanatory series are optimized, minimizing the residuals, which is the difference between the model result and the observations. The analysis was executed using Menyanthes software (von Asmuth et al. 2012). In this research, a time series model is considered valid when the explained variance percentage (EVP) exceeds 70 %, the evaporation factor (EVAP) is between 0.5 and 1.5, and the root mean squared error (RMSE) is smaller than 0.15 (Von Asmuth et al. 2011).
The response of the groundwater in the Veluwe system with the thick unsaturated zone is slow, and therefore changes such as abstraction and infiltration result in a long transition period from one state to another. Because this transition period must be omitted from the time series analysis (Von Asmuth et al. 2008), the length of the time series is important for the result of the analysis. Therefore a selection was made from the available groundwater head time series in the area, selecting all monitoring wells that were measured for at least 10 years during the infiltration period. The wetland area Wisselse Veen is a fast responding system on the edge of the Veluwe (van Engelenburg et al. 2017) and monitoring in this area only started in 2010. Because the aim of the infiltration is to restore seepage flows in this area, and because it is a fast responding system with a short transition period, the researchers chose to include these wells in the research despite the short data period. This led to the selection of 101 time series for 86 monitoring wells (Fig. 1b).
The precipitation data time series from the Royal Dutch Meteorological Institute (KNMI) meteorological station at Epe was used. Although there are longer time series of precipitation available, this local station was used because the Veluwe precipitation volume is affected by topography and land use (ter Maat et al. 2013). For evaporation, the Makkink reference crop evaporation data are used that were calculated by KNMI (Hiemstra and Sluiter 2011) using the data from the meteorological station at Deelen. Both stations (Epe and Deelen) are located in the Veluwe area and near to the Epe abstraction/ infiltration location (Fig. 1).
The Menyanthes analysis was executed in three steps with different explanatory series: first only precipitation and evaporation, second with additionally net abstraction, and third with gross abstraction and infiltration instead of net abstraction. When the addition of an explanatory series results in a valid model with an EVP improvement of more than 5% compared to the previous model steps, the impact of the addition is considered a significant improvement.
Menyanthes calculates the response of the groundwater head (in meters) to an abstraction of 1 m 3 /day (parameter: M0). As the indicator to compare the impact of net abstraction, (gross) groundwater abstraction and surface-water infiltration of 1 million m 3 /year (approximately 2,800 m 3 /day), the stationary impact of 1 million m 3 /year was calculated by multiplying the value of M0 with 2,800. A summary of the valid models is given in Appendix 1. The electronic supplementary material (ESM1) provides the Menyanthes results for all 101 time series.

Water quality
At the Epe infiltration site, untreated surface water from local brooks is infiltrated. The aim of the water quality analysis is to understand how the infiltration water influences the groundwater quality, by analyzing trends in macro parameters. First, a comparison was made of the trends in the infiltration water, the individual abstraction wells of Epe, and a monitoring well positioned between the infiltration ponds and the abstraction wells. Second, the process of cation exchange during infiltration was modelled. PHREEQC is a hydrogeochemical modeling code used to describe the interaction of several hydrogeochemical processes during transport of groundwater through the aquifer (Parkhurst and Appelo 1999). The aim of the PHREEQC modeling is to simulate the trends of the different anions and cations in the abstracted groundwater after startup of infiltration, to understand which hydrogeochemical processes take place. The quality of the intake water was used as input. The value for dispersion was calibrated on the trend for sulphate. For calcium and other cations, a 1D PHREEQC transport model was set up to calibrate cation exchange on the breakthrough curves of potassium, calcium and magnesium.

Results for groundwater heads
In the first step of the Menyanthes modeling, all 101 time series were modeled using precipitation (P) and evaporation (E) as explanatory series, resulting in 47 valid models (Appendix 1). The results of this first step confirm that in areas with a fast responding water system, such as the Wisselse Veen area at the edge of the Veluwe, groundwater-head response is explained well by precipitation and evaporation time series (Fig. 3).
In the second step, the net abstraction (Qn) was added as an explanatory series to all 101 time series. For 11 time series this resulted in significantly improved model results (Fig. 4). The results indicate that stress caused by net abstraction is decreasing with the distance between monitoring well and abstraction/ infiltration site, but is also influenced by the tilted clay layers: six of the monitoring wells with significantly improved model results are situated just north of the area that has surface water, within the same compartment between tilted clay layers as the abstraction/infiltration. The impact of the net abstraction reaches south as far as monitoring well B27D0517, situated at the north side of the fast-responding Wisselse Veen area (Fig.  4). The monitoring wells B27D0143 and B27D0045 are close to the abstraction/infiltration, but apparently east of the tilted clay layer: the impact from the net abstraction is comparable to the impact in B27D0125, which is further away but within the same compartment of the abstraction/infiltration (Fig. 4). This indicates that there is a limited groundwater flow through the tilted clay layers, causing a preferential flow direction parallel with the clay layers.
The results of step 2 also show that the addition of net abstraction as an explanatory factor does not improve the results for the monitoring wells in the fast-responding area (Figs. 3 and 4). The change in groundwater flow resulting from the reduction of the net abstraction volume may cause groundwater to be discharged through the Tongerense Beek surface-water system instead of increasing the groundwater levels, which would also prevent seepage flow restoration in the Wisselse Veen area.
In step 3, net abstraction was replaced by (gross) groundwater abstraction and surface-water infiltration separately. Eight time series showed a significant improvement compared to the models with net abstraction (Fig. 5). These eight time series included data for groundwater heads measured at four deep monitoring screens (sometimes referred to as filters) that are situated in the underlying Maassluis aquifer (Fig. 6). These results were unexpected and called for further analysis of this result.
While the Menyanthes analysis did not result in valid models for the shallow screens of these four deep monitoring wells, the groundwater heads time series were analyzed visually (Fig. 6). The differences in average groundwater heads between the monitoring wells confirm the discontinuity in the groundwater flow in the west-east cross section caused by the tilted clay layers. Although surface levels decline in the eastward direction, the groundwater heads indicate that both deep and shallow groundwater flow also is partially directed westward. The groundwater head in the shallow screen B27B0155-1 has the lowest average groundwater head of the shallow screens. This may be explained by the location of the monitoring well west of a tilted clay layer system, but also indicates that the groundwater flow is directed westward in this area. B27B0255-1 is situated in the direct proximity of the abstraction and infiltration and therefore shows strong annual groundwater head fluctuations, following the net abstracted volume. The fluctuations due to the abstraction and infiltration volume decrease with the distance to the well field, but additionally there are significant groundwater head jumps between B27B0155-1, B27B0238-1 and B27B0156-1, probably caused by the tilted clay layers (Fig. 6).
The deep screens are situated in the Maassluis aquifer, a marine formation below the fluviatile Peize Formation. The average groundwater heads in the deep screens exceed the elevation of the top of the aquifer, which indicates that the Maassluis aquifer is a confined aquifer where the groundwater is under pressure. The groundwater heads in the deep screens differ from the groundwater heads in the screens situated in the overlying Peize aquifer and the pushed moraine complex. and are within a much smaller range compared to the groundwater heads in the shallow screens (Fig. 6). This indicates a limited permeability of the interface between the push moraine complex/Peize Formation aquifer, which is confirmed by the presence of a clay layer at an elevation of -123 to -128 m asl (see Appendix 2). However, the Menyanthes model results show that the abstraction and infiltration volume time series partially explain the groundwater fluctuations in the deep aquifer, indicating that there must be some groundwater interaction between this aquifer and the Peize aquifer. The groundwater fluctuation patterns of the four deep screens also reflect the same annual pattern as the shallow screen B27B0156-1 (Fig. 6). This could indicate that there is a hydraulic connection between the compartment of the push moraine complex in between the adjacent tilted clay layers where B27B0156 is situated, and the Peize and Maassluis aquifers.
The infiltration takes place at the surface level in the cone of depression caused by the groundwater abstraction. Due to the limited hydraulic conductivity through the tilted clay layers, the largest impact of the infiltration is found in the compartment between the adjacent clay layers. The groundwater is abstracted at a depth of 30-80 m below the surface, and originates partially from the deep aquifers below the push moraine complex. The abstraction from the aquifer affects the groundwater flow directly; the major part of the infiltration fills the cone of depression and is abstracted within 1-3 years after infiltration. The remaining infiltrated water will only reach the saturated zone after slowly passing through the thick unsaturated zone, causing a slower and smaller response to the infiltration than to the abstraction.

Discussion on groundwater heads
The results of the analysis of the groundwater heads confirm the indicative location of the identified tilted clay layers near the infiltration site. The difference in modeling results between analysis with net abstraction (Qn) and with gross The results confirmed that the elevated part of the Veluwe is a slow responding groundwater system due to the presence of a thick unsaturated zone, where precipitation and evaporation do not translate directly into the groundwater head (van Engelenburg et al. 2017). In the slow-responding elevated part of the Veluwe, the simulated models with precipitation and evaporation as explanatory data are not valid, in contrast to the results in the fast-responding system in the Tongerense Beek and the Wisselse Veen area.
The horizontal groundwater exchange in the push moraine complex is limited due to the presence of tilted clay layers, and there is a hydraulic interaction between the push moraine complex and the Peize aquifer, which seems to function as a combined aquifer. The results indicate that although there is a clay barrier between the Peize aquifer and the Maassluis aquifer, there is groundwater interaction between both aquifers. Furthermore, the impact of the abstraction and infiltration to the groundwater levels and thus to the seepage flows in the Wisselse Veen wetland area is limited. This can be explained by the assumption that groundwater flow that is restored by the infiltration discharges through the surface-water system of the Tongerense Beek before reaching the Wisselse Veen area. These local hydrological characteristics limit the effectiveness of the infiltration towards the Wisselse Veen area.
To further study this interaction, the gained knowledge on tilted clay layers and the interaction between the push moraine complex and the Peize and Maassluis aquifers should be added to the regional hydrological model of the area. An optimization study could be performed by modeling the impact of the tilted clay layers under different slopes, and the interaction between the deep aquifers and the push moraine complex. This will enhance the understanding of the hydrogeology further, which will support an adequate design of additional MAR in the Veluwe area.
Time series analysis of monitoring data provides historical information on groundwater head trends, which provides knowledge on groundwater dynamics in a hydrogeologically heterogeneous groundwater system, that can help to build or improve a groundwater model (van de Vliet and Tiebosch 2001). The time series analysis did not predict all groundwater-head data series well. Only large  Fig. 1a), based on REGIS II v2.2 (TNO 2018) combined with groundwater monitoring data and approximate locations and directions of tilted clay layers; b Groundwater-head time series in shallow-screen and deep-screen monitoring wells in cross section changes in the groundwater dynamics, such as abstractions and infiltration, can be explained through this statistical analysis and optimization. The impact of smaller, or less measurable, influences such as changes in the surfacewater systems that affect groundwater recharge or discharge, cannot be predicted with this method. However, because time series models are easy to construct (von Asmuth et al. 2012), they can be used as a quick-scan of an area to gain hydrological knowledge before designing a MAR system, and thus can also contribute to the development of a spatial hydrological model for the area.

Results for water quality analyses
The groundwater quality of the Veluwe area is gradually changing due to acid deposition , which is shown by the increase of sulphate in the abstraction wells ( Table 2). The monitoring well 27B-0356 and the abstraction well B5 are the wells nearest to the infiltration ponds (Fig. 1c). On average, the water quality in both wells after 15 years of infiltration matches the water quality of the infiltration water. This indicates that the original groundwater in these wells was fully replaced by infiltrated surface water in this period of time ( Table 2). The water quality in other abstraction wells is also affected, but to a lesser extent. This is shown by the differences in the average concentrations of macro parameters such as pH, HCO 3 -, and the cations Al 3+ , K + , Ca 2+ , and Mg 2+ , just before the infiltration started (1997)(1998)(1999), compared to more recent concentrations (2014-2016; Table 2). Abstraction well B5 and well group A9/B7/B9 are situated near the infiltration. Here the effects of infiltration (dominant) and acidification (minor) are combined, by mixing of native groundwater and infiltration water. The average quality of the infiltration water almost equals the groundwater quality in the monitoring well 27B-0356.2, and the quality of abstracted groundwater in well B5. In the monitoring well only pH shows a slight decline from 7.4 to 7.0, probably due to precipitation of iron during transport and infiltration. Well B5 shows the strongest effect of the infiltration (Table 2), whereas pH rises to an almost neutral value of 6.8, and aluminium concentration declines to a very low level of 6 μg/L, and the concentrations of calcium, bicarbonate and potassium nearly double. The increasing calcium concentration in well B5 from 10 to 24 mg/L is the most notable change in water quality (Fig. 7). Well group A9/B7/B9 shows a less pronounced effect. This can be due either to delay by longer travel time, or to a different mixing ratio between infiltration water and native groundwater. In well group A6/A7/A8/B8 the water quality changes only slightly, mainly due to acidification. The decreasing pH in these wells from 6.3 to 6.0 and the increasing SO 4 2concentration show the moderate impact of acidification. Acidification causes an increase of dissolution of lime, here resulting in a small increase of the concentrations of HCO 3 -, Ca 2+ and Mg 2+ .
The trend of calcium (Fig. 7) illustrates the length of the transition period to a new equilibrium. In well B5 it has taken about 3 years to reach the calcium concentration of the infiltration water of approximately 25 mg/L (Fig. 7). In well group A9/B7/B9 a more gradual transition occurs, where it is unclear if an equilibrium already has been reached in 2017. In well group A6/A76/A8/B8 the calcium concentration is stable around 10 mg/L, confirming the assumption that these wells seem uninfluenced by infiltration water.
The modelled trend for sulphate indicates the impact of dispersion, because sulphate can be considered conservative in this aerobic environment. The trend (breakthrough curves) of potassium, calcium and magnesium indicates the impact of cation exchange. The 1D PHREEQC model simulates the quality changes over the period 1993-2018 for cations and sulphate adequately, and shows that the breakthrough of surface water can be described by the main hydrochemical processes dispersion and cation exchange (Fig. 8).

Discussion on water quality
The first step in geochemical modeling is to start with basic processes like transport, dispersion and ion exchange. If the simulation is not satisfactory, more complex processes can be added to the model, like (kinetic) solution of minerals or redox processes. In the studied case with an oxic, sandy aquifer containing only low contents of reactive minerals, it can be assumed that dispersion and ion exchange are the dominant processes (Appelo 1982). This assumption is confirmed by the presented model results, which show a good simulation by only modeling transport, dispersion and ion exchange.
Based on the modeling results with PHREEQC, the travel time is estimated at approximately 2 years from infiltration pond to abstraction well B5. Therefore, it is expected that the effect of the extension of the infiltration volume that started in 2015 will protrude in the abstraction wells gradually from 2017 onwards. The years 2014-2016 are used in the comparison, because they represent the period just before the recent extension of the infiltration volume in 2015.
Managed aquifer recharge is frequently used as a treatment step for reclaimed water, because the passage through the unsaturated zone potentially improves the quality of the infiltrated water (Bekele et al. 2011). In this research, the PHREEQC modeling results in an average cation exchange capacity for the aquifer. The water quality analysis shows that the influence of the infiltration to the quality of the abstracted groundwater is relatively small, and the rise of hardness is the most remarkable change in the quality of the abstracted groundwater.
The impact of the infiltration on the water quality in the adjacent area is uncertain, because the available groundwater quality data were limited to the abstraction wells and one nearby monitoring well. To avoid the risk of spreading of infiltration water to the surrounding groundwater system, the net yearly abstraction will be kept positive, by controlling the groundwater abstraction volume as well as the infiltration volume. Although in winter periods the infiltration rate can exceed the abstraction rate, the winter surplus will remain within the cone of depression caused by the groundwater abstraction, due to the slow response of the Veluwe groundwater system. To prove this assumption, more data from monitoring wells combined with nonstationary modeling of the surrounding groundwater system is required.

Conclusions and recommendations
In this research, monitoring data associated with over 20 years of infiltration at Epe, the Netherlands, was analyzed with time series analysis and 1D modeling of groundwater quality. The aim of the research was to evaluate the impact of 20 years of MAR as compensation for the hydrological impact of the Epe drinking-water abstraction in the Veluwe glacial moraine complex, to enhance the understanding of hydrological processes in a glacial moraine complex to support effective MAR design.
The results of the research suggest that MAR in a glacial moraine complex can be an effective strategy for storage of surplus surface water, compensation of groundwater abstraction, or water quality improvement. The results also show that the complexity of the hydrogeology in a glacial moraine complex may affect the possibility to reach a specific goal, and may cause unintended side effects such as problems with water quality changes or increased groundwater levels in other areas. A precondition for an effective MAR system is an adequate intake as well as infiltration facility. To design effective MAR systems, a thorough knowledge of the hydrology and hydrochemistry of the groundwater, as well as the surface-water system that is used as the water resource for infiltration, is essential.
In the Epe infiltration case, long time series of monitoring data were available. The time series analysis confirmed that the glacial moraine complex of the Veluwe area, with a thick coarse-grained, sandy and gravelous aquifer, is suitable for MAR, but also showed that its hydrogeological heterogeneities may reduce the effectiveness of MAR. The infiltration has provided a positive impact to nearby groundwater levels; however, from the research it can be concluded that the infiltration does not contribute significantly to the restoration of the water system in the Wisselse Veen area, because of the groundwater discharge in the deep aquifers below the glacial moraine complex and the surface-water discharge in the Tongerense Beek system. The water quality data and modeling showed that due to the good quality of the infiltrated water, 20 years of infiltration mainly resulted in a change in acidity and hardness of the abstracted groundwater. The 1D geochemical modeling indicated that cation exchange and dispersion are the main geochemical processes in the aquifer between infiltration and abstraction. From the research it is concluded that time series analysis and 1D geochemical modeling may contribute valuable knowledge of the water system, provided that there are sufficient monitoring data available. A well-functioning monitoring system of groundwater heads and quality is essential, preferably starting before the MAR system is designed.
Finally, from this research it is concluded that the Veluwe glacial moraine complex in the Netherlands is fit for MAR by infiltration. To predict the impact of MAR on the local and regional groundwater system of the Veluwe area, a further understanding of the hydrogeology of the area is necessary. The findings of this research can be elaborated to gain additional knowledge on the complex hydrogeology of the push moraine with tilted clay layers.
Data availability statement The summarized results of the Menyanthes time series analysis, the used groundwater quality data and the PHREEQC modeling input file are available in the electronic supplementary material (ESM1, ESM2, ESM3). The complete Menyanthes results are available on request from the corresponding author (JvE). The hydrogeological and groundwater data are publicly available at Dinoloket.nl.

Compliance with ethical standards
Conflict of interests There is no conflict of interests. Appendix 1: summary results of Menyanthes time series analysis A summary is given here of the valid models. A model is considered valid with an explained variance percentage (EVP) > 70%, an evaporation factor between 0.5 and 1.5, and the root mean squared error (RMSE) smaller than 0.15. An overview of the Menyanthes results for all 101 time series is available in ESM1.  Table 4 Menyanthes time series analysis X --47 - Table 5 X X -57 11 Table 6 X -X 6 2 8 Table 7 Analysis deep monitoring wells hydro(geo)logical data -----    Appendix 3: statistics water quality data In the following, a summary of the statistical data for Table 2 can be found.  Table 9 Summary of water quality statistics. The full source data are available in ESM2. n number of samples; min minimum concentration; max maximum concentration; average average concentration Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.