Numerical simulation of a managed aquifer recharge system designed to supply drinking water to the city of Amsterdam, The Netherlands

Managed aquifer recharge (MAR) is increasingly used to secure drinking water supply worldwide. The city of Amsterdam (The Netherlands) depends largely on the MAR in coastal dunes for water supply. A new MAR scheme is proposed for the production of 10 × 106 m3/year, as required in the next decade. The designed MAR system consists of 10 infiltration ponds in an artificially created sandbank, and 25 recovery wells placed beneath the ponds in a productive aquifer. Several criteria were met for the design, such as a minimum residence time of 60 days and maximum drawdown of 5 cm. Steady-state and transient flow models were calibrated. The flow model computed the infiltration capacity of the ponds and drawdowns caused by the MAR. A hypothetical tracer transport model was used to compute the travel times from the ponds to the wells and recovery efficiency of the wells. The results demonstrated that 98% of the infiltrated water was captured by the recovery wells which accounted for 65.3% of the total abstraction. Other sources include recharge from precipitation (6.7%), leakages from surface water (13.1%), and natural groundwater reserve (14.9%). Sensitivity analysis indicated that the pond conductance and hydraulic conductivity of the sand aquifer in between the ponds and wells are important for the infiltration capacity. The temperature simulation showed that the recovered water in the wells has a stable temperature of 9.8–12.5 °C which is beneficial for post-treatment processes. The numerical modelling approach is useful and helps to gain insights for implementation of the MAR.


Introduction
With the increase of groundwater demand, the managed aquifer recharge (MAR) technique is increasingly being applied worldwide to augment drinking-water sources, restore depleted storage, prevent saltwater intrusion, and improve water quality . A properly designed MAR system improves water quality through physical, chemical and biological processes, facilitating the removal of viruses, bacteria, organics and many more (Massmann et al. 2006).
A considerable number of review papers have been published, including those focused on the use of MAR for water recycling ), recharge and recovery using wells (Pyne 2005), technologies and engineering for artificial recharge of water (Bouwer 2002), riverbank filtration (Stuyfzand et al. 2006), MAR case risk assessment (Page et al. 2011), sustainable groundwater management through artificial recharge (Alqahtani et al. 2021;Jha et al. 2008), development of MAR in China (Wang et al. 2014), development of MAR across Europe (Sprenger et al. 2017) and history of artificial recharge in the coastal dunes of the Netherlands (Olsthoorn and Mosch 2020). Most of the MAR systems are established for economic benefit while also alleviating consequences from extensive groundwater extraction on the environment and ecology. Large-scale MAR systems are operated in many countries including Australia, India, The Netherlands, Germany, South Africa, and the USA. Some examples of operational MAR systems are the large-scale infiltration basins in Burdekin Delta in Queensland, Australia (Dillon 2009), the subsurface infiltration galleries of the aquifer recharge system of Geneva in Switzerland (De los Cobos 2002), the aquifer storage and recovery system in the Upper Floridan Aquifer System in the USA (Pyne 1995), and the Amsterdam Water Supply Dunes (AWD) system in the Netherlands (Olsthoorn and Mosch 2020). Amsterdam Water Supply Dunes is one of the sustainable MAR systems, consisting of a large number of small infiltration ponds operated for more than 60 years and can be considered a good example of MAR design in the Netherlands. The dune system is used for securing drinking water and maintaining groundwater level for natural vegetation and animals by recharging the pretreated water from the river Rhine (Geelen et al. 2016;Olsthoorn and Mosch 2020;Van Breukelen et al. 1998). The minimum residence time of 60 days and the average residence time of about 90 days in the aquifer brings improvement in the quality of water (Olsthoorn and Mosch 2020).
Numerical models have been widely used for MAR design and impact assessment. Ringleb et al. (2016) published a review of model tools for MAR investigation based on a survey of published studies over 30 years. The widely used models are MODFLOW, MT3DMS and SEAWAT, and most of the applications were for feasibility studies, system optimisation, computing residence time, and simulating temperature variations. Some examples of such studies are Jonoski et al. (1997) for optimising infiltration and pumping system, Chenini and Ben Mammou (2010) for feasibility assessment, Vandenbohede and Van Houtte (2012) for heat transport simulation, Hashemi et al. (2015) for estimating recharge rate of floodwater spreading, and Liu et al. (2019) for simulating the effect of temperature variation on infiltration rate and residence time.
Waternet, the drinking water supply company of the city of Amsterdam and surroundings, is the responsible organisation for the operation and maintenance of AWD. Waternet supplies about 75% of the total production from AWD (also called the western production site) and the remaining 25% from the drinking-water production plant in Loenderveen (also called the eastern production site) by treating water pumped from the Bethune polder. The prognosis, based on the internal analysis, suggests that the drinking water demand is expected to rise beyond the current supply capacity. There is great need for expansion of the water production capacity to meet the rise in demand by 10 × 10 6 m 3 /year. Alternative sources of water supply were investigated in the Netherlands to cope, and some of the research output includes Van Hateren (2019), Antoniou et al. (2017), Stuijfzand et al. (2013) and Faber and Grimm (1996). As climate variability has caused higher fluctuations in the available surface discharges and water supply companies have limits on groundwater pumping (overseen by the government to prevent greater drawdown of the water tables), use of managed aquifer recharge and recovery techniques has become an interesting option.
This study aims to propose a new MAR system, consisting of 10 infiltration ponds and 25 recovery wells, in the eastern production site for a production capacity of 10 × 10 6 m 3 /year. The identified raw water sources include surface water from Bethune polder and Amsterdam Rhine Canal (ARC). The source water will be pretreated, stored in the reservoir, and used for infiltration in the MAR ponds. The quality of water is expected to improve further by aquifer passage with sufficiently long residence time.
This study used an integrated modelling approach for the design and simulation of the proposed MAR system. The conceptual design of the MAR system includes: (1) a layout of the infiltration ponds system, including all the ponds and their codes/names, shapes, and dimensions; (2) the codes/ names, locations, and pumping rates of the recovery wells. The main criteria for the design include: (1) minimum of aquifer residence time of 60 days; (2) maximum drawdown of 5 cm in the residential areas. Groundwater flow and transport models were used to: (1) compute drawdowns caused by the MAR system; (2) compute infiltration rate of the ponds; (3) compute residence time of infiltrated water; (4) determine fractions of various sources of water contributions to the wells; and (5) predict temperature variations of the recovered water in the wells. Furthermore, sensitivity analysis was carried out to analyse any effects on the infiltration capacity from uncertainties in pond conductance and low hydraulic conductivity of the sand layer. The results show that the proposed MAR system is robust enough to produce the desired abstraction with acceptable constraints.

