Abstract
Thermal perturbation produced in the subsurface by open-loop groundwater heat pump systems (GWHPs) must be predicted and constantly controlled, especially in the shallow aquifers of more densely urbanized areas, in order to guarantee plants’ long-term sustainable use and to avoid adverse effects on adjacent geothermal systems. Transient conditions in the flow dynamic can be successfully modelled by means of numerical modelling tools. However, for small plants in suitable hydrogeological systems, an alternative tool for predicting the thermally affected zone (TAZ) around the injection well can be found in analytical solutions for steady advective transport in a shallow aquifer. The validity of using steady analytical solutions to predict the TAZ development at the end of two different cooling seasons (2010 and 2016) was tested in the Politecnico di Torino GWHP system (NW Italy). When fixing the constant thermal difference () between the injection and abstraction wells at 5°C, results revealed that a rather reliable assessment of the TAZ of Politecnico di Torino GWHPs, in Turin shallow aquifer, can be performed by plotting the cumulative distribution function of the injected discharge rate () and setting 63% as a steady value.
1. Introduction
Among low-enthalpy geothermal systems, open-loop groundwater heat pumps (GWHPs) currently represent one of the major technologies applied in the heating and cooling of buildings. These systems were designed to take advantage of the heat available in the shallow subsurface (<400 m) by withdrawing water from a well or a surface water source, passing it through a heat exchanger and redischarging water into an injection well or nearby river [1].
However, the continuous increase in the implementation of open-loop GWHPs and their reinjection of warmer water could potentially cause, even in the short term, a significant subsurface environmental impact associated with local variations in groundwater temperature within the thermally affected zone (TAZ), especially in the shallow aquifers of more densely urbanized areas. As thermal plumes could adversely affect adjacent and neighbouring geothermal systems, the TAZ extension must be well predicted and constantly controlled in order to guarantee the systems’ long-term sustainable use. For this purpose, numerical models have been widely applied in the past years. Simulation approaches have made it possible to better understand and predict subsurface heat transport mechanisms [2–6]; however, these simulations require the use of complex, expensive, and time-consuming numerical calculation tools [7].
Therefore, many authors have attempted to develop analytical methods to analyse heat transport in the subsurface for small plants in simplified conduction–advection-dominated systems [8–10].
Among others, Banks [10] proposed a robust and simplified analytical solution to solve the problem of the assessment of the expected extension of TAZ associated to simple configuration plants in a homogeneous shallow aquifer.
The use of analytical solutions in place of more complex numerical modelling creates several advantages in terms of the reduction of the time and technical complexity required for plant design; therefore, the applicability of such solutions to real plants needs to be verified.
To allow for correct analytical solution application, many assumptions should in fact be validated concerning site-specific characteristics. Moreover, analytical solutions have been strictly designed for steady abstraction and reinjection flow conditions both in terms of discharge rate, temperature variation in the heat pump, and temperature of the reinjected groundwater.
The described constraints on the current analytical solutions limit their applicability because, even for the cases in which the hydrogeological conditions are similar to those that are assumed for the validity of the equations, steady operating arrangements are usually far from reality.
The ordinary functioning conditions of open-loop GWHPs are dynamic and time variable.
In modern plants, the inverter technologies tend also to modify the operating conditions of the heat pump together with the pumps in the abstraction wells, following the real energy demand of the building’s heat needs.
Daily and hourly cycles of the outdoor temperatures variations, along with variabilities in the usage of the building, usually induce a dynamic running functioning of the systems including the heat pump, abstraction wells, and reinjection wells. As a result, the pumping rate in the abstraction wells (with the correspondent reinjection rate) and the temperature of the reinjected groundwater are time variable.
In order to provide reliable modelling, transient flow conditions should therefore be considered both for the discharge rate and for the reinjection temperature.
While each of these main conditions allow to enhance the field of applicability of numerical modelling to solve the problem of the TAZ assessment, at the same time, they reduce the potential of the analytical solutions.
In this paper, we address the main problem of the applicability of analytical solutions to a real hydrogeological and plant context by evaluating the accuracy of their results for the estimation of TAZ dimensions at the end of different system-functioning seasons.
In detail, we verified the reliability of the analytical solution proposed by Banks [10] to predict the thermal perturbation induced by a real open-loop GWHP plant serving the Politecnico di Torino (Turin, northern Italy). We initially validated the transient condition numerical modelling results provided by the finite-element code FEFLOW® 6.2 [11] through a comparison with the groundwater temperature data collected during two different cooling seasons (2010 and 2016).
Then, following a statistical approach, we individuated the steady discharge rate and the constant reinjection groundwater temperature that should be used in Banks’s [10] analytical solutions for assessing the TAZ dimensions of Politecnico di Torino GWHP plant, in order to provide reliable results considering a system where flow rate and temperature are both time-variable parameters.
2. Materials and Methods
2.1. Test Site
The analysed test site (Politecnico di Torino open-loop GWHPs) is located in Turin, Piedmont region in NW Italy (geographical coordinates: 45°0345N, 7°3941E, elevation: 248 m a.s.l).
As shown in Figure 1, the Politecnico di Torino open-loop geothermal system (CF1) is composed of a 40 m deep pumping well (P2), a 47 m deep injection well (P4), and a 35 m deep monitoring piezometer (S2), located downgradient of P4. Distances among the described components, measured along flux lines, are, respectively: , , and .
The described well-doublet system works only during the spring and summer seasons, to cool some university buildings.
All CF1 system components affect the Turin shallow unconfined aquifer [12] (Figure 2). From a geological and hydrogeological point of view, the Politecnico di Torino test site is already well known, as it has been studied in previous works [5, 13].
2.2. Geological Setting
The Turin plain extends from the Rivoli-Avigliana Morainic Amphitheatre (RAMA) on its western extreme to the Torino Hill on its eastern border (Figure 2). Several rivers delimit the Turin urban area: the Stura di Lanzo River (northern boundary), the Sangone River (southern boundary), and the Po River (eastern boundary).
Piemonte [14] documented the hydrogeology of the plain area with a high degree of confidence. Information from drilled wells and downhole log tests has confirmed the presence of two different lithological units with distinct hydraulic properties, stratigraphic units 1 and 2 (Figure 3).
Stratigraphic unit 1 (Middle Pleistocene–Holocene) consists of a continental alluvial cover, primarily composed of coarse gravel and sandy sediments with local subordinate clay lens inclusions (≤1.5 m thick). Stratigraphic unit 2 (early Pliocene–Middle Pleistocene) includes shallow marine environment deposits (Sabbie di Asti and Argille di Lugagnano) composed of fossiliferous sandy-clayey layers with subordinate fine gravelly and coarse sandy marine layers as well as quartz-micaceous sands [15, 16]. The outwash plain substrate (unit 3 in Figure 2) consists of a Cenozoic terrigenous marine succession composed of conglomerates, sandstones, and marls (Piemont Tertiary Basin) [17].
2.3. Hydrogeological Setting
The main surface water drainage network of the Turin plain (Stura di Lanzo, Sangone, Dora Riparia, and Po rivers) communicates hydraulically with an unconfined aquifer extending over the entire urban plain (NNW–SSE gradient of 0.29% towards the Po river).
By means of a step drawdown test on the abstraction well (47 m deep) of the GWHP system, an unconfined aquifer hydraulic conductivity value (K1) of was evaluated. Based on the aquifer lithology determined by examining well-drilling logs, an effective porosity of 0.20 was estimated. In addition, for stratigraphic unit 1, a volumetric heat capacity value equal to J·m-3K-1 was calculated as a composition-weighted mean based on the lithological records [5].
2.4. Heat Transport in Porous Shallow Aquifers
Heat transport theories have been studied for several decades [8, 9, 19], and attention on underground heat transport has been continuously increasing over the last years due to the growth in the exploitation of shallow geothermal energy.
Focusing on shallow unconfined aquifers, thermal energy movement is predominately realized by advection caused by the fluid motion inside the medium and conduction driven by the temperature gradient through the solid phase and the heat exchange between the solid and fluid phases.
Additionally, the differential groundwater flow paths at the pore scale and the heterogeneity of the permeability field at macroscopic scale create the thermal dispersion phenomenon.
These mechanisms are described by the equation for heat transport in porous media with the assumption of a groundwater flow aligned with the -axis (2D form, - horizontal plane), as follows [19]:
(for nomenclature, see Table 1).
The first term in Equation (1) represents the energy variation over time, where can be computed as the weighted arithmetic mean of the solids of the aquifer and water [19]:
The second term in Equation (1) is the advection term, which is a function of Darcy velocity ():
The third and fourth terms in Equation (1) are the conduction terms, which, respectively, depend on the longitudinal () and transverse () thermal conductivity of the porous medium [], defined by the bulk thermal conductivity of the porous media () and the thermal dispersion coefficient ():
Referring to the first term in Equation (1), instantaneous thermal equilibration between the groundwater and soil can be assumed, which means that a single point in an aquifer system (both the matrix and the mobile groundwater) is characterized by a single temperature () at any given time () [20].
In addition, groundwater flow advection usually has the principal impact on heat transport in the porous medium. If there is a lack of groundwater flow, heat transport occurs only due to diffusion. Due to the fact that natural groundwater flow is present under most conditions, thermal conduction can be considered insignificant and be neglected [12].
Lo Russo and Taddia [21], by analysing the groundwater monitoring data from the surroundings of the injection wells of open-loop geothermal plants, have highlighted the prevalence of the advective heat transport component, enabling us to correctly refer to the thermal front velocity associated with the Politecnico di Torino open-loop system as an advective-type velocity.
The value of thermal plume advective velocity differs from the natural groundwater velocity by a factor (thermal retardation factor) (Equation (5)), derived by factoring out the second term of heat (Equation (1)):
By assuming, for the Turin saturated porous shallow aquifer, an effective porosity value of 0.20, Lo Russo et al. [12] estimated an factor equal to 2.4: natural groundwater velocity would therefore have 2.4 times the migration velocity of the thermal plume (approximately 1.30 mday−1).
2.5. FEFLOW Model Setup
Heat transport mechanisms in the unconfined aquifer have been modelled by using the finite-element code FEFLOW® 6.2 [11], to evaluate the thermal perturbation induced by the water reinjection in Politecnico di Torino open-loop well doublets.
Starting with hydrodynamic and thermal parameters assigned to unit 1 and unit 2 by Lo Russo et al. [12] (see Table 2) and available hourly datasets of [m3day-1] and [°C] measured continuously by probes fixed in Politecnico di Torino GWHP’s injection well (P4), we performed a two-unit conceptual model simulation wherein hydrostratigraphic unit 1 represents the exploited unconfined alluvial aquifer (Figure 4).
A complete list of the hydrodynamic and thermal parameters assigned to Unit 1 and Unit 2 is provided in Table 2.
The dimensions of the plan-view grid in the model were set equal to 1914 m (NW–SE) and 1525 m (SW–NE): the model area was designed to be larger than the site under investigation in order to reduce any impact of the assumed boundary conditions on model outcomes. The average mesh spacing domain was 15 m, refined to 3 m in areas near the wells to improve thermal plume resolution. The selected grid spacing was defined after a suitable trial test.
Due to lack of infiltration data and the surface characteristics of the model area, which is mostly covered by buildings and roads and therefore considered impermeable, rainfall infiltration rate was not taken into account.
In addition, based on groundwater level continuous monitoring, the unperturbed groundwater flow can be considered stable throughout the year. Constant heads (Dirichlet conditions) are set on the western (230 m a.s.l.) and eastern (220 m a.s.l.) boundaries in accordance with available onsite potentiometric surface measurements [14].
Using the results from different thermal logs performed in the extraction and injection wells (P2-P4) before pumping and in the piezometer (S2), which showed negligible vertical temperature variation, average natural groundwater temperature was set to 15.0°C throughout the aquifer.
Due to the unsaturated zone thickness and the lack of detailed information about its hydraulic and thermal characteristics, adiabatic conditions were established for the ground-surface over the unsaturated zone.
Eventually, as the Politecnico di Torino GWHPs operate at variable loads , depending on the building cooling demands, heat transport mechanisms in the subsoil were simulated by considering transient flow conditions. By means of multiparameter probes installed in the P2 and P4 wells and in the S2 piezometer, potentiometric levels [m], temperature [°C], and pumping and reinjection rates [m3day-1] were constantly recorded.
In detail, the performed simulation covered the hourly cooling period dataset from 2010 (from 25th May to 30th September—operating time of 129 days) and the hourly cooling period dataset from 2016 (from 6th July to 30th September—operating time of 90 days). To assess the reliability of the generated numerical models, calibration for both the 2010 and 2016 datasets was performed by comparing different simulated temperature values with the available temperature measurements from the piezometer (S2).
Following model calibration, a second modelling scenario was created by moving the reinjection well location (P4) downstream with respect to the abstraction well (P2). This condition was necessary in order to create a good comparison between the results obtained by the simulation model and Banks’s analytical solution. The presence of a reinjection well located downstream of the abstraction well along groundwater flow direction is indeed one of the primary assumptions for the analytical solution’s applicability.
2.6. Banks’s (2011) Analytical Solutions
Thermal plume size results from the simulation model were compared to those obtained using the 2D analytical solutions for thermal plume evolution under steady state conditions provided by Banks [10].
According to Banks [10], the thermal plume length () can indeed be approximated as shown in Equation (6), while the maximum downstream width () is given by Equation (7) (see Figure 5). where is a constant fraction of maximum injected flow rate [m3day-1] which is not recirculated to the abstraction wells.
In addition, using Equation (7), the Darcy velocity can be derived as follows:
By using Equation (8), it is then possible to also rewrite the 2D analytical solution for thermal plume length () as a function of (Equation (9)):
The applicability of the described Banks’s [10] analytical solutions is however strongly bound to the validity of the following main assumptions: (1)A homogeneous, isotropic aquifer of constant thickness(2)The presence of confined or unconfined aquifers where head variations are small with respect to total thickness(3)The presence of fully penetrating wells(4)Instantaneous thermal equilibrium between the groundwater and matrix(5)An injection well immediately downgradient of the abstraction well(6)Lack of vertical heat conduction (thermally confined aquifer)(7)Neglected dispersion phenomena(8)Constant pumping rate and injection temperature
Considering the analysed Politecnico di Torino test site, each of these assumptions can be considered valid except for the presence of a constant pumping rate and injection temperature, as both pumping rate and injection temperature are parameters that continuously change over time following the real energy demands associated with building needs.
We therefore provided a reconstruction of an equivalent steady state by means of estimating the temperature [°C] and average pumping rates [m3day-1] over long-term monitored periods (summer seasons 2010 and 2016) by analysing the statistical distribution of available hourly flow rate values [m3day-1].
From a statistical point of view, hourly reinjection flow rates () introduced into the aquifer during the plant’s operating period represent a real function defined on a sample space where for every real value which can take the random variable , there exists a probability [].
Every random variable is characterised by its probability-density-function (PDF) and its cumulative distribution function (CDF), described by its statistical parameters: mean value and variance (or standard deviation) [22].
Mean value () is the centre of gravity of the probability density function (Equation (10)) while the standard deviation (Equation (11)) is a measure of the dispersion around the mean value. where is the population size.
Starting from available hourly values of reinjected [m3day-1] and [°C], for both analysed seasons, we first estimated, by means of Equation (12), the hourly and total thermal load value []:
(for nomenclature, see Table 1).
Second, in order to make the available datasets of the two different analysed seasons (2010 and 2016) statistically comparable, we recalculated the hourly flow rate values () by using the available hourly values of the thermal load () and fixing a constant positive thermal difference () between the injection and abstraction wells () at 5°C (Figure 5).
Recalculated hourly flow rate data [m3day-1] ( values: 3672 for 2010; 2172 for 2016) were then plotted in cumulative distribution functions (CDFs) to determine a reinjection steady discharge rate value () which could be correctly employable in Banks’s [10] analytical solution (Equation (7) and Equation (9)) to provide reliable results of TAZ dimensioning comparable to those obtained by simulations.
3. Results and Discussion
3.1. Calibration Models
To assess the reliability of the generated numerical models, calibration for both the 2010 and 2016 datasets was performed by comparing the different simulated temperature values with the available piezometer measurements (S2).
The 2010 calibration simulation returned a root-mean-square error () of 0.15°C; the simulation model performed by considering the 2016 dataset returned a root-mean-square error () of 0.65°C. The match between the simulated and measured temperature is a confirmation of the validity of the modelling assumptions and of the reliability of the TAZ extension obtained by FEFLOW simulation.
Although the RMSE value for 2016 season is higher than that calculated for 2010, its value can be considered acceptable likewise. Both RMSE values turned in fact to be lower than the reference unit (1°C) chosen for the sizing of the thermal plume (16°C with undisturbed temperature of the aquifer equal fixed equal to 15°C).
3.2. Numerical Simulation Results
By working on the 2010 and 2016 heat transport simulations, total cooling thermal load and maximum values of thermal plume and were calculated at the end of each analysed season.
The values of thermal load obtained at the end of the simulations (Table 3) are comparable with those analytically estimated using Equation (12), which confirms the reliability of the elaborated model. Thermal plume dimensions were determined considering the isotherm of alteration equal to 1°C (isotherm 16°C) compared to the initial temperature of 15°C (natural groundwater temperature): was calculated as the downstream distance (m) from the rejection well; was defined as the maximum extension of the 16°C isotherm perpendicular to the groundwater flow direction (Figure 6).
Detailed results of the maximum dimensions of the thermally affected zone after seasons 2010 and 2016 are represented in Figure 7 and reported in Table 4.
As can be also seen in Figures 8 and 9, the TAZ reached its maximum dimensions at a depth between 30 and 32 meters for both analysed seasons, with a maximum simulated length of the plume () along the groundwater flow direction of 166.69 and 141.79 meters, respectively.
3.3. Comparison with Banks (2011) Analytical Solutions
Subsequently, using the recalculated hourly flow rate data distribution [m3day-1], the cumulative distribution function (CDF) was performed (Figure 10).
(a)
(b)
For both simulation periods, thermal plume dimensions were compared to those obtained using the reinjection steady discharge rate value () with values associated to a CDF of 63% (63% of chance that a randomly selected data point from the hourly dataset is lower than the corresponding selected value) (Figure 10).
By means of error analysis, a CDF of 63% was identified as a value, valid for the combined analysis of both datasets, which allows for the estimation of the thermal plume width by minimizing the dimensioning errors (for detailed results, see Table 5).
As reported in Tables 5 and 6, Banks’s [10] analytical formula (Equation (7) and Equation (9)) was able to provide values relating to TAZ dimensioning at the end of the 2010 plant functioning season with a good degree of accuracy: both TAZ dimensions are in fact conservatively overestimated by fixing the CDF value equal to 63%.
In contrast, for the 2016 cooling season, despite the fact that thermal plume width and length turned to be slightly underestimated, dimensioning error was limited to about 20 meters for both dimensions (for detailed results, see Tables 5 and 6).
4. Conclusions
Thermal plumes associated with open-loop geothermal systems must be well predicted and constantly controlled, especially in the shallow aquifers of more densely urbanized areas, not only to guarantee long-term sustainable use of plants but also to avoid adverse effects on adjacent and neighbouring geothermal systems.
For the small geothermal plant of the Politecnico di Torino, the analytical solutions proposed by Banks [10] was found to be a possible alternative to the use of more complex numerical simulation software for TAZ sizing at the end of each operating season.
In this paper, to verify the reliability of the analytical solutions proposed by Banks [10] to predict the thermal perturbation, we reconstructed an equivalent stationary state, starting from parameters that, for the Politecnico di Torino open-loop GWHP plant, are time variable.
In detail, we recalculated equivalent hourly flow rate data [m3day-1] by using measured thermal load values () and fixing a constant temperature difference () between the injection and abstraction wells at 5°C. Obtained values were then analysed by plotting cumulative distribution functions (CDFs).
For this particular case study, results of TAZ dimensions of the Politecnico di Torino GWHP system, obtained by using a value [m3day-1] associated to a cumulative frequency of 63% in Banks’s [10] analytical solutions, were found to be overall comparable with values modelled in transient mode by using the finite-element code FEFLOW® 6.2 [11].
If correctly validated for (1) hydrogeological systems characterised by geological and hydrogeological features comparable to those of the case study analysed and (2) geothermal plants with characteristics similar to the one considered, the proposed approach could represent a useful tool to preliminarily predict TAZ dimensions for cases in which hourly flow rate data are available, but it is not possible to develop any simulation models.
Data Availability
The hourly temperature and flow rate data used to support the findings of this study have not been made directly available because they are ownership of Politecnico di Torino. However, they are reported as graphs.
Conflicts of Interest
The authors declare that they have no conflicts of interest.