Thermal Impact Assessment of Groundwater Heat Pumps (GWHPs): Rigorous vs. Simplified Models

Groundwater Heat Pumps (GWHPs) are increasingly adopted for air conditioning in urban areas, thus reducing CO2 emissions, and this growth needs to be managed to ensure the sustainability of the thermal alteration of aquifers. However, few studies have addressed the propagation of thermal plumes from open-loop geothermal systems from a long-term perspective. We provide a comprehensive sensitivity analysis, performed with numerical finite-element simulations, to assess how the size of the thermally affected zone is driven by hydrodynamic and thermal subsurface properties, the vadose zone and aquifer thickness, and plant setup. In particular, we focus the analysis on the length and width of thermal plumes, and on their time evolution. Numerical simulations are compared with two simplified methods, namely (i) replacing the time-varying thermal load with its yearly average and (ii) analytical formulae for advective heat transport in the aquifer. The former proves acceptable for the assessment of plume length, while the latter can be used to estimate the width of the thermally affected zone. The results highlight the strong influence of groundwater velocity on the plume size and, especially for its long-term evolution, of ground thermal properties and of subsurface geometrical parameters.


Introduction
Ground Source Heat Pumps (GSHPs) utilize the ground as a heat source or sink, respectively, for the heating or cooling of buildings. At moderate depths (5 to 20 m), the ground has an almost constant temperature, with a value close to the yearly average air temperature [1]; hence, GSHPs have a higher energy efficiency (coefficient of performance, COP) compared to Air Source Heat Pumps (ASHPs) [2], and their potential for the reduction of the carbon intensity of building air conditioning is widely acknowledged [3].
The energy efficiency of a GSHP strongly depends on site-specific properties of the ground, for both Borehole Heat Exchangers (BHEs) [4][5][6][7][8] and Ground Water Heat Pumps (GWHPs) [9][10][11][12]. GWHPs exchange heat with groundwater extracted through one or multiple water wells, and are increasingly employed as an air conditioning system, especially for large commercial and public buildings [13][14][15][16]. Groundwater is usually reinjected through wells into the same aquifer to avoid its depletion [17], but this leads to the formation of thermally altered zones, called thermal plumes. Thermal plumes pose two sustainability issues: the upstream propagation of the plume, which can reach the abstraction well (thermal short-circuit) [9,18,19], thus reducing the COP of the heat pump, and the downstream propagation of the plume, which can impair drinking water wells or other GWHPs. Assessing the subsurface thermal impact of open-loop geothermal systems is therefore essential for their design.

Flow and Heat-Transport in Porous Media
Heat transport in saturated porous media occurs by conduction (driven by temperature gradient), advection (due to the flow of the fluid phase), and dispersion (induced by the heterogeneity of the groundwater velocity field).
These mechanisms are described by the heat conservation equation in a two-phase (solid and liquid) medium [34], with the assumption of a groundwater flow aligned with the x-axis: where ρc and ρ w c w are the thermal capacities (J·m −3 ·K −1 ) of the porous medium and of the liquid phase, respectively; T is the temperature (K), considered equal for the solid and the liquid phases (i.e., the thermal equilibrium is considered instantaneous); v D is the Darcy velocity (m·s −1 ); α L and α T are, respectively, the longitudinal and transverse dispersivities (m) relative to the groundwater flow direction; λ is the thermal conductivity of the porous medium (W·m −1 ·K −1 ); and H is the heat source/sink (W·m −3 ). The first term of Equation (1) is the temperature variation over time, which depends on the volume-averaged (bulk) heat capacity of the porous medium (ρc): ρc = n e ρ w c w + (1 − n e )ρ s c s (2) where ρ s c s is the thermal capacity (J·m −3 ·K −1 ) of the solid phase. The second term describes the advection, which is function of the Darcy velocity (v D ): where K is the hydraulic conductivity (m·s −1 ), h is the hydraulic head (m), and i is the hydraulic gradient along the x direction (dimensionless). Conduction depends on the bulk thermal conductivity in the third term of Equation (1): where λ s and λ w are the thermal conductivities (W·m −1 ·K −1 ) of, respectively, the solid matrix and the groundwater, and n e is the effective porosity (dimensionless). Thermal dispersion is described in Equation (1) by the corresponding coefficients α L and α T . Dividing Equation (1) by ρc, we get: where v e is the effective velocity of groundwater through the pores, calculated as: Since heat is exchanged between the fluid and the solid phase, the advective velocity of the thermal plume is lower than the groundwater flow velocity [19]. This phenomenon is similar to solute sorption in porous media with linear kinetics. By analogy, it is therefore possible to define the thermal retardation coefficient of Equation (5) as the ratio between the effective velocity of groundwater (v e ) and heat front (v th ): (1 − n e )ρ s c s + n e ρ w c w n e ρ w c w = ρc n e ρ w c w (7) D x , D y , and D z are the longitudinal and transversal thermal dispersion coefficients (m 2 /s):

