Can Borehole Heat Exchangers Trigger Cross-Contamination between Aquifers?

Borehole heat exchangers (BHEs) commonly reach depths of several tens of meters and cross different aquifers. Concerns have been raised about the possibility of boreholes to act as preferential pathways for contaminant transport among aquifers (cross-contamination). This article employs numerical modelling of contaminant transport in the subsurface to address these concerns. A common hydrogeological setup is simulated, composed of three layers: A shallow contaminated and a deep uncontaminated aquifer separated by an aquitard, all crossed by a permeable borehole. The hydraulic conductivity of the borehole and, to a lesser extent, the vertical hydraulic gradient between the aquifers are the key factors of cross-contamination. Results of the numerical simulations highlight that, despite the severe conditions hypothesized in our modelling study, the cross-contamination due to the borehole is negligible when filled with a slightly permeable material such as a geothermal grout properly mixed and injected. A good agreement was found with analytical formulas used for estimating the flow rate leaking through the borehole and for studying the propagation of leaked contaminant into the deep aquifer.


Introduction
Ground source heat pumps (GSHPs) gained increasing popularity in recent years due to the economic and environmental benefits they achieve as they are used to replace conventional techniques for heating and cooling of buildings and for domestic hot water production [1][2][3][4][5][6][7]. GSHPs are composed of a heat pump, which uses the underground as a heat source or sink, based on two different operating principles: open-loop systems, exchanging heat with groundwater abstracted by a well, and closed-loop systems, exchanging heat with the ground through the circulation of a heat carrier fluid in a pipe loop buried into the ground [8,9]. Closed-loop geothermal systems equipped with borehole heat exchangers (BHEs) are the most diffused type of GSHP thanks to the lower maintenance requirements and the possibility to install them even in the absence of an aquifer. BHEs are composed of one or two U-pipes installed in a borehole with a depth usually ranging between 50 and 200 m. At such depths, BHEs generally intercept different aquifers and, if incorrectly grouted during the completion, they tend to become preferential pathways through the less permeable layers (aquitards), thus making deep aquifer more vulnerable to contamination [10]. This concern led to restriction and bans of BHE installation in different protection areas of some European countries, as reported in a recent report by Prestor et al. [11]. sources of cross-contamination. Natural aquitard discontinuities are permeable lenses, whereas manmade discontinuities are generally identified with wells, both improperly cemented and abandoned (and, hence, likely to collapse [13]). Another major issue induced by the presence of wells is related to the leakage, through the permeable well annulus, between two adjacent aquifers with different hydraulic heads. Most studies available in the literature [14][15][16][17][18][19] addressed wells as the main anthropogenic discontinuity in aquitards. The study by Bucci et al. [20] is, to our knowledge, the only one addressing BHEs, providing a qualitative description of voids formed along an experimental grout-filled borehole, drilled and filled on the escarpment of a gravel pit and subsequently uncovered by excavation.
The aforementioned works highlight the key role of the hydraulic conductivity of the borehole filling in preventing the leakage among aquifers. In the last 30 years, specific cement-bentonite mixtures (geothermal grouts) were developed to seek a compromise between hydraulic insulation, strength, and thermal conductivity [21]. The ongoing research on geothermal grouts is focused on increasing the thermal conductivity, e.g., with the addition of graphite flakes [22][23][24], and reducing the hydraulic conductivity, generally with the addition of bentonite. Figure 1 from Casasso and Sethi [10] compares the experimental results of Allan and Philippacopoulos [25], Park et al. [26] and Indacoechea-Vega et al. [27], who tested grout specimens in different conditions. Hydraulic conductivities range between 10 −12 m/s, far below the conductivities of aquitards and of nearly all aquicludes, and 10 −6 m/s, i.e., a value typical of fine sand. Therefore, geothermal grouts are generally less permeable than the majority of clayey and silty aquicludes. Figure 1 also highlights that tests conducted on grout-pipe specimens [25][26][27] showed higher hydraulic conductivities compared to grout-only specimens. This is due to the fact that the interfaces between pipes and grout are prone to form discontinuities which act as preferential pathways for water [25,26]. In addition, the hydraulic conductivity also increases when the water/cement ratio is too high [27] and/or the grout is subjected to extreme thermal (freeze-thaw) [25,27] or hydraulic (wet-dry) [25] alternated stresses. Previous studies on aquifer cross-contamination only focused on the quantification of flow among aquifers. To our knowledge, the contaminant propagation due to this flow was only addressed by Lacombe et al. [28]. The authors performed numerical solute transport simulations modelling the borehole crossing an aquitard, with two different approaches. The borehole was modelled with two approaches, namely as a column with the same hydraulic conductivity of the aquifer (i.e., with the Darcy model) to reproduce the case of a collapsed well and as a conduit (i.e., with the Hagen-Poiseuille model) to simulate the presence of fissures along the borehole. However, Previous studies on aquifer cross-contamination only focused on the quantification of flow among aquifers. To our knowledge, the contaminant propagation due to this flow was only addressed by Lacombe et al. [28]. The authors performed numerical solute transport simulations modelling the borehole crossing an aquitard, with two different approaches. The borehole was modelled with two approaches, namely as a column with the same hydraulic conductivity of the aquifer (i.e., with the Darcy model) to reproduce the case of a collapsed well and as a conduit (i.e., with the Hagen-Poiseuille model) to simulate the presence of fissures along the borehole. However, Lacombe et al. [28] simulated operating conditions that were far more critical than those induced by the presence of BHEs.
The purpose of this paper is to quantitatively evaluate the cross-contamination induced by the presence of an unproperly sealed BHE completed into a dual aquifer system separated by an aquitard. The case of a shallow contaminated and a deep clean aquifer is considered, as it is the most common one [28]. A numerical sensitivity analysis is performed with the finite element numerical software FEFLOW to assess how the aquitard thickness and the backfilling hydraulic conductivity affect the contaminant leakage from the upper aquifer to the lower one.
This paper is organized as follows. Section 2 delineates the conceptual model and describes the most relevant parameters involved in the study. Section 3 presents and discusses the results: The leaking flow rates, the concentration distributions in the deep aquifer and the possible use of analytical solutions to approximate them. Conclusions (Section 4) summarize methods and findings and identify possible fields of future research activity.

