Impact of river reconstruction on groundwater flow during bank filtration assessed by transient three-dimensional modelling of flow and heat transportImpact de la reconstruction d’une rivière sur l’écoulement des eaux souterraines via la filtration sur berge évalué par un modèle tridimensionnel

Bank filtration (BF) is an established indirect water-treatment technology. The quality of water gained via BF depends on the subsurface capture zone, the mixing ratio (river water versus ambient groundwater), spatial and temporal distribution of subsurface travel times, and subsurface temperature patterns. Surface-water infiltration into the adjacent aquifer is determined by the local hydraulic gradient and riverbed permeability, which could be altered by natural clogging, scouring and artificial decolmation processes. The seasonal behaviour of a BF system in Germany, and its development during and about 6 months after decolmation (canal reconstruction), was observed with a long-term monitoring programme. To quantify the spatial and temporal variation in the BF system, a transient flow and heat transport model was implemented and two model scenarios, ‘with’ and ‘without’ canal reconstruction, were generated. Overall, the simulated water heads and temperaturesmatched those observed. Increased hydraulic connection between the canal and aquifer caused by the canal reconstruction led to an increase of ~23% in the already high share of BF water abstracted by the nearby waterworks. Subsurface travel-time distribution substantially shifted towards shorter travel times. Flow paths with travel times <200 days increased by ~10% and those with <300 days by 15%. Generally, the periodic temperature signal, and the summer and winter temperature extrema, increased and penetrated deeper into the aquifer. The joint hydrological and thermal effects caused by the canal reconstruction might increase the potential of biodegradable compounds to further penetrate into the aquifer, also by potentially affecting the redox zonation in the aquifer.


