**Anthropic and Meteorological Controls on the Origin and Quality of Water at a Bank Filtration Site in Canada**

**Janie Masse-Dufresne 1,\*, Paul Baudron 1,2, Florent Barbecot 3, Marc Patenaude 1, Coralie Pontoreau 1, Francis Proteau-Bédard 1, Matthieu Menou 4, Philippe Pasquier 1, Sabine Veuille <sup>1</sup> and Benoit Barbeau <sup>1</sup>**


Received: 2 October 2019; Accepted: 20 November 2019; Published: 28 November 2019

**Abstract:** At many bank filtration (BF) sites, mixing ratios between the contributing sources of water are typically regarded as values with no temporal variation, even though hydraulic conditions and pumping regimes can be transient. This study illustrates how anthropic and meteorological forcings influence the origin of the water of a BF system that interacts with two lakes (named A and B). The development of a time-varying binary mixing model based on electrical conductivity (EC) allowed the estimation of mixing ratios over a year. A sensitivity analysis quantified the importance of considering the temporal variability of the end-members for reliable results. The model revealed that the contribution from Lake A may vary from 0% to 100%. At the wells that were operated continuously at >1000 m3/day, the contribution from Lake A stabilized between 54% and 78%. On the other hand, intermittent and occasional pumping regimes caused the mixing ratios to be controlled by indirect anthropic and/or meteorological forcing. The flow conditions have implications for the quality of the bank filtrate, as highlighted via the spatiotemporal variability of total Fe and Mn concentrations. We therefore propose guidelines for rapid decision-making regarding the origin and quality of the pumped drinking water.

**Keywords:** anthropic forcing; meteorological forcing; lake bank filtration; mixing ratios; environmental tracer; time-varying mixing model; sensitivity analysis

#### **1. Introduction**

Bank filtration (BF) is known as a cost-effective treatment step to produce drinking water [1,2]. This natural or artificially induced process occurs as surface water infiltrates into the aquifer from the banks and/or bed of a lake or a river and is subsequently intercepted by a pumping well [3]. During subsurface passage, water is exposed to physical, chemical, and biological processes, which may attenuate contaminants initially present in the surface water but also release unwanted minerals [4,5]. BF systems have proven to be efficient for the removal of turbidity [6–8], pathogens [9–11], and organic compounds [12–15]. The efficiency of BF systems to attenuate contaminants is strongly controlled by travel times [13,16] and redox conditions [17–19], which in turn depend on numerous site-specific natural and engineered parameters. Natural parameters include the hydrological and hydrogeological conditions, surface and groundwater quality, and prevailing physico–chemical

conditions [20]. Engineered parameters refer to the number of wells, the distance between wells and surface water, well spacing, the well type, depth, radius, location, and screen length [21,22].

Most BF systems are in the vicinity of rivers, where the bank filtrate is a mixture of surface water and ambient groundwater [4,23]. Numerous studies have shown that the dilution of contaminants by high-quality groundwater can also help to attenuate contaminants, enhancing the efficiency of a BF system. For instance, Kvitsand et al. [24] reported that dilution with ambient groundwater was significant enough to lower concentrations of natural organic matter. Derx et al. [25] numerically studied the effects of flooding on virus removal by bank filtration. They reported that a rapid decrease in river water level can lead to a hydraulic gradient towards the river and a dilution of virus concentrations by regional groundwater. In addition, some BF systems are placed in hydrogeological contexts with low-quality groundwater but can still achieve high-quality raw water with adequate regulation of mixing ratios [26–29]. Hence, when assessing the performance of a BF system, estimating mixing ratios is crucial to: (1) correctly differentiate between dilution and removal mechanisms and (2) control the occurrence of groundwater-borne contaminants [30]. BF systems typically show spatial variability of mixing ratios at the pumping wells, since they are affected by the distance to the surface water body [22]. Another factor governing the mixing ratios is the drawdown at the pumping wells [27]. The latter is subject to spatial and temporal variations, since BF systems are rarely operated under steady-state hydraulic conditions (e.g., river stage) and/or pumping regimes. However, when calculating mixing ratios, authors rarely discuss the temporal variations and the factors controlling this variability, even though erroneous estimation of the mixing ratios can lead to misinterpretation of the performance of the BF system.

This study aims to provide a better understanding of the relationship between anthropic (i.e., pumping regimes) and meteorological (i.e., hydraulic gradients) effects on the origin of bank filtrate. To this end, we investigated the spatiotemporal variability of flow patterns and mixing ratios at a two-lake BF site, where two surface water types (Lake A and Lake B) contribute to seven pumping wells. A time-varying mixing model based on electrical conductivity (EC) was developed in order to quantify the contributions of Lake A and Lake B (i.e., two water sources and further referred to as end-members) over a one-year period. A sensitivity analysis was conducted in order to test the assumptions concerning the definition of the end-members.

#### **2. Site Description**

#### *2.1. Hydrogeological Context*

#### 2.1.1. Description of the Bank Filtration and Aquifer System

The studied BF system supplies drinking water to more than 18,000 people in a town near Montreal, Canada (Figure 1a). A total of eight pumping wells are located between two artificial lakes (Figure 1a,b), which were created by sand dredging activities. The exploitation stopped a few decades ago at Lake B, while Lake A is still in operation. As described by Ageos [31], the aquifer is a buried valley embedded in the Champlain Sea clays (Figure 1b,d). The aquifer is mainly composed of alluvial fine to medium sands. A small lens (≤3.45 m thick) of alluvial gravel (with a sandy matrix in places) lies between the Champlain Sea clays and the alluvial sands near pumping wells P4 and P5 (see Figure 1d). The aquifer is fully unconfined. Hydraulic conductivity was estimated as 2.7 <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>m</sup>/s [31].

**Figure 1.** Study site location maps (**a**–**c**) and schematic lithological cut along pumping wells (**d**) (adapted with permission from Ageos, 2010 [31]).

The maximum thickness of the aquifer is 26 m and the static water level is about 4 m below the ground surface. The sandy bank is 100 m to 120 m wide and approximately 500 m long. All the pumping wells are screened at the base of the aquifer over an 8 m long section, except for pumping well P5, which only has a 4 m long screen due the shallower depth of the aquifer at this location. The distance between Lake A and the well cluster is 70 m to 80 m, whereas a distance of 30–35 m separates the wells from Lake B. Finally, the wells are spaced 30–60 m from one another.

#### 2.1.2. Lake A and Lake B

Lake A (2.8 <sup>×</sup> 105 <sup>m</sup>2) is fed by a stream named S1, which discharges from the North with a mean annual rate of 0.32 m3/s [31]. It drains a small watershed (14.4 km2), where land use is mostly industrial and agricultural. A 1 km long channeled stream named S2, located at the southeastern bank, allows water to exit Lake A and flow towards Lake C. The flow direction between Lake A and Lake C can be temporally reversed (Figure 1b) when the surface water level of Lake C exceeds both the elevation of Lake A and a topographic threshold at 22.12 m.a.s.l. [31]. Under these hydraulic conditions, Lake A receives surface water inputs from Lake C. This process typically occurs during spring (from April to May) and more occasionally during autumn (from October to December) due to snowpack melting and/or abundant precipitations. Ageos [31] reported that surface water input into Lake A seems to control its geochemistry, as it features a Ca-HCO3 water type.

Lake B (7.6 <sup>×</sup> 104 m2) is a groundwater-fed lake without any inlet stream. An artificial outlet channel can drain Lake B water towards the town's stormwater collection system (when Lake B elevation is above approximately 21.8 m.a.s.l.). A NaCl water type is found in Lake B [31]. Pazouki et al. [32] stated that the salinity of Lake B originates from de-icing road salts that are applied during wintertime. This is supported by the fact that a regional and widely used road is located less than 100 m from the study site. Precipitations are approximately 1000 mm/year and contribute to the water mass balance of both lakes. Runoff is likely a negligible contribution to the water mass balances of the lakes, considering the nearly flat topography. The maximum observed depths at Lake A and Lake B are 20 m and 19 m,

respectively (at LA-P2 and LB-P2). Based on the lithological cross sections at the pumping wells and observation wells [31], it is believed that lake bottoms roughly correspond to the elevation of the marine clay sediments. In this geological context, no or only minor groundwater flow could occur beneath the lake bottom. The sediments at the bottom of the lakes were not sampled and no quantitative information concerning clogging is available. However, while sampling for surface water, relatively high turbidity (denoted by the color and the milky appearance of water) was observed at Lake A, which indicates that the sediments at the lake–aquifer interface are susceptible to clogging [33,34]. Sampling of the sediments would be needed to confirm this hypothesis.

#### *2.2. Hydraulics of the Two-Lake BF System*

A water table monitoring program was performed by Ageos [35] from 2012 to 2015. This study reported that, prior to the activation of the BF system in October 2012, surface water levels of Lake B were higher than in Lake A. Such conditions forced surface water to infiltrate and flow naturally through the sandy bank from Lake B to Lake A (Figure 2a). For instance, during summer 2012, the water level difference was about 0.1 m, which created a natural hydraulic gradient of approximately 0.001 between the lakes. Based on Darcy's law, the mean residence time of the water in the bank was approximately one month.

**Figure 2.** Schematic representation of the flow patterns and directions at the study site when (**a**) elevation of Lake B > elevation of Lake A, (**b**) elevation of Lake A > elevation of Lake B, and (**c**) the pumping wells are in operation. Black, blue, and red arrows refer to regional groundwater and water originating from Lake A and Lake B, respectively. Theoretical elevation difference between Lake A and Lake B in (**d**). Numbers 1 to 3 correspond to different hydraulic conditions, namely high, moderate, and low hydraulic gradients between Lake A and Lake B.

The above-mentioned water table monitoring program also demonstrated that the water level in Lake A was significantly higher than in Lake B during springtime from 2012 to 2016 (i.e., up to 1 m water level difference). This is due to the intermittent hydraulic connection between Lake C and Lake A and supporting surface water inputs into Lake A (see Section 2.1). Under such hydraulic conditions, the direction of groundwater flow into the bank is reversed, i.e., from Lake A to Lake B (Figure 2b).

Since the implementation of the BF system (on 3 October 2012), the relative surface water elevations of Lake A and Lake B have not been the only controlling factors on the direction and intensity of groundwater flow through the sandy bank. Drawdown of the water table in the vicinity of the active pumping wells induces an artificial hydraulic gradient, forcing surface water from both lakes (A and B) to infiltrate into the sandy bank and travel toward the pumping wells (Figure 2c). A schematic representation of the theoretical elevation difference between Lake A and Lake B is shown in Figure 2d. When analyzing the data from the monitoring program conducted by Ageos [35], we depicted three typical hydraulic conditions recurring each year. First, a high hydraulic gradient between Lake A and Lake B develops in response to the hydraulic connection between Lake A and Lake C (as explained above). Second, in summertime, the hydraulic connection between Lake A and Lake C stops and water demand increases. This leads to a moderate hydraulic gradient between the lakes. Finally, in wintertime, a low hydraulic gradient is expected, as surface water inputs into Lake A are very limited and municipal water demands are reduced. In sum, the lake dynamics and the pumping regimes both influence the relative surface water elevations of Lake A and Lake B and allow for a gradual transition from high (during springtime) to low (during wintertime) hydraulic gradient between the lakes.

#### **3. Materials and Methods**

#### *3.1. Surface and Groundwater Sampling*

Monitoring of surface water and groundwater was conducted on a monthly basis and included measurements of physico–chemical parameters and water sampling for geochemical analyses. Surface water sampling was performed near the shore (see location of lake sampling points LA-S and LB-S in Figure 1c). Additional sampling campaigns were conducted at Lake A (on 15 February 2017 at LA-P1 and LA-P2) and Lake B (on 9 September 2016 at LB-P1 and LB-P2 and on 3 March 2017 at LB-P3 and LB-P4) to assess for vertical heterogeneity. Physico–chemical parameters were measured along vertical profiles at 1 to 2 m intervals and water was sampled at multiple depths (e.g., 3 m, 7 m, and 12 m) with a submersible pump. Groundwater sampling was conducted at the pumping wells via a bypass faucet, as submersible pumps permanently regulate flow rate at each well. Water sampling was conducted at least 30 min after pumping started, allowing the stagnant water to be purged. In the case of observation wells, a submersible pump (WSP-12V-5 Tornado, Proactive Environmental Products, Bradenton, FL, USA) with a 30 m long polyvinyl chloride (PVC) tube was used and sampling was conducted after purging at least three well volumes and stabilizing the physico–chemical parameters.

Measurements of temperature, pH, electrical conductivity (EC), and redox potential (Eh) were performed with a multiparameter probe (YSI Pro Plus 6051030 and Pro Series pH/ORP/ISE and Conductivity Field Cable 6051030-1, YSI Incorporated, Yellow Springs, OH, USA) installed in an airtight cell connected to the pump. Samples for major ions and alkalinity were collected in 50 mL low-density polyethylene (LDPE) containers and were filtered through a 0.45 μm hydrophilic polyvinylidene fluoride (PVDF) membrane (Millex-HV, Millipore, Burlington, MA, USA) prior to analysis. Water samples were transported and stored at 4 ◦C. The same sampling and transport procedures were applied for total and dissolved metals analysis (Fe and Mn). Following on-site filtration, acidification with HNO3 (in order to lower pH < 2) was performed in the laboratory within a 24 h delay.

#### *3.2. Analytical Techniques*