Study area
The study area is a lake area located in the province of North Holland, about 32 km southeast of the city of Amsterdam (Fig. 1). The focus area is on the lake Loenderveen Plas where an artificial sandbank will be created for implementing a MAR scheme with a production capacity by 10 × 10 6 m 3 /year. Loenderveen Plas is located adjacent to Waterleidingplas which is currently used as a reservoir for the eastern drinking water production plant. The lakes herein are termed as 'water supply lake 1' (WL1) for Waterleidingplas and 'water supply lake 2' (WL2) for Loenderveen Plas in this research, for simplicity ( Fig. 1).
The study area is characterised by the presence of several lakes and polders (Fig. 1). The regional flow in the study area is controlled by recharge at an ice-pushed low sandy hill ridge in the east called the Utrechtse Heuvelrug, flowing into deep polder Groot Mijdrecht in the west (Van Hateren 2019). Utrechtse Heuvelrug is also called Utrecht Hill Ridge with a maximum elevation of 68 m above the Amsterdam mean sea level (called Normaal Amsterdams Peil, NAP). The local flow is influenced by the water bodies with higher water levels, like Amsterdam Rhine Canal (ARC), the river Vecht, and the polders with low water levels like Horstermeer and Bethune polder (Fig. 1). The water bodies are connected through culverts, weirs, and pumps, and water levels are measured, controlled and fixed with a targeted minimum and maximum level by the water authority Waternet, in order to maintain the groundwater levels throughout the year and especially to prevent damage in the surrounding area used for settlement and agriculture. ARC is a canal with high traffic which is more polluted and flows through a lined canal; thus, the canal sides are considered impermeable. River Vecht is a natural river and has higher bottom conductance than ARC. The target water levels in the lakes and polders maintained by the water authority are given in Table 1.
Because of intensive drainage by underground drains in polders and the tight control of surface water levels, the water table is very shallow, less than 0.5 m below polders and lakes, and about 1.5 m underneath the roads and canals/channels.
The climate in the Netherlands is affected by the North Sea and is quite uniform in most parts of the country. The precipitation measured from 2012-2021 at the station De Bilt (station No. 260, about 15 km south of the study area) indicates that the average annual precipitation over the research area is 845 mm with a minimum of 582 mm and a maximum of 948 mm. Precipitation in the Netherlands is classified as a relatively wet season of autumn and winter, and a drier season of spring and summer. Average annual temperature in the research area is about 10.4 °C as recorded in the same station for the same time period. The temperature in the winter can drop below 0 °C (retrieved from the Koninklijk Nedarlands Meteorologisch Instituut KNMI 2022).

Data collection
The initial model for the research area was prepared by Van Hateren (2019) for the study of deep infiltration wells. The data for the model input files were collected from the Waternet database. The hydrogeological units from the regional model REGIS II (DINOloket 2020) were used to define the aquifer and aquitard layers of the model. Topographical, hydrological, and climatological data required for setting up the model were acquired from different national websites. For the groundwater heads, 23 observation wells with screens at different depths and with multiannual observation records were taken to ensure the well-rounded average head in the observations. Polder discharge measurements were directly taken from pumping station records when available. Some discharges were also calculated using water balance or using the weir formula if any measuring weirs were present. Surface-water levels, discharges and groundwater head were also collected from the Waternet database. Table 2 presents the data collected and their sources.