Conceptual Model
The hydrogeological scenario hypothesized for this work is illustrated in Figure 2: A shallow aquifer (with a hydraulic conductivity K 1 and contaminated with a concentration C 0 of a certain solute) is separated from the deep aquifer (with a hydraulic conductivity K 2 and a zero-value concentration of the solute) by an aquitard (i.e., a layer with conductivity K K 1 , K 2 ). A borehole heat exchanger was added, modeled as a cylinder with a diameter d, hydraulic conductivity (K f ill ), crossing all model layers. If the borehole is permeable than the aquitard (K f ill > K ), a preferential flow is generated, that depends on the hydraulic head difference (∆h = h 1 − h 2 ) between the aquifers crossed. As highlighted in Santi's paper [12], an unfavorable situation occurs if the hydraulic head of the polluted aquifer is higher than the clean aquifer one (h 1 > h 2 ). In this case, contaminants flow vertically through the aquitard, through the permeable borehole and propagate in the deep aquifer. This work assesses the deep aquifer incremental contribution to contamination due to the preferential flow through the borehole, i.e., the additional concentration in the deep aquifer, compared to the case in which the borehole has the same hydraulic conductivity of the aquitard (K f ill = K ). A sensitivity analysis is then performed to assess the role of K f ill on the leakage. The aquitard thickness (b ) is also considered, since it influences both the preferential leakage through the borehole and the contaminant transport through the rest of the aquitard. Lacombe et al. [28] simulated operating conditions that were far more critical than those induced by the presence of BHEs. The purpose of this paper is to quantitatively evaluate the cross-contamination induced by the presence of an unproperly sealed BHE completed into a dual aquifer system separated by an aquitard. The case of a shallow contaminated and a deep clean aquifer is considered, as it is the most common one [28]. A numerical sensitivity analysis is performed with the finite element numerical software FEFLOW to assess how the aquitard thickness and the backfilling hydraulic conductivity affect the contaminant leakage from the upper aquifer to the lower one.
This paper is organized as follows. Section 2 delineates the conceptual model and describes the most relevant parameters involved in the study. Section 3 presents and discusses the results: The leaking flow rates, the concentration distributions in the deep aquifer and the possible use of analytical solutions to approximate them. Conclusions (Section 4) summarize methods and findings and identify possible fields of future research activity.