Introduction
Bank filtration refers to a process in which river water infiltrates at the riverbed, migrates through the aquifer and is finally extracted at pumping wells that are located adjacent to the river (Hiscock and Grischek 2002;Kuehn and Mueller 2000). It is a pretreatment option often used for conventional drinking water supply. In Europe, bank filtration is widely used, and contributes 16% of the potable water supply in the whole of Germany; it may even dominate as a treatment option locally (Hiscock and Grischek 2002;Massmann et al. 2008a). Along the subsurface flow paths from the river towards the pumping wells, raw water quality can be improved through several beneficial attenuation processes such as equilibration of the temperature variation, pathogen filtration, mixing, biodegradation and sorption; those processes act on suspended solids, biodegradable compounds, other dissolved organic carbon (DOC) or bacteria, viruses, parasites, etc. (Hiscock and Grischek 2002;Sharma et al. 2012). The quality of water gained via bank filtration strongly depends on the subsurface capture zone, and the mixing ratio between river water versus ambient groundwater, as well as the spatial and temporal distribution of subsurface travel times (Hiscock and Grischek 2002). For all properties mentioned in the preceding, the hydraulic conductivity of riverbed materials plays a decisive role (Hiscock 2005;Massmann et al. 2008a;Vogt et al. 2009). The formation of a clogging layer is almost inevitable during the long-lasting operation of bank filtration systems (Schubert 2006). It leads to a decrease of the riverbed hydraulic conductivity, limiting the infiltration rate of surface water and thus to a decreased recoverable ratio of surface water at the production wells (Goldschneider et al. 2007;Tufenkji et al. 2002;Ulrich et al. 2015); however, it would not decrease infinitely but rather reach a state of equilibrium (Gutiérrez et al. 2018). The formation of a clogging layer is commonly caused by physical (sedimentation of fine-grained particles) and biogeochemical (growth of microbes, alga and biofilms) factors, while the dynamic hydrological condition also plays an important role (Ghodeif et al. 2016;Gutiérrez et al. 2018;Rinck-Pfeiffer et al. 2000;Schubert 2002). Natural declogging processes (decolmation) involve an increase in flow velocity giving rise to removal and restorage of river sediments and, in part, by digging activities of the animal stock at the river bottom. In managed aquifer recharge schemes, the clogging layer is regularly abraded to reestablish a higher hydraulic conductivity (corresponding recharge rates may increase by the factor of 10; Greskowiak et al. 2005;Escalante et al. 2015), although artificial excavation in surface-water beds at bank filtration sites has rarely been studied.
Furthermore, the quality of water gained via bank filtration depends on subsurface temperature patterns (Munz et al. 2019). Temperature is a key factor controlling redox reactions, enzyme activity and biodegradation of organic substances (Maeng et al. 2010;Miettinen et al. 1996;Umar et al. 2017). Besides the short-term fluctuations of the daily temperature signal in the riverbed itself (e.g. , seasonal temperature changes penetrate to greater aquifer depth (Kalbus et al. 2006), strongly affecting the bioactivity of the microorganisms (Massmann et al. 2008b). If the subsurface temperature falls below 14°C, biodegradation has been shown to be strongly constrained (Massmann et al. 2006;Prommer and Stuyfzand 2005). At an artificial recharge and bank filtration site at Lake Tegel in Berlin and at the site used for this study in Potsdam, Germany, an aerobic condition develops between the river and the pumping well in winter, while in summertime the oxygen is consumed within the first meter below the riverbed followed by the prevailing anaerobic condition (Munz et al. 2019;Gross-Wittke et al. 2010;Massmann et al. 2006).
Numerical models are widely adopted in water-resource studies (Umar et al. 2017) and are the most common tool for understanding the flow regime at a bank filtration site. By recalibrating a three-dimensional (3-D) groundwater flow model against pumping-test results, Schafer (2006) confirmed riverbed colmation at a bank filtration site at the Ohio River (USA) several months after commissioning. Shankar et al. (2009) quantified the contributing ratio between infiltrating water and natural groundwater to the production well, and also delineated the catchment area. May and Mazlan (2014) evaluated the change of the groundwater flow field under the influence of heavy abstraction along Langat River, Malaysia, by comparing the scenarios with and without pumping activities in a 3-D numerical model. By considering seasonal temperature variations in reactive transport studies, more accurate simulation results can be achieved for water quality changes (Prommer and Stuyfzand 2005) as well as for variation of O 2 and NO 3 − during bank filtration (Henzler et al. 2016).
Changes in the flow field (e.g. by colmation or decolmation) will influence the heat transport through the riverbed and further influence the temperature distribution along the subsurface flow path to pumping wells, but this issue has not been investigated quantitatively in the field so far. The objective of this research was to determine the hydraulic connection between the river and groundwater at a bank filtration site with substantial geological heterogeneity and further investigate the temporal changes caused by a river reconstruction (artificial deepening by removal of till sediments and river bank renewal). A 3-D flow and heat transport model has been calibrated first against water heads and then against temperature data to (1) investigate the groundwater flow regime, (2) quantify the changes of infiltration rates, the mixing ratio between river water versus ambient groundwater, and the spatial and temporal distribution of subsurface travel times, and (3) quantify the changes in the subsurface temperature distribution brought on by the reconstruction. The coupled water flow and heat transport simulations have been performed as a transient 3-D model for a period of 4.5 years, delivering spatiotemporal distributions of the hydraulic and thermal properties.

Site description
The study site is located in Potsdam, Germany (Fig. 1a), where the pumping wells of a local waterworks are operated south of a manmade canal named 'Nedlitzer Durchstich' (ND) between two lakes "Weisser See" [Lake Weisser] and "Jungfernsee" [Lake Jungfern] (Fig. 1b). The study site covers a total area of around 6.9 km 2 , with a maximum extension of 3,484 and 3,256 m in the north-south and east-west directions, respectively (Fig. 1c). Average surface elevation at the river is about 30 m a.s.l. (meters above sea level), while hilly   Fig. 1 a Location of Potsdam in Germany. b Location of observation wells and pumping wells beside the ND canal (yellow area, divided based on the effective schedule of reconstruction in the ND). c Study area with groundwater equipotential lines (m a.s.l.) of the uppermost aquifer, location of pumping and remediation wells, and the location of all borehole profiles used to design the geological model structure. Background map of DEM was adopted from NASA's Land Processes Distributed Active Archive Center (LP DAAC) (2019) accessed on 24 September 2019 areas in the northern part of the study site reach 80 m a.s.l. The reconstruction work is illustrated in Fig. S1a,b of the ESM.
The bank filtration takes place mainly through Quaternary aquifers, which are semi-confined aquifers stratified by glacial till layers from the later Pleistocene. Featured as rather impermeable material, the glacial tills, from the bottom to the top, were formed in the Saale I glacial period, Saale II glacial period and Weichsel I glacial period, but structure and extent are substantially different between them. The glacial tills of Saale-I (SI-GT) and Saale-II (SII-GT) cover most of the study area and join at some places forming thicker glacial blocks. The glacial till formed in the Weichsel I glaciation (WI-GT) is the youngest and uppermost glacial till in the study area, and it is not as continuous as the two former glacial deposits due to periglacial processes and erosion during the Holocene. Directly in the area below the ND, there is a lens of WI-GT, which covers half of the area of the canal bottom. Sand exists above the WI-GT, while mud and silt show up at the lake bottom. Stratified by semi-permeable glacial tills, the shallow aquifers could be divided into four sandy layers. At the bottom of the deepest of these Quaternary aquifers, a Miocene layer consisting of silt, fine sand, lignite, and clay exists, which separates these shallow aquifers from the deeper aquifers. Sapropel existed in the top few meters of Weisser See as well as Jungfernsee, which are located at the western and eastern end of the ND canal, respectively. Water exchange might occur at some locations due to the existence of hydrogeological windows possibly formed by erosion during the Pleistocene.

Observation wells and data collection
The observational network consists of 104 groundwater observation wells which were grouped into three categories according to their sampling frequency and sampling procedure: Group 1. Water heads were measured in 89 observation wells distributed within the study site- Fig. S2a of the electronic supplementary material (ESM)-in mid-April 2012. The filter screens are mostly 2 m long, while the bottom elevation varies from 19 to 34 m a.s.l. These water heads were used as the calibration targets for a numerical steady-state flow model. Group 2. Water heads and temperatures were measured continuously at 15 observation wells using automatic data loggers with time intervals varying between 5 min and 6 h (Rugged TROLL 100, In-Situ, USA) with an accuracy of ±0.01 m for water head and ±0.3°C for temperature. Measurements took place from June 2009 to December 2014. Most of these observation wells (except U100) are located between the canal and the pumping well group, and some of them are situated along the canal banks (Fig. 1b). Some of the observation wells (W03, W04, W07, and W09) comprise two adjacent wells with different filter screen elevations, which is denoted by HIGH (upper well screen) and LOW (lower well screen; Table 1). For the observation wells W03, W07 and W09, the data loggers in HIGH and LOW were set at the same elevation. These water heads were used as calibration targets for a transient flow model. Group 3. Representative groundwater temperatures were measured in 12 observation wells weekly or biweekly, over a period of 18 months from June 2013 to December 2014. Except for two additional measurements in W05 and W06, the sampling points are identical with group 2. These temperatures were used as calibration targets for the heat transport model.
The river-water temperature for the simulation period was retrieved from two stations. A local data logger was put into operation in the canal from October 2012. Temperature data before that date were derived from a gauging station around 6 km upstream, in a parallel branch of the River Havel flowing through the city of Potsdam. The two temperature time series showed good data consistency and thus were used in combination to define the canal temperature for the entire simulation period. Temperature amplitudes were calculated as half of the difference between maximum and minimum temperature.
Daily precipitation and air temperature (2 m above ground) were measured at the nearby Leibniz Institute for Agricultural Engineering and Bioeconomy (about 2 km west of the study site). Highest rainfall magnitudes and intensities were observed mainly in July/August each year. Maximum rainfall intensity was 42.8 mm/day. Also for high rainfall intensities, no surface runoff was noticed at the study site and thus also highly intense precipitation was assumed to recharge into the aquifer without substantial losses at the surface. Groundwater recharge was estimated based on a large-scale conceptual model (Landesamt für Umwelt Brandenburg 2019). In the study area, the evaluated ratio of (Precipitation minus Actual evapotranspiration)/ Precipitation (1991-2010) varied from 18.7 to 19.7%, and a rounded-up value of 20% was taken as a fixed recharge ratio in the modelling work. Correspondingly, the estimated yearly groundwater recharge ranged between 76.1 mm/year in 2011, as driest year, and 138.2 mm/year in 2010, as wettest year. The pumping rates at the waterworks slightly varied between 6,000 and 7,000 m 3 /day. In July 2014 a remediation pumping (~2,090 m 3 /day) located east of the pumping wells of the waterworks (Fig. 1c) was started, but was not related to the reconstruction measures. Overall, groundwater recharge in the catchment area of the waterworks plays only a minor role due to the large contribution of river bank filtration to the water balance.

Model structure
Groundwater equal potential lines were interpolated from synoptic head measurements of mid-April, 2011 (adapted from the environmental agency, MLUL Brandenburg) and used to define the model domain (Fig. 1c). About 1 km south of the pumping well group, a groundwater divide at a flat hilltop extends as a flowline at the hillslope towards the southern bank of Jungfernsee, forming a no-flow boundary at the south-eastern boundary. From the western end of the water divide, a flowline to the north, ending at the southern part of "Fahrländer See" [Lake Fahrländer], is taken as no-flow boundary of the model domain. Both are relatively far from the main area of interest between ND and the waterworks wells (Fig. 1b). The water levels of Jungfernsee and Fahrländer See were included in the model as fixed, but temporally variable water heads. The small northern area between these lakes is characterized by a hilly topography (Fig. 1c) formed by a glacial end-moraine. Here, local groundwater flow paths tend to directly end at the nearest surface water body and the groundwater divide approximately follows the surface-water divide delineating the northern groundwater subcatchment. Together, these boundaries define a model domain that includes more than the complete groundwater catchment area of the pumping wells and embeds the main area of interest between the pumping wells and the ND canal.
A total of 145 borehole profiles were available to construct the geological model (Fig. 1c). Most of the drilling locations are placed between the channel bank and pumping wells, while comparatively little information is available for the southern and western parts of the study site. Vertically, 67 of the drillings near the canal reached shallow depths only (mainly above 10 m a.s.l.), while 78 reached deeper layers, though only a few of them reached down to the Miocene layer (5.16 m to −12,8 m a.s.l).
The geological condition was implemented in a model based on the software FEFLOW 6.2 (Diersch 2014). Using this hydrogeological 3-D model (Fig. 2), steady-state and transient simulations of combined water and heat flow were performed also in FEFLOW. A 2-D triangular finite element mesh of 30,384 elements was generated for the top of the domain (Fig. S2b of the ESM). Local grid refinement was applied around the area of the ND, the observation wells and the pumping wells with the final element size varying between 5 and 10 m, aiming to accurately depict the relatively strong hydraulic gradient caused by pumping. With increasing distance to the ND the element size was enlarged to 300 m in the southern part. The generated mesh was smoothed (to meet the Delaunay criterion) and superimposed onto a 3-D finite element mesh consisting of nine model layers to represent the main hydrogeological units (Table 2). Eighteen extra grid node layers were implemented in between to improve the numerical accuracy and to achieve a better vertical resolution. The unsaturated layer above the groundwater was represented by a phreatic top layer in the model. Three cross-section of U100-W03-W04, Canal bank-W07 and Canal bank-W09 could be seen from Fig. S5 of the ESM.
In respect to heat transport, in FEFLOW, the volumetric heat capacity (c) and heat conductivity (λ) are defined separately for solid and water phase, while the total porosity (ɛ) and longitudinal and transversal thermal dispersivities (λ disp ) are set for the bulk porous medium. Total porosity mainly controls the volumetric ratio between the solid and liquid phases, as both phases contribute to heat transport, which is a difference to non-sorbing solute transport. The total porosity of the glacial tills and sapropel sludge were inferred from saturated gravimetric moisture content and grain density of sediment samples taken in the river section at the south end of the Fahrländer See. The average value of total porosity was 0.26. The thermal parameters for water are well known, while parameters of the solid phase could vary between sediments and have been subjected to calibration. The unsaturated zone on the top was set as 'phreatic mode' in FEFLOW, implying that its hydraulic conductivity adopted is taken as the saturated K f -value times the ratio between the saturated thickness and full element height.

Steady-state model
A steady-state model was run and calibrated with the aim of providing a consistent representation of the geological structure and their corresponding hydraulic parameters. The model calibration was based on the observed water heads at 104 locations distributed within the study site (data groups 1 and 2). A Dirichlet boundary of 29.5 m was chosen for the surface-water bodies. The pumping wells were set as multi-layer wells with a total pumping rate of 6,556 m 3 /day. The distribution of the total pumping rate to each layer in contact with the well screen was calculated based on the particular K f -values and the hydraulic head distribution in this layer. In the steadystate model, water level and pumping rate represent the mean of values measured between mid-March and mid-April 2012. Groundwater recharge was set as 0.3 mm/day, which represents 20% of the mean precipitation of 2011 and 2012. Average daily groundwater recharge volume at the land surface of the model area was about 1,780 m 3 / day and thus only about 21% of the average volume pumped daily (sum of production wells of the waterworks and nearby groundwater remediation wells), reflecting a high fraction of bank filtration at this site. All outer boundaries not representing surface waters and the bottom of the model were set as no-flow boundaries according to the configuration of the model domain. Saturated   (Doherty 2015) framework. Parameter zones of homogeneous hydraulic conductivity were defined based on the existing geological units. The parameter ranges used in the calibration are given in Table S1 of the ESM). The objective function was based on the observed water heads and the anisotropy ratio of hydraulic conductivities, ensuring vertical conductivities (K v ) being linked to the horizontal ones (K h , where K y = K x ). A maximum limit was set so that K h /K v must remain larger than 2 in all geological units, according to Todd and Mays (2005). With respect to data quality and spatial relevance, different weighting factors were assigned to the hydraulic heads contributing to the objective function. The 12 observation points of group 2 that were directly located between the canal and the pumping well group were given a weighting factor of 10, while for the remaining observation points of group 1 it was set as 1.