Hydrogeologic setting of the study area
The hydrogeological properties of the aquifers depend on their sedimentary origin. The hydrogeological information about the study area is based on the national subsurface model REGIS II and the BRO GeoTOP models on the website of DINOloket-Data and Information of the Dutch Subsurface (DINOloket 2020). The hydrogeological information is schematically visualised in Fig. 2.
The hydrogeological units were taken from borehole B31F0183 for lower aquifers and borehole B31F0231 for the details of upper aquifer layers in the database of the BRO REGIS II model. The bottom of the aquifer system consists of the Maassluis Formation (with alternating clayey and sandy units), which is considered the hydrogeological base for the study area. This layer is located between NAP -144 and -204 m, coarsening upward, with the bottommost part consisting of tough impermeable marine, light-to-dark-grey clay (either silty or sandy), calcareous materials and micas. The fluvial Peize-Waalre Formation is a sandy complex deposited alternately by the processes associated with the Peize or Waalre Formation and lies above Maassluis Formation. This layer is located at the depth from NAP of -60 to -144 m and is the most permeable layer of the overall regional system. The Waalre Formation consists of clay with a mixture of fine to coarse sand and a hint of gravel. These are river deposits with fine layers at the top and becoming coarser as moving downwards. This layer is about 8 m thick exceeding from NAP of -60 to -48 m. From this layer to the ground surface, an aquifer system is present. The Sterksel Formation has a clayey and a sandy unit which is a coarse-grained river deposit, running from NAP -34 to -48 m; this layer contains some sandy clay deposits, but only local in some areas. The Urk Formation is present from NAP -17 to -34 m and consists of moderately fine to coarse sand particles with a high mineral content. The Drente Formation extends from NAP -13 to -17 m; it originated from the glacial and preglacial age and consists of clay and loam deposits with sand and gravel, and the deposits are layered ranging to meters in thickness. The second layer from the top is the aeolian Boxtel Formation between NAP -4 and -13 m; this contains sand, ranging from fine to coarse but locally it may also contain silt and loam, and these deposits are finely sorted and contain fossils. The Holocene top layer (NIHO Formation) originates from the late and middle Holocene periods and consists of alternate sand and clay ranging from fine to medium; locally, it contains a large amount of peat and clay at the surface. Since the thickness of the Holocene cover layer differs from one location to another, the transmissibility and permeability of the top layer show a large spatial variation. For the research area, the cover layer is thinner or nonexistent in the east and gets thicker in the west.

Identification of the water sources for recharge
With the intended MAR objective, careful consideration of raw water for recharge is essential. The surrounding surface-water bodies were evaluated as possible sources for recharge in the ponds. Water pumped from lakes and polders is favored as the source of raw water because it has better quality. A water budget study carried out for the water bodies in the surrounding areas suggests that when the target levels in the lakes have been maintained, the excess water that is pumped from the Loosdrechtse Plassen (LP) to the river Vecht is one of the possible sources for recharge. Calculations were made on the daily volume of recorded discharge from 1984 to 2019. In order to produce 10 × 10 6 m 3 /year, a discharge rate of 833,333 m 3 /month can be supplied to the recharge lakes. The calculations suggest that excess water from LP can, at maximum, support 43% to the total required recharge, which is 359,440 m 3 /month during the wet period, while at minimum, they provide only 1.2% of the total, which is 10,310 m 3 /month during the dry period. In the absence of another source, ARC can supply raw water which has a high discharge of more than 25 m 3 /s and is stable all year round in terms of quantity. However, the quality of raw water is compromised after mixing with the water from ARC; therefore, raw water from two sources is to be pretreated by coagulation and sedimentation, a system similar to the existing pretreatment facility for seepage from the Bethune polder that currently feeds lake WL1. The monthly discharge variability in the excess water that can be used as a source from LP, along with the remaining water to be pumped from ARC, is given in Fig. 3. Figure 4a shows a cross-section of the site before constructing a MAR system. The exact location for the cross-section relative to the study area is added in Fig. 1. The Holocene top clay layer (NIHO Formation) at the site should be removed before a sandbank is placed. Removal of the Holocene clay layer increases the permeability for effective infiltration. The sandbank will then be placed 100 m east of the existing dike between the two lakes WL1 and WL2 (Fig. 4b). A small lake thus created in between the preexisting dike and the sandbank preserves the beautiful ecology in the banks of the dike. The aeolian sand deposits of the Boxtel Formation underneath the Holocene layer in WL1 are to be excavated to build the sandbank. This sand of the aeolian deposits has similar properties to the sand in the AWD system (Liu et al. 2019). The sand layer acts as a filter when the water is recharged to the unconfined aquifer underneath. The designed MAR system would consist of a number of recharge ponds for infiltration and pumping wells placed at suitable depths for recovery. Multiple smaller infiltration ponds were preferred over one or two large ones, considering the maintenance needs of the MAR system when the ponds are clogged. The use of underground drains to recover the recharged water was also tested. However, despite the feasibility and advantages of using drains to maintain consistent drawdown, wells were selected in preference, when considering the existing landscape, available construction area and risks associated with installation, breakdown and maintenance of drains as a recovery system. The effectiveness of wells is already verified in the AWD as they are used to recover water from the deep aquifer system (Olsthoorn and Mosch 2020).