Conceptual Model
The hydrogeological scenario hypothesized for this work is illustrated in Figure 2: A shallow aquifer (with a hydraulic conductivity 1 and contaminated with a concentration 0 of a certain solute) is separated from the deep aquifer (with a hydraulic conductivity 2 and a zero-value concentration of the solute) by an aquitard (i.e., a layer with conductivity ′ ≪ 1 , 2 ). A borehole heat exchanger was added, modeled as a cylinder with a diameter , hydraulic conductivity ( ), crossing all model layers. If the borehole is permeable than the aquitard ( > ′), a preferential flow is generated, that depends on the hydraulic head difference (Δℎ = ℎ 1 − ℎ 2 ) between the aquifers crossed. As highlighted in Santi's paper [12], an unfavorable situation occurs if the hydraulic head of the polluted aquifer is higher than the clean aquifer one (ℎ 1 > ℎ 2 ). In this case, contaminants flow vertically through the aquitard, through the permeable borehole and propagate in the deep aquifer. This work assesses the deep aquifer incremental contribution to contamination due to the preferential flow through the borehole, i.e., the additional concentration in the deep aquifer, compared to the case in which the borehole has the same hydraulic conductivity of the aquitard ( = ′). A sensitivity analysis is then performed to assess the role of on the leakage. The aquitard thickness ( ′) is also considered, since it influences both the preferential leakage through the borehole and the contaminant transport through the rest of the aquitard.

Numerical Model
The scenario described in Section 2.1 was modelled with the software FEFLOW 7.1 by DHI WASY [29]. FEFLOW solves the flow and transport equations in saturated porous media with a finite-