Transient flow model
The transient flow model was set up to simulate the temporal variation in the groundwater flow system caused by the transient condition in surface water heads, pumping rates and groundwater recharge, and by changes in riverbed hydraulic conductivity induced by the river reconstruction (decolmation). The transient simulation started mid-June 2010 and covered a period of 4.5 years until the end of December 2014. The hydraulic heads of the steady-state model were used as initial condition for the transient simulation. The model boundaries were consistent in type with the steady-state model, but now applying daily average values obtained from the measured time series. The specific yield was set as 6% for glacial tills, 8% for silt and 23% for sand (Domenico and Schwartz 1998) and specific storage was taken as 1 × 10 −4 m −1 (Freeze and Cherry 1979). Limited by the computation time, an automated calibration of the hydraulic conductivities was not feasible and therefore a manual calibration by trial and error was performed.
In the model, the river reconstruction was implemented in a way that all the grid elements located in the river area with a top elevation higher than 24.5 m were assumed to be influenced by the excavation. Since it was not possible to remove or replace inner grid cells during a transient simulation, the change in hydraulic resistance was instead implemented in an indirect way, i.e. all elements affected by the excavation kept their shape but obtained an increased hydraulic conductivity. Its value was fitted manually to reproduce the increasing water head after the start of the canal reconstruction and the finally adopted values were K h = 1.16 × 10 −4 and K x,y = 2.31 × 10 −5 m/s.