Conceptual design of the managed aquifer recharge scheme
The shape of the sandbank, the size and arrangements of ponds, and the water level in the ponds largely determine the infiltration capacity of the MAR system. The hydrogeological conditions, and the locations and abstraction rates of the wells also influence the infiltration rate and recovery efficiency. In the current location, the infiltration ponds have to be positioned such that there is sufficient thickness for the unsaturated zone and the bottom level of the ponds should be kept above the surface-water level of surrounding lakes, thereby facilitating the dewatering of ponds during the time of maintenance. The criteria used for the selection of the conceptual design are given in Table 3.
The constraints were carefully chosen considering the quality of recovered water, the safety of neighbouring assets due to changes in water level, and the experience of artificial recharge and recovery in AWD. For the design, it was preferred that the well screens be placed in an aquifer layer with large hydraulic conductivity and native groundwater of good quality. Large hydraulic conductivity ensures better yield and smaller drawdown; therefore, recovery wells are to be placed in the Urk Formation.
The layout of the MAR system, shape and size of the ponds, the distance between the infiltration ponds and the wells, the number of wells and the water level in the ponds were designed by an iterative process to meet the design criteria. The final layout is shown in Fig. 5. For each pair of infiltration ponds, five wells were designed to recover the infiltrated water. The pumping rate of each well was calculated to be 1,096 m 3 /day, corresponding to the production of 10 × 10 6 m 3 /year from 25 wells in the MAR system. The conductance value of 0.65 m 2 /day per unit area of polygon for the infiltration ponds was calculated by using the average hydraulic conductivity of the Boxtel sand used for constructing the sandbank.
The required infiltration area of the ponds was then determined by the infiltration amount and infiltration rate to be maintained in the infiltration ponds. Lower values of infiltration rate were used for infiltration ponds to decrease the likelihood of clogging. The AWD system is used as a reference where 40 infiltration ponds are used to infiltrate at an average rate of 20 cm/day. The area of the infiltration ponds was iteratively increased to enlarge the infiltration capacity. Along with an increase in area, the water level in the ponds was also increased to meet the required rate and total infiltration capacity. For infiltration capacity to be 10 × 10 6 m 3 /year in the ponds with an area of 9,105 m 2 , the water level in the ponds has to be maintained at NAP +2.10 m, which was too high considering that the sandbank had to be elevated at NAP +3.00 m. The height of the sandbank is restricted by cost and changes in the landscape, which is a major concern for surrounding residents. An increase in the levels in the ponds will also cause infiltrated water to drain back to the lakes due to a steeper hydraulic gradient. Therefore, the area of the ponds was enlarged and the final adopted area of the infiltration ponds was approximately 11,432 m 2 per pond with the water level in the ponds maintained at NAP +0.50 m. The thickness of the pond was set at 0.8 m with the bottom of the infiltration ponds at NAP -0.30 m, which is higher than the maximum water level of the surrounding lakes, i.e. NAP -1.15 m for draining water during maintenance. To limit the capture zone of the recovery wells to more or less within the infiltration zone, the ponds were made curved at the ends.
The water arriving at the well should maintain a minimum residence time of 60 days, to ensure acceptable Fig. 3 Average monthly percentage of two water sources associated with the MAR scheme: Loosdrechtse Plassen (LP) and Amsterdam Rhine Canal (ARC). The scheme is aiming for production of 10 × 10 6 m 3 /year purification of the water (Olsthoorn and Mosch 2020). Multiple tests were carried out for the selected layout using tracer in the MT3DMS model to find the suitable aquifer layer and the horizontal distance of the wells from the infiltration ponds. Results suggested that, while the wells were placed in Urk Formation between NAP -22 and -36 m, and at 40 m horizontal distance from the ponds, the minimum residence time of 60 days could be maintained to abstract 10 × 10 6 m 3 /year.

Modelling objectives
The objectives of modelling in this study are multiple. First, the flow model was used to select a configuration of recharge ponds and recovery wells so that the production of 10 × 10 6 m 3 /year is achieved while the constraints are met. Second, the tracer transport model Fig. 4 Conceptual crosssections depicting the creation of recharge ponds and a recovery system using wells. a Hydrogeological cross-section of the existing situation, and b Conceptual layout of the proposed ponds and a recovery well. Information for the hydrogeological cross-section was derived using the BRO REGIS II v2.2 model (DINOloket 2020) and modified for additional features Preferably within the infiltration ponds Pond water level 0.50 m above NAP, considering the surrounding water level and preservation of the landscape was used to compute fractions of source waters contributing to the wells and the recovery efficiency. Third, the heat transport model was used to check temperature variations in the recovery wells so that post-treatment can be optimized.