Numerical Model
The scenario described in Section 2.1 was modelled with the software FEFLOW 7.1 by DHI WASY [29]. FEFLOW solves the flow and transport equations in saturated porous media with Water 2020, 12, 1174 4 of 11 a finite-element approach including fluid density and viscosity changes (fully coupled flow and transport modelling). Figure 3 reports the plan view (A) and the cross section (B) of the modelling domain. A large domain size was chosen (500 m × 500 m × 150 m of depth) to avoid boundary effects. The domain was divided into 1,870,498 elements and 959,176 nodes with a strong refinement along the flow direction (aligned with the x axis) in correspondence of the borehole ( Figure 3A). The vertical discretization presents 46 layers (47 slices) with a variable vertical resolution represented in Figure 3B and hereby described. The upper aquifer is 35-m thick, and it is composed of 7 layers with a constant thickness of 5 m. The aquitard is composed of 1-m thick layers in a number that depends on the prescribed thickness b . The lower aquifer is divided into 5-m thick layers up to a depth of 150 m from ground surface, with a further refinement of 1-m deep layers between the top of the aquitard (35 m from ground surface) and a depth of 55 m from ground surface, i.e., the maximum depth reached by the aquitard bottom in the case b = 20 m. In this way, although the thickness of the deep aquifer can vary from 95 m (b = 20 m) to 111 m (b = 4 m), the domain mesh is exactly the same, i.e., with the same number of layers and nodes, for all the simulations run. An additional simulation was run with a double vertical resolution (i.e., with layers 2.5 m thick for the aquifers and 0.5 m thick for the aquitard) to assess whether a mesh refinement would change results. As reported in Section S3 of the Supplementary Materials, only minor changes are observed for the time series of concentrations in the deep aquifer.
The regional groundwater flow in both aquifers was simulated by assigning constant hydraulic head values on the western and eastern domain borders, with a 1% gradient on both aquifers as shown in Figure 3. A strong hydraulic head difference (∆h = 10 m) was assigned between the shallow contaminated aquifer and the deep uncontaminated one. From a cross-contamination point of view, this represents a severe, yet possible scenario as reported in Refs. [12,[30][31][32][33] and in Section S4 of the Supplementary Materials. Values of the aquitard thickness (b ) between 4 m and 20 m were assigned in the simulations since they are representative of the stratigraphy of alluvial plains [30,33].
Water 2019, 11, x FOR PEER REVIEW 4 of 11 element approach including fluid density and viscosity changes (fully coupled flow and transport modelling). Figure 3 reports the plan view (A) and the cross section (B) of the modelling domain. A large domain size was chosen (500 m × 500 m × 150 m of depth) to avoid boundary effects. The domain was divided into 1,870,498 elements and 959,176 nodes with a strong refinement along the flow direction (aligned with the x axis) in correspondence of the borehole ( Figure 3A). The vertical discretization presents 46 layers (47 slices) with a variable vertical resolution represented in Figure 3B and hereby described. The upper aquifer is 35-m thick, and it is composed of 7 layers with a constant thickness of 5 m. The aquitard is composed of 1-m thick layers in a number that depends on the prescribed thickness ′. The lower aquifer is divided into 5-m thick layers up to a depth of 150 m from ground surface, with a further refinement of 1-m deep layers between the top of the aquitard (35 m from ground surface) and a depth of 55 m from ground surface, i.e., the maximum depth reached by the aquitard bottom in the case ′ = 20 m. In this way, although the thickness of the deep aquifer can vary from 95 m ( ′ = 20 m) to 111 m ( ′ = 4 m), the domain mesh is exactly the same, i.e., with the same number of layers and nodes, for all the simulations run. An additional simulation was run with a double vertical resolution (i.e., with layers 2.5 m thick for the aquifers and 0.5 m thick for the aquitard) to assess whether a mesh refinement would change results. As reported in Section S3 of the Supplementary Materials, only minor changes are observed for the time series of concentrations in the deep aquifer.
The regional groundwater flow in both aquifers was simulated by assigning constant hydraulic head values on the western and eastern domain borders, with a 1% gradient on both aquifers as shown in Figure 3. A strong hydraulic head difference (Δℎ = 10 m) was assigned between the shallow contaminated aquifer and the deep uncontaminated one. From a cross-contamination point of view, this represents a severe, yet possible scenario as reported in Refs. [12,[30][31][32][33] and in Section 4 of the Supplementary Materials. Values of the aquitard thickness ( ′) between 4 m and 20 m were assigned in the simulations since they are representative of the stratigraphy of alluvial plains [30,33]. The borehole crosses all the layers and has a radius = 0.075 m, i.e., typical values of borehole heat exchangers. The cross-section of the borehole was considered as homogeneous, and a unique value of hydraulic conductivity ( ) ranging from 10 −9 m/s to 10 −2 m/s was assigned to it. This range also includes values much higher than those of geothermal grouts and typical of sand and gravel ( Figure 1). Indeed, this work also aims at assessing the order of magnitude of values from which an appraisable cross-contamination may be induced by a BHE. The BHE cross section was The borehole crosses all the layers and has a radius r a = 0.075 m, i.e., typical values of borehole heat exchangers. The cross-section of the borehole was considered as homogeneous, and a unique value of hydraulic conductivity (K f ill ) ranging from 10 −9 m/s to 10 −2 m/s was assigned to it. This range also includes K f ill values much higher than those of geothermal grouts and typical of sand and gravel ( Figure 1). Indeed, this work also aims at assessing the order of magnitude of K f ill values from which an appraisable cross-contamination may be induced by a BHE. The BHE cross section was treated as a homogeneous zone, i.e., without differentiating between grout and pipes. This choice is justified by the fact that, as mentioned in Section 1, tests conducted on geothermal grout specimens are generally performed on representative BHE sections (i.e., including pipes), rather than on grout-only samples. In addition, including BHE pipes would have required a much denser mesh with a consequent increase of computational effort.
Homogeneous porosity values ε = 0.3 (total) and n e = 0.2 (effective) were assigned all over the modelling domain, whereas the hydraulic conductivity was set to K 1 = 10 −3 m/s in the shallow aquifer, K 2 = 10 −4 m/s in the deep aquifer, and K = 10 −9 m/s in the aquitard. Two additional simulations were run setting higher values of the hydraulic conductivity of the deep aquifer (K 2 = 10 −3 m/s and K 2 = 10 −2 m/s). Results are reported in Section S1 of the Supplementary Materials.
The initial contaminant concentration was imposed at 100 mg/L in all nodes of the unconfined aquifer, thus assuming an entirely and uniformly polluted upper aquifer. Similarly, the initial and the boundary concentrations in the lower aquifer were forced to 0 mg/L, thus simulating that the only contamination occurring in the deep aquifer comes from the shallow aquifer. The contaminant was set as conservative, with no sorption nor decay.
The longitudinal dispersivity was set to α L = 5 m for all simulations reported in this paper, and the transverse dispersivity was set to α T = 0.1α L = 0.5 m. Three simulations were run with lower dispersivity values (α L = 0.5 m, 1 m, 2 m and α T = 0.1α L ) and results are reported in Section S2 of the Supplementary Materials. Table 1 reports the complete parameter set, including initial and boundary conditions of the numerical model.

Results
A total number of 40 simulations, combining 8 hydraulic conductivities K f ill and 5 different aquitard depths b , were performed using the modelling scheme described in Section 2. Results of the modelling activities are presented in this Section and are divided as follows. The leakage flow rate through the borehole is analyzed in Section 3.1, and a correspondence is found with available analytical formulas. Section 3.2 presents the time trends of concentrations in the deep aquifer downstream the borehole. Section 3.3 analyzes the difference between the concentrations of each scenario and reference cases with a borehole as permeable as the aquitard (K f ill = K ): predictably, the concentration increments in the deep aquifer are proportional to the leakage flow rates. Therefore, in Section 3.4 the numerical results are compared to the analytical solution of a continuous point source release in 3D geometry (Hunt,[34]).

Evaluation of Borehole Leakage Rates
As reported in Section 1, most studies on aquifer cross-contamination have focused on the quantification of the leakage flow rate (Q l ) of water crossing an aquitard permeable zone (well, borehole, permeable lens). Analytical formulas were developed to calculate the flow rate with different degrees of complexity [15][16][17][18][19], the simplest one being an application of the Darcy law: where r a is the radius of the borehole (m), and b is the aquitard thickness (m). The Darcy formula assumes a constant and uniform gradient i = ∆h/b along the borehole section across the aquitard, with ∆h remaining equal to the hydraulic head difference of the two aquifers. On the other hand, Bonte et al. [19] recently proposed an analytical formula that takes into account the local alteration of hydraulic heads around the borehole, which changes the value of ∆h between the inlet and the outlet of the borehole crossing the aquitard. The flow rate value calculated with the Darcy formula (Q l,Darcy ) is modified by a correction coefficient: where K 1 and K 2 are respectively the hydraulic conductivities (m/s) of the shallow and deep aquifer.
Since the parameters K f ill , D, r a , K 1 and K 2 are positive, it follows that Q l,Bonte < Q l,Darcy . Figure 4A reports a comparison between these analytical formulas with the setup adopted in the numerical models (K 1 = 10 −3 m/s, K 2 = 10 −4 m/s, D = 10 m, K f ill = 10 −9 to 10 −2 m/s). An appraisable difference between the two models is observed only for K f ill > 10 −4 m/s, i.e., with a very high leakage flow rate. In Figure 4B, leakage flow rates resulting from FEFLOW simulations are compared with those calculated with the analytical formulas above described, i.e., Equations (1) and (2). The Darcy equation has an excellent agreement with the numerical model results, whereas the Bonte equation has a slightly worse fit for higher leakage rates as it underestimates the leakage flow.
Water 2019, 11, x FOR PEER REVIEW 6 of 11

Evaluation of Borehole Leakage Rates
As reported in Section 1, most studies on aquifer cross-contamination have focused on the quantification of the leakage flow rate ( ) of water crossing an aquitard permeable zone (well, borehole, permeable lens). Analytical formulas were developed to calculate the flow rate with different degrees of complexity [15][16][17][18][19], the simplest one being an application of the Darcy law: where is the radius of the borehole (m), and ′ is the aquitard thickness (m). The Darcy formula assumes a constant and uniform gradient = ℎ/ ′ along the borehole section across the aquitard, with Δℎ remaining equal to the hydraulic head difference of the two aquifers. On the other hand, Bonte et al. [19] recently proposed an analytical formula that takes into account the local alteration of hydraulic heads around the borehole, which changes the value of Δℎ between the inlet and the outlet of the borehole crossing the aquitard. The flow rate value calculated with the Darcy formula ( , ) is modified by a correction coefficient: where 1 and 2 are respectively the hydraulic conductivities (m/s) of the shallow and deep aquifer. Since the parameters , , , 1 and 2 are positive, it follows that , < , . Figure 4A reports a comparison between these analytical formulas with the setup adopted in the numerical models ( 1 = 10 −3 m/s , 2 = 10 −4 m/s , = 10 m , = 10 −9 10 −2 m/s ). An appraisable difference between the two models is observed only for > 10 −4 m/s, i.e., with a very high leakage flow rate. In Figure 4B, leakage flow rates resulting from FEFLOW simulations are compared with those calculated with the analytical formulas above described, i.e., Equations (1) and (2). The Darcy equation has an excellent agreement with the numerical model results, whereas the Bonte equation has a slightly worse fit for higher leakage rates as it underestimates the leakage flow.

Propagation of Contaminants in the Deep Aquifer
The contaminant propagation from the shallow to the deep aquifer occurs through the aquitard through two different pathways: the aquitard itself and the borehole which, in the case it has a higher hydraulic conductivity ( > ′), represents a preferential contamination pathway. The preferential contaminant migration through the borehole was assessed by simulating the effect of different values of ranging from 10 −9 m/s to 10 −2 m/s. The reference case with = ′ = 10 −9 m/s provides as a result the contaminant distribution in the deep aquifer, observable even in the absence of any

Propagation of Contaminants in the Deep Aquifer
The contaminant propagation from the shallow to the deep aquifer occurs through the aquitard through two different pathways: the aquitard itself and the borehole which, in the case it has a higher hydraulic conductivity (K f ill > K ), represents a preferential contamination pathway. The preferential Water 2020, 12, 1174 7 of 11 contaminant migration through the borehole was assessed by simulating the effect of different values of K f ill ranging from 10 −9 m/s to 10 −2 m/s. The reference case with K f ill = K = 10 −9 m/s provides as a result the contaminant distribution in the deep aquifer, observable even in the absence of any borehole. Figure 5A shows the influence of the aquitard thickness b on the time trends of the "base contamination" at the top of the aquifer. As expected, concentrations propagate more slowly as the aquitard thickness increases, due to both the increase of the preferential pathway length (i.e., b ) and the reduction of the hydraulic gradient (i.e., ∆h/b ). Figure 5B shows the contaminant concentration trends 200 m downstream the borehole outlet for different values of K f ill . As K f ill increases, the deep aquifer reaches higher concentrations and the arrival time of the contaminant is faster, due to its rapid propagation through the permeable borehole. borehole. Figure 5A shows the influence of the aquitard thickness ′ on the time trends of the "base contamination" at the top of the aquifer. As expected, concentrations propagate more slowly as the aquitard thickness increases, due to both the increase of the preferential pathway length (i.e., ′) and the reduction of the hydraulic gradient (i.e., Δℎ/ ′). Figure 5B shows the contaminant concentration trends 200 m downstream the borehole outlet for different values of . As increases, the deep aquifer reaches higher concentrations and the arrival time of the contaminant is faster, due to its rapid propagation through the permeable borehole.

Cross-Contamination Effect of the Borehole
The contribution of a permeable borehole to the contamination of the deep aquifer was identified as the increment of concentration (Δ ) in the cases with > ′, compared to the respective reference cases with = ′ and with the same aquitard thickness ′: Δ ( , , , , ) = ( , , , , > ′) − ( , , , , = ′ ) (3) Figure 6 shows the variation of the maximum values of Δ / 0 at 200 m from the borehole, depending on the aquitard thickness ′ and the borehole conductivity . The range of from 10 −6 m/s to 10 −2 m/s was selected to show the cases in which an appraisable impact is exerted by the contaminant transport through the borehole. In the severe scenario hypothesized, the contamination brought by the preferential pathway of the borehole falls below 1% of the concentration in the shallow aquifer (i.e., Δ < 0.01 0 ) for a value of = 10 −4 m/s , i.e., over 100 times higher than the maximum measured for geothermal grouts (see Figure 1). These results confirm that geothermal grouts available in commerce provide enough protection against cross-contamination.

Cross-Contamination Effect of the Borehole
The contribution of a permeable borehole to the contamination of the deep aquifer was identified as the increment of concentration (∆C) in the cases with K f ill > K , compared to the respective reference cases with K f ill = K and with the same aquitard thickness b : ∆C x, y, z, t, K f ill = C x, y, z, t, K f ill > K − C x, y, z, t, K f ill = K (3) Figure 6 shows the variation of the maximum values of ∆C/C 0 at 200 m from the borehole, depending on the aquitard thickness b and the borehole conductivity K f ill . The range of K f ill from 10 −6 m/s to 10 −2 m/s was selected to show the cases in which an appraisable impact is exerted by the contaminant transport through the borehole. In the severe scenario hypothesized, the contamination brought by the preferential pathway of the borehole falls below 1% of the concentration in the shallow aquifer (i.e., ∆C < 0.01C 0 ) for a value of K f ill = 10 −4 m/s, i.e., over 100 times higher than the maximum measured for geothermal grouts (see Figure 1). These results confirm that geothermal grouts available in commerce provide enough protection against cross-contamination. contaminant transport through the borehole. In the severe scenario hypothesized, the contamination brought by the preferential pathway of the borehole falls below 1% of the concentration in the shallow aquifer (i.e., Δ < 0.01 0 ) for a value of = 10 −4 m/s , i.e., over 100 times higher than the maximum measured for geothermal grouts (see Figure 1). These results confirm that geothermal grouts available in commerce provide enough protection against cross-contamination.

Analytical Modelling of the Propagation of Contaminants in the Deep Aquifer
The borehole section at the bottom of the aquitard acts as a continuous point source for the release of contaminants in the deep aquifer. After leaking in the aquifer, the pollutant migrates along the three directions generating a contamination plume. However, since the aquitard acts as a virtually impermeable limit on top of the domain, the pollutant propagation along the z axis is allowed only in the downward direction. A closed form analytical solution for contaminant transport in a 3D semi-infinite aquifer was provided by Hunt [34]: with B (m) and U (unitless): where Q l (m 3 s −1 ) is the leakage flow rate through the borehole, here calculated using the Darcy law; C 0 (kg m −3 ) is the contaminant concentration in the leaking water; D x = α L v e and D y = D z = α T v e are the longitudinal and transverse dispersion coefficients (m 2 s −1 ), which are proportional to the respective values of dispersivity α L and α T (m); R is the retardation coefficient (unitless), which is equal to 1 as no sorption was considered; and λ (s −1 ) is the degradation coefficient set equal to 0 s −1 (no degradation). Hunt's solution was compared with the concentrations resulting at the end of FEFLOW simulations, performed for different combinations of the aquitard hydraulic conductivity and thickness. A scatterplot of the results is shown in Figure 7 in terms of potential concentration increment at the top of the deep aquifer and at different distances downstream the borehole outlet. A good agreement was found for all simulations, except for a few ones at very low values of ∆C/C 0 (due to precision issues of the numerical model) and at very high values (due to the local perturbation of the flow field, which is not considered by Hunt's model). Such result suggests that Hunt's solution represents a useful mean to perform a preliminary estimation of the potential environmental impact of borehole installation in multi-layer aquifers. thickness. A scatterplot of the results is shown in Figure 7 in terms of potential concentration increment at the top of the deep aquifer and at different distances downstream the borehole outlet. A good agreement was found for all simulations, except for a few ones at very low values of Δ / 0 (due to precision issues of the numerical model) and at very high values (due to the local perturbation of the flow field, which is not considered by Hunt's model). Such result suggests that Hunt's solution represents a useful mean to perform a preliminary estimation of the potential environmental impact of borehole installation in multi-layer aquifers.