Major ion quantification was performed via either atomic absorption (Aanalyst 200 Atomic Absorption Spectrometer, Perkin Elmer, Waltham, MA, USA) or ion chromatography (ICS 5000 AS-DP DIONEX Thermo Fisher Scientific, Saint-Laurent, QC, Canada) for all surface water samples and groundwater samples collected at observation wells, depending on the availability of the equipment. Total Fe and Mn concentrations were measured via atomic absorption for all surface water and observation wells samples. Inductively coupled plasma mass spectrometry was used for the quantification of major ions and total and dissolved Fe and Mn concentrations for the water samples collected at the pumping wells. The limit of detection (LOD) was 0.2 mg/L for all major ions and 0.01 mg/L or 0.05 mg/L for total and dissolved Fe and Mn, depending on the quantification method. For subsequent calculations and interpretations, all results ≤ LOD will be considered equal to LOD/2. Duplicates were analyzed to confirm the repeatability of the quantification methods. Bicarbonate concentrations were derived from alkalinity, which was measured manually in the laboratory according to the Gran method [36]. On samples with measured alkalinity (*n* = 98), the ionic balance errors were all below 10%. The mean and median ionic balance errors were 1% and the standard deviation was 3%.

#### *3.3. Estimating Mixing Ratios*

The mixing between two end-members can be quantified via a binary mixing model which can be described by the following equations:

$$f\_{\mathcal{A}} + f\_{\mathcal{B}} = 1,\tag{1}$$

$$X\_{\!\!\!A}f\_{\!\!A} + X\_{\!\!\!B}f\_{\!\!B} = \ \!X\_{\!\!W\_{\!\!\!A}} \tag{2}$$

where *f* represents the fraction of the different sources and *X* the concentration (or value) of the tracer. *A* and *B* correspond to the two water sources, whereas *W* represents the water sampled at the well.

Tracer-based approaches can be used to estimate mixing ratios and travel times, as long as the tracer presents conservative or predictable behavior [37,38]. Various natural tracers, such as chloride (Cl−), electrical conductivity (EC), and stable isotopes of waters (δ18O-δ2H), have been applied in numerous BF or alluvial aquifer contexts [39–44]. In this paper, we used EC values as a quantitative mass balance tracer for the application of the mixing model, with the assumption that it behaves conservatively. Violation of this assumption was unlikely at the study site, considering that the aquifer matrix is alluvial sands (mainly siliceous with no calcite). Good correlation (R<sup>2</sup> = 0.95) between EC values and Cl− (a conservative tracer) was also observed. The advantages of using EC instead of Cl− are that measurements can be done at a low-cost, as well as remotely and continuously.

#### **4. Results and Discussion**

#### *4.1. Highly Transient Pumping Schemes*

In this section, we (1) identify typical pumping schemes and (2) depict the seasonal variability of the total pumped volume.

Figure 3a–d shows the pumping rates for P1, P3, P5, and P6 over a typical one-week period (from 16 January 2017 to 23 January 2017). P1 was mainly active during daytime for 1–12 h (Figure 3a). A similar pumping scheme was applied to P2, P7, and P8 during summertime (data not shown). P3 and P6 were operated at rates ranging from 1000 m3/day to 3000 m3/day. Both were typically active on a daily basis, although P3 was turned off during night time (for less than 6 h) as water demand diminished (Figure 3b,d). P5 and P4 were typically activated on a monthly basis for monitoring and sampling procedures (Figure 3c). Three general pumping schemes emerged from this analysis of pumping rates and made it possible to distinguish three groups: (1) wells operated at nearly continuous rates (P3 and P6); (2) wells operated intermittently (P1, P2, P7 and P8); and (3) wells operated occasionally (P4 and P5).

**Figure 3.** Pumping rates for wells (**a**) P1, (**b**) P3, (**c**) P5, and (**d**) P6 during a typical one-week period (from 16 January 2017 to 23 January 2016). Monitoring and water sampling were conducted on 17 January 2017 at all the pumping wells.

Figure 4 illustrates the monthly mean total pumping rate for all wells from March 2016 to March 2017. The mean pumping rate was about 4400 m3/day, excluding summertime (May 2016 to September 2016), during which it was approximately 7000 m3/day. Throughout most of the year, with the exception of summer months, 71% to 83% of the total daily pumped volume was provided by the continuously pumping wells. The intermittently pumping wells provided 16% to 29% of the pumped volume. The remaining volume (<1%) was supplied by the occasionally pumping wells. In summertime, pumping rates increased at all wells, except for P5. Continuously pumping wells were operated at mean rates of approximately 2000 m3/day, representing from 52% to 63% of the total pumping rate. The intermittently pumping wells together supported 36% to 46% of the total pumped rate and the occasionally pumping wells supplied together the remaining 3%.

Over the study period, the total pumped volume fluctuated daily and seasonally to accommodate the municipal water demand. Indeed, higher pumping rates prevailed during (1) mornings and evenings, (2) weekends, and (3) summertime. This well field is typically operated with a hierarchical system, giving priority to the continuously pumping wells. If the water demand increases, intermittently pumping wells are subsequently activated. Lastly, the occasionally pumping wells can be solicited. This implies that anywhere from one to eight pumping wells were solicited to fulfill the water demand and accommodate for the daily and seasonal water demand fluctuations.

**Figure 4.** Monthly mean total pumping rate from March 2016 to March 2017. Above each bar are the proportions of the total pumped volume supplied by the continuously (in black) and intermittently (in grey) pumping wells. The occasionally pumping wells supply only <1–2% of the total pumped volume.

#### *4.2. Geochemistry as a Proxy of the Hydrosystem Dynamics*

The objective of this section was to examine the geochemistry of Lake A, Lake B, regional groundwater, and the bank filtrate in order to identify the contributing water sources to the pumping wells.

Box plots of the temperature, EC, pH, and Eh at Lake A (<1 m depth), Lake B (<1 m depth), the pumping wells, and the observation wells Z12, Z15, and Z16 are illustrated in Figure 5. Concerning Lake A and Lake B, note that the presented data correspond to measurements at the surface of the lakes (i.e., <1 m depth). Hence, the medians and the 25th and 75th quartiles values may not be representative of the entire water column. Observed temperatures at Lake A and Lake B ranged from 1.3 ◦C to 27.5 ◦C and from 3.9 ◦C to 27.5 ◦C, respectively. For the pumping wells, box plots are spatially sorted (P6; P1; P8; P2; P7; P3; P4; P5) from the northwest to the southeast ends of the well field (see location of the pumping wells in Figure 1c). Temperatures ranged from 3.4 ◦C to 16.2 ◦C, with minimum and maximum values being observed in occasionally and continuously pumping wells, respectively. EC values at the pumping wells ranged from 491 μS/cm (at P5) to 895 μS/cm (at P8), which is in between observed EC values in Lake A and Lake B. Note that the EC values for Lake A and Lake B in Figure 5 are associated with water sampled at <1 m depth. Higher EC values were measured in situ in Lake B at >12 m deep (further details in Section 4.3). Observed EC values at Z12 were similar to those in Lake B, whereas Z15 showed lower values, similar to Lake A. The highest EC values were observed at Z16, which is representative of regional groundwater. Measured pH values at the pumping wells tended to increase spatially from NW to SE (P6 to P5). Redox conditions also varied spatially and decreased from NW to SE. A H2S odor was noticed when sampling at P4, P5, Z15, and Z16, which is consistent with Eh measurements that indicate more reduced conditions.

**Figure 5.** Boxplots of (**a**) temperature, (**b**) electrical conductivity (EC), (**c**) pH, and (**d**) redox potential (Eh) at the lakes, pumping wells (PW), and observation wells (OW). Blue and red boxes are associated with Lake A (<1 m depth) and Lake B (<1 m depth), whereas dark, medium, and light grey boxes correspond to continuously, intermittently, and occasionally pumping wells, respectively. Numeric values above each box correspond to the median.

Figure 6 shows the spatial variability of total Fe and Mn concentrations at Lake A (<1 m depth), Lake B (<1 m depth), the pumping wells, and the observation wells Z12, Z15, and Z16. Concentrations in total Fe ranged from <0.01 mg/L to 1.28 mg/L at the pumping wells, with median concentrations increasing from NW to SE. Median total Fe concentrations at P4 and P5 were high relative to Canada's aesthetic objective for total Fe in drinking water (i.e., 0.3 mg/L) [45]. Analyses also reported high concentrations (from 0.05 mg/L to 2.12 mg/L) at Z15 (near P5). The highest total Fe concentrations were observed at Z16. Total Mn concentrations ranged from 0.1 mg/L to 1.3 mg/L at the pumping wells, which exceeded the aesthetic objective for total Mn in drinking water in Canada (i.e., 0.02 mg/L) [46]. The highest concentrations were measured at the intermittently pumping wells. Total Mn concentrations at the surface of Lake A and Lake B were relatively low (i.e., typically ≤0.03 mg/L). However, it is important to note that 1.06 mg/L was observed at 6 m depth in Lake B (see red circle in Figure 6b). This result highlights that total Mn concentrations may be more important at greater depths in Lake B. Release of Fe and Mn to the water column in Lake A may potentially occur from sand dredging activities, as this lake is still actively mined for sand. However, no data are available to discuss the evolution of Fe and Mn concentrations in relation to these anthropic activities. Dissolved Fe concentrations were all <0.05 mg/L, while dissolved Mn concentrations were similar to the total ones.

**Figure 6.** Boxplots of (**a**) total Fe, and (**b**) total Mn concentrations at the lakes, pumping wells, and piezometers. Blue and red boxes are associated with Lake A (<1 m depth) and Lake B (<1 m depth), whereas dark, medium, and light grey boxes correspond to continuously, intermittently, and occasionally pumping wells, respectively. Numeric values above each box correspond to the median. Red circle represents the maximal observed total Mn in Lake B (at 6 m depth).

Figure 7 shows the relationship between (Ca2<sup>+</sup> + Mg2<sup>+</sup>)/Na<sup>+</sup> and cationic content (i.e., sum of major cations) for Lake A, Lake B, pumping wells, and regional groundwater. Potassium (K+) was excluded from these calculations since only a few samples were analysed for K<sup>+</sup> and concentrations in K<sup>+</sup> only represent a small fraction of the total cation content (i.e., approximately 1.5%). Lake A and Lake B samples are plotted in opposing regions of the graph, Lake A having a high (Ca2<sup>+</sup> + Mg2+)/Na<sup>+</sup> ratio and low cationic content and Lake B having a low (Ca2<sup>+</sup> + Mg2+)/Na<sup>+</sup> ratio and high cationic content. Concerning the samples from the pumping wells, they are mostly plotted in the area extending from the Lake A to Lake B regions. Occasionally pumping wells had a geochemical signature similar to Lake A, whereas continuously and intermittently pumping wells spread between both lake signatures. Regional groundwater samples were sampled from one observation well, namely Z16, located on the NE side of Lake B (see Figure 1c). These samples were characterized by the lowest (Ca2<sup>+</sup> + Mg2+)/Na<sup>+</sup> ratios and highest cationic content. It is believed that direct contribution to the pumping wells from regional groundwater is not likely at this site, due to the hydrogeological context (see Section 2.1). Hence, we hypothesize that the spreading of pumping well samples relative to the Lake A–Lake B mixing line is potentially due to an indirect contribution from regional groundwater, which discharged into Lake B. Only three wells (i.e., P2, P7, and P8) were affected from November 2016 to February 2017. During this period, the three wells together supplied <10% of the total pumped volume. Based on these observations, we propose that the mixing between Lake A and Lake B is the dominant process governing the geochemical facies of the pumping wells.

#### *4.3. EC Time-Varying Mixing Model*

It was discussed in Section 4.2 that the geochemical facies at the pumping wells are controlled by mixing between Lake A and Lake B. Hence, we used the binary mixing model of Equations (1) and (2) to estimate the relative contributions of each lake to the pumping wells. In this section, we first present the temporal and vertical variability of EC at Lake A and Lake B in order to define the end-member values. Then, estimations of the mixing ratios are evaluated with respect to a reference scenario, and spatiotemporal evolution is discussed. We also provide a sensitivity analysis, which helps strengthen the conclusions of the model.

**Figure 7.** Comparison of the geochemical facies of Lake A, Lake B, pumping wells, and the regional groundwater (GW). The solid black line represents the mixing line between Lake A and Lake B mean values.

#### 4.3.1. Temporal and Vertical EC Variability at Lake A and Lake B

Temporal variability in EC values at the surface (<1 m depth) of Lake A and Lake B are illustrated in Figure 8a,b. At the surface of Lake A, minimal and maximal EC values were observed in springtime and wintertime, respectively. Low EC values were expected for springtime, since it corresponds to the period of hydraulic connection between Lake A and Lake C. During this period, surface water with low EC is discharged into Lake A from streams S1 and S2 with inverted flow direction (see Figure 1b). In Lake B, EC values at the surface (<1 m depth) were also found to be variable over time. During springtime and summertime, EC values were relatively constant. A significant increase in EC values was observed in autumn–winter. Figure 8b also depicts the EC time series at an observation well (namely Z12) which was located between pumping well P1 and Lake B. It is screened at the bottom of the aquifer over a 9.14 m long section. EC measurements at Z12 are thus representative of the mixing between multiple flow lines originating from various depths in Lake B. The mean EC value at Z12 was 848 μS/cm and values were typically higher than at the surface of Lake B. These results reveal that (1) the EC measurements at the surface of Lake B are not representative of the infiltrating water and (2) considering the vertical variability in EC in Lake B is important.

Figure 9 shows vertical EC profiles for Lake A and Lake B. EC was measured at depth in winter (on 15 February 2017) at Lake A and in summer and winter (on 9 September 2016 and 2 March 2017) at Lake B. For each campaign, at least two vertical profiles were conducted in order to assess the horizontal variability (see Figure 1c for location of the vertical profiles). Maximal EC differences (at the same depth) were 6 μS/cm and 25 μS/cm for Lake A and Lake B, respectively, suggesting no significant horizontal variability at both lakes.

**Figure 8.** Time series of electrical conductivity (EC) at (**a**) Lake A (<1 m depth), and (**b**) Lake B (<1 m depth) and Z12. The grey shaded area represents the timing for hydraulic connection between Lake A and Lake C. Error bars (±20 μS/cm) represent the maximal expected measurement error on EC. The yellow dashed line corresponds to the mean EC value at Z12.