Conceptual hydrogeological model
A large area of about 120 km 2 was chosen to set up the regional flow model, and the site area for MAR construction is ~6 km 2 located in the centre area. The model boundaries were chosen far away from the site so that the influence of boundary uncertainty on the site is negligible. Since the regional groundwater flow is from the eastern high-elevation area to the western deep polder, the boundaries on the east and west were defined as specified head boundaries. The northern and southern boundaries that divide Horstermeer and Bethune polders into half are considered groundwater divides, and thus are defined as no-flow boundaries. Hydrogeological formations were identified from the regional hydrogeological model of the Netherlands and were simplified to form eight layers based on their hydraulic conductivities. The first layer and seventh layer consist of mainly clay and silt and the rest are sandy layers (Table 4). Groundwater inflow components include: recharge from precipitation; leakage from the canal, river, and lakes; and lateral inflow from the east boundary. Groundwater discharges mainly to subsurface drains in the polders and lateral outflow through the west boundary.
The Groundwater Modelling System, GMS, software version 10.4 (Aquaevo 2021), was used for constructing the hydrogeological conceptual model. In GMS, boundary coverage for boundary conditions, layer property coverage for hydraulic parameters, recharge coverage for precipitation infiltration, river coverage for the canal, river, lakes and drains, and observation coverage for groundwater level measurements, were constructed. The data in the conceptual model were mapped and transferred to the numerical model automatically. -2005- (Harbaugh 2005) was used to construct the groundwater flow model. The initial regional groundwater flow model was developed for simulating injections to the deep aquifer in the Loendeerveen area (Van Haterent 2019). This model was modified in this study, whereby the regional model consists of a uniform grid size of 100 × 100 m. The grid was refined to 10 × 10 m in the MAR site (Fig. 6). The model consists of eight model layers corresponding to eight hydrogeological formations. The top of the model Fig. 5 Proposed configuration of the recharge ponds and recovery wells in GMS model. The area per pond is 11,431.6 m 2 and the total area of 10 ponds is 114,316 m 2 . The 25 recovery wells are located in the layer 4 of the model, which is the Urk Formation represents the surface elevation and the bottom is bounded by the Maassluis Formation which is considered the hydrogeological base for the study area. The general properties of the eight model layers are presented in Table 4.

Model inputs
First, a steady-state groundwater flow model was constructed and calibrated with MODFLOW. Since surfacewater levels are maintained at target levels and no groundwater abstractions take place, the groundwater system is in an equilibrium state. The average annual precipitation in the last 10 years has been 845 mm, as recorded at the De Bilt weather station. A filtration coefficient of 0.432 is usually used to compute an average net recharge of 1 mm/day (Gehrels 1999). This value was assigned in the recharge package (RCH). Recharge is important in the eastern part where precipitation infiltrates into the aquifer from elevated hill ridge, while in the western part, due to the presence of surface-water bodies, surface-water level controls the groundwater head. The surface-water bodies including the polders, lakes, river and main canals were simulated using the River package (RIV). Surface-water levels were taken as average values (Table 1). Average groundwater levels were assigned to the specified head cells on the east and west boundaries. The aquifers are heterogeneous and anisotropic in the vertical direction. The top layer consists of clays in most areas with very low hydraulic conductivity values (Table 4). The Boxtel sand is rather homogeneous with a low hydraulic conductivity value. Layer 4 consists of Urk coarse sand and gravels with very high hydraulic conductivity values where the recovery wells are to be placed. For layer 7, most of the area is covered by clay with very low hydraulic conductivity, and a high hydraulic conductivity value is assigned to the area with clay absent.
Second, a transient groundwater flow model was constructed. The simulation time covers a period of 10 years from 2012 to 2021 with monthly stress periods. The computed heads from the steady-state model were assigned as initial heads. Monthly groundwater recharge was computed from monthly precipitation series using the same infiltration coefficient. Monthly water levels were assigned to the lakes and polders.
Once the model was calibrated, the proposed MAR system (Fig. 5) was implemented in the model. The artificial sandbank for the MAR system was introduced by changing the hydrogeological properties of the area to be covered/replaced by sand. The infiltration ponds were simulated using the General Head Boundary (GHB) package. Recovery wells were introduced in model layer 4 and simulated with the Well (WEL) package. The MODLOW model was able to compute infiltration rate from the ponds, effects on water budget, and changes in the groundwater head caused by the proposed MAR system.
Furthermore, MODPATH was used to track the particles backwards from the abstraction wells to delineate the capture zone of the wells. The MT3DMS model (Zheng and Wang 1999) was used to identify the minimum residence time, the mixing of different sources of water in the well, and the recovery efficiency of the MAR system. A hypothetical tracer of arbitrary concentration was added to the cells of the infiltration ponds and the advective transport was simulated with the MT3DMS model. The arrival time of the tracer at the abstraction wells was considered the minimum residence time of the infiltrated water from the ponds. The advective transport package was selected which simulates the transport of the tracers at the same velocity as groundwater. The third-order TVD scheme (ULTIMATE) was used for the advective transport simulation in order to produce smooth breakthrough curves.

Model calibration
The groundwater flow model was calibrated under steady-state and transient conditions. The observed heads and discharges were used for the model calibration. For the steady-state model calibration, average groundwater levels from the 42 piezometers were computed and used for the model calibration. In order to improve model reliability, measured discharges from lakes and polders were used to compare computed discharges from the same polders and lakes. For the transient model calibration, monthly groundwater levels measured in 2012-2021 were used to compare the computed heads.
The parameters that were targeted for calibration included hydraulic conductivity, specific yield and storage, and riverbed conductance. A trial-and-error method was used for the calibration. The parameter values were changed to achieve the calibration target. Statistics for evaluating the calibration included the mean residual error, root mean squared error (RMSE), and coefficient of determination (R 2 ).