Conclusions
This article addressed the concerns about BHEs cross-contamination by performing a numerical modelling study with the finite-element flow and solute transport code FEFLOW. The setup simulated for this study is a shallow aquifer with a contaminant concentration C 0 , separated from a deep aquifer by an aquitard with thickness b . The substantial hydraulic head difference (∆h = 10 m) between the two aquifers makes the hypothesized setup a very severe simulation scenario. The three layers are crossed by a borehole with a hydraulic conductivity K f ill , ranging from 10 −9 m/s (equal to the aquitard) to 10 −2 m/s (typical of a coarse sand). The sensitivity analysis focused on the influence of the parameters b and K f ill on the leakage from the shallow to the deep aquifer. Two contributions were assessed, i.e., the preferential migration through the borehole-which occurs if the borehole is more permeable than the aquitard crossed (K f ill > K )-and the transport through the aquitard, which would occur even in the absence of the borehole.
The main findings of our study are hereby summarized: The flow rate crossing the borehole (Q l ) is well described by the Darcy's law (Equation (1)), i.e., Q l is proportional to the hydraulic head difference between the aquifers crossed (∆h) and to the hydraulic conductivity of the borehole filling (K f ill ); The concentrations in the deep aquifer at a certain point are well correlated with the flow rate leaking through the borehole (Q l ); The spatial distribution of concentration at the top of the deep aquifer is well reproduced by Hunt's analytical formula [34] reported in Equation (4); In the severe scenario hypothesized in our modelling study, the incremental concentration in the deep aquifer due to the preferential flow through a borehole is below 1% of the contaminant concentration in the shallow aquifer (C 0 ) if the borehole filling is as permeable as fine sand (K f ill = 10 −4 m/s). Laboratory experiments on grout specimens reported values of K f ill between 10 −11 m/s and 10 −6 m/s (see Figure 1).
The results of the study presented in this paper provide insights on the possible cross-contamination induced by borehole heat exchangers. Although, to our knowledge, this occurrence is not reported in the literature, this possibility is worth being evaluated in the proximity of water wells. Based on our results, appraisable cross-contamination effects are observed only if a very permeable filling material is used, e.g., coarse sand or gravel. Conversely, pre-mixed geothermal grouts available in the market, if correctly mixed and injected, provide enough guarantees against cross-contamination even in severe conditions. Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/4/1174/s1. Figure S1. Effect of the hydraulic conductivity of the deep aquifer on the contaminant concentrations at the top of the deep aquifer, 200 m downstream the borehole outlet; Figure S2. Effect of the contaminant longitudinal dispersivity value α L . The transverse dispersivity is set to to α T = 0.1α L for each case; Figure S3. Effect of vertical resolution: comparison of the time series of contaminant concentrations at different distances downstream the borehole outlet (40, 80, 120, 160 and 200 m) with the "standard" resolution (layer thickness of 5 m for aquifers and 1 m for the aquitard) and the doubled resolution (layer thickness of 2.5 m for aquifers and 0.5 m for the aquitard); Figure S4. Position of Ciriè, the location of the two monitored wells examined, in North-Western Italy; Figure S5. Monitoring wells in Ciriè, screened in the shallow aquifer (red line) and in the deep aquifer (light blue line).
Author Contributions: Conceptualization, methodology, software, writing-original draft preparation, writing-review and editing, A.C., N.F., S.P., C.B., R.S. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.