**Figure 9.** Electrical conductivity (EC) measurements against depth at Lake A and Lake B in winter and summer. Solid lines represent depth-average value at the lakes, while the yellow dashed line is associated with mean EC values at observation well Z12. Blue and red shaded areas illustrate the variability in EC at the surface (<1 m depth) of Lake A and Lake B, respectively.

Concerning Lake A, no significant vertical variability in EC values was observed. This suggests that Lake A was vertically well mixed in wintertime. Given that Lake A receives surface water from a stream and that some industrial activity (i.e., sand dredging) takes place in the lake during the ice-free period (typically form early May to late October), it is likely that some currents in Lake A are stimulating mixing of the water column. Hence, we assumed that Lake A is fully mixed and does not develop any significant vertical EC stratification over a hydrological year. Temporal variability in EC at the surface of Lake A (at <1 m depth) is presumably representative of the evolution of the entire water body.

At Lake B, observed EC values increased with depth for all vertical profiles. Higher EC at greater depth could be induced by regional groundwater inputs into Lake B. In Canada, groundwater inputs are typically found at the bottom of lakes, due to thermal (and density) contrast [47]. Smaller vertical variability was observed in wintertime (in comparison to summertime). However, both summertime and wintertime depth-averaged values were similar (861 μS/cm and 884 μS/cm; a difference of 23 μS/cm being barely significant). It is important to note that the wintertime vertical profiles were conducted in a shallower zone of the lake and could explain the discrepancy between the depth-averaged values. Additionally, it is interesting to note that the depth-averaged EC values were similar to the Z12 mean EC value. This suggests that the depth-averaged EC value was adequate to depict the EC signal originating from Lake B.

#### 4.3.2. Reference Scenario

Based on the assumption that Lake A is well mixed, we considered that EC measurements at <1 m depth were representative of the Lake A end-member. Hence, the Lake A end-member is a time-varying EC signal. Temporal interpolation between discrete measurements was done with the cubic spline method (using the *spline* function in MATLAB). The result of this calculation is represented in Figure 8a by the solid blue line and gives the best available estimate of the Lake A end-member from 27 April 2016 to 9 February 2017. No temporal shifting was considered, since travel times were expected to be much smaller than the observed changes in EC in Lake A. Concerning Lake B, a constant value of 873 μS/cm (i.e., mean of the wintertime and summertime depth-average values) was considered to correctly represent this second end-member. There was no need to consider temporal variation for Lake B end-member as both depth-averaged values were found to be similar.

The results of the mixing model are shown in Figure 10. By considering the relative pumping rates and the estimated contributions from Lake A at each well, we calculated that 62% of the annual pumped volume originated from Lake A. The continuously pumping wells are characterized by 54% to 78% of water originating from Lake A, with the highest contributions from Lake A occurring from April to July, i.e., during the highest hydraulic gradient period. The lowest contributions from Lake A were observed from July to September, i.e., during the moderate hydraulic gradient period. This is likely related to an increase in the total pumped volume during summertime (see Section 4.1).

The intermittently pumping wells showed the widest distribution in mixing ratios, with contributions from Lake A ranging from 0% to 87%. Similar to the continuously pumping wells, the highest contributions from Lake A were observed during the high hydraulic gradient period, while its contribution decreased as the hydraulic forcing became less important. It was estimated that during the low hydraulic gradient period, the fraction of Lake B can reach up to 100% for the intermittently pumping wells. In fact, the mixing model yields such mixing ratios when the measured EC at the pumping wells is greater or equal to the Lake B end-member (i.e., 873 μS/cm). This condition was observed four times and EC measurements at the concerned pumping wells ranged from 879 μS/cm to 895 μS/cm. However, expectations were that Lake A and Lake B would always contribute to the pumping wells, since a radial depression cone normally develops in the vicinity of an active pumping well, forcing water to infiltrate from both sides of the sandy bank. Hence, we considered that a calculated contribution of 100% of Lake B depicts a limit of the developed mixing model. In reality, this result could be an indication that, during winter, water preferentially infiltrates from the bottommost zone (≥12 m) of Lake B (where EC >873 μS/cm), leading to higher EC values at the intermittently pumping wells. Controls on the development of such preferential flow paths are not within the scope of this paper and, thus, will not be further discussed. Future work concerning the spatiotemporal variability of the hydraulic conductivity is still needed to draw any conclusions on this topic. The combination of various environmental tracers would help to differentiate between contributions from the surface and bottommost zones of the lakes.

**Figure 10.** Estimated contribution from Lake A to the pumping wells according to the reference scenario. Lake A end-member is defined as a time-varying electrical conductivity (EC) signal, which was derived from the observed EC values at <1 m depth. Lake B end-member is a fixed EC value which corresponds to the mean of the wintertime and summertime depth-average value.

Contribution from Lake A at the occasionally pumping wells is typically >90%. However, in April and May, it was estimated that the former wells were receiving a relatively smaller contribution from Lake A (from 74% to 86%). This result possibly reflects that pore water with EC >500 μS/cm and originating from Lake A could have been stored during winter within the sandy bank and pumped only in April and May. Hence, mean residence time of the water in the sandy bank could reach months when the pumping wells are not active. Greater attenuation of the surface water temperature signal, at the occasionally pumping wells, also testifies to longer residence times. We thus highlighted the need for considering the variability of the residence time of water into the sediments when applying a time-variant mixing model.

In sum, contribution from Lake A is typically greater than Lake B throughout the year. However, the mixing ratios are temporally and spatially variable. Strong variability was found especially during the period of low hydraulic gradient at the intermittently and occasionally pumping wells. In such a context, the pumping regime seems to have a decisive impact on the mixing ratios and, ultimately, on water quality of the bank filtrate. Under high hydraulic gradient the mixing ratios tend to be more similar, regardless of the pumping regimes.

#### 4.3.3. Sensitivity Analysis

A sensitivity analysis was conducted in order to investigate (1) the representativity of EC end-members values and (2) the uncertainties related to the EC measurements. Mixing ratios were therefore recalculated according to various scenarios where Lake A and Lake B end-member values varied from 458 μS/cm (i.e., minimal observed value) to 576 μS/cm (i.e., maximal observed value) and from 824 μS/cm to 924 μS/cm (i.e., reference scenario ±50 μS/cm), respectively. Also, a variation of ±40 μS/cm was applied to all the EC measurements at the pumping wells. Differences between the results of the scenarios were typically <10%, except for the ones concerning the Lake A end-member. When considering fixed EC values for the Lake A end-member, the estimation of the mixing ratios diverged up to 30% compared to the reference scenario. This result helped to quantify the importance

of considering the temporal variability of the end-members to obtain reliable results when estimating mixing ratios. Despite the sensitivity of the model to Lake A end-member variability, general trends for mixing ratios were conserved for all the scenarios. This result was expected, as the mixing model was linear. Overall, the sensitivity analysis revealed that the relative estimations of mixing ratios were acceptable and that measurement errors were not likely to influence our conclusions.

The temporal resolution of the applied monitoring program did not allow for discussion of the short-term (i.e., hourly to daily) EC variability. Hence, it is not clear whether hourly variations in pumping rates could influence the observed EC values at the pumping wells.

#### *4.4. Dominant Controls on the Origin of the Bank Filtrate*

The time-variant binary mixing model highlights that the contribution of Lake A to the bank filtrate can vary from 0% to 100%. This section aims to understand the competing roles of anthropic and meteorological forcings on the origin of the bank filtrate. We define anthropic forcing as a process via which the origin of the bank filtrate at a given well is affected by its own pumping scheme and/or rate. Anthropic forcing can also occur indirectly, as drawdown of the water table in the vicinity of a given pumping well can influence the origin of water at less active adjacent pumping wells. Meteorological forcing is considered a natural process. Concerning our study site, the surface elevations of Lake A and Lake B showed seasonal variations, which are mainly controlled by meteorological conditions allowing or limiting surface water inputs into Lake A. Hence, we considered that the hydraulic gradient between Lake A and Lake B is meteorological forcing acting on the BF system.

Figure 11a,b shows EC against the one-month average pumping rate prior to the sampling date. Distinction between the pumping regimes (i.e., continuous, intermittent and occasional) is illustrated in a, while hydraulic gradients between Lake A and Lake B (i.e., high, moderate, low) are represented in b. Figure 11c is a schematic representation of the dominant forcing in relation to the different hydraulic contexts. First, in Figure 11a,b, we observed little variability in EC measurements if the pumping rate was >1000 m3/day. In fact, most samples associated with high pumping rates (Figure 11a) showed EC values ranging from 583 μS/cm to 689 μS/cm, despite the variability of the hydraulic gradient (Figure 11b). Hence, for high pumping rates (i.e., >1000 m3/day), it appears that anthropic forcing is dominant over the meteorological forcing (Figure 11c). However, two samples (see downward arrow in Figure 11c) showed lower EC while being operated at >1500 m3/day. These exceptions were observed exclusively when the hydraulic gradient between Lake A and Lake B was maximal (i.e., in May and June). Under such hydraulic conditions, the anthropic forcing cannot counteract the meteorological forcing, resulting in an increase in the contribution from Lake A during springtime. Second, higher EC values (>750 μS/cm) are associated with the intermittent pumping regime, while low EC values (<600 μS/cm) are mostly related to the occasionally pumping regime (Figure 11a). For these two pumping regimes, meteorological forcing was clearly dominant over the anthropic one (Figure 11c). In fact, high EC values were strictly observed during the low hydraulic gradient period (Figure 11b). In short, this revealed that when a pumping rate of approximately 1000 m3/day is applied continuously, the mixing ratios are less variable due to direct anthropic forcing. When wells are operated only intermittently or occasionally, indirect anthropic and/or meteorological forcings control the mixing between Lake A and Lake B waters.

**Figure 11.** Relationship between electrical conductivity (EC) and the one-month average pumping rate prior to the sampling date, according to (**a**) the pumping regime and (**b**) the hydraulic gradient between Lake A and Lake B. A schematic representation of the dominant forcing is illustrated in (**c**), where solid and dashed lines represent the range of observed values and the delimitation between regimes where meteorological and anthropic forcings are dominant, respectively.

#### *4.5. Implications for the Quality of the Bank Filtrate*

In Section 4.2, we highlighted that geochemical analyses showed spatial variability in both total Fe and Mn concentrations at the pumping wells. This section aims to discuss the relationship between total Fe and Mn concentrations and the origin of water.

High total Fe concentrations were found at the occasionally pumping wells (i.e., at P5 and, to a lesser extent, at P4) and were associated with the highest contributions from Lake A (see Figure 10) and more reduced conditions (see Figure 5d). In comparison to the more anthropized section of the BF system, the residence times of the infiltrating water in the vicinity of the occasionally pumping wells are likely to be longer, because meteorological forcing alone is controlling groundwater flows (see Section 4.4). Since relatively low temperatures were observed at P4 and P5, it is also likely that higher viscosity, resulting in lower hydraulic conductivity, was responsible for longer residence times of the bank filtrate in the vicinity of these wells [48–50]. The longer residence times are potentially responsible for the high total Fe concentrations at P4 and P5. Evolution of redox conditions (from oxic to anoxic) is typically observed at BF systems [51] and can result in the dissolution of iron and/or manganese along flow paths [52]. However, as dissolved Fe concentrations are very low (i.e., generally <LOD), total Fe is controlled predominantly by the particulate fraction. Hence, it is more likely that the high rate and/or occasional pumping are causing the mobilization and resuspension of particulate Fe at P4 and P5. In fact, when activated for monthly sampling and monitoring, P4 and P5 typically operate at 150 m3/h, while the other wells operate at lower rates (i.e., from 40 m3/h to 125 m3/h). Moreover, P5 was the only pumping well equipped with a 4 m long screened section (i.e., half of those of the other wells). The mean effective velocity of water entering P4 and P5 screens was from 2 to 4 times greater than at the other wells. The total Fe concentration at P4 and P5 could thus potentially be reduced by lowering the hourly mean pumping rates and operating on a daily basis. However, such engineered operational strategy would not help to lower the total Fe concentration to <0.2 mg/L (see Figure 12a).

The highest total Mn concentrations were concomitant with the highest fraction of Lake B water at the intermittently pumping wells. The presence of total Mn in the raw water could also be explained by the evolution of redox conditions along the flow path. Besides this, an elevated concentration in total Mn (1.06 mg/L) was measured at a 6 m depth in Lake B. As the latter was found to be geochemically stratified, relatively reduced conditions can develop in the epilimnion and promote the solubilization of Mn. Hence, it is also likely that Mn reaches the pumping wells by advective transport with water originating from the deeper zones of Lake B. Further investigation is needed to better understand the site-specific drivers of the Mn occurrence in the bank filtrate, since its mobility is controlled by numerous factors, such as travel times, temperature, pH, microbial activity, the extent of a clogging layer, and the degree of oxygen consumption [52]. Figure 12b illustrates the relationship between total Mn concentrations and the one-month average pumping rate prior to the sampling date. Total Mn concentrations decrease with higher pumping rates (for intermittently and continuously pumping wells). This suggests that pumping rate can be used as an operational tool to control the total Mn concentration in the pumped water.

In sum, high total Fe and Mn concentrations in the pumped water are governed by two distinct processes. Total Fe seems to originate from particulate iron mobilization and resuspension when effective velocities of water entering the screens of the pumping wells are high, whereas high total Mn concentrations seem to be associated with an increase in the contribution from the bottom of Lake B. Total Fe and Mn concentrations could potentially be regulated by lowering the mean effective velocity of water entering the screens and adjusting mixing ratios (i.e., by operating at adequate pumping rates).

**Figure 12.** Relationship between (**a**) total Fe and (**b**) total Mn concentrations and the 1-month average pumping rate prior to the sampling date.

#### **5. Conclusions**