Determination of sources of water contributing to the recovery well
The advective transport in the MT3DMS model was used to compute fractions of source waters contributing to the recovery wells. There are potentially four sources of water contributing to the recovery well: (1) infiltration pond water; (2) surface water; (3) recharge from precipitation; and (4) natural groundwater reserve. The fractions can be solved with the following equations: where f a = Q a /Q recovery , f b = Q b /Q recovery , f c = Q c /Q recovery an d f d = Q d /Q recovery are the fractions of sources a, b, c and d to the recovery well when the break-through curve approaches steady concentration (C 1recovery , C 2recovery , C 2recovery ). The coefficients, C 1a …, C 3d , are hypothetical tracer concentrations at the four sources, respectively.
The fractions of the four source waters can be obtained by solving the equations as: In this study, three tracer transport models were run for 100 years until the break-through curves approached a constant value. For each model, a tracer concentration of 100 mg/L was assigned to one source water, while concentration of other sources was kept zero. The fractions of the four source waters contributing to the wells were then calculated by using the mixing mass balance (Eq. 5).
The recovery efficiency refers to the effectiveness of the MAR system in recovering infiltrated water from the ponds. The total infiltration rate from the ponds can be obtained using the MODFLOW water budget. The amount of infiltrated water contributing to the recovery wells can be determined by tracer mixing (Eq. 5); therefore, the recovery efficiency of the MAR system can be calculated as: The amount of infiltrated water from the ponds recovered by the wells The total infiltration rate from the ponds

Simulation of temperature variations
MT3DMS was also used for heat transport simulation by considering the temperature as one of the species. MT3DMS does not consider density and viscosity effects. Since the density effect can be neglected for a fresh groundwater system, ignorance of the viscosity effect will compute the fastest possible travel times, since the groundwater temperature in the area is below 25 °C (Liu et al. 2019). Required inputs for the simulation of heat transport using MT3DMS were calculated as follows. Thermal distribution factor: Heat retardation factor: where C p solid and C p fluid are the specific heat capacity of the solid and fluid respectively; ρ is the bulk density of an aquifer; n is the porosity of the aquifer; K d _ temp is the thermal distributive factor. The heat exchange between solids and water was simulated with the sorption reaction package in MT3DMS. The heat conduction was simulated with equivalent diffusive transport in MT3DMS. Bulk thermal diffusivity was calculated as: where D m _ temp is the bulk thermal diffusivity; K T bulk is the bulk thermal conductivity calculated using K T fluid (the thermal conductivity of the fluid) and K T solid (the thermal conductivity of the solid); θ is the porosity; ρ is the bulk density of the aquifer; C p fluid is the specific heat capacity of the fluid. The parameter values for the heat transport model are given in Table 5. References used for selecting parameters are found in DINOloket (2020), Witte et al. (2002), Langevin et al. (2008) and Caljé (2010).
The retardation factor R d is calculated with the equation: For heat simulation, regular temperature measurements are important for the source water, the surrounding surface water as external heat sources and the groundwater temperature of the aquifer layers. Measured temperature data in lake WL1, which could represent the temperature of surrounding surface water, was measured only once a month. Similarly, the temperature measurement at ARC was limited as well. The measurements were not sufficient for heat transport modelling. Liu et al. (2019) found that in AWD, the infiltration ponds are usually shallow and variation is mainly caused by atmospheric temperature changes. Therefore, for this research, temperature input in the infiltration ponds was represented by the measured temperature value in the ponds of the AWD due to their common shallow nature and the influence of atmospheric temperature, which was similar at both pond locations. The temperature measured in the surface water of the AWD between 2014 and 2020 shows that the maximum temperature can be as high as 29 °C during summer and drop below 0 °C during winter. For the aquifer, the temperature was taken from the two data loggers recording hourly from years 2004-2020 at two different depths, of NAP -31.9 and -51.1 m. Average monthly temperature for both the loggers f luctuated between 9 and 15 °C; thus, an average groundwater temperature value of 11.6 °C was used in the model. The transient input of temperature for the infiltration ponds was specified in the GHB package. Figure 7a shows a scatter plot of the computed versus measured heads indicating that most points fall on the diagonal line. The mean residual error is 5 cm, the RMSE is 9.87 cm, and R 2 is 0.98. All statistics show a good fit of the computed heads to the observed heads in the 42 piezometers. Figure 7b shows a scatter plot of the computed versus measured discharges. The computed discharge values also fit very well with the measured discharge values. The RMSE is 953 m 3 /day, and R 2 is 0.99. The good fit to the measured discharges provides confidence that the model can be used for the MAR simulation. Table 6 presents the groundwater budget under the current conditions before the MAR implementation. Groundwater recharge comes from surface water leakage (46%), precipitation infiltration (43%), and subsurface inflow from the east boundary. Groundwater discharges mainly to the polders and lakes (80%), and a small amount of water flows out through the west boundary (20%). Figure 8 plots two fitting curves of computed monthly groundwater-level series to observed series from 2012 to 2021. It can be seen that computed heads match well with the observed heads. Both of the observation wells shown in Fig. 8 are located in the proposed MAR site. Observation well 31F0234 is located in the model layer 4 (where recovery wells are to be placed). The mean residual error is 3 cm and the RMSE is 5 cm. The range of the seasonal variation is about 15 cm. Observation well 31F0231 is also located in the model layer 4. The mean residual error is 1 cm and the RMSE is 5 cm. The range of the seasonal variation is about 10 cm. Both curves show stable seasonal fluctuations of groundwater levels over the past 10 years. Figure 9 plots monthly groundwater budget components in the last 10 years. Leakage from surface water is rather stable because the water levels are tightly controlled. Recharge from precipitation shows seasonal fluctuations: more recharge in autumn and winter wet seasons and less recharge in spring and summer dry seasons. Groundwater discharge to polders and lakes fluctuates in response to groundwater level changes. Net discharge through boundaries is small and stable. The change of storage as in Fig. 10 shows an increase in wet seasons and a decrease in dry seasons, and the mean of storage changes is almost zero, indicating a dynamic equilibrium of the groundwater system. Therefore, the steady-state model based on the average values of time-dependent inputs over the last 10 years is good enough for the MAR simulations.