Model Setup
The objective of the first part of this work is to evaluate and compare the influence of thermal and hydrogeological subsurface properties (i.e., hydraulic conductivity of the aquifer, K; hydraulic gradient, i; effective porosity, n e ; aquifer and vadose zone thickness, b and d, respectively; thermal capacity ρc), as well as plant parameters (i.e., well spacing L), on the time evolution and spatial distribution of the thermal plume originated by a GWHP.
A simple but representative plant configuration was therefore adopted for the simulations, i.e., a well doublet in an unconfined aquifer, where the wells are aligned with the groundwater flow direction and the reinjection well is downstream of the abstraction well. The open-loop geothermal system was modelled with the finite-element code FEFLOW ® 6.2 (DHI-WASY, Berlin, Germany) [34], performing transient coupled flow and heat-transport simulations over an operating lifetime of 50 years. The default dimensions of the 3D modelling domain are 6000 m × 3000 m × 75 m, and variations were introduced only in the mesh thickness. The default vertical discretization is of 15 horizontal layers (16 slices), yielding a triangular mesh of 461,850 elements and 252,144 nodes for the default well configuration ( Figure S1 in the Supplementary Materials), with slight variations for the others. A mesh density convergence study was performed ( Figure S2 in the Supplementary Materials. The layers crossing the vadose zone, the aquifer, and the upper part of the underlying aquitard (layers 1 to 13) have a constant thickness of 2.5 m. The bottom two layers are thicker (10 m and 30 m, respectively), since this portion is deemed to be negligibly affected by the GWHP.
A simplified 2D cross-section of the model is shown in Figure 1.

Model Setup
The objective of the first part of this work is to evaluate and compare the influence of thermal and hydrogeological subsurface properties (i.e., hydraulic conductivity of the aquifer, ; hydraulic gradient, ; effective porosity, ; aquifer and vadose zone thickness, and , respectively; thermal capacity ), as well as plant parameters (i.e., well spacing ), on the time evolution and spatial distribution of the thermal plume originated by a GWHP.
A simple but representative plant configuration was therefore adopted for the simulations, i.e., a well doublet in an unconfined aquifer, where the wells are aligned with the groundwater flow direction and the reinjection well is downstream of the abstraction well. The open-loop geothermal system was modelled with the finite-element code FEFLOW ® 6.2 (DHI-WASY, Berlin, Germany) [34], performing transient coupled flow and heat-transport simulations over an operating lifetime of 50 years. The default dimensions of the 3D modelling domain are 6000 m × 3000 m × 75 m, and variations were introduced only in the mesh thickness. The default vertical discretization is of 15 horizontal layers (16 slices), yielding a triangular mesh of 461,850 elements and 252,144 nodes for the default well configuration ( Figure S1 in the Supplementary Materials), with slight variations for the others. A mesh density convergence study was performed ( Figure S2 in the Supplementary Materials. The layers crossing the vadose zone, the aquifer, and the upper part of the underlying aquitard (layers 1 to 13) have a constant thickness of 2.5 m. The bottom two layers are thicker (10 m and 30 m, respectively), since this portion is deemed to be negligibly affected by the GWHP.
A simplified 2D cross-section of the model is shown in Figure 1. The regional groundwater flow was reproduced by assigning constant hydraulic head values (hydraulic boundary condition (BC) of the first kind) on the western (upgradient) and eastern (downgradient) domain borders, based on both the aquifer thickness and the hydraulic gradient. Hydraulic head initial conditions were then assigned by running the model in steady state conditions and flow-only mode, without pumping wells. Two water wells, one for abstraction and  The regional groundwater flow was reproduced by assigning constant hydraulic head values (hydraulic boundary condition (BC) of the first kind) on the western (upgradient) and eastern (downgradient) domain borders, based on both the aquifer thickness and the hydraulic gradient. Hydraulic head initial conditions were then assigned by running the model in steady state conditions and flow-only mode, without pumping wells. Two water wells, one for abstraction and one for reinjection, were implemented in the model as a Multi-layer Well (flow BC of the fourth kind, according to [34]) screened over the total aquifer thickness. The abstraction well was placed at coordinates x, y (1000; 1500) m, in order to allow the thermal plume to propagate entirely inside the domain, while the reinjection well was located downgradient at a selected distance from the abstraction wells. The well radius was set equal to 0.25 m, a typical value for wells operating at the flow rates considered in this work.
The thermal balance of the aquifer was reproduced by imposing the undisturbed aquifer temperature (assumed equal to 10 • C) on the upstream side of the domain (first kind BC). On the top slice, the heat exchange with the external air was reproduced imposing a Cauchy boundary condition of the third kind [35], i.e., a heat flux proportional to the difference relative to the above-mentioned reference temperature (10 • C). The geothermal temperature gradient was not taken into account, since its effect is usually negligible at depths between 30 m and 50 m [36][37][38][39]; hence, a null heat flux (second kind BC) was imposed on the lowest slice. The same BC was applied to the remaining domain borders. Consistently with the heat boundary conditions, an initial temperature of 10 • C was assigned to the whole domain.
The operation of the GWHP was simulated by imposing a constant temperature difference (∆T = T inj − T abs ) of ±4 • C between the injection and abstraction wells (negative when heating and positive when cooling). The flow rate Q(t) varies with the thermal load P(t) exchanged with the groundwater according to the following equation: The yearly evolution of P(t) follows a sinusoidal trend, as shown in Figure 2. The thermal load is negative when heat is extracted from groundwater, and positive when heat is injected. Details on the model validation are available in the Supplementary Materials. abstraction wells. The well radius was set equal to 0.25 m, a typical value for wells operating at the flow rates considered in this work. The thermal balance of the aquifer was reproduced by imposing the undisturbed aquifer temperature (assumed equal to 10 °C) on the upstream side of the domain (first kind BC). On the top slice, the heat exchange with the external air was reproduced imposing a Cauchy boundary condition of the third kind [35], i.e., a heat flux proportional to the difference relative to the above-mentioned reference temperature (10 °C). The geothermal temperature gradient was not taken into account, since its effect is usually negligible at depths between 30 m and 50 m [36][37][38][39]; hence, a null heat flux (second kind BC) was imposed on the lowest slice. The same BC was applied to the remaining domain borders. Consistently with the heat boundary conditions, an initial temperature of 10 °C was assigned to the whole domain.
The operation of the GWHP was simulated by imposing a constant temperature difference (∆ = − ) of ±4 °C between the injection and abstraction wells (negative when heating and positive when cooling). The flow rate ( ) varies with the thermal load ( ) exchanged with the groundwater according to the following equation: The yearly evolution of ( ) follows a sinusoidal trend, as shown in Figure 2. The thermal load is negative when heat is extracted from groundwater, and positive when heat is injected. Details on the model validation are available in the Supplementary Materials. Using this model setup, a parametric sweep study was performed starting from a default setup and varying one parameter value per simulation. Table 1 reports the complete parameter set, including default values and the ranges considered in this analysis. These parameters are representative of aquifers in which the application of the geothermal system hypothesised in this work could be considered [9][10][11].
The relative weight of each parameter was assessed by calculating the "sensitivity index" (SI), as suggested by Heiselberg et al. [40], defined as: Using this model setup, a parametric sweep study was performed starting from a default setup and varying one parameter value per simulation. Table 1 reports the complete parameter set, including default values and the ranges considered in this analysis. These parameters are representative of aquifers in which the application of the geothermal system hypothesised in this work could be considered [9][10][11]. The relative weight of each parameter was assessed by calculating the "sensitivity index" (SI), as suggested by Heiselberg et al. [40], defined as: where Y max and Y min are the maximum and the minimum values, respectively, of the output variable considered in the analysis, obtained by varying the parameter over the defined range. The output variables considered in our analysis are the length and the width of the thermal plume, as explained later in this section. We assigned a homogenous horizontal hydraulic conductivity value: values adopted for the aquifer are shown in Table 1, while the conductivity of the aquitard was set to K xx = K yy = 10 −7 m·s −1 . The vertical hydraulic conductivity is usually lower than the horizontal one, and K zz = 0.1K xx was used, as suggested in [12].
The thermal conductivity of groundwater was kept equal to λ w = 0.6 W·m −1 ·K −1 , while four different values were considered for the solid phase conductivity λ s . Such values were calculated with Equation (4), in order to obtain bulk values of thermal conductivity that are consistent with those reported in the VDI 4640 guidelines for saturated sedimentary lithologies [41].
A similar approach was adopted for the volumetric heat capacity (ρc) of the aquifer, which was assigned using Equation (2) and literature reference values for water-saturated sand and/or gravel [11]. As for the thermal conductivity, a fixed value was set for the fluid phase (ρ w c w = 4.2 M·J·m −3 ·K −1 ), while three different values were adopted for the solid phase.
Thermal dispersivity in aquifers has scarcely been studied in the literature. According to De Marsily [41], its value depends on the spatial scale of the observed phenomenon, and several empirical formulae are available in the literature [42]. Therefore, a wide range of longitudinal (α L ) and transverse (α T ) dispersivities were employed in the analysis, assuming α T = 0.1α L .
Besides thermal and hydrogeological parameters, the plant setup also significantly influences the propagation of thermal plumes. The effect of the distance between the two wells (L) and of the maximum injected flow rate (Q max ) were investigated, considering three values of well spacing and four values of the maximum annual water flow rate. The analysed range of total thermal power exchanged through the heat pump was 80-800 kW, thus covering most of the building heating and cooling applications for GWHPs. Subsequently, further simulations were run to understand how the propagation of thermal plumes differs if time-varying thermal loads (P(t)) are considered, compared to their yearly average P avg , calculated as follows: where t y is the time span considered (365 days). Different levels of unbalance between the heating and cooling demand were studied, ranging from heating-only to almost perfectly balanced, with the default demand being heating-dominated ( Table 1). The ratio of the yearly heating need E H over the total heat exchanged with the aquifer each year (E TOT = 1820 MWh·y −1 ) ranges between 0.6 and 1, with a default value of 0.75, as depicted in Figure 2. This means that, as a default configuration, the heat extracted in the heating season is 75% of the total heat exchanged; hence, three times the heat injected in the cooling season, which is the remaining 25%. For each of these combinations, a simulation imposing the constant yearly average of the thermal load (i.e., a net heat extraction) was run for comparison. This approximation, which is adopted by software such as Groundwater Energy Designer (GED) [33], allows one to drastically reduce the computational time required for the simulations. A perfectly balanced thermal demand (E H = 0.5 E TOT ) was not considered because the equivalent yearly average load is null, and would therefore lead to a null thermal impact on the ground.
Finally, the previous results of plume size for different configurations were compared with those obtained using the analytical formulae for 2D plume evolution problems provided by Banks [22]. Plume length for long distances can be approximated as: where x is the downstream distance (m) from the reinjection well; v th is the thermal advective velocity, as defined in Equation (7); v e is the groundwater effective flow velocity (m·s −1 ); and t is the simulation time (s). The width of the plume can be calculated with [22]: where y is the asymptotic plume width (m) perpendicular to the groundwater flow direction, i.e., the maximum width of the well release front, which depends on the maximum flow rate escaping down-gradient Q pl (m 3 ·s −1 ), the aquifer thickness b (m), and the Darcy velocity v D (m·s −1 ). Q pl is the fraction of the maximum annual injected flow rate (Q max ) which is not recirculated to the abstraction well and which, therefore, travels down-gradient; it depends on the intensity of thermal recycling between the wells: The parameter X of Equation (15) allows one to assess whether thermal recycling occurs (if X > 1) and, in this case, how strong it is: If X ≤ 1, no thermal recycling occurs and Equation (15) can be simplified to: In our simulations, the plume dimensions were determined considering the isotherm of alteration equal to 1 • C compared to the initial temperature (T 0 = 10 • C), i.e., the isotherm at 9 • C, since heating-dominated thermal loads are used. The plume length was calculated as the downstream distance (m) from the reinjection well, while the width was defined as the maximum extension of the 9 • C isotherm perpendicular to the groundwater flow direction ( Figure 3).
If ≤ 1, no thermal recycling occurs and Equation (15) can be simplified to: In our simulations, the plume dimensions were determined considering the isotherm of alteration equal to 1 °C compared to the initial temperature ( = 10 °C), i.e., the isotherm at 9 °C, since heating-dominated thermal loads are used. The plume length was calculated as the downstream distance (m) from the reinjection well, while the width was defined as the maximum extension of the 9 °C isotherm perpendicular to the groundwater flow direction (Figure 3).  Both the length and the width of the thermal plume were calculated at half the depth of the aquifer, where the strongest thermal alterations are observed, as verified for each simulation. Since the flow rate of injected water varies on an annual basis, and the plume size is subjected to oscillations throughout the simulation time, the maximum annual values of plume width and length are plotted to show the time evolution of the plume dimensions.

Sensitivity Analysis
The results of the sensitivity analysis aim at identifying the parameters that have the strongest influence on the dimensions of a thermal plume. Starting from the hydrodynamic parameters, hydraulic conductivity (K) and hydraulic gradient (i), we observe that plume propagation only depends on the thermal velocity (v th = Ki/(n e R th ), from Equations (3), (6) and (7)). Conversely, the plume size is almost insensitive to the single values of K and i ( Figure 4a). Indeed, the hydraulic conductivity only influences the spatial distribution of the hydraulic head and, consequently, of the flow field close to the abstraction and injection wells. This effect is negligible far downstream of the well doublet. Both the length and the width of the thermal plume were calculated at half the depth of the aquifer, where the strongest thermal alterations are observed, as verified for each simulation. Since the flow rate of injected water varies on an annual basis, and the plume size is subjected to oscillations throughout the simulation time, the maximum annual values of plume width and length are plotted to show the time evolution of the plume dimensions.

Sensitivity Analysis
The results of the sensitivity analysis aim at identifying the parameters that have the strongest influence on the dimensions of a thermal plume. Starting from the hydrodynamic parameters, hydraulic conductivity ( ) and hydraulic gradient ( ), we observe that plume propagation only depends on the thermal velocity ( = /( ), from Equations (3), (6) and (7)). Conversely, the plume size is almost insensitive to the single values of and ( Figure 4a). Indeed, the hydraulic conductivity only influences the spatial distribution of the hydraulic head and, consequently, of the flow field close to the abstraction and injection wells. This effect is negligible far downstream of the well doublet.
(a) (b) Based on this result, the successive parametric sweep focused on variations of the Darcy velocity . As expected, higher values of lead to a considerable increase in the plume length and a dramatic decrease in its width (Figure 4b). A one order of magnitude increase in the regional groundwater velocity reduces the width of the thermal plume by more than 60%.  Based on this result, the successive parametric sweep focused on variations of the Darcy velocity v D . As expected, higher values of v D lead to a considerable increase in the plume length and a dramatic decrease in its width (Figure 4b). A one order of magnitude increase in the regional groundwater velocity reduces the width of the thermal plume by more than 60%.
The plume extension also reaches a stationary value in a shorter time when the Darcy velocity increases (Figure 5a), with a trend similar to the Moving Infinite Line Source analytical solution [43]. It is also interesting to note that, as the regional groundwater flow velocity increases beyond a certain value (v D = 1 × 10 −5 m/s), the plume length suddenly drops (Figures 4b and 5a). If the Darcy velocity value increases further, the plume length starts to increase again. Indeed, if the groundwater velocity is high enough, the thermal plumes produced in each heating and cooling season are separated into two or more smaller plumes which travel downgradient, while cold and warm plumes are merged at lower velocities ( Figure 6). In this case, the farthest plume was considered for length determination. conductivity is the parameter most subject to uncertainties, depending on the method of determination [12,48]; therefore, it strongly affects the thermal impact quantification of GWHPs.
(a) (b) Unlike the Darcy velocity, the effective porosity has little effect on the plume length, exerted only in the long term (Figure 5b), and does not affect the plume width. Understanding how effective porosity affects the plume dimensions is not straightforward, since different variations of heat transport parameters occur as varies. When the effective porosity increases, the thermal capacity of the porous medium generally increases, since the thermal capacity is usually higher for groundwater than for the solid phase ( > ). An increment of results in a decrease in the effective velocity of groundwater ( ); on the other hand, the thermal retardation factor ( ) also decreases, according to Equation (7). This eventually results in a slight decrease in the thermal velocity = / . Therefore, if we consider only advection, a slight decrease in the plume length would be expected. However, an increase in the effective porosity also causes a reduction of the total thermal conductivity of the porous medium (Equation (4)), since groundwater is usually less conductive than the solid phase ( < ), and this makes the plume length increase, as explained later in this section. These contrasting effects are reflected in plume propagation. As long as the plume is expanding in the aquifer and the propagation towards the unsaturated layers is negligible (i.e., at the onset of GWHP operation), the length of the thermally affected zone diminishes as the effective porosity increases. On the other hand, as the plume also starts to expand in the unsaturated layers, the reduction of the thermal conductivity is the parameter most subject to uncertainties, depending on the method of determination [12,48]; therefore, it strongly affects the thermal impact quantification of GWHPs.
(a) (b) Unlike the Darcy velocity, the effective porosity has little effect on the plume length, exerted only in the long term (Figure 5b), and does not affect the plume width. Understanding how effective porosity affects the plume dimensions is not straightforward, since different variations of heat transport parameters occur as varies. When the effective porosity increases, the thermal capacity of the porous medium generally increases, since the thermal capacity is usually higher for groundwater than for the solid phase ( > ). An increment of results in a decrease in the effective velocity of groundwater ( ); on the other hand, the thermal retardation factor ( ) also decreases, according to Equation (7). This eventually results in a slight decrease in the thermal velocity = / . Therefore, if we consider only advection, a slight decrease in the plume length would be expected. However, an increase in the effective porosity also causes a reduction of the total thermal conductivity of the porous medium (Equation (4)), since groundwater is usually less conductive than the solid phase ( < ), and this makes the plume length increase, as explained later in this section. These contrasting effects are reflected in plume propagation. As long as the plume is expanding in the aquifer and the propagation towards the unsaturated layers is negligible (i.e., at the onset of GWHP operation), the length of the thermally affected zone diminishes as the effective porosity increases. On the other hand, as the plume also starts to expand in the unsaturated layers, the reduction of the thermal The drop of the plume length at high groundwater velocities can be explained by the fact that the separation of cold and hot plumes results in a more efficient heat exchange with the neighbouring layers (the vadose zone and aquitard), and, hence, the total length of the plume is reduced. Our results confirm the conclusion drawn in previous studies that the hydrodynamic parameters of the aquifer have a very strong influence on thermal plume size [44][45][46][47]. Hydraulic conductivity is the parameter most subject to uncertainties, depending on the method of determination [12,48]; therefore, it strongly affects the thermal impact quantification of GWHPs.
Unlike the Darcy velocity, the effective porosity has little effect on the plume length, exerted only in the long term (Figure 5b), and does not affect the plume width. Understanding how effective porosity affects the plume dimensions is not straightforward, since different variations of heat transport parameters occur as n e varies. When the effective porosity increases, the thermal capacity of the porous medium ρc generally increases, since the thermal capacity is usually higher for groundwater than for the solid phase (ρ w c w > ρ s c s ).
An increment of n e results in a decrease in the effective velocity of groundwater (v e ); on the other hand, the thermal retardation factor (R th ) also decreases, according to Equation (7). This eventually results in a slight decrease in the thermal velocity v th = v e /R th . Therefore, if we consider only advection, a slight decrease in the plume length would be expected. However, an increase in the effective porosity also causes a reduction of the total thermal conductivity of the porous medium (Equation (4)), since groundwater is usually less conductive than the solid phase (λ w < λ s ), and this makes the plume length increase, as explained later in this section. These contrasting effects are reflected in plume propagation. As long as the plume is expanding in the aquifer and the propagation towards the unsaturated layers is negligible (i.e., at the onset of GWHP operation), the length of the thermally affected zone diminishes as the effective porosity increases. On the other hand, as the plume also starts to expand in the unsaturated layers, the reduction of the thermal conductivity of the porous medium above and below the saturated zone due to a higher porosity value results in a slight long-term increase in the thermal plume length (Figure 5b).
Thermal dispersivity is one of the most influential aquifer properties for the propagation of thermal plumes; however, the scarcity of information and reliable experimental data available on this topic causes significant design issues [49][50][51]. As is evident from Figure 7a, the thermal dispersivity has a considerable effect on plume length, while its effect on plume width is negligible, especially for low dispersivity values. conductivity of the porous medium above and below the saturated zone due to a higher porosity value results in a slight long-term increase in the thermal plume length (Figure 5b). Thermal dispersivity is one of the most influential aquifer properties for the propagation of thermal plumes; however, the scarcity of information and reliable experimental data available on this topic causes significant design issues [49][50][51]. As is evident from Figure 7a, the thermal dispersivity has a considerable effect on plume length, while its effect on plume width is negligible, especially for low dispersivity values. When the value of dispersivity increases, energy dissipation in the aquifer increases, so the plume propagates slower, along the groundwater flow direction. Furthermore, longitudinal thermal dispersivities below 2 m do not sensibly influence the extension of the plume: we observe a negligible difference between = 0.5 m and = 2 m . This is an interesting finding for addressing the numerical modelling issue of oscillatory results, which arise when null or very low dispersivity values are adopted in the presence of strong advection [34]. Usually, in these cases, upwinding schemes are used to introduce a numerical dispersion, thus avoiding the oscillation of model results, but also altering them [34]. It is therefore possible to avoid using upwinding schemes by assigning thermal dispersivities of < 2 m and < 0.2 m; this prevents both oscillatory results and the underestimation of the thermal impact on the aquifer. Values of dispersivity higher than 2 m should be set carefully, since they could lead to a strong reduction of the estimated plume When the value of dispersivity increases, energy dissipation in the aquifer increases, so the plume propagates slower, along the groundwater flow direction. Furthermore, longitudinal thermal dispersivities below 2 m do not sensibly influence the extension of the plume: we observe a negligible difference between α L = 0.5 m and α L = 2 m. This is an interesting finding for addressing the numerical modelling issue of oscillatory results, which arise when null or very low dispersivity values are adopted in the presence of strong advection [34]. Usually, in these cases, upwinding schemes are used to introduce a numerical dispersion, thus avoiding the oscillation of model results, but also altering them [34]. It is therefore possible to avoid using upwinding schemes by assigning thermal dispersivities of α L < 2 m and α T < 0.2 m; this prevents both oscillatory results and the underestimation of the thermal impact on the aquifer. Values of dispersivity higher than 2 m should be set carefully, since they could lead to a strong reduction of the estimated plume length (e.g., −45% for α L = 50 m relative to α L = 2 m), as shown in Figure 7a. Notably, variations in the longitudinal and transverse dispersivity do not affect the maximum plume width, which mainly depends on the hydrodynamic aquifer properties.
A few numerical simulation studies of open-loop plants in the long-term (in the order of tens of years) have been published. In previous studies, limited to the saturated layer [44] and usually considering a time span up to one year [32,44], the heat conduction mechanism is negligible compared to advection and dispersion. However, when a sufficiently long time is considered, the conductive propagation of heat in the ground also involves the vadose zone above the aquifer and the underlying layers [52]. Figure 7d depicts the effect of such a mechanism on the longitudinal temperature distribution, comparing the assumptions of a no-flux boundary condition (second kind) applied to the ground surface and of the aforementioned reference temperature (third kind BC). The effect of heat exchange with neighbouring layers gets stronger over time (Figure 7b,d) and is negligible in the short term; this explains why previous studies, which considered short operating periods (less than one year [32,44,46]), concluded that heat conduction is negligible in GWHPs.
Conversely, the plume width is mainly affected by the hydrodynamic parameters of the aquifer, while the heat exchange has a negligible effect.
The volumetric heat capacity of the solid matrix has a minor influence on the plume size in the long-term, as shown in Figure 7c. Both the plume width and length slightly decrease as the thermal capacity of the porous medium (ρc) increases. This is due to the capability of a portion of the aquifer to store a greater amount of heat, given the same thermal alteration.
Long-term propagation of thermal plumes is also influenced by the thickness of the saturated and vadose zones, which both play a substantial role in the expansion of the thermal plume. As the saturated thickness increases (Figure 8a), the plume takes longer to reach the top of the aquifer and exchange heat with the unsaturated zone. Hence, the conductive exchange with neighbouring layers occurs later and heat remains confined for a longer period of time in the aquifer. This implies that the plume propagates considerably further along the groundwater flow direction, as confirmed by the analytical solutions of a 2D plane-symmetric transient heat transport problem reported by Tan et al. [53]. In previous studies, limited to the saturated layer [44] and usually considering a time span up to one year [32,44], the heat conduction mechanism is negligible compared to advection and dispersion. However, when a sufficiently long time is considered, the conductive propagation of heat in the ground also involves the vadose zone above the aquifer and the underlying layers [52]. Figure 7d depicts the effect of such a mechanism on the longitudinal temperature distribution, comparing the assumptions of a no-flux boundary condition (second kind) applied to the ground surface and of the aforementioned reference temperature (third kind BC). The effect of heat exchange with neighbouring layers gets stronger over time (Figure 7b,d) and is negligible in the short term; this explains why previous studies, which considered short operating periods (less than one year [32,44,46]), concluded that heat conduction is negligible in GWHPs.
Conversely, the plume width is mainly affected by the hydrodynamic parameters of the aquifer, while the heat exchange has a negligible effect.
The volumetric heat capacity of the solid matrix has a minor influence on the plume size in the long-term, as shown in Figure 7c. Both the plume width and length slightly decrease as the thermal capacity of the porous medium ( ) increases. This is due to the capability of a portion of the aquifer to store a greater amount of heat, given the same thermal alteration.
Long-term propagation of thermal plumes is also influenced by the thickness of the saturated and vadose zones, which both play a substantial role in the expansion of the thermal plume. As the saturated thickness increases (Figure 8a), the plume takes longer to reach the top of the aquifer and exchange heat with the unsaturated zone. Hence, the conductive exchange with neighbouring layers occurs later and heat remains confined for a longer period of time in the aquifer. This implies that the plume propagates considerably further along the groundwater flow direction, as confirmed by the analytical solutions of a 2D plane-symmetric transient heat transport problem reported by Tan et al. [53]. On the other hand, a higher saturated thickness ( ) leads to larger aquifer transmissivity, which reduces the width of the plume. Therefore, within the range of saturated thicknesses commonly found in shallow alluvial aquifers (where most of the GWHPs are installed), plume size can vary substantially in the long term.
Similarly, an increase in the unsaturated thickness (d) above the aquifer causes a significant increase in the plume length, while the plume width is insensitive to the depth of the saturated zone (Figure 8b). Indeed, heat transport across the unsaturated layers towards the atmosphere is conductive, and a thicker vadose zone has a higher thermal resistance, which limits this heat transfer mechanism. Figure 8b shows that the effect of vadose zone thickness is visible as the thermal On the other hand, a higher saturated thickness (b) leads to larger aquifer transmissivity, which reduces the width of the plume. Therefore, within the range of saturated thicknesses commonly found in shallow alluvial aquifers (where most of the GWHPs are installed), plume size can vary substantially in the long term.
Similarly, an increase in the unsaturated thickness (d) above the aquifer causes a significant increase in the plume length, while the plume width is insensitive to the depth of the saturated zone ( Figure 8b). Indeed, heat transport across the unsaturated layers towards the atmosphere is conductive, and a thicker vadose zone has a higher thermal resistance, which limits this heat transfer mechanism. Figure 8b shows that the effect of vadose zone thickness is visible as the thermal alteration approaches the ground surface. For example, after five years, the plume length curve for d = 10 m starts to diverge from the others, and the same occurs after 10 years for the curve with d = 20 m. This study, therefore, confirms the importance of depth to water table and aquifer thickness on the migration of the thermal plume, as was previously reported [11,19].
To conclude the sensitivity analysis, some plant settings also have an appraisable influence on the propagation of thermal plumes in GWHPs.
The thermal alteration at the injection well depends on the temperature of the re-injected water, which in turn strongly depends on the distance L between the injection and abstraction wells. As shown in Figure 9a, an increase in the distance L results in a decrease in the plume length and an increase in its width. This is due to a reduction of the recycled flow rate (see Equations (15) and (16)), which results in a reduction of the thermal alteration at the reinjection well (which influences the length of the plume), and an increase in the flow rate released downstream of the well doublet (Q pl ), thus increasing the plume width (see Equation (14)). For the considered plant configuration, no thermal recycling occurs when the injection well is located 100 m downstream of the abstraction well: hence, no variations in the plume dimensions are expected over this value of L. = 10 m starts to diverge from the others, and the same occurs after 10 years for the curve with = 20 m. This study, therefore, confirms the importance of depth to water table and aquifer thickness on the migration of the thermal plume, as was previously reported [11,19].
To conclude the sensitivity analysis, some plant settings also have an appraisable influence on the propagation of thermal plumes in GWHPs.
The thermal alteration at the injection well depends on the temperature of the re-injected water, which in turn strongly depends on the distance between the injection and abstraction wells. As shown in Figure 9a, an increase in the distance results in a decrease in the plume length and an increase in its width. This is due to a reduction of the recycled flow rate (see Equations (15) and (16)), which results in a reduction of the thermal alteration at the reinjection well (which influences the length of the plume), and an increase in the flow rate released downstream of the well doublet ( ), thus increasing the plume width (see Equation (14)). For the considered plant configuration, no thermal recycling occurs when the injection well is located 100 m downstream of the abstraction well: hence, no variations in the plume dimensions are expected over this value of . The injected-water flow rate and, therefore, the thermal power exchanged with the ground, has a very strong influence on thermal plume propagation (Figure 9d). The plume length has a roughly logarithmic correlation with the flow rate, and an increase in by a factor of ten results in a fourfold increase in the plume length. On the other hand, the width of the thermally affected area increases almost linearly with the injected flow rate (Figure 9c), as expected from Equation (14). The thermal alteration induced downstream of the injection well during each heating season interacts with the plume generated during the previous cooling season, determining the temperature The injected-water flow rate and, therefore, the thermal power exchanged with the ground, has a very strong influence on thermal plume propagation (Figure 9d). The plume length has a roughly logarithmic correlation with the flow rate, and an increase in Q max by a factor of ten results in a fourfold increase in the plume length. On the other hand, the width of the thermally affected area increases almost linearly with the injected flow rate (Figure 9c), as expected from Equation (14). The thermal alteration induced downstream of the injection well during each heating season interacts with the plume generated during the previous cooling season, determining the temperature distribution in the ground. For this reason, keeping the same total heat exchanged (E TOT = E H + E C ) and varying the proportion between the heating (E H ) and cooling (E C ) demands results in different thermal plume dimensions. As expected, the higher the ratio E H /E TOT , the longer and the wider the cold thermal plume originated by the well doublet (Figure 9d).
Finally, the relative importance of the analysed parameters for thermal plume size after 50 years of plant operation was assessed by means of the Sensitivity Index described by Equation (11). Among the hydrogeological and thermal parameters of the ground, the Darcy velocity v D is the most influential for plume propagation, both in the longitudinal and the transversal direction ( Figure 10). A strong impact on the plume length is also exerted by the dispersivity α L and the thermal conductivity of the solid matrix λ s , while the porosity n e and the volumetric heat capacity of the solid matrix ρ s c s have a negligible effect on thermal plume size. Strong variations in the plume length and width were observed while varying geometry parameters, such as the aquifer thickness b and the well spacing L, whereas the depth to the water table d is only relevant for plume length. The most important parameter for plume propagation in the long term is the injected water flow rate Q max , followed by the ratio of the yearly heating demand over the total demand of the heat pump E H /E TOT . and varying the proportion between the heating ( ) and cooling ( ) demands results in different thermal plume dimensions. As expected, the higher the ratio / , the longer and the wider the cold thermal plume originated by the well doublet (Figure 9d).
Finally, the relative importance of the analysed parameters for thermal plume size after 50 years of plant operation was assessed by means of the Sensitivity Index described by Equation (11). Among the hydrogeological and thermal parameters of the ground, the Darcy velocity is the most influential for plume propagation, both in the longitudinal and the transversal direction ( Figure 10). A strong impact on the plume length is also exerted by the dispersivity and the thermal conductivity of the solid matrix , while the porosity and the volumetric heat capacity of the solid matrix have a negligible effect on thermal plume size. Strong variations in the plume length and width were observed while varying geometry parameters, such as the aquifer thickness and the well spacing , whereas the depth to the water table is only relevant for plume length. The most important parameter for plume propagation in the long term is the injected water flow rate , followed by the ratio of the yearly heating demand over the total demand of the heat pump ⁄ .
Detailed results of this parametric study for plume length and width are reported in Table S1 and S2 of the Supplementary Materials.

Comparison with Simplified Models
Transient numerical models are the most precise and rigorous tool for the simulation of the subsurface thermal alteration induced by a GWHP. However, the expense and the time required to perform such simulations are usually very large. Two alternative simplified methods were therefore assessed: numerical simulation with a constant thermal load (equal to the yearly average value, as from Equation (12)) and analytical formulae available for thermal plume length and width evaluation (Equations (13) and (14)) [22]. By imposing a constant thermal load, the computational time for the simulations can be reduced from several hours to a few minutes. Analytical formulae allow fast and inexpensive evaluations of the thermal impact of GWHPs, although they neglect thermal dispersion and conduction.
We used these two approaches to calculate plume length and width at time = 50 years, comparing the results with those of the numerical model reported in the previous paragraph in order to assess whether the approximations introduced by these methods are acceptable or not. As shown in Figure 11, simulations performed with a constant thermal load introduce an error in the values of both the plume length and width.
A remarkable result is that, for heating-only systems (i.e., / = 1), the plume length is practically the same but, as the heating-cooling imbalance of the thermal load diminishes, the Detailed results of this parametric study for plume length and width are reported in Tables S1 and S2 of the Supplementary Materials.

Comparison with Simplified Models
Transient numerical models are the most precise and rigorous tool for the simulation of the subsurface thermal alteration induced by a GWHP. However, the expense and the time required to perform such simulations are usually very large. Two alternative simplified methods were therefore assessed: numerical simulation with a constant thermal load (equal to the yearly average value, as from Equation (12)) and analytical formulae available for thermal plume length and width evaluation (Equations (13) and (14)) [22]. By imposing a constant thermal load, the computational time for the simulations can be reduced from several hours to a few minutes. Analytical formulae allow fast and inexpensive evaluations of the thermal impact of GWHPs, although they neglect thermal dispersion and conduction.
We used these two approaches to calculate plume length and width at time t = 50 years, comparing the results with those of the numerical model reported in the previous paragraph in order to assess whether the approximations introduced by these methods are acceptable or not. As shown in Figure 11, simulations performed with a constant thermal load introduce an error in the values of both the plume length and width. discrepancy increases ( Figure 12) by up to about 30% for / = 70% (which is a fairly acceptable error). As the thermal load approaches the perfect balance between heating and cooling, i.e., / = 0.5, the error introduced by this approximation rapidly increases to unacceptable levels.
(a) (b) Figure 11. Relative error on plume length (a) and width (b), evaluated with a constant thermal load simulation and analytical formulae, compared to a variable thermal load simulation, for the considered range of parameters; isotherm ∆T = 1 °C, t = 50 years. The yearly averaged load usually leads to an overestimation of the plume length, thus proving a conservative approach, but also to an underestimation of the plume width (Figure 11b). This is due to the fact that the effect of peak flow-rates on the width of the injection well release front is not considered. The tested analytical formula largely overestimates the long-term plume length ( Figure  11a), since the advective model described by Equation (13) neglects thermal dispersion and conduction. On the other hand, the plume width can conservatively be estimated with Equation (14), with an acceptable overestimation of 30-60%. Underestimation of the plume width occurs only for A remarkable result is that, for heating-only systems (i.e., E H /E TOT = 1), the plume length is practically the same but, as the heating-cooling imbalance of the thermal load diminishes, the discrepancy increases ( Figure 12) by up to about 30% for E H /E TOT = 70% (which is a fairly acceptable error). As the thermal load approaches the perfect balance between heating and cooling, i.e., E H /E TOT = 0.5, the error introduced by this approximation rapidly increases to unacceptable levels. discrepancy increases ( Figure 12) by up to about 30% for / = 70% (which is a fairly acceptable error). As the thermal load approaches the perfect balance between heating and cooling, i.e., / = 0.5, the error introduced by this approximation rapidly increases to unacceptable levels.
(a) (b) Figure 11. Relative error on plume length (a) and width (b), evaluated with a constant thermal load simulation and analytical formulae, compared to a variable thermal load simulation, for the considered range of parameters; isotherm ∆T = 1 °C, t = 50 years. The yearly averaged load usually leads to an overestimation of the plume length, thus proving a conservative approach, but also to an underestimation of the plume width (Figure 11b). This is due to the fact that the effect of peak flow-rates on the width of the injection well release front is not considered. The tested analytical formula largely overestimates the long-term plume length ( Figure  11a), since the advective model described by Equation (13) neglects thermal dispersion and conduction. On the other hand, the plume width can conservatively be estimated with Equation (14), The yearly averaged load usually leads to an overestimation of the plume length, thus proving a conservative approach, but also to an underestimation of the plume width (Figure 11b). This is due to the fact that the effect of peak flow-rates on the width of the injection well release front is not considered. The tested analytical formula largely overestimates the long-term plume length (Figure 11a), since the advective model described by Equation (13) neglects thermal dispersion and conduction. On the other hand, the plume width can conservatively be estimated with Equation (14), with an acceptable overestimation of 30-60%. Underestimation of the plume width occurs only for the highest considered value of v D (2 × 10 −5 m/s), because of the plume separation effect and the increased heat dispersion. Hence, the equations proposed by Banks ([22], from Equation (13) to (16)) provide a good and conservative estimation of the plume width, while they cannot be used to assess the longitudinal extension of the thermally affected area.
For detailed results of the methods comparison study, for simulations at 10 years and 50 years, we refer the reader to Tables S1 and S2 of the Supplementary Materials.

Conclusions
Assessing the long-term subsurface thermal impact of open-loop geothermal systems is of crucial importance for their design and their proper spatial planning. Various numerical models and analytical solutions can be used to estimate the thermal plume produced by a GWHP, each characterized by different costs, degrees of complexity, and precision.
This work presents a parametric study through numerical flow and heat transport simulations carried out to identify key parameters and compare their relative influence within their typical ranges of variation. Moreover, simplified numerical modelling with a constant thermal load and analytical solutions for 2D plume propagation in an aquifer were investigated to quantify the discrepancy compared to transient numerical modelling.
Among the hydraulic and thermal subsurface parameters, the most influential for both plume length and width is the Darcy velocity (v D ), i.e., the product of hydraulic conductivity (K) and gradient (i). In this light, it is very important to determine the hydraulic conductivity with reliable in-situ tests, in order to perform realistic simulations and avoid major errors in plume size estimations.
The thermal properties of the ground affect the plume length but have a negligible influence on its width. In particular, thermal conductivity (λ) and dispersivity (α) of the ground play an important role in the long-term, while variations of heat capacity (ρc) have minor effects on plume size.
Depth to water table (d) and aquifer thickness (b) are also key factors in the propagation of thermal plumes, thus confirming the necessity of a 3D modelling geometry. Plume length dramatically increases with the thickness of the saturated and unsaturated zones, as thermal exchange with air is reduced, and propagation of the thermal plume is confined to the saturated zone. On the other hand, an increase in the saturated thickness results in higher transmissivity and, consequently, in a reduction of the plume width with an almost linear inverse trend.
The distance between abstraction and injection wells (L) noticeably influences the width of the thermal plume, while its impact on its longitudinal extension is lower, though not negligible.
Regarding thermal loads, a logarithmic increase in the plume length with the abstracted/injected flow rate (and consequently with the thermal power) is observed, while the correlation is almost linear for the plume width. The balancing of the heating and cooling demand of the building sensibly affects the size of the thermal plume, which reaches its maximum for a fully heating-dominated (or cooling-dominated) demand and its minimum for a perfectly balanced load.
The assessment of the thermal plume adopting a constant thermal load (equal to its yearly average) is not reliable at high groundwater velocities, due to the separation of plumes originated during heating and cooling seasons. In such cases, a simulation with a time-varying thermal load is required.
On the other hand, for heating or cooling-dominated loads, constant load simulations provide a good approximation of the thermal plume lengths. In this way, the computational time is drastically reduced, and simpler and cheaper software can be used. For balanced or slightly unbalanced thermal loads, however, the length of the thermal plume is largely overestimated using the constant load approach.
Finally, analytical formulae for 2D advective propagation provided by Banks [22] can be used to quickly assess the width of thermal plumes, while they provide unrealistically high estimates of the plume length.
We provide useful insights for the evaluation of the thermal impact of GWHPs on aquifers, and hence for their planning in densely inhabited areas. The key parameters are identified, along with the error margins deriving from a poor characterisation of the site.
The effects of simplifying assumptions, such as using the yearly average of the thermal load or adopting available analytical solutions, are also analysed, thus defining under what conditions such methods can be used with acceptable accuracy.
Supplementary Materials: The following are available online at www.mdpi.com/1996-1073/10/9/1385/s1, Figure S1: Mesh adopted in the analysis, abstraction and injection wells in red (252k nodes), Figure S2: Mesh convergence study results: thermal plume length (a) and width (b) after 50 years, for different levels of mesh refinement, Table S1: Length of the thermal plume after 10 years and 50 years of operation, calculated with different methods-parametric study, Table S2: Width of the thermal plume after 10 years and 50 years of operation, calculated with different methods-parametric study.