In this study, we demonstrated the controls of variable meteorological conditions and pumping schemes on the origin and quality of bank filtrate. Through a pumping rate analysis, the pumping schemes could be separated into three categories, namely the continuously, intermittently, and occasionally pumping wells. The continuously pumping wells (i.e., P3 and P6) supported 71% to 83% of the total pumping rate, except in summertime, when they contributed from 52% to 63%, of the total pumping rate as it increased from approximately 4000 m3/day to 7500 m3/day. An investigation of the geochemical facies of Lake A, Lake B, regional groundwater, and the bank filtrate revealed that the geochemistry of the pumped water is governed by mixing of Lake A and Lake B. Therefore, a two end-member mixing model was developed to estimate the contribution from both lakes to the pumping wells over a one-year period. To this end, EC measurements were used as a quantitative environmental tracer. A time-varying EC signal was considered for the Lake A end-member, whereas

a fixed EC value was used to depict the Lake B end-member. This simple mixing model revealed the following:


This study highlights how understanding the competition between anthropic and meteorological forcings can help to recommend guidelines for rapid decision-making regarding the quality of the pumped water. For instance, by identifying contexts for which the anthropic forcing is dominant, one can control the origin of the bank filtrate. Moreover, predicting periods under which the meteorological forcing is governing the flow patterns can help to adjust post-BF treatment in order to secure high quality of the distributed drinking water.

**Author Contributions:** Conceptualization, J.M.-D., P.B. and F.B.; Data curation, J.M.-D.; Formal analysis, J.M.-D.; Funding acquisition, P.B.; Investigation, J.M.-D., M.P., C.P., F.P.-B., M.M. and S.V.; Methodology, J.M.-D., P.B. and F.B.; Project administration, J.M.-D. and P.B.; Supervision, P.B., F.B., P.P. and B.B.; Visualization, J.M.-D.; Writing—original draft, J.M.-D.; Writing—review and editing, P.B., F.B., P.P. and B.B.

**Funding:** This research was funded by NSERC, grant numbers CRSNG-RDCPJ: 523095-17 and CRSNG-RGPIN-2016-06780.

**Acknowledgments:** The authors gratefully acknowledge the Town and M. Rybicki to allow access and water sampling on their property. The authors would also like to thank the students (A. Siméon, T. B. Quézel, T. Crouzal, D. Dufresne, J.-S. Grenier, just to name a few) who participated in the fieldwork. For the laboratory work, we would like to thank M. Leduc and J. Leroy from the Laboratoire de géochimie de Polytechnique Montréal.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Evaluating Bank-Filtration Occurrence in the Province of Quebec (Canada) with a GIS Approach**

**Marc Patenaude 1,2,3,\*, Paul Baudron 1,2,3,4,\*, Laurence Labelle 1,2,3 and Janie Masse-Dufresne 1,2,3**


Received: 20 November 2019; Accepted: 20 February 2020; Published: 1 March 2020

**Abstract:** Due to the abundance of surface water in the province of Quebec, Canada, it is suspected that many groundwater wells are pumping a mixture of groundwater and surface water via induced bank filtration (IBF). The regulatory framework in Quebec provides comprehensive guidelines for the development and monitoring of surface water and groundwater drinking water production systems. However, the regulations do not specifically address hybrid groundwater-surface water production systems such as IBF sites. More knowledge on the use of IBF in the province is needed to adjust the regulations with respect to the particularities of these systems. In order to provide a first evaluation of municipal wells potentially using IBF and the corresponding population served by these wells, a Geographic Information Science framework (GISc) was used to implement an IBF spatial database and calculate the distance from each well to the nearest surface water body. GISc is based on open source GIS programs and openly available data, to facilitate the reproducibility of the work. From this provincial scale approach, we show that nearly one million people are supplied by groundwater from municipal wells located <500 m from a surface water body, and half a million have a significant probability to be supplied by IBF wells. A more focused look at the watershed scale distribution of wells allows us to improve our interpretations by considering the aquifer type and other regional factors. This approach reveals strong spatial variability in the distribution of wells in proximity to surface water. Of the three selected regions, one has a high potential for IBF (Laurentides), one requires additional information do draw precise conclusions (Nicolet), and the third region (Vaudreuil-Soulanges) is unlikely to have widespread use of IBF. With this study, we demonstrate that extensive use of IBF is likely and that there is a need for improved understanding and management of these sites in order to properly protect the drinking water supply.

**Keywords:** managed aquifer recharge (MAR); induced bank filtration (IBF); geographic information science (GISc); geographic information systems (GIS); drinking water supply; guidelines

#### **1. Introduction**

Induced bank filtration (IBF) is a widely used method of managed aquifer recharge (MAR) [1]. In an IBF system, surface water is drawn through the banks and bed of a lake or river towards a pumping well by an induced hydraulic gradient. This results in a pumped mixture of both groundwater and infiltrated surface water. This process reduces the risk of intensive use [2] of the aquifer and improves water quality relative to surface water [3]. The drawback of this method, however, is that there is a higher risk of contamination from certain sources in comparison to standard groundwater exploitations [4–7],

which makes it important to develop a specific regulatory framework adapted to this type of drinking water supply system.

In recent reviews of the international use of bank filtration [8–11], it is clear that IBF has been used worldwide and, in particular, in Europe for more than 100 years. In the last decades, the use of IBF also became relatively popular in the USA [3 and references therein]. More recently, few IBF systems have also been implemented in developing countries [3 and references therein]). The use of this technology is, however, completely overlooked in Canada and the province of Quebec in particular. As of this point in time, no inventory of IBF sites exists for the province, and the extent of its use throughout the province remains unknown. There is, however, a high probability that IBF is widely used. Indeed, Quebec is a province rich in both surface water and groundwater, with 22% of its surface covered by water either in the form of lakes, rivers, or wetlands [12]. It is estimated that Quebec contains 3% of the Earth's renewable freshwater resources [13]. In addition, Quebec was settled through its waterways, which has had long-term effects on its population distribution. The abundance of water and the distribution of the population mean that there is a high probability that pumping wells are located in close proximity to surface water. The process of unintentional IBF is thus likely to be occurring in many municipalities throughout the province.

In addition to the number of IBF sites throughout Quebec being unknown, the current regulations and guidelines in Quebec do not specifically address wells that pump a mixture of surface water and groundwater. In Quebec, the *Règlement sur le prélèvement des eaux et leur protection* (RPEP) provides comprehensive guidelines for protecting surface water and groundwater extractions from contamination [14]. This regulation provides a framework for monitoring contaminants that allows for appropriate intervention when groundwater resources contain surface microbiological contaminants, which ensures that the population is not likely to be affected by changes in water quality. Additionally, drinking water production systems that are considered Groundwater under the Direct Influence of Surface Water (GWUDI) are subjected to the same regulations as surface water, which are more stringent than for all groundwater wells. For instance, in such a case, the raw pumped water is tested weekly for microbiological parameters, whereas groundwater can be tested monthly [15]. Despite its name, the GWUDI classification does not aim at characterizing infiltration from surface water bodies. Rather, what is referred to as "surface water" within this classification system is in reference to any source of contamination from the surface (including septic tanks) that might provide recurrent microbial or viral contamination to a well. There is, therefore, no correlation with surface water bodies that are hydraulically connected to an aquifer unless said surface water bodies is considered a potential source of contamination or if there is persistent bacteriologic or virologic presence in the pumping well during initial characterization. Existing non-GWUDI IBF sites are therefore necessarily treated as standard groundwater extractions. This means that even if the groundwater wells are located within a few meters of a water body, and there is a hydraulic connection between them, there is a lack of protection guidelines for the nearby surface water.

Due to the lack of regulations specific to hybrid groundwater-surface water systems, many problems common to IBF sites may be overlooked or unanticipated during the planning and operation of these sites. IBF sites are sensitive to changes in surface water quality [6,16], changing redox conditions [17], and changes in the hydraulic conditions of the site [18]. Additionally, there are many undesirable chemical components and contaminants that can be found at pumping wells in IBF sites including Mn [19–21], Fe [22,23], NO−3<sup>−</sup> [24], organic micropollutants [25–28], cyanobacteria [29–31], coliforms [32]. In order to have a more resilient water supply, the risk posed by these contaminants should be identified as early on in the development of the site as possible in order to minimize unforeseen costs and develop plans to reduce the risk. By identifying the potential contaminants at the sites related to changing hydraulic and chemical conditions, the risks associated with them would be reduced by allowing for strategic planning to avoid the conditions associated with poorer water quality. For example, well configuration and pumping schemes could be modified to increase transit times from surface water to pumping wells producing a more consistent water quality.

Geographic Information Systems (GIS) and Multi-Criteria Decision Analysis (MCDA) approaches are being developed for site suitability mapping for the development of MAR ([33] and references therein) and IBF [34,35]. As an example, Jamarillo Uribe [36] has studied the impact of stream morphology on bank filtration sites. Nonetheless, to our knowledge, this is the first study that uses a GISc approach to identify the number of existing IBF sites. The objective of this study is to provide a first overview of the potential extent of IBF in the province of Quebec. The framework is based on the concept of Geographic Information Science (GISc), described in [37] and [38]. GISc is based on the use of open-source GIS programs and openly available data. This concept facilitates reproducible research to other areas of study and data sets and could improve the transparency of research using GIS. First, we have processed and homogenized the sources of information from three different agencies. Secondly, we have carried out a pre-quantification of municipal wells with a higher likelihood of providing drinking water through IBF using easily available government data. Following this pre-selection, zones with varying characteristics are considered in greater detail, and the likelihood of IBF taking place in these areas is discussed.

#### **2. Materials and Methods**

#### *2.1. Data Set*

#### 2.1.1. Well Data

The well data in the province of Quebec is contained in a variety of databases. All private wells require that various information is reported following drilling, such as coordinates, ownership, well design (i.e., type of tubing and depth of well), and stratigraphic sequence. This data is compiled in the *Système d'Information Hydrogeologique* (SIH) database [39], which has previously been described by Sterckx [40]. The SIH database consists of information largely provided by well-drilling companies. This can lead to inconsistencies in the geological descriptions and the precision of the coordinates. This database is, however, considered a reliable source of information concerning the depth of the contact between overburden and bedrock [41].

The database with the most extensive information is available through a series of studies entitled "*Programme d'Acquisition de Connaissances sur les Eaux Souterraines*" (PACES) [42]. These projects were initiated in 2008 and aimed at improving knowledge of groundwater resources in the southern regions of the province of Quebec in order to protect them and ensure their sustainability. These studies have led to a number of subsequent publications [43–46]. As of today, the PACES studies have been completed on a total of 13 regions throughout the province and led to the compilation of varied information (i.e., depth of well, depth of screen, depth of water table, type of aquifer, type of well and coordinates) on a total of roughly 180,000 wells. Supplemental information, including geochemistry, geology, and hydrogeological data (e.g., hydraulic conductivity), was also compiled for a small subset of wells (n = 15,162). The PACES database contains significant overlap with the SIH database, and therefore, some issues with the reliability of the coordinates and geological descriptions also affect this database.

The third database is comprised of municipal wells and surface water extraction points [47]. This database is typically used by decision-makers for public health and land-use planning. The version of the database used in this study was acquired in 2017 (i.e., prior to the most recent update in November 2018). The information contained in this database is centered on geographic coordinates, population served, and types of potabilization treatment. It contains information on 2116 individual wells and surface water extraction points. The geographic coordinates are generally more recently acquired and derived from official declaration documents, resulting in a better precision of the localization of the wells than the two other databases.

The municipal database formed the basis of this study with complementary information pulled from the other sources. The decision to focus on the municipal well database was made since (i) IBF requires sufficiently high pumping rates to induce a hydraulic gradient from the surface water to the well, (ii) the localization of the wells is precise and (iii) all these wells are supplying drinking water to the population. The PACES database was used to extract complementary information and draw some general conclusions about the distribution of wells within the province. Table 1 summarizes the data sets and the variables used in this study.



\* Databases used for subsequent calculations; (X) not systematically compiled.

#### 2.1.2. Surface Water Bodies and Other Data

Natural Resources Canada's (NRCAN) website contains a variety of vector data in the Canvec portion of the site [48]. The surface water files are subdivided into two categories, "watercourses" and "water bodies". Many of the features in the "watercourses" files are drainage ditches and ephemeral streams that are not likely to be supplying bank filtrate to wells throughout the year. The features in the "water bodies" files correspond to water bodies of larger size and more permanent nature, which are more likely to contain a sufficient volume of water to support a municipal water supply. Considering the above, it was decided to conduct our province-wide study with the "waterbodies" files only. The named rivers from the "watercourses" files were added to the regional studies. These files contain two types of data, i.e., polygons and polylines, respectively representing the water bodies and shorelines.

#### *2.2. Description of the GISc Framework*

In this study, we used a GISc framework in order to calculate the minimal distance between each municipal pumping well and the nearest surface water body (i.e., lake or river). Throughout this study, the distance to surface water is calculated with respect to the "waterbodies" files retrieved from the NRCAN website, as mentioned above. The work process has been carried out with the Quantum GIS program (QGIS) [49]. In this subsection, we list a number of steps that were taken in order to perform the spatial analysis.

2.2.1. Processing and Homogenization of Spatial Data Sets

• Removal of duplicate wells

The work process is initialized with spatial data debugging. The database of municipal wells contained a number of duplicated geometries, i.e., wells with identical coordinates, which caused issues when performing the calculations. The duplicate data most often resulted when wells were supplying water to multiple municipalities. In order to remedy this, duplicate geometries were automatically identified and removed from the spatial database. In cases where discrepancies were found in the attribute table, the correct points were identified and retained.

• Homogenization of the spatial reference system

Coordinates needed to be converted into a common projected coordinate system in order to do subsequent distances calculations. Spatial data were reprojected from the original coordinate system (i.e., NAD 83) to Lambert EPSG: 32,198. The Quebec Lambert was selected since it is suitable for use in Quebec (Canada) and for applications with an accuracy of <2m[50].

• Conversion of geometries

First, the geometries were converted from multipart to single part. This step was necessary in order to ensure that the conversion to lines in the following step inserted all the necessary points. Without this step, only the existing nodes of the polyline file were converted to points, which led to an uneven and wide-ranging spread of points along the shorelines.