Performance of the proposed MAR system
The steady-state flow model was used to assess the performance of the proposed MAR scheme (Fig. 5).
The total abstraction rate from 25 recovery wells was assigned to the target production of 10 × 10 6 m 3 /year (27,400 m 3 /day). A local water budget zone was defined at the grid refinement area to compute the infiltration capacity and water exchanges. The local water budget zone covers the MAR scheme. The computed water budget is presented in Table 7. Under the designed production target (27,400 m 3 /day), the MAR recharge from infiltration ponds amounts to 18,146 m 3 /day, which is about 66% of the total production. The other main sources of water contributing to the wells are leakage of surface water and recharge from precipitation. Water exchanges with the regional model are very small. The corresponding infiltration rate of the ponds was computed to be 16 cm/day, which is 4 cm/day smaller than the infiltration rate in the AWD area.
The drawdown in the unconfined aquifer influenced by the proposed MAR scheme does not exceed 5 cm outside the lakes (Fig. 11a). The maximum drawdown of 1.97 m was found at the location of the wells in model layer 4. All pathlines starting from the recovery wells end at the edges of the infiltration ponds indicating that the capture zones of the recovery wells are within the MAR site (Fig. 11b).
The residence time varies depending on the length of pathlines from the filtration ponds to the wells. A hypothetical tracer of 100 mg/L was assigned in all infiltration ponds and was moved with the advective transport for 10 years. Break-through curves at 25 recovery wells were analysed. The results indicated that on average, the minimum residence time, measured by the arrival of the tracer at the wells, was around 65-70 days. The minimum residence time associated with each well exceeded the criterion of 60 days. The average break-through curve for the 25 wells is plotted in Fig. 12, where it can be seen that the average tracer concentration from the 25 wells approaches a constant value of 65.3 mg/L after 4 years, indicating that the fraction of infiltrated water contributing to the wells is ~65.3%.
The fractions of the different sources of water contributing to the wells are given in Table 8. The contribution from infiltration ponds is about 65.3% of the total abstracted water. Leakage from lakes accounts for 13.1%, whereby ~14.9% comes from groundwater reserves, and a small amount of ~6.7% comes from precipitation recharge.
Comparing Tables 7 and 8, one can see the water inflow from infiltration ponds into the aquifer is at the rate of 18,146 m 3 /day, out of which 17,892 m 3 /day was captured by the wells. The recovery efficiency of the MAR system can then be calculated using Eq. (6) as 98%.

Sensitivity analysis of the MAR system
From the calibration of the steady-state and transient models, it was found that the recharge from precipitation and discharges to polders and lakes were quite accurately estimated and they had small influences on the MAR scheme. However, Boxtel sand determines the efficiency of the MAR scheme. From the cross-sections in Fig. 4 and Table 4, one can see that the hydraulic conductivity of the Boxtel sand below the ponds is low, which controls the vertical groundwater flow and infiltration rate of the ponds. Since the ponds will be created over the excavated Boxtel In the simulation model, the hydraulic conductivity of the Boxtel sand was given a value between 4 and 5 m/day and the corresponding conductance value was 0.65 m 2 /day/m 2 . A 50% increase in the hydraulic conductivity, i.e., 6-7.5 m/day, results in conductance of 0.975 m 2 /day/m 2 . These values are similar to values used in AWD area and were simulated as an optimistic scenario. Furthermore, with the operation of the ponds for some years, pond conductance might decrease. Therefore, a 50% decrease of pond conductance was simulated as a worst-case scenario. The results are summarised in Table 9.
The sensitivity results show that the infiltration rate depends on the hydraulic conductivity; the larger the hydraulic conductivity value is, the higher the infiltration rate. For all cases, the fraction of water from the infiltration ponds contributing to the wells is almost the same. The recovery efficiency decreases with the increase of infiltration rate. Furthermore, the percentage change in the infiltration rate is smaller than the percentage change in the hydraulic conductivity, demonstrating that the proposed MAR scheme is robust.