Transient heat transport model
The transient heat transport simulation was carried out from mid-June 2010 to the end of December 2014 corresponding to the transient flow model. Dirichlet boundary conditions were assigned for the surface water bodies (according to measured river temperatures) and land surfaces (according to measured air temperatures), while the remaining ones were defined as no-flux boundaries. All temperatures were applied as time series of daily-averaged values. The initial temperature in the model was set to 11.5°C, which is the temporal average of all observed temperatures.
The thermal parameters, heat conductivity, volumetric heat capacity and thermal dispersivity (longitudinal and transversal) of the existing geological units were adapted during manual calibration. The parameter range used varied between 1.6 and 6.8 W m −1 /K for thermal conductivity, between 1.8 and 3.2 MJ m −3 /K for the volumetric heat capacity, and between 0.1 and 1 m for the longitudinal thermal dispersivity following established literature values (VDI 4640/1 2010;Abu-Hamdeh 2003;Stonestrom and Blasch 2003). Gelhar et al. (1992) summarized the studies existing at that time and have shown that longitudinal dispersivity is bigger than transversal dispersivity, in the horizontal as well as vertical direction. In this model, the transversal dispersivity was set to 10% of the longitudinal value, a common assumption and experimentally supported by the study of Huang et al. (2003). To illustrate the uncertainties in simulated temperatures resulting from the given range in thermal parameters, two extreme scenarios were adopted representing good ability (minimum c and maximum λ) and limited ability (minimum λ and maximum c) of heat transport through the saturated porous media.
The accuracy of temperature simulations was quantified by the deviation in timing of the maximum summer temperature or the minimum winter temperature (ΔΦ) by the availability and difference in temperature amplitude (ΔA) between simulated and observed temperatures. For these analyses, the measured groundwater temperatures of group 3 were used as they are assumed to be representative of the groundwater temperature and directly linked to the position of the filter screen of each observation well. Only for W05 and W06 no sampling data were available and instead logger temperatures were used for calibration.

Field observation of water heads and temperatures
The study site is located in a lowland catchment, where several lakes of high storage capacity are connected along the surface waters. During the sampling period, the seasonal variation of canal stage fluctuation was rather small (29.17-29.86 m a.s.l.). Along the groundwater flow path, the water heads dropped to below 25 m a.s.l. when reaching the pumping wells. Showing a similarity in both the distance to the river and the screen elevation, three near-bank wells, W07 LOW , W08, and W09 LOW , were selected and the recorded water head is shown in Fig. 3a. The resulting hydraulic gradient remained relatively stable due to low water level variations in the canal and a stable pumping rate at the waterworks. Onwards from the start of reconstruction, an increase of around 0.5 m in water heads could be observed in W07 LOW , W08, and W09 LOW, while no increase was occurring in the surface-water levels and it did not correlate with the observed daily precipitation (Fig. 3a).
As an annually periodic signal, the temperature time series take the form of a sine function. This original signal gets smoothed, damped and delayed on its way from the infiltration at the canal bottom and banks to the observation wells (Fig. 3b). The distance from wells W07 LOW, W08, W09 LOW to the canal is comparable, as well as the elevation of the data logger installed inside these wells (Table 1), but the temperature patterns in those wells feature significant differences in both amplitude and peak time. The average temperature amplitude in W07 LOW was~7.5 K, while the average amplitudes at W08 and W09 LOW were much smaller (2.9 and 1.0 K, respectively). The peak time for W07 LOW was postponed by about one month, but for W08 and W09 LOW , the signals were shifted by about half a year.