Second, waterbody files were converted from lines to points with 2 m spacing with the function *Convert lines to points* in the SAGA [51] toolbox. This function was run as a batch process for all of the 1:50,000 National Topographic System (NTS) zones individually.

2.2.2. Performing Distance Calculation

• Calculate the distance from each well to the nearest point

Using the smaller 1:50,000 zones allowed for the direct calculation from each well to the nearest point in each of the NTS zones. This calculation was completed using the *Distance to nearest hub* function. This algorithm identifies the nearest feature to each point and the Euclidean distance between all points. This, however, led to the creation of roughly 220 different files, each containing one distance value for each well. This avoided some problems that occurred when the calculation was done in the reverse order, but it created extremely large files that were not easy to work with in QGIS, especially for the larger PACES files.

• Calculate the minimum value for each well

The subsequent step was to merge all of these roughly 220 files and calculate the minimum value for each well. This was done in QGIS using the *merge* function, and the minimum selection was made using the *Select by expression* dialogue box and the following line of code HubDist = minimum (HubDist, Well ID). This could only be done in QGIS for the smaller municipal well files. For the larger PACES files, the merge and distance calculations were completed using R [52] programming language.

Once the distance was calculated, manipulation of the data was done manually using filters in QGIS, and R. Additional data was also joined to the well files either manually, or using common fields or spatial relationships using a variety of functions in QGIS.

This process was repeated for the regional studies with named rivers from the watercourse file and the water bodies used in the province-wide study.

#### **3. Results and Discussion**

There are two main conditions that will determine whether a well is pumping a mixture of surface water and groundwater. First, a hydraulic connection between the surface water and the aquifer is necessary. Second, assuming groundwater typically discharges into surface water bodies in the study region [53], it is necessary to pump a sufficient volume to reverse that natural gradient and draw surface water toward the pumping well. In order to conduct a study at the scale of the province, we used the municipal well database (description available in Section 2.1.1), which reports the geographical coordinates and the population served by each of the 2075 municipal wells and surface water extraction points. Since the type of aquifer is not reported within this database, the distance from the surface water bodies was deemed to be the most important criterion for estimating the number of IBF sites available for a province-wide scale. Then, we selected three areas with distinct characteristics in order to conduct regional analyses, including integrating the type of aquifer and other hydrogeological information. The results of both province and regional scale approaches are presented and discussed in the next subsections.

Bank filtration sites can be located at various distances from surface water bodies [32,54] from a few meters to greater than one kilometer. In this study, as a starting point, distances >500 m are considered less likely to be performing IBF, as wells were located <500 m from surface water in several studies on IBF [32,54,55]. Groundwater will discharge from aquifers towards rivers and lakes in most of Quebec for the majority of the year, especially during drier months [53]; therefore, wells must be able to reverse that natural gradient in order to perform IBF. In addition, most of the sites have been operational for at most a few decades, are not likely to have experienced intensive use [2], and therefore, large scale flow regime reversals, as seen in the Netherlands [54], are not likely. These combined factors make it unlikely that wells at a distance >500 m are performing IBF since reversal of a gradient over those distances would be more difficult. This choice for the cut-off could potentially omit some bank filtration sites from the province; however, this is considered to have minimal impact on the overall conclusions of this study.

#### *3.1. Overview of the Province*

#### 3.1.1. Municipal Drinking Water Supply: Surface Water vs. Groundwater

Municipal drinking water sources can be classified into three categories, i.e., groundwater, surface water, and groundwater considered surface water, as shown in Figure 1. The first category (in red) includes all municipal drinking water sources which rely on one or multiple groundwater wells. Contrastingly, the second category (in blue) refers to municipal drinking water sources supplied by surface water. The third category (in green) corresponds to the few municipal drinking water sources relying on groundwater wells, which are documented as GWUDI according to the protocol detailed in the *Guide de conception des installations de production d'eau potable* described in Section 1 [15]. Of the 2075 municipal extraction points, 87% (n = 1799) are groundwater pumping wells, whereas 13% (n = 276) are directly extracting surface waters. This results in approximately 15% of the population of Quebec relying on groundwater resources for drinking water, representing roughly 1,260,000 citizens. A similar estimate (i.e., 20%) was also reported by the *Ministère de l'environnement et lutte contre les changements climatiques* [42]. The discrepancy between these two estimates is likely due to the proportion of the population supplied by private wells rather than a municipal distribution network.

In the more densely populated areas, surface water is the main water source for drinking water supply systems. This is likely a consequence of the greater drinking water demand in the cities, and the larger population's ability to support the more costly water treatment required for surface water. In smaller municipalities, groundwater sources are preferred as they are less costly to operate. In fact, in rural areas of the province, 90% of the population is served by groundwater sources [56]. As explained above, the important number of surface water bodies and the population's distribution results in a high likelihood of having groundwater wells near a surface water body, suggesting that many municipal drinking water systems may be benefiting from IBF processes.

**Figure 1.** Spatial distribution of municipal drinking water sources. The water sources are classified into three categories, i.e., groundwater (in blue), surface water (in red), and groundwater considered surface water (in green). (**a**) overview of the study area; (**b**) view of the distribution of all wells throughout the province; (**c**) view of the most densely populated area of the province along the Saint-Lawrence River

#### 3.1.2. Water Bodies Distribution around Wells

Among the municipal wells, almost all (97% of the cases; n = 1749) are located at <2000 m from a surface water body. As illustrated in Figure 2a, it is evident that the closest water body to most municipal wells is lakes (72%; n = 1262). This is likely due to the extremely high number of lakes in Quebec, which means that there is likely always a lake within a reasonable distance from any given point in the province, resulting in a geographically homogeneous distribution of municipal wells near lakes. The distribution of wells in close proximity to rivers (28%; n = 487) is less homogenously distributed throughout the province compared to those in proximity to lakes (see Figure 2). Many municipalities are located in close proximity to a river due to the settlement of the province through its waterways. Secondly, the area around rivers can often have more favorable properties for larger-scale water extraction due to the higher likelihood of containing sandy granular deposits compared to the more common glacial-marine deposits that cover the rest of the province [57].

Figure 2b also reveals the same trend of homogenous distribution of wells in proximity to lakes for each bin of 10 m width. In fact, 14% (n = 245) of the municipal wells are located at <200 m from a river. The overall distribution of wells shows a significant decrease in the density of wells at around the 120 m marks. Since the number of wells in proximity to lakes remains relatively constant in all 10 m bins, and the distribution of wells near rivers decreases around the 120 m marks, the overall trend in this graph of more wells in the first 120 m is controlled mostly by the greater number of wells in proximity to rivers.

**Figure 2.** The distance of municipal wells from water surfaces for (**a**) 0 to 2000 m and (**b**) 0 to 500 m.

3.1.3. Insights into the Potential Population Supplied by IBF

Figure 3 illustrates the cumulative population served by municipal wells that are located at a distance of 0 m to 2000 m from a surface water body. It shows that approximately 1,200,000 people are served by those municipal wells. It also reveals that outside of the major cities, wells located at less than 500 m from a surface water body account for 74% of the population (n = 920,000) whose drinking water is supplied by groundwater from municipal wells. In fact, there is a rapid increase in the population served by wells within the first few hundred meters of a water body. Wells within 250 m and 100 m account for 57% (n = 720,000) and 34% (n = 410,000) of the population served by municipal wells respectively. Although many of these sites may not be using IBF due to various factors, these results reveal a significant probability that half of the groundwater-fed population connected to a municipal drinking water supply in Quebec, i.e., more than half a million citizens, might depend on IBF. This initial overview demonstrates the great opportunity of a better assessment of IBF occurrence in the province of Quebec in order to better protect our resources. This high proportion of close-to surface water wells is likely an indication that planners and developers are, in fact, selecting pumping sites in close proximity to surface water. Still, as aforementioned, the existing knowledge of IBF is not fully integrated into the planning process when developing and managing water abstraction plants.

**Figure 3.** Cumulative population supplied by municipal groundwater wells with respect to the distance from the surface water body.

#### *3.2. Insights on Selected Sub-Basins*