Temperature variations in the recovered water in wells
In the heat transport model, weekly average temperature measurements in the AWD ponds were used as inputs for recharge, surface water and infiltration ponds. The model was run from January 2016 to March 2020 with a weekly stress period. Figure 13 shows a comparison between input temperature series and computed average temperature series from the 25 recovery wells. The input temperature has a mean of 10.9 °C with a minimum of -3 °C and a maximum of 22.9 °C. The average temperature of water captured in the wells is ~11.1 °C, fluctuating between 9.8 and 12.5 °C. The results show that the peak temperature in the wells was delayed with a lag time which represents the slow transport of the heat. The lag time of the peak was computed with the cross-correlation analysis, resulting in 210 days. With a retardation factor of 2.2, the average residence time of sources of water contributing to the wells is ~95 days, which is sufficiently long to improve the water quality. The abstracted groundwater has a more stable temperature with small variations due to heat conduction and storage in the aquifer. The stable temperature helps in improving the efficiency of the slow sand filter used for post-treatment in the drinking-water production plant. Especially in wintertime, the mixed water sources will have a higher temperature which will increase hydraulic conductivity of the slow sand filters and thus maintain productivity. Fig. 12 Average breakthrough curve for a tracer moving from the infiltration ponds and arriving at the 25 abstraction wells over more than 4 years

Considerations for implementation
For the successful implementation of the proposed MAR, the following issues must be considered: 1. The hydraulic conductivity of Boxtel sand, used to create the sandbank in the model, determines the infiltration rate. The sensitivity analysis shows that the infiltration capacity can increase by 18% with an increase of the hydraulic conductivity of the Boxtel sand by 50%. Insitu tests should be carried to ascertain the Boxtel sand hydraulic conductivity values. 2. Maintenance of the ponds is also important. The sensitivity analysis shows that a decrease of the pond conductance by 50% due to clogging will reduce the infiltration capacity by 9%. 3. Beaches between the old dike and the sandbank should be carefully constructed to create an ecological environment for aquatic species and birds. This ecological environment is important for nature conservation. 4. Systematic measurements of temperature should be carried out in the lakes, river/canal, and observation wells with daily frequency. Temperature measurements in a year can be used to calibrate heat transport models and further improve model results including viscosity effects (Liu et al. 2019).

Conclusion
A unique MAR system is proposed which consists of 10 infiltration ponds created on a built-up sandbank and 25 recovery wells which are located in the productive aquifer. The shapes and dimensions of the ponds and locations of the wells were determined with an integrated modelling approach according to the design criteria. The integrated modelling approach was found to be effective for the MAR assessment and included the following five steps. First, steady-state and transient natural groundwater flow models with local grid refinement in the MAR site were constructed and calibrated with measured heads and discharges. Second, the steady-state model with the MAR system was constructed. The infiltration ponds were built in the MAR model with the General Head Boundary package of  . 13 Comparison of the input temperature in the infiltration ponds and the output temperature at the abstraction wells MODFLOW. 25 pumping wells were located in the productive aquifer with large hydraulic conductivities and simulated with the MODFLOW Well package. The MAR flow model computed the infiltration rate from the ponds and the drawdown caused by abstraction. Third, MODPATH was used to delineate the capture zones of the 25 wells to ensure the capture zones ended at the boundaries of the ponds. Fourth, an advective transport model was constructed with MT3DMS and used to simulate the movement of a hypothetical tracer from the ponds to the wells to determine the travel time. The first appearance of the tracer at the well indicated minimum travel time. The transport model was also used to simulate three different tracers from four sources of water to arrive at the 25 wells. Tracer mixing mass balance equations were created and solved to determine the fractions of sources of water contributing to the wells. The results were also used to compute the recovery efficiency. Fifth, a heat transport model was constructed with MT3DMS using equivalent heat transport parameters. The large seasonal variations of temperature in the infiltration ponds and surface-water bodies were smoothed by the aquifer and resulted in a stable temperature in the abstracted water.
The results showed that the proposed MAR system is capable of producing 10 × 10 6 m 3 of water per year. The major source of water contributing to the wells comes from the infiltration ponds (65.3%). Other sources of water include leakage from lakes (13.1%), natural groundwater reserves (14.9%), and recharge from precipitation (6.7%). The recovery efficiency was computed as 98%, which means 98% of the infiltrated water from the ponds is recovered by the wells. The computed travel times from the ponds to the wells are more than 60 days meeting the quality requirement. The drawdown outside the MAR site is less than 5 cm, satisfying the environmental impact requirement. The temperature variation in the abstracted water varies between 9.8 and 12.5 °C which is beneficial for post-treatment. The conductance of the ponds and the hydraulic conductivity of the Boxtel sand should be determined with pumping tests during the well construction. Furthermore, the constructed sandbank and designed infiltration ponds will create a unique nature area for the ecosystem where birds will find a resting place on the beaches surrounding the sandbank and unique flora and fauna may emerge.