Steady-state flow modelling
The scatter plot of water heads between observation and simulation after model calibration is shown in Fig. 4. The evaluation of the goodness of fit for the steady-state model showed good agreement between observed and simulated values (groups 1 and 2), with a mean absolute error (MAE) of 0.71 m and a root mean square error (RMSE) of 1.04 m (Fig. 4a). Observed water heads lower than 25 m a.s.l. are mainly from the observation points of the pumping wells (green points in Fig. 4a). A systematic deviation could be observed in this region, which is to be expected when applying a steady-state model while pumping was varying on a dayby-day basis. Points with high water head are mainly located in the southeast and southwest of the model, where a larger degree of uncertainty exists in the geological conditions, and  Table 1 thus also in the simulated heads, due to sparser drilling information. Generally, the geological condition at the study site is known to be discontinuous and heterogeneous, because of the unconsolidated composition of the Quaternary deposits that have likely been modified by erosion or formed naturally by depositional fluvial system dynamics.
For the observation points between the ND and pumping well gallery the deviation between observed and simulated water heads (group 2 only) was even smaller, with a MAE of 0.18 m and RMSE of 0.28 m (Fig. 4b). Almost all of the points are located within the interval of ±0.5 m, except for W03 HIGH (absolute error of 0.7 m, highlighted by red circle in Fig. 4b). Generally, the model did show a good overall fit and is assumed to represent the complex hydraulic conditions at the study site moderately well. The hydraulic conductivity for sand was largest at the bottom aquifer, with a value of 3.5 × 10 −4 m/s for the horizontal and the K h /K v of 8.3, respectively. The hydraulic conductivity of the remaining aquifers ranged from 2.9 × 10 −5 to 1.0 × 10 −4 m/s in the horizontal direction and from 5.1 × 10 −6 to 4.2 × 10 −5 m/s in the vertical direction. The K h of the aquitards (glacial tills and silt) was much lower, varying between 2.0 × 10 −7 and 2.0 × 10 −6 m/s. Details of the calibrated hydraulic conductivities can be found in Table 3. The large anisotropy ratio of SI-GT was set in order to fit the hydraulic heads in the southern part of the model domain which were otherwise underestimated. The same result could have been achieved by a lower horizontal hydraulic conductivity but that would also influence the hydraulic conductivities in the near-shore area. Therefore, the large anisotropy ratios of SI-GT mimic the impermeable second aquifer bottom without affecting the good calibration results of the first aquifer.

Simulations of transient flow over 4 years
Transient recalibration of the hydraulic conductivities yielded minor modifications of only a few sandy aquifer and glacial till units, mainly for the bottom aquifer (K transient > K steady state ) as well as for the SI-GT and WI-GT lenses (K transient < K steady state ; Table 3). The simulation results of the recalibrated model were grouped into 'before' and 'after' the start of the reconstruction period and were evaluated by both Pearson's correlation coefficient (r) and RMSE (Table 4). The resulting time series are exemplified in Fig. 5 (and more in Fig. S3 of the ESM).
The RMSE before reconstruction ranged between 0.09 and 0.33 m, with an average value of 0.24 m. Before the start of the reconstruction, three observation points (W03 LOW , W07 HIGH and W07 LOW ) were characterized by a moderate correlation level (<0.4). For W03 LOW (Fig. S3f of the ESM), the mismatch is overemphasized by the uneven distribution of the sampling frequency. The observed dynamic changes before Oct 2012 were actually in good accordance with the simulation results; however, it takes much smaller weight on evaluating the simulating accuracy, due to the sparsity of the data points. For W07 HIGH (Fig. 5b) as well as W07 LOW (Fig. S3a of the ESM), the mismatch resulted from a period with an observed decreasing trend in October 2012 before an increasing trend started in April 2013. This deviation was also observed at W05, W06 (FS3c,d of the ESM), and, to a lesser    5d); that is, at almost all observation points close to the canal except W08. This could indicate that part of the main driving force influencing the water heads in that winter period might be missing. The winter 2012/2013 was characterized by a long period in time with minimum daily temperatures below 0°C, and average daily temperatures of only 2.4°C. Potentially, freezing processes at the study site which were not represented in the numerical model might have caused this deviation between nominally observed and simulated water levels.
Besides, for all simulations the hydraulic conductivity (also fluid viscosity and density) was constant in time and did not depend on temperature. This simplification may contribute to the general deviation of the simulation results compared to the observations. However, the temperature dependent variations in fluid viscosity and density for the given temperature range were rather small and thus the uncertainties caused by temperature dependent variations in hydraulic conductivity were much smaller than the general uncertainty in heterogeneous hydraulic conductivity itself. Generally, the performance of the model including the reconstruction improved in respect to the scenario without reconstruction as expected. The increasing groundwater head after the start of reconstruction was captured at all observation points by the simulation including the hydraulic effect of reconstruction ( Fig. 5 and Fig. S3 of the ESM). During the riverbed reconstruction, low-permeability riverbed sediments (including parts of the glacial tills) were dredged out, while the Sand zone of WI-GT layers 1.0 × 10 −4 10.0 1.0 × 10 −4 10.0 Sand layers between WI-GT and SII-GT 5.6 × 10 −5 2.8 8.1 × 10 −5 2.3 Sand zone of SII-GT layers 8.4 × 10 −5 5.9 8.4 × 10 −5 5.9 Sand layers between SII-GT and SI-GT 2.9 × 10 −5 5.6 1.9 × 10 −4 5.3 Sand zone of SI-GT layers 3.5 × 10 −4 8.3 1.2 × 10 −4 2.0 Sapropel of Weisser See 8.9 × 10 −6 2.9 8.9 × 10 −6 2.9 Sapropel of Jungfernsee  bank reconstruction led to the removal of lateral hydraulic resistances. This increased the hydraulic connection between the canal and aquifer. To maintain the consistent pumping rates after the reconstruction, a lower hydraulic gradient between the canal and pumping wells was sufficient due to the increased hydraulic connection. Therefore, the water heads in the pumping wells increased as well as in the area between the canal and the waterworks. Nevertheless, the RMSE on average increased for the period of reconstruction. This reflects that the system change is less well represented than the stable period previously, which could be attributed to fact that the real local changes in the riverbed are not known and that this period was a forward prediction and not subject to recalibration.
In the steady-state model, W03 HIGH has shown the largest error among all the wells in group 2 and it continues this tendency in transient state (Fig. S3e of the ESM). After reconstruction, although the rising trend could be reflected in the scenario of including the reconstruction, the absolute error was smaller in the scenario without reconstruction. The deviation in W03 HIGH could be caused by the existence of lowconductivity material along the groundwater flow path from the canal not being captured in the geological model. During and after reconstruction, W07 HIGH and W07 LOW were characterized by a moderate correlation level. For W07 HIGH the short-term dynamics in the observed water level including a substantial drop of about 0.7 m on 12 January 2014, an increasing trend of around 0.2 m before a slow recovery of about 0.5 m until the end of February; after that, step changes were not captured by the model simulations (Fig. 5b). During this period the southern canal bank adjacent to W07 was renewed (see Fig. S1c of the ESM), i.e. a new geotextile was placed and covered with armour stones, potentially affecting the observed water levels at the nearby observation well W07 (distance between W07 and channel is only 13.2 m paired with a good hydraulic connection at this spot). From July 2014 onwards, the observed water level at W07 HIGH continuously decreased until the end of the observation period (Fig. 5b). In December 2014 the water level at W07 HIGH was close to the simulated water level of the model scenario without canal reconstruction, whereas the water level of the model scenario including canal reconstruction was about 0.6 m higher. This indicates that the better hydraulic connection between the canal and groundwater caused by the artificial decolmation of the canal bed at this canal section lasted for about 6 months only. The hydraulic connection there continuously decreased to its initial level before river reconstruction, probably by colmation. Comparable effects were observed for W07 LOW (Fig. S3a of the ESM). These results agree with Schafer (2006) who confirmed that substantial riverbed colmation was observed at a bank filtration site of the Ohio River several months after commissioning. However, at this study site it was observed just for a local spot.
The increasing riverbed hydraulic conductivity corresponds to the rising infiltration rate. The ratio of surfacewater recharge to groundwater (infiltration rate) at the ND in both the scenarios with and without reconstruction is shown in Fig. 6a. It is evaluated as a net volume exchange from surface water to the groundwater by bank filtration during each day. The infiltration rate increases promptly together with reconstruction, which started by removing clogged riverbed and bank fortifications, and established a higher water level soon thereafter. On average, the bank filtration volume rate from the river increased by 23%, i.e. 521 m 3 /day, compared to the scenario without reconstruction (Fig. 6a). This increase is equivalent to around 9% of the total pumping rate, and will have led to some spatial rearrangement of well catchments and an increase in the already high share of bank filtration water abstracted by the waterworks.
Influenced by the increasing infiltration rate, the travelling time distribution must also change. The results of forward particle tracking prove that the travel time distribution is shifted towards shorter travel times (shorter than 550 days); i.e. subsurface residence time of water in the aquifer is reduced (Fig. 6). Thus, flow paths with travel times between 100 and 550 days occurred more often in the model scenario including reconstruction compared to the scenario without reconstruction-for example, 55% of flow paths had travel times shorter than 300 days, whereas without reconstruction this would have been only 40% of the flow paths. Also an according decrease of long residence time (longer than 550 days) resulted, which has to be expected because pumping as the general driver was continuing on the same level as before. Steady-state backward particle tracking from the pumping wells indicates the general shape of the capture zones and corresponding subsurface travel times/residence times. Based on the boundary conditions of 30 June 2014, flowlines were extracted for the model scenario without and with reconstruction (Fig. 7). For the scenario without reconstruction (Fig.  7a), the capture zone of the north-eastern pumping wells extends toward the adjacent canal with dominant subsurface travel times ranging between 50 and 150 days; and some flowlines could be tracked back to start at the northern lake "Krampnitzsee" with subsurface travel times up to several years. North-western pumping wells mainly received their water from the lake 'Weisser See'. Along these flow paths, the travel times ranged between 150 and 600 days. Only the southern pumping wells turned out to cover the southwestern, larger part of the catchment.
The reconstruction work caused some major changes in subsurface flow regime. The northern capture zone shrinks towards the canal, i.e. no flow lines were traced back anymore to the northern lake "Krampnitzsee" [Lake Krampnitz] (Fig.  7b). The main reason for this substantial change could be attributed to removal of impermeable glacial tills below the canal bottom, which had so far impeded the vertical infiltration from the ND, and had made the nearby lake of same water level to deliver the water instead. After the reconstruction, part of the glacial tills has been removed or replaced by more conductive material and the direct infiltration from the canal became the origin of flow paths. In contrast, no main differences in subsurface capture zone and corresponding travel times were observed in the southern catchment. Overall, the increased hydraulic connection between canal and aquifer caused a reduction in subsurface capture zone in the area north of the canal and shortened subsurface travel times from canal to pumping wells.

Simulation of heat transport
The thermal parameters of the transient heat transport model were assigned similarly to each sedimentary material (e.g. glacial tills, sand, sapropel) and changed individually in each geological unit by trial and error. Due to this moderate level of calibration the heat transport simulations still rather have the character of a forward modelling. The calibrated thermal conductivities were 1.7 W and 2.7 W m −1 /K for silt/glacial tills and for sand, respectively. The volumetric heat capacity slightly varied depending on geological composition between 2.5 MJ m −3 /K (silt/glacial tills) and 2.7 MJ m −3 /K (sand). The agreement between the simulated temperature and the observed ones varies in different wells (Fig. 8, and more in Fig. S4 of the ESM). Mean differences in time shift of the maximum and minimum temperatures (ΔΦ) and temperature amplitudes (ΔA) between simulated and observed temperatures for the time period after river reconstruction for all observation wells with substantial temperature variations are presented in Table 5. The model scenario including the reconstruction characteristically shows that after the start of canal reconstruction the temperature amplitude of the annual cycle increased and the arrival time of the peak became earlier (Fig.  8b-e, cyan lines) compared to the simulation scenario without reconstruction (Fig. 8b-e, dark blue lines). This general change in the temperature signal directly corresponds to an increased groundwater flow rate (i.e. increased advective heat transport), which is generated by the better hydraulic connection between the river and aquifer after the canal reconstruction. Furthermore, for W07 ( Fig. 8b and Fig. S4a of the ESM), the match of observations on different levels supports a reasonable vertical representation of the groundwater flow distribution. In contrast, the temperature dynamics in W05 and W06 ( Fig. 8a and Fig. S4c of the ESM) did not differ substantially between both model scenarios. This indicates that the groundwater flow field in the eastern catchment is not substantially affected by the canal reconstruction.
The observation wells W04 HIGH , W04 LOW and U100 were characterized by almost constant temperatures throughout the year, which is well represented by both simulation scenarios ( Fig. S4e-g of the ESM). The constant signal was not expected to be influenced by the canal reconstruction, as the average yearly temperature (equal to long-term average in air/river temperature) is independent of the groundwater flow velocity. Mean simulated temperatures were 1.17, 1.12 and 0.70 K lower than observed temperatures in W04 HIGH , W04 LOW and U100 respectively.
The greatest deviation between observed and simulated temperatures of the model scenario including canal reconstruction occurred in W09 HIGH (Fig. S4b of the ESM). The overestimated amplitude and underestimated arrival time of the temperature signal were caused by overestimated groundwater flow rates (advective heat transport). The model is not reproducing the slow flow regime in W09 HIGH , indicating that the hydraulic connection between the canal and the shallow aquifer is weaker than incorporated in the current geological model. Thus, the heat transport simulations here could indicate and rectify a local refinement of the geological model to be appropriate.
The temperature distribution is shown for a depth of around 10 m (±2 m) in Fig. 9a-d. The zone adjacent to the river was characterized by strong seasonal temperature oscillations, whereby the wells which are not located in the infiltration zone (U100, northern bank) or are far away from the river (W04, distance to the channel~126 m) were characterized by yearly constant temperatures. Driven by the strong hydraulic gradient from the canal towards the pumping wells, the seasonal temperature signal was transported into the saturated sediment via advection. Both maximum summer temperature and minimum winter temperatures infiltrated into the canal bed, and further penetrated into the aquifer, but are retarded to a certain degree depending on heat capacity and thermal conductivity of the sediment (i.e. the thermal front moves slower than the water particle or conservative solutes itself; the retardation factor of heat is about 2.15). Minimum temperatures in the canal occurred at the end of January (~1°C). For the scenario without reconstruction, after 6 months (summer, 30th June) the cold front rose to 6°C after traveling about 50 m from the canal, whereas the temperatures in the river already increased to about 20°C (Fig. 9a). Variations in the transport distance (penetration depth) along the shoreline of the channel correspond to local differences in hydraulic conductivity; i.e. at W07 the cold front already passed at the presented summer snapshot (Fig. 9a). Analogously, the winter temperature pattern (Fig. 9c) could be interpreted in the same way. The rest of the aquifer tended to constant temperatures of~11.4°C (Fig. 9a,c). The seasonal temperature amplitude was completely reduced after a transport distance larger than~100 m in the southern aquifer. In the northern aquifer, the penetration depth of the seasonal temperature signal (thermal front) was limited by heat conduction to roughly 1-2 m (Fig. 9a,c).
In comparison, the reconstruction scenario (Fig. 9b,d) shows a more continuous area of larger temperature variation range at the southern river bank induced by the annually cycling canal temperatures which represents a more direct hydraulic contact caused by the canal bed reconstruction. The snapshot of the vertical temperature distribution at 30th June 2014 (Fig. 9e,f) illustrates that in the reconstruction scenario, the warm water from the canal could penetrate much deeper into the aquifer, almost reaching the model bottom, with less damping during its propagation, and a similar situation would be observed in winter.
Generally, the periodic temperature signal travelled further before it became flat, both in the vertical and horizontal  Fig. 9 Spatial temperature distribution (color-coded according to legends) and water head isolines (white) of the aquifer below WI-GT for a characteristic summer (30th June 2014) snapshot (a-b) and winter (31st December 2014) snapshot (c-d). Also shown is fluid flux and spatial temperature distribution of a cross-section between U100 and W03 for the characteristic summer snap shot (e, f). The fluxes are indicated by arrows representing fluid flow directions. All temperature patterns are shown for the model scenario without (a, c, e) and for the model scenario including canal reconstruction (b, d, f) directions. It is the evident that the reconstruction has driven substantial changes to water flow and heat transport at the surface-water/groundwater interface.

Discussion
Insights from heat transport simulation into the geological model Featured as a naturally occurring periodic signal with easy and rapid measurement, the time series of temperature shows greater sensitivity to variations during transport, above potential environmental noise, compared to conventional or pulse tracers (Constantz 2008;Cox et al. 2007;Mahinthakumar et al. 2005).
Research by  illustrated that by taking heat transport into consideration, the model non-uniqueness could be better constrained and the accuracy of flux estimation across the interface between surface water and groundwater could be improved. Delsman et al. (2016) showed that conditioning the model to temperature measurements could worsen other model results such as electrical conductivity and salinity but also recognized that this could be mainly caused by errors in the model structure itself. However, few studies so far have used heat as a natural tracer for the numerical modelling of bank filtration. One of them is Henzler et al. (2016), who especially investigated its details with 3-D modelling of flow and heat transport in areas with geological heterogeneity. By employing a 3-D numerical model for transient water flow and heat transport incorporating the heterogeneous geological structure, it is only partly possible to match both the seasonal temperature patterns and water heads at the same time. However, substantial deviations can be used to trace back shortcomings of the model structure, and heat transport can help as an additional test to evaluate model performance and realizing model limitations. The simulation results at W09 showed the lowest RMSE and the highest Pearson's correlation coefficient (see section 'Steady-state flow modelling'), but also substantial deviation between observed and simulated temperatures (see section 'Simulations of transient flow over 4 years'). The poor match in the temperature pattern of W09 indicated an overestimation of water flux, which could not be compensated by reasonably adapting the thermal parameters (Fig. 9d); i.e. the hydraulic connections in the geological structural model, especially between the canal and the two aquifer layers leading to W09, are not representing the real situation well, and could be adapted based on this insight resulting from temperature observations and heat transport modelling.

Influence of river reconstruction on redox reaction
The quality of water gained via bank filtration strongly depends on the concentration of surface-water-born substances and pollutants, the subsurface capture zone, the mixing ratio between river water versus ambient groundwater, and the spatial and temporal distribution of subsurface travel times (Hiscock and Grischek 2002). Both the increased infiltration rate along the canal and shorter subsurface travel times may have implications for the occurrence, fate and distribution of biodegradable compounds, other DOC or even pathogens etc. in the adjacent aquifer. Generally, removal efficiency (purification capacity) of the underground passage will be diminished by shorter path lengths, higher groundwater flow velocities, and accompanying shorter travel times. Resulting concentrations of surface-water-borne substances and pollutants within the aquifer, as well as the mobility of these substances, might increase, especially for very quick flow paths. However, a sufficiently long setback distance could compensate this by supplying enough time for the adaption of biosystems and other slow processes (Schoenheinz and Grischek 2011). The results of this study indicate that a change in the heat transport by reconstruction can very well also be relevant for pollutant transport to the waterworks, especially for the canal bottom which acts as a sensitive interface between the groundwater and surface water. Generally, the periodic temperature signal, and finally the summer and winter temperature extrema in groundwater increased and penetrated deeper into the aquifer. The potential for attenuation and degradation (removal efficiency) of biodegradable compounds and other DOCs depends on the subsurface temperature, redox state, and the evolving microbial community; all of them are site-specific and strongly vary in time. Especially, temperature is recognized to control microbial processes in saturated sediments (Munz et al. 2019;Acuña et al. 2008;Thamdrup and Fleischer 1998). The more pronounced temperature variation within the aquifer may thus lead to changes in the redox zonation. The oxic zone may reach deeper into the aquifer by a reduced microbial activity in winter (colder temperature); vice versa, the penetration depth of the oxic zone may be reduced by the increased microbial activity in summer (higher temperatures). However, the aforementioned variation might be offset by the increasing flow velocity in summer, and intensified in winter. With short travelling distance (13~16 m), the change may be observed in the near-bank wells of W07, W08 and W09 but to different degrees, due to influence from geological heterogeneity.

Conclusions
With an intensive and long-term monitoring programme, it was possible to observe the seasonal behaviour of a bank filtration system, as well as its development, during and shortly after the reconstruction of the ND canal as the main supply of bank filtration water towards the pumping wells of the waterworks. Evidenced by the change of the groundwater heads in the aquifer close to the canal, the riverbed reconstruction by its artificial decolmation had substantially increased the hydraulic connection between the canal and the aquifer. The change started to take effect directly with the beginning of reconstruction, but intensified during the reconstruction period and had a clear effect after it was accomplished. Generally, this effect lasted for at least a year, until the end of this observation period, except for W07. Additionally, to these qualitative conclusions, a much more detailed and quantitative interpretation was enabled by employing a 3-D numerical model for transient flow. Initially, a steady-state model supplied a calibrated catchment-scale flow field and checked the feasibility and consistency of the geological model. Based on the set of calibrated parameters, a transient model was implemented that allowed for the development of the bank filtration groundwater system to be tracked over several years, including the reconstruction period. It showed to be a useful tool to quantify the effect of the reconstruction by directly comparing two model scenarios 'with' and 'without' reconstruction, which went far beyond just interpreting the measured time series. A key benefit is that this modelling approach eliminates the effects of changing hydrologic boundary conditions that otherwise would add a lot of uncertainty to the interpretation of the situation after reconstruction with inherently varying drivers hampering hard conclusions based merely on experimental data.
The increased hydraulic connection between the canal and the aquifer led to some spatial rearrangement of well catchments and an increase in the already high share of bank filtration water abstracted by the waterworks. The bank filtration rate along the canal has increased on average by 23%. Correspondingly, the increasing flow rate substantially changed the subsurface travel time distribution towards shorter travel times. Flow paths with travel times below 200 days increased by around 10% and those below 300 days by 15%. The potential of biodegradable pollutants and other DOCs to penetrate into the aquifer may thus be increased, especially in winter, due to the joint hydrological and thermal effects caused by the river reconstruction. Nevertheless, a significant deterioration of water quality due to the canal reconstruction was not observed so far in the waterworks but could well be the case at other locations.
Calibration performance limited the model accuracy on reflecting the flow field and thermal field. It could be attributed to both (1) the insufficient or incomplete interpretation of the geological conditions and (2) the calibration accuracy. Heterogeneity brought great difficulty and uncertainty in concluding on the geological condition, and further detailed investigation would be required, especially on the area where the modelling result deviated largely from the observation. Autocalibration may also reach a better parameter composition for the model, however, due to the large calculation load, it could be applied only for the steady-state model. Other simplifications and assumptions might also contribute to the inaccuracy in the model simulation, e.g. the fixed groundwater recharge, or not considering the influence from temperature on the viscosity, which would require another set of underlying equations and thus modelling environment. Pumping well settings in the model might also contribute to the possible error in calibration. As the current pumping rate distribution in a well is decided by the K f -value of different layers, a clogging of the well screen may modify the real situation to show a different local flow field. Logger-recorded water head and biweekly sampling frequency for temperature was sufficient for describing the temporal dynamic. The number of wells and their distribution limited the spatial resolution of the study, but this will be the rule rather than the exception at other sites, too. Multi-level sampling could at least improve the vertical sampling resolution of both flow and thermal field in a given sparse set of groundwater observation wells.
The investigations of bank filtration and particularly the numerical modelling allowed for a quantitative assessment of the hydraulic impacts of river reconstruction within a dynamic, seasonally and annually changing hydrological system. Furthermore, the modelling was able to show how far the heat patterns penetrate into the aquifers during bank filtration and leave a seasonal imprint, even stronger after the canal reconstruction and a resulting higher hydraulic connection. Last but not least, such quantitative information is also required as the basis to interpret the fate of pollutants, e.g. micropollutants (Munz et al. 2019), during bank filtration, including their sensitivity to temperature and travel time distributions.