In light of the results presented in the previous subsection, it was determined that making widespread generalizations about the province would not be straightforward. We opted to focus on a selection of areas with diverse population sizes, diverse geological settings, a variety of distances from surface water, a variety of surface water body types, and a variety of land cover types. By cross-checking the ID and localization of wells in the available databases, it was possible to manually assign the type of aquifer (fractured bedrock, unconfined granular, confined granular) to 101 municipal wells within the three areas, namely Laurentides (area #1), Nicolet (area #2) and Vaudreuil-Soulanges (area #3). These areas correspond to watersheds "du Nord", "Nicolet", and "Vaudreuil-Soulanges"respectively. It is important to note that the "Vaudreuil–Soulanges" watershed also extends into the neighboring province of Ontario to the West, but the investigation is limited to the portion within the Quebec border. The localization and the extent of each area are illustrated in Figure 4. The main geological contexts are also reported in this figure. The Grenville Province (area #1) is mainly composed of Archean autochthonous rocks dominated by highly metamorphosed gneissic complexes. The Appalachian Province (area #2) is composed of various types of strongly deformed rocks (i.e., sedimentary, volcanic, and ophiolitic rocks), whereas the St-Lawrence Platform (area #3) corresponds to sedimentary rocks. The criteria used to evaluate each region are the nature of quaternary deposits, the type of nearest water body in proximity, the population of the municipalities, the distance of wells, as well as certain qualitative factors that are mentioned in the PACES reports for each region [58–61].

**Figure 4.** Location of the three selected watersheds (area#1: Laurentides, area #2 Nicolet, and area #3 Vaudreuil–Soulanges).

#### 3.2.1. Area #1: Laurentides

The area of Laurentides has a mixture of forested areas and urban areas and is underlain by the hilly Grenville Province. It has not been covered by regional PACES studies at this time, and therefore, a detailed description of the surficial deposits underlying the region is somewhat dated. Portions of the area along the Saint Lawrence and Ottawa Rivers were described by Lajoie [58]. Bedrock in this area is overlain by till that varies in thickness from 7 to 12 m. Deposits of sand and gravel can also be found within this region, that were deposited by glacial rivers. Along existing rivers, including the Rivière-du-Nord, recent fluvial deposits composed mainly of loam are present [59]. Supplemental Materials containing the map of the quaternary geology and well distribution for each area is provided with this paper (Figure S1, Tabel S1).

In this area, 31 municipalities (of a total of 58) rely on groundwater to produce drinking water from a total of 136 municipal wells, and 91% (n = 124) are located <500 m from a surface water body. A total of 63 municipal wells (51%) are located at <500 m from a lake, the remaining part (49%; n = 61) being at <500 m from a river. As shown in Figure 5, a significant number (31%; n = 39) of municipal wells are found along the main river (i.e., Rivière-du-Nord). These wells are often serving municipalities with a population of >1000 people and are typically found in unconfined granular aquifers. There are three other rivers and a lake with wells in close proximity, namely, Red River with four wells serving two municipalities, Ottawa River with four wells serving three municipalities, Lac de Deux Montagnes with 12 wells serving one municipality, Achigan River with two wells serving one municipality. Further from the Rivière-du-Nord, the population generally decreases, and the wells are more frequently found in proximity to lakes and in fractured bedrock aquifers. The type of aquifer could be compiled for a total of 56 wells. Of these wells, 73% (n = 41) are found in an unconfined granular aquifer, with the majority located <100 m from a surface water body. Since the largest municipalities are more likely to be pumping sufficiently to induce a hydraulic gradient forcing the surface water to infiltrate the aquifer, the municipalities with the highest likelihood of performing IBF are the more populated ones in proximity to the Rivière-du-Nord. Meanwhile, wells located further from surface water and installed in fractured bedrock present a lower confidence level in the probability that IBF is taking place. It is important to note that the above-mentioned results include some wells outside of the selected watershed, as they shared a similar geological setting.

**Figure 5.** (**a**) The spatial distribution of municipal wells located at less than 500 m from lakes and rivers in area #1 and the distribution of wells according to (**b**) the type of surface water and (**c**) the type of aquifer.

#### 3.2.2. Area #2: Nicolet

The Nicolet area straddles two geologic provinces, the Saint-Lawrence Platform in the lower altitude portion to the North and the Appalachian Orogen in the higher altitudes to the South. This area is principally covered by agricultural land. A total of 84 municipal wells are actively producing drinking water from groundwater resources to serve a total of 35 municipalities (of a total of 39). These municipalities are generally smaller than those in other regions, with many serving <500 people. As illustrated by Figure 6, in area #2, a high number of wells are in the vicinity of different rivers and tributaries, similar to area #1. In fact, a total of 48 municipal wells (57%) are in the 0–500 m range from a surface water body. Of these wells, most are located close to rivers (71%; n = 34), and the remaining 29% (n = 14) are located near lakes. However, the distance of these municipal wells varies more widely than in area #1. The Nicolet River and its tributaries are the rivers with the largest number of wells in close proximity. There are five rivers in the Nicolet river system with wells in close proximity. In addition, the Saint-Francois River has three wells in proximity serving two municipalities. Within this region, there are many more rivers whose quality could impact drinking water quality than in the other regions.

**Figure 6.** (**a**) The spatial distribution of municipal wells in the Nicolet area that are located at less than 500 m from lakes and rivers in area #2 and the distribution of wells according to (**b**) the type of surface water and (**c**) the type of aquifer.

The Nicolet region has a distinctly different sequence of quaternary deposits when compared to the other two zones. As reported in the PACES study [60], thick quaternary deposits exceeding 100 m in certain areas overlie the bedrock in certain parts of the region. A rough 20 km wide band along the Saint–Lawrence River is underlain by a significant thickness of marine clays, which can be partially or completely overlain by marine and lacustrine sands. The central portion of the basin located between elevations of 80 m and 120 m is dominated by aeolian and coastal sands, underlain by impermeable till around which peatlands can form. Thick glacial-fluvial sand and gravel deposits can lie directly on the bedrock in certain areas. Superficial deposit thicknesses at the municipal well locations calculated from the raster database that accompanies the PACES study [60] ranged from <5 m up to 30 m. Supplemental Materials containing the map of the quaternary geology and well distribution for each area is available (Figure S2, Table S1).

Geological information could be compiled for a total of 26 municipal wells located <500 m from a water body. In this region, 50% (n = 13) of the municipal wells procure water from unconfined granular aquifers. Of these wells, 77% (n = 10) are located <200 m from a surface water body. These results suggest that there are favorable aquifers for pumping near rivers, while the regions more distant from rivers are less favorable. This region, although it consists mostly of smaller municipalities, has a high probability of IBF taking place, especially in wells within the first few hundreds of meters from a river. More information relative to the pumping rates for these wells, combined with a geochemical and isotopic approach, would be needed to better estimate the potential use of IBF in this region.

#### 3.2.3. Area #3: Vaudreuil–Soulanges

The area of Vaudreuil–Soulanges is located near Montreal and within the St-Lawrence Platform. The main types of land uses are agricultural and residential. This area, unlike the Nicolet area, is underlain solely by more recent quaternary deposits that were deposited uniquely during the last glacial cycle. The region is predominantly covered by marine clays (64% of the surface area). Only a small number of areas, along the Ottawa River and higher relief areas (i.e., till deposits, Mount Rigaud, "butte Saint-Lazare", and "butte de Hudson"), remain uncovered by these clays. Certain uncovered areas are composed of glacial-fluvial deposits that host productive granular aquifers. Examples of these deposits can be found along the Ottawa River in the northern portion of the zone and also in the Saint-Lazare region. Supplemental Materials containing the map of the quaternary geology and well distribution for each area is available (Figure S3, Table S1).

In this area, 36 wells are in operation and produce drinking water from groundwater resources for eight municipalities (of a total of 13) with populations typically >1000 people. As illustrated in Figure 7, the area contains a distribution of wells mostly in proximity to lakes. In fact, 66% (n = 24) of the municipal wells are located <500 m from a lake. Among these wells, 92% (n = 22) are located along a 6 km-long North-South transect near the region of Saint–Lazare. The PACES report highlights the presence of a thick and unconfined sandy deposit with a very productive aquifer in this zone [60]. Wells, in proximity to rivers, for the most part, are located along the periphery of this zone near the more major water bodies. Another zone of concentrated wells (20%; n = 10) is found in the region of Mount Rigaud, the majority of these wells are located in fractured rock or clay deposits [60], making the probability of IBF unlikely. It can also be seen in Figure 7 that there is a limited number of wells (33%; n = 16) located within the first 200 m of a water body, which sharply contrasts with the other two regions. In this region, wells are located in proximity to two rivers only. The Viviry River with four wells serving one municipality and the Saint–Lawrence River with four wells serving one municipality are the two most prominent ones with some smaller streams making up the difference. The type of aquifer could be compiled for 22 municipal wells in area #3. From these wells, the majority is found in confined granular aquifers (31%; n = 7) and fractured bedrock aquifers (40%; n = 9). The fractured bedrock aquifers are known to be confined [60], except in the regions with the highest concentrations of wells (i.e., St-Lazare, Hudson, and Rigaud).

Overall, the greater distance from surface water in this area seems to indicate that IBF in most of the municipalities in this area is not likely. We identified only two municipalities with a total of 3 wells in close proximity to surface water. However, these wells are in confined granular or fractured bedrock aquifers and, thus have limited potential for IBF. Moreover, there is a strong possibility that the permeable areas of the Saint–Lazare region are influenced by infiltrating water from the spring thaw and precipitations.

**Figure 7.** (**a**) The spatial distribution of municipal wells located at less than 500 m from lakes and rivers in area #3 and the distribution of wells according to (**b**) the type of surface water and (**c**) the type of aquifer.

#### **4. Conclusions**

This research has provided a reliable starting point for determining the impact of IBF on the population of Quebec. A simple process based on GISc was used by incorporating several open-source programs (QGIS, SAGA, and R) and openly available data. The initial use of distance from a surface water body as well as additional information extracted from the government database revealed that nearly one million people in the province of Quebec are supplied drinking water via wells in close proximity (i.e., <500 m) to a river or a lake. This first overview has also demonstrated that there is a high degree of regional variability with regard to the probability that IBF is being performed. One investigated area has the greatest potential use of IBF (area #1: Laurentides), a second one requires more information to draw precise conclusions (area #2: Nicolet), and the third one has the lowest probability of IBF (area #3: Vaudreuil–Soulanges).

There are a number of shortcomings that are evident in the current development, management, and understanding of hybrid groundwater-surface water extraction sites in the province of Quebec. The current regulations do not contain general insights into the protection or management of these sites. The existing GWUDI classification is not designed to assess risk at an IBF site as it does not expressly consider the contribution of surface water to a pumping well, and only under certain circumstances are IBF wells considered GWUDI. With mounting pressures on our water resources from climate change and population growth, we hope that this work demonstrates the need to better understand and improve the regulatory framework to specifically address hybrid groundwater-surface systems such as IBF in order anticipate problems common to IBF sites and ensure that our water resources are well protected and exploited sustainably.

It will be important to precisely develop guidelines for determining what municipalities should be targeted within future regulations. The policy changes would require the definition of a new category of wells for IBF sites that requires short, medium, and long-term assessment of risk related to the permanent or intermittent contribution of surface water bodies to wells.

In order to determine more precisely which wells are using IBF, a possible approach would be based on a dedicated sampling program for specific environmental tracers. These tracers would be used to quantify the spatio-temporal evolution of transit times (e.g., stable isotopes of water) and relative proportions of groundwater and surface water (e.g., electrical conductivity). A longer characterization period would allow for the minimization of unforeseen costs over the life of the water treatment plant.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/3/662/s1, Figure S1: Map of Areas #1 (Laurentides) with quaternary geology and well distribution, Figure S2: Map of Areas #2 (Nicolet) with quaternary geology and well distribution, Figure S3: Map of Areas #3 (Vaudreuil Soulanges) with quaternary geology and well distribution, Table S1: Codes for the quarternary geology units in Figures S1–S3.

**Author Contributions:** Conceptualization, M.P. and P.B.; Data curation, M.P.; Funding acquisition, P.B.; Investigation, M.P.; Methodology, M.P. and P.B.; Project administration, P.B.; Resources, P.B.; Supervision, P.B.; Validation, L.L. and J.M.-D.; Visualization, L.L. and J.M.-D.; Writing—original draft, M.P.; Writing—review & editing, M.P., P.B., L.L. and J.M.-D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by FRQNT through grant number FQRNT 2018-NC-205703 and NSERC through grants number CRSNG—RDCPJ: 523095-17 and CRSNG-RGPIN-2016-06780.

**Acknowledgments:** We would like to thank the "Ministère de l'Environnement et Lutte contre les changements climatiques" for the data used in this study, Marc-André Bourgault for help during conceptualization, Guillaume Meyzonnat and Sylvain Gagné for discussions on various zones and insights into the reliability of data sources, Geneviève Boisjoly for help with the GIS calculations, and Gabriel Dion for the initial study that evolved into this one. Great thanks to Francisco Gomariz Castillo for his constructive review from a GIS expert point of view.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Reactive Barriers for Renaturalization of Reclaimed Water during Soil Aquifer Treatment**

**Cristina Valhondo 1,2,\*, Jesús Carrera 1,2, Lurdes Martínez-Landa 2,3, Jingjing Wang 1,2,3, Stefano Amalfitano 4, Caterina Levantesi <sup>4</sup> and M. Silvia Diaz-Cruz <sup>1</sup>**


Received: 20 November 2019; Accepted: 27 March 2020; Published: 2 April 2020

**Abstract:** Managed aquifer recharge (MAR) is known to increase available water quantity and to improve water quality. However, its implementation is hindered by the concern of polluting aquifers, which might lead to onerous treatment and regulatory requirements for the source water. These requirements might make MAR unsustainable both economically and energetically. To address these concerns, we tested reactive barriers laid at the bottom of infiltration basins to enhance water quality improvement during soil passage. The goal of the barriers was to (1) provide a range of sorption sites to favor the retention of chemical contaminants and pathogens; (2) favor the development of a sequence of redox states to promote the degradation of the most recalcitrant chemical contaminants; and (3) promote the growth of plants both to reduce clogging, and to supply organic carbon and sorption sites. We summarized our experience to show that the barriers did enhance the removal of organic pollutants of concern (e.g., pharmaceuticals and personal care products). However, the barriers did not increase the removal of pathogens beyond traditional MAR systems. We reviewed the literature to suggest improvements on the design of the system to improve pathogen attenuation and to address antibiotic resistance gene transfer.

**Keywords:** organic amendments; managed aquifer recharge; contaminants of emerging concern (CECs); pathogens; new water challenges

#### **1. Introduction**

Climate change and the expansion of urban areas is a major worldwide threat to sustainable and safe drinking water supplies [1]. Managed aquifer recharge (MAR) is a technique that allows groundwater-dependent ecosystems, including rivers, to be maintained, enhanced, and/or protected with limited consumption of energy and chemicals [2,3]. MAR systems based on water filtration during soil passage have been proven to retain suspended particles and colloids, including microorganisms [4], and to favor biodegradation of chemical contaminants, resulting in significant water quality improvement [5–7]. The processes affecting pathogen transport in these aquifers are retention and inactivation, and an extensive number of factors influence them [8]. However, periodic detection of pathogens in groundwater, some with severe human health impacts [9–13], has led to strict quality requirements that effectively impede the use of lesser quality water for MAR. For instance, rainfall fails to meet Spanish regulations for reuse (too low pH and too high suspended solids), which are the regulations adopted in practice for MAR [14]. This is paradoxical because potable water treatment during the 19th century simply consisted of sand filtering to remove pathogens and resulted in a life expectancy increase of some 20 years [15,16]. This paradox is well reflected in the ongoing debate about quality requirements for artificial recharge. Health protection authorities recommend strict controls on the water used for MAR but, at the same time, several major cities have shown that recharge using wastewater can be safe [17]. As a result, the European Commission's Joint Research Center (JRC) failed to reach a consensus on MAR water quality recommendations [18]. The situation is inadequate. Prudence demands regulation, while fear hinders the actual implementation of MAR, which impedes the restoration of ecosystem services of groundwater-dependent water bodies.

Overcoming resistance requires the addressing of not only old problems (e.g., water scarcity, recovery of groundwater-dependent water bodies), but also emerging concerns [19]. Among these, we include chemicals of emerging concern (CECs), antibiotic-resistant bacteria (ARBs), and antibiotic resistance genes (ARGs). The term CECs encompasses a wide range of substances, including pharmaceuticals, personal care products, and nano- and micro-plastics, among others, which are characterized by their continuous release into the environment and their potential to impact aquatic ecosystems and eventually human health [20]. Several studies have demonstrated that even after extensive treatment, such as advanced oxidation processes and reverse osmosis, some recalcitrant CECs are still detectable in reclaimed water [21–23]. Until the turn of the millennium, it was unknown that these chemicals presented a hazard to the environment, as they generally occur at trace levels, and pharmaceuticals in particular were always found at concentrations far below the therapeutic doses prescribed for humans [24,25]. However, studies carried out since then have provided evidence that even sub-therapeutic concentrations of certain pharmaceuticals affect microbes, plants, fishes, and insects [26–28]. Consequently, the concentrations of CECs measured in reclaimed water can be biologically relevant or can increase to such levels in the unavoidable co-occurrence with other chemicals that may increase their biological activity [29]. Under the certainty that the reclaimed waters still contain CEC residues, the use of these waters as source waters in MAR may pose a risk to human and environmental health.

Biodegradation and sorption appear to be the main processes involved in water quality improvement during MAR, especially regarding CECs' behavior [30–32]. The biomass and biodiversity of the microbial community is relevant for CEC degradation [33,34] Therefore, parameters controlling microbial community such as temperature, and the amount of organic substrate available and its quality (which controls the redox conditions), have a direct effect on biodegradation rates [4,35,36]. Changing parameters within an aquifer could lead to an increase in the microbial biodiversity, and a continuous source of organic substrate should allow the biomass to increase [37].

Sorption might be relevant as well since it retards contaminants [38], thus increasing the time available for the microbial community to degrade them. The discussion on whether retardation is favorable remains open. On the one hand, increasing the residence time would increase the ability of microorganisms to degrade them. On the other hand, some authors argue that contaminants are not biologically available when adsorbed, and therefore they are not potentially biodegradable [39,40].

To favor these two processes, we proposed adding a reactive barrier at the bottom of the infiltration basin in a MAR system. The barrier provides a reactive surface and diverse sorption sites, and adds organic carbon to yield a range of redox states. Ideally, this should allow diverse microbial communities to develop, thus increasing CECs' removal.

In this context, the goal of this paper was two-fold. First, we summarized what we have learned in two experiences of MAR using a reactive barrier [6,7,40,41]. Second, based on this adquired knowledgeand the results of others, we have discussed how to improve the system design and operation to enhance not only the removal of CECs but also the attenuation of pathogens, to minimize the transport of ARGs.

#### **2. Materials and Methods**

#### *2.1. The Concept of the Reactive Barrier*

We designed a reactive barrier to be installed at the bottom of soil aquifer treatment (SAT, the specific term for intermittent infiltration of reclaimed water) infiltration basins. The barrier consisted on an organic substrate able to release dissolved organic carbon (DOC) to the infiltrated water and to provide potential sorption surfaces. The purpose of the reactive barrier is to favor biodegradation by generating a redox zonation and enhanced adsorption for the widest possible range of CECs.

Figure 1 displays redox zonation during infiltration periods in a conventional SAT system and in a system with a reactive barrier. The source water should contain a labile DOC concentration higher than 6–9 mg/L in order to consume the oxygen and start to consume nitrate as the next electron acceptor. The implementation of the reactive barrier increases the concentration of DOC in the recharged water, so that available electron acceptors are consumed and redox zonation is developed further, reaching Fe- and Mn-reducing conditions, and hopefully SO4-reducing conditions. Conditions should return to aerobic during drying periods in both cases.

**Figure 1.** Schematic description of redox zonation during infiltration without and with a reactive barrier. The barrier adds dissolved organic carbon (DOC) to the recharge water, thus promoting highly reducing conditions. Ideally, the vadose zone becomes aerobic during drying periods in both cases.

This approach was tested at two sites and based on two organic substrates as organic carbon source, i.e., compost and woodchips. The first site was a pilot-scale MAR system located at Sant Vicenç dels Horts (close to Barcelona, Spain) where a reactive barrier based on compost was installed. The system operated with the barrier for four consecutive years. The second site consisted of six MAR systems with small variations in the configuration, located in Palamós (close to Gerona, Spain). One of the systems, the reference one, operated without a reactive barrier, four systems operated with reactive barriers based on compost, and the remaining system operated with a reactive barrier based on woodchips. To date, it appears that the implementation of these barriers has favored the infiltration capacity. The characteristics of each site and the performance of the tested reactive barriers were described in previous publications [6,7,41,42] and are summarized later.

Compost and woodchips were selected as organic substrates due to their capability to release DOC, their low cost, and the ease of their handling and transportation. The amount of DOC released by the organic substrate is expected to decrease with time (or with the volume of water infiltrated), but it may be compensated by the release from biomass growing in the basin, including plant roots. Still, after a period of operation, the barrier may have to be replaced.

#### *2.2. Site Description*

#### 2.2.1. Sant Vicenç Dels Horts

The Sant Vicenç site is a complex of two basins (settlement and recharge, each ~5000 m2) constructed at the side of the Llobregat River, some 15 km upstream of Barcelona (Figure 2A). The MAR system was constructed over the lower Llobregat valley sedimentary aquifer, formed mainly by gravels, sand, and a small fraction of clay [43]. The saturated aquifer thickness and vadose zone ranged from 12 to 14 m, and from 5 to 9 m, respectively, during the recharge experience (Figure 2C).

**Figure 2.** (**A**) Sant Vicenç dels Horts infiltration basin, and (**C**) the cross-section with the monitoring points and their screened sections; (**B**) six Palamós replicate MAR systems named T1 to T6, and (**D**) the cross-section of one generic replicate. Flow direction is from left to right in both cross-sections.

The MAR system was fed with the Llobregat River water, which is heavily impacted by wastewater treatment plants' (WWTPs') effluents [44]. The river water was diverted to the settlement basin, where it remained for 2 to 4 days. From there, water flowed to the recharge basin. Flow rate was measured hourly into the connecting pipe. The average infiltration rate was 1 m/d.

We installed a 65 cm thick reactive barrier on the bottom of the infiltration basin. This barrier consisted of vegetal compost and aquifer sediments in equal volumetric portions and a small quantity of clay and iron oxide. The role of the vegetable compost was to release degradable organic matter to the infiltrating water to favor changes in redox conditions underneath the basin, promoting microbiological diversity to enhance the removal of chemical contaminants [6,7,37,45], and to provide surfaces for neutral organic compound adsorption. Clay increased the sorption of cationic compounds and iron oxide facilitated the sorption of the anionic ones.

#### 2.2.2. Palamós Site

The pilot MAR system was constructed in a municipal WWTP facility on the northeastern Spanish Mediterranean coast. This facility collects wastewater from several municipalities. The population served increases to include some 90,000 inhabitants during the summer, reaching the maximum treatment capacity of the plant. As a consequence, effluent water quality varies throughout the year.

We constructed six pilot recharge systems (15 <sup>×</sup> 15 m<sup>2</sup> excavated structures, divided into six 2.38 × 15 m channels; Figure 2B), to test the effect of the reactive barrier's composition and the role of plants on the fate of CECs and pathogens. The system was fed with the secondary treatment effluent of the WWTP, which infiltrated from the basin through the barrier and further flowed along the 15 m simulated aquifer, to finally discharge at the base of the 1.5 m thick aquifer. Indeed, the pilot MAR operated as a tertiary treatment (Figure 2D). In this case, two organic carbon sources were tested: compost and wood chips. We assessed their performance by comparing the removal of more than 50 CEC and pathogen indicators to that of a reference system (infiltration without reactive barrier).

#### *2.3. Analytical Methods*

Pressure, temperature, and electrical conductivity were continuously recorded using conductivity, temperature and depth submersible dataloggers (CTD-Divers, Schlumberger water services, Delft, The Netherlands) in the source water and several monitoring points at San Vicenç dels Horts and Palamós sites (Figure 2C). Additionally, samples for chemical analysis were collected during several recharge events in both sites.

Target CECs were selected based on the frequency of their detection in the aquatic environment, and since the analytical methodology for each site was different, the final list of CECs analyzed in each case was defined according to the methodology requirements and the source water type (urban, hospital effluents, agricultural, or industrial).

At Sant Vicenç dels Horts, 51 CECs were analyzed in the collected samples following the method described by Nödler et al. [46] (Table S1). Briefly, the samples were allowed to settle overnight at 4 ◦C and the supernatant was recovered. A 500 mL aliquot of the supernatant was spiked with 10 μL of an internal standard solution and with 5 mL of a buffer solution before solid-phase extraction (SPE). The extraction and purification was performed using OASIS HLB (6 mL, 500 mg; Waters, Eschborn, Germany) cartridges. The analytes were eluted from the cartridges, and the extracts were evaporated with a stream of nitrogen and reconstituted with ammonium acetate solution before its transference to an auto-sampler LC-vial. The analyses were performed by high-performance liquid chromatography tandem-mass spectrometry (HPLC-MS/MS, Varian Inc., Palo Alto, CA, USA).

At the Palamós site, 58 CECs were determined by online solid-phase extraction coupled to high-performance liquid chromatography–tandem mass spectrometry (online-SPE-HPLC-MS/MS) in accordance with Gago-Ferrero et al. [47] (Table S2). In this method, water samples previously spiked with an isotopically labeled surrogate standard solution were isolated, pre-concentrated, and purified using an automated SymbiosisTM Pico online SPE-(Spark Holland; Emmen, the Netherlands). The online SPE of all samples, calibration standard solutions, and methodological blanks were performed by loading 5 mL of the water samples through PLRP-s cartridges. The trapped compounds were eluted from the cartridge to the HPLC column by the chromatographic mobile phase. The chromatographic separation was achieved with a HibarPurospher® STAR® HR R-18 ec. (50 mm <sup>×</sup> 2.0 mm, 5 <sup>μ</sup>m) column from Merck using a mobile phase consisting of HPLC-grade water and acetonitrile, both with 0.1% formic acid for positive electrospray ionization, and with 5 mM ammonium acetate buffer (pH 6.8) for the negative ionization mode. MS/MS detection was performed on a 4000 Q TRAPTM MS/MS hybrid mass spectrometer from Applied Biosystems-Sciex (Foster City, CA, USA). Selected reaction monitoring (SRM) mode was applied for improved selectivity and sensitivity. Four identification points were considered, in compliance with the European Council Directive 2002/657/EC [48].

Additionally, a microbiological analysis was carried out at the Palamós site. Gram-positive and Gram-negative fecal bacteria indicator analysis was done via most probable number (MPN) detection tests following manufacturer's instructions (Colilert IDEXX, US). The log removal values were calculated considering the average of total coliforms and *Escherichia coli* concentrations in the source water, in the monitoring point immediately below the reactive barrier (O-points), and effluents (E-points) for each MAR system of Palamós site (Figure 2D). We analyzed also other substances, namely DOC, cations and anions, but they are out of the focus of the present study.

#### *2.4. Assesing the Reactive Barrier E*ffi*ciency*

We estimated first-order degradation rates (λ) and retardation coefficient (R) for 10 CECs at Sant Vicenç dels Horts site [41,49] and compared them with those reported in the literature from other experiments. While λs can be highly uncertain, they can be considered "relative measures for comparison" of results [50]. To this end, we measured pressure, electrical conductivity, and temperature at different distances from the infiltration basin. We also performed a pulse injection tracer test to obtain the residence time distribution of the recharged water at six monitoring points. The heads and breakthrough curves were used to estimate the flow and conservative transport parameters of the aquifer using a quasi-3D numerical model [49]. The model was built using the finite element code Transdens [51–53]. Secondly, CEC concentrations measured in the source water and nine monitoring points were used to estimate λ and R (in three and two defined subdomains, respectively) for 10 CECs [41,49]. The subdomains for λ estimation were BARR (reactive barrier), UZ (unsaturated zone), and AQF (aquifer). The subdomains for R estimation were BAR (reactive barrier) and UZ+AQF (unsaturated zone and aquifer). Pathogen indicators were not analyzed at the Sant Vicenç dels Horts site.

At the Palamós site, we compared the reduction in the concentration of 58 CECs, classified into four groups namely UV filters, paraben preservatives, pharmaceuticals, and total contaminants' load Table S2). Besides, the pathogen indicators (total coliforms and *E. coli*) were also measured along the reference system and compared to those obtained for the systems operating with the two reactive barriers to assess the efficiencies of these designs.

#### **3. Results and Discussion**

#### *3.1. CECs Behavior*

Figure 3 shows the comparison among the estimated degradation rates, λ (A), and retardation coefficients, R (B), at Sant Vicenç dels Horts and those reported from other studies carried out both in laboratory experiments and at field sites [54–62]. Considering the typical maximum duration for columns (1–2 years) and field experiments (10–20 years), the minimum value for λ was allowed to be 10<sup>−</sup>3, and 10−<sup>4</sup> d<sup>−</sup>1, respectively, following the approach from Greskowiak et al. [63].

Estimated λs for the Sant Vicenç dels Horts site, operating with the reactive barrier, were similar to or higher than those reported in the literature. Indeed, the λs estimated for the reactive barrier subdomain (BARR) tended to be much larger than literature values, whereas the λs estimated for the aquifer domain (AQU) were comparable, suggesting the proper performance of the reactive barrier. Estimated Rs were also much higher than literature values for the barrier (BARR) and comparable for the aquifer (UZ + AQU). At this site, we did not have the opportunity to operate the system with and without reactive barrier simultaneously while keeping the remaining variables identical to assess the performance of the barrier. However, the proper assessment of the barriers' efficiency was performed at the Palamós site. There, we estimated removal efficiencies of the analyzed CECs for the reference system (T2) and the systems operating with the two reactive barriers (T4 and T5). For the sake of clarity, we compiled the information on the analyzed CECs grouped into four categories, as aforementioned i.e., UV filters (ΣUVF), paraben preservatives (ΣPBs), pharmaceuticals (ΣPhAC), and total load of the analyzed CECs (ΣTOTAL) (Table S2).

Figure 4 displays the removal efficiencies for the target contaminants (Table S2), estimated for the reference system (T2) and the systems operating with reactive barriers (T4 and T5) in two recharge periods, i.e., January and March 2018. Overall, significant removal, from 40% to 100%, was observed in all three systems, supporting the robustness of MAR in improving recharged water quality. Additionally, the systems operating with the reactive barrier performed equally well to or better than the reference system, indicating that the reactive barrier was successful at enhancing CEC removal. However, differences between the two barrier types' efficiencies were compound-dependent. This finding was in agreement with the results reported from a laboratory study carried out by Bertelkamp et al. [64], who investigated the sorption and biodegradation behavior of 14 CECs in soil columns under oxic conditions. They concluded that the presence of ethers and carbonyl groups increased biodegradability, whereas ring structures, amines, aliphatic ethers, and sulfides hindered degradation [64].

**Figure 3.** (**A**) First-order degradation rates (λ) (estimated for the three domains) and (**B**) retention factors (R) (estimated for two domains) at the Sant Vicenç dels Horts site for carbamazepine (CBZ), gemfibrozil (GMZ), ibuprofen (IBU), primidone (PRM), sulfamethoxazole (SMX), tolyltriazole (TOL), iohexol (IOHe), iomeprol (IOMe), iopamidole (IOPa), and iopromide (IOPr) and comparison with data from the literature obtained in field and laboratory column experiments.

**Figure 4.** Removal efficiencies for organic UV filters (ΣUVF), paraben preservatives (ΣPB), pharmaceuticals (ΣPhAC), and total CECs (ΣTOTAL) estimated at the effluent of the reference system (T2), the system operating with the reactive barrier based on compost (T4), and the system operating with the barrier based on woodchips (T5) during (**A**) January 2018 and (**B**) March 2018 recharge episodes.

Our results demonstrate that a reactive barrier improved the removal of CECs, but other options are possible. Several authors have reported that oxic conditions and low biodegradable dissolved organic carbon (BDOC) favor CECs' degradation [57,65], whereas others suggest that anoxic conditions favor the degradation of a broader range of CECs [6,7,66–68]. Regnery et al. [69] proposed a method to reduce BDOC while boosting aerobic conditions during MAR by coupling two MAR systems, riverbank filtration followed by an aeration step prior to soil aquifer treatment. The goal was to reduce the BDOC during the riverbank filtration and to induce aerobic conditions with the aeration, providing oxic and low-BDOC conditions to the second MAR system [69,70]. In our experiences, we achieved reducing conditions during the recharge, assuming that aerobic conditions would be reached in the aquifer after recharged water mixed with native groundwater. This approach reduces the demands on the technique, since only one MAR system is needed. The issue is not settled, but it is clear that the optimal design of a MAR system is driven by several factors, including the water source, the hydrological characteristics of the aquifer, the geographical situation, and land availability, among others. Therefore, it is desirable to have a large set of MAR configuration options to select the most suitable for each particular case.

#### *3.2. Pathogen Removal*

The barriers' efficiency at reducing total coliforms and *E. coli* at Palamós during the recharge episode in March 2018 is shown in Figure 5. We compared their concentrations in the source water (INF) to those immediately below the vadose zone (O) and in the effluent (E) for the reference system (T2) and the systems operating with the reactive barriers based on compost and woodchips (T4 and T5, respectively). Figure 5 displays reductions between 2.5 and 5 log units for both, total coliforms and *Escherichia coli*. (*E. coli*). Similar results were observed during other recharge episodes. Overall, few or no differences among the systems were observed.

**Figure 5.** Concentrations (log MPN/100 mL) of total coliforms and *Escherichia coli* (*E. coli*)measured at the inflow (INF), immediately below the vadose zone (O-points), and effluent (E-points) of the reference system (T2, control) and the systems operating with reactive barriers based on compost (T4) and woodchips (T5) during the March 2018 recharge episode.

These results suggest that the barrier did not increase the attenuation of pathogens beyond traditional MAR systems, and, therefore, the system needs to be improved. Unfortunately, pathogen removal has traditionally been taken for granted in aquifers and, as a consequence, conceptual understanding is limited so far. Still, much research is available in the sand filtration literature, both rapid and slow sand filtration (more relevant for MAR) [71,72]. Materials other than sand have also been studied. Perez-Mercado et al. [73] explored the performance of biochar in reducing bacterial indicators from wastewater, with varying success.

While the extent of the supporting evidence is highly variable, the consensus is that pathogen fate is governed by two processes: retention and inactivation.

Retention refers to the immobilization of pathogens by straining or adsorption [8]. Straining refers to physical blocking of particles at small pores, and it is often assumed to occur at pores smaller than the bacteria size [74]. This contradicts what is known about colloidal straining, where filtration is largely caused by the lumping of particles, but it is hard to falsify. Adsorption refers to the retention of particles by electrostatic forces (actually, the set of mechanisms is much broader, ranging from diffusion into immobile water pores to the formation of surface complexes). Adsorption is usually explained using the double-, sometimes triple-layer theory.

Regardless of the actual mechanism (for quantification and upscaling purposes, identifying the retention mechanism is important), several factors affect pathogen retention, including size and accessible surface (in practice, they are hard to separate, since a small grain size leads to a high specific surface), hydraulic gradient (rather than velocity), pH (which controls the surface charge of both microorganisms and mineral surfaces), temperature, ionic strength (which affects surface potential, so that both increased and reduced retentions have been reported), and biofilm [75]. While a unified approach needs to be synthesized from the extensive literature on this topic, it is clear that a broad range of surface types and chemical states (ionic strength and pH appear to be the most controlling parameters) favor retention [76].

Inactivation is more difficult to ascertain. Pathogens tend to die outside the human body, which provides optimal conditions for their survival. The question, therefore, is how long it takes and whether it can be confirmed. Inactivation may occur in the liquid phase or, after adsorption, in the solid phase. Numerous factors contribute, including pH (in general, acidic conditions facilitate removal), temperature (survival has been observed to decrease with increasing temperature), and the presence of predators [77]. The conclusion from these studies is that most pathogens die off after a few weeks of residence in the soil. However, a few form spores or adopt spore-like forms that become inactive under unfavorable conditions but "resuscitate" when these conditions are favorable. For these, ascertaining elimination at the solid phase is critical [78].

It is precisely this variability in responses that underscores the need to observe a range of microorganisms. Given the impossibility to analyze all of them, it is common to use "indicator microorganisms" (bacteriophages, *E. coli*, *Cryptosporidium*, *Clostridium*).

The presence of metals favors inactivation. Urfer [79] showed that the addition of aluminum to (slow) sand filters enhanced the removal of bacteria. This is relevant because reclaimed water is often rich in aluminum (generally used for promoting flocculation during primary wastewater treatment). This aluminum will not precipitate as bauxite, but as gibbsite [Al(OH)3] at neutral pH. Still, its trivalent state should favor retention and would explain why biofilm ageing increases pathogen retention. Park et al. [76] conducted column experiments with *Cryptosporidium parvum* oocysts (a frequent indicator microorganism) and found that retention was greatly enhanced by the presence of iron coatings on the sand medium and that suspended illite clay drastically enhanced oocyst deposition. Increasing ionic strength (up to a certain value) and decreasing pH also enhanced attachment efficiency.

Pathogenic behavior has been studied in numerous artificial recharge sites. Weiss et al. [80] calculated reductions from 2 to 4 log units for aerobic spores and 5.5 log units for total coliforms after travel distances between 24 and 177 m in three different riverbank filtration areas in the USA.

Bekele et al. [81] conducted a 39 month study on changes in water quality during infiltration through an unsaturated zone of 9 m and concluded that the elimination of microbial species was efficient: it detected adenovirus in only 6% and enteric viruses in 4% of the samples after 4.2 days.

Betancourt et al. [82] evaluated the elimination of enteric viruses in three artificial recharge facilities in the USA: one induced recharge (Colorado), an infiltration basin in Tucson, and another one in California. They concluded that enteric viruses group was below the bioanalytical method limit of detection after 5 days of transit time, and that residence time played a key role in the elimination of pathogens. The only infectious virus detected in the study was a reovirus. As it is difficult to associate reoviruses with a specific disease, they have not been paid much attention. However, they were found at higher concentrations than enteroviruses in treated and untreated wastewater, proving to be more resistant to UV disinfection than enteroviruses. Moreover, it appears that they also survive in water for long periods of time. According to these outcomes, reoviruses should be monitored in MAR system studies.

Elkayam et al. [83] studied several indicators at Shafdan (Israel), which recharges secondary effluent from a WWTP through 54 infiltration basins, covering an area of 270 ha in the coastal aquifer of Israel. The aquifer is formed of calcareous sandstones with intercalations of conglomerates, silts, and clay layers. Aquifer thickness is 180–200 m. The unsaturated area under the basins is 30–40 m deep. The wells are placed in two rings around the basins. The first ring extracts only recharged water from the rafts, whereas the second extracts also between 15% and 30% of water from the aquifer. The system has been operating on a large scale for more than 30 years. The study focused on indicators of bacteria, pathogenic viruses, coliphages, microbial source-tracking indicators (MST), and ARGs. The results showed a complete elimination of pathogenic viruses (enteroviruses, adenoviruses, noroviruses, parechoviruses), coliphages, and indicators of total and fecal bacteria and coliforms, fecal streptococcus, and bacteroides (MST) in the unsaturated area under the rafts. ARGs were detected in several wells, but they were also found in wells not impacted by the effluent, suggesting that these genes were related to the native microbial communities of the aquifer.

Beyond the actual mechanism controlling the fate of pathogens, it is important to take into account that (1) the soil is a living organism in itself, and (2) a broad range of pathogens, with different properties, may be present in treated wastewater or other source waters. The former implies that the soil system will be sensitive to external perturbations and its behavior will evolve in time. This, together with the specificities of every microorganism type, may explain the broad range of often contradictory sets of results reported in the literature.

The role of redox state has hardly been analyzed [84], which may reflect that most work is motivated by sand filters. Since we argue that a sequence of redox states improves the removal of CECs, it would also be desirable to understand the effect of redox state on the fate of pathogens. However, the indirect impact of redox variability becomes apparent. Bringing the water back into aerobic conditions favors the oxidation of iron (ferrous iron is mobile, while ferric iron tends to precipitate as goethite or its precursors). The positively charged surfaces of ferric oxides should retain pathogens that are mobile in other environments [84].

Thick unsaturated zones consistently lead to excellent removal rates of pathogens, and specifically viruses. While we attribute this success to air–water interfaces, which tend to retain colloids, other factors may be at play. The thickness of the unsaturated zone can only be controlled by pumping or, at the design stage, by site selection. Unfortunately, we did not have a thick unsaturated zone at our Palamós pilot site.

#### *3.3. Antibiotic-Resistant Bacteria (ARB) and Antibiotic-Resistant Genes (ARGs)*

Antibiotic resistance is a growing issue. Resistance is one of a number of consequences of the misuse and overuse of human and veterinary antibiotic drugs [85–88]. Antibiotics are substances that prevent the growth of bacteria. The term encompasses a wide range of pharmaceuticals with quite different physicochemical properties that are used in the treatment of bacterial infections, and also as prophylaxis for cattle and poultry. However, over the past decade, bacteria have been found to resist the drugs developed to suppress their growth and biological activity. In other words, bacterial strains are developing antibiotic resistance.

According to a recent report on the occurrence of antibiotic resistance in the USA [89], almost 3 million antibiotic-resistant infections are diagnosed every year in the USA, resulting in the deaths of more than 35,000 people. Globally, 10 million deaths per year caused by antibiotic-resistant bacterial infections are expected by 2050. Antibiotic-resistant bacteria, along with their resistance genes, are spread globally among foodstuffs, animals, plants, people, and the environment [90]. In this context, the World Health Organization (WHO) proposed a holistic action plan on microbial resistance involving humans and the environment [91], and the European Commission launched "One Health" as the European Action Plan to fight against antimicrobial resistance [92]. Additionally, the Watchlist for European Union monitoring defined in the Decisions 2018/840/EU (5 June 2018) included five antibiotics to be used to gather occurrence data to estimate the associated environmental risk of certain potentially hazardous compounds [93].

One of the reasons why ARGs cause great concern is because they are related to mobile genetic elements, and can therefore be easily transferred among microorganisms by horizontal gene transfer. This transfer can occur from bacteria, phages, free DNA, and dead cells to living cells [94].

WWTPs, in particular those treating urban sewage, have been recognized as one of the major receptors/sources of antibiotic resistance in the environment [88]. Even worse, ARGs may be enhanced in WWTPs [88]. The relationship between residual antibiotics in WWTPs and ARGs remains unclear. According to some studies, the presence of antibiotic residues during wastewater treatment may influence antibiotic resistance [87,95,96]. In contrast, other authors have pointed out that no correlation exists between the load of antibiotic residues in WWTPs and ARG abundance [97]. Nevertheless, a strong relationship between clinical and environmental antibiotic resistance has been reported. Temperature and humidity have been identified as two key factors controlling antibiotic resistance.

Since WWTPs are implicated as hotspots for the dissemination of antibiotic resistance into the environment and secondary effluents display higher relative abundance than the influents [98], it is important to assess whether ARGs could be reduced during soil aquifer treatment. Lack of a solid conceptual model for the fate of ARGs makes any proposal highly conjectural. However, some pieces of side evidence appear hopeful. The generation of anaerobic areas should help because the activity of microorganisms is reduced and the transfer of ARG is inhibited [99]. The presence of plants has proven efficient at removing ARGs in constructed wetlands [100]. The plant species selection deserves further research. In a recent study, the addition of biochar to soil resulted in notable changes in the microbial community, and these changes were different depending on the type of biochar used [101]. Changes in bacterial phylogenetic compositions can result in a change of ARGs. Therefore, the use of biochar as a component of reactive barriers might reduce ARGs. Nanotechnology may also be tested in MAR systems by including new nanomaterials in the infiltration pathway. To date, a few studies have pointed out the capability of selected nanomaterials to eliminate both ARB and ARGs [102].

#### *3.4. Public Acceptance of MAR*

A positive perception of MAR by the public is essential for its smooth implementation as a feasible and effective solution to increase water resources. The results of a public consultation about the use of treated wastewater in MAR operations were published recently [103]. Among the opposed respondents, the major concern was the lack of confidence in the wastewater treatment effectiveness before recharge. It is feared that if chemical and biological contaminants are recalcitrant to the treatments, recharge with WWTP effluents will lead to worsening groundwater quality. This kind of concern has led to broad negative public opinion on water reuse in general [104]. Negative perceptions may lead to failure [105–107].

In this regard, the implementation of operational strategies or designs to increase pathogen attenuation, chemical contaminant removal, and ARB and ARGs mitigation in MAR systems would facilitate public acceptance. If pilot studies are undertaken and results are successful, the MAR benefits supported by reliable scientific information should be stated and reasonably presented to the public. It is likely that the availability and understanding of this information will facilitate their positive perception and ultimately achieve public support.

Our concern, however, is that this may not suffice in times of "post-truth" and "fake news". Traditional approaches are needed: working with community organizations, promoting positive local media coverage of projects involving MAR, giving messages in clear non-technical language emphasizing its benefits and safety, offering public visits to the facilities, etc. Additional avenues of action that are specific to MAR include:


#### **4. Conclusions and Current**/**Future Challenges in MAR**

Our work supports the extensive literature body on water quality improvement during soil passage. Specifically, adding a reactive barrier improved the removal of chemical contaminants. Still, pathogen attenuation is significant (2–5 log units in our case), but was not particularly improved by the addition of the reactive barrier. Therefore, further improvement in the design of the reactive barriers and the operation of the system is needed.

We have discussed several options to enhance the degradation of recalcitrant chemical contaminants and the mitigation of pathogens. These include new compositions of the reactive barrier to broaden the types of sorption surfaces (biochar, zeolite, etc.), addition of metals to promote pathogen inactivation, implementation of thicker unsaturated zones to increase pathogen retention, and changes in the system operation to favor ferric oxide precipitation to create positively charged surfaces for further pathogen attenuation. We will test these approaches in the coming years at the Palamós pilot site and will assess the performance of the optimized system through the monitoring not only of CECs and pathogens, but also the development of antibiotic resistance, a serious emerging concern nowadays.

Public support must be achieved for the broad success of MAR. In the current context of climate change, where events of water scarcity and floods are occurring daily, improving water quality and increasing its quantity deserve determined action.

**Supplementary Materials:** Supplementary Materials are available online at http://www.mdpi.com/2073-4441/12/ 4/1012/s1.

**Author Contributions:** Conceptualization, C.V., J.C., L.M.-L., and M.S.D.-C.; Investigation and Methodology, C.V., L.M.-L., S.A., C.L., J.W. and M.S.D.-C.; Supervision, J.C. and M.S.D.-C. Writing—original draft, C.V., J.C. and M.S.D.-C.; Writing—review and editing, C.V., J.C., L.M.-L., S.A., C.L., J.W., and M.S.D.-C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Spanish Ministry of Science and Innovation CEX2018-000794-S), Water JPI (MARadentro-PCI2019-103603) and Catalan Water Agency (RESTORA-CA210/18/00040).

**Acknowledgments:** The authors acknowledge the Consorci Costa Brava for its constant support and the staff working at Palamós site for their unconditional help.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
