Thermal modeling of a Swiss urban aquifer and implications for geothermal heat pump systems

The progressive electrification of the building conditioning sector in recent years has greatly contributed to reducing greenhouse gas emissions by using renewable energy sources, particularly shallow geothermal energy. This energy can be exploited through open and closed shallow geothermal systems (SGS), and their performances greatly depend on the ground/groundwater temperature, which can be affected by both natural and anthropogenic phenomena. The present study proposes an approach to characterize aquifers affected by high SGS exploitation (not simulated in this work). Characterization of the potential hydro/thermogeological natural state is necessary to understand the regional flow and heat transport, and to identify local thermal anomalies. Passive microseismic and groundwater monitoring were used to assess the shape and thermal status of the aquifer; numerical modeling in both steady-state and transient conditions allowed understanding of the flow and heat transport patterns. Two significant thermal anomalies were detected in a fluvio-glacial aquifer in southern Switzerland, one created by river water exfiltration and one of anthropogenic nature. A favorable time lag of 110 days between river and groundwater temperature and an urban hot plume produced by underground structures were observed. These thermal anomalies greatly affect the local thermal status of the aquifer and consequently the design and efficiency of current and future SGS. Results show that the correct characterization of the natural thermo-hydrogeological status of an aquifer is a fundamental basis for determining the impact of boundary conditions and to provide initial conditions required to perform reliable local thermal sustainability assessments, especially where high SGS exploitation occurs.


Introduction
Shallow geothermal systems (SGS) are a type of technology that is becoming increasingly used to provide building conditioning across Europe, by exploiting solar heat stored in the shallow portion of the subsurface: from approximately 10 m depth, ground temperature is stable during the year and is mainly affected by elevation and latitude (Banks 2012). Solar heat is stored both in the ground and in groundwater and SGS can exploit this heat by using different methods, depending on the use of predominant conduction or advection processes. The heat can be tapped by ground source heat pump systems (GSHP) through conduction by commonly using 100 m deep borehole heat exchangers (BHEs), which are boreholes with an installed and grouted heat-exchanger where a thermo-vector fluid absorbs or release heat in the ground in a closed-loop (conduction only) and transfers it to the heat pump. Heat can also be exploited by advection through the use of groundwater heat pump systems (GWHP), which are constituted of a shallow production well in a productive aquifer where a hydraulic submersible pump is installed. Groundwater is therefore tapped from the aquifer and the heat is exchanged with the building through a heat pump: thermally altered groundwater is then reinjected or infiltrated back to the aquifer, commonly through a restitution well. Both GSHP and GWHP systems use a small portion of electricity, which in turn provides a large amount of thermal energy, in the range of 3-4 times that of the consumed electricity for GSHP systems and more than 4 times that for GWHP systems (Banks 2012). In the last few decades SGS has proved to be a reliable and efficient technology to provide heat/cooling to both residential (Qian et al. 2020) and nonresidential buildings (Spitler and Gehlin 2019) or even to industry/research facilities-e.g. ELI-NP "extreme light infrastructure" Magurele plant, Romania, with 1080 BHEs which fully cover heating and cooling demand (Romanian Geoexchange Society 2019)-with great thermal efficiency, low to nearly zero climate-change related emissions (Ahmadi et al. 2018). The installation of these systems has been proficient both in cold (Bakirci 2010;Zhen et al. 2017) and hot climates (Beckers et al. 2018;Roy et al. 2020), in both rocky (Han and Yu 2016;Liebel et al. 2012) and unconsolidated (Perego et al. 2016) geological frameworks and it was even coupled with other renewable technologies in order to provide the electric portion of energy without relying on the use of fossil fuels (Emmi et al. 2017), assuring a zero direct emission geothermal system. SGS advantages and the progressive electrification of the building conditioning sector explain why these systems are increasingly being installed all across the world. The phenomenon of high SGS density due to large deployment, especially in urban areas, is currently being observed in different locations across Europe (García-Gil et al. 2014;Attard et al. 2020;Vienken et al. 2019;Bakr et al. 2013). The topic is growing in importance, since the correct balance between the management of both the shallow geothermal resource and groundwater quantity/quality is becoming fundamental to ensure SGS sustainable operation in the long-term (García-Gil et al. 2020a). In particular, understanding heat transport in groundwater is essential for the preemptive design, performance analysis and impact assessment of SGS (Fujii et al. 2005) and therefore is extremely important to characterize the aquifer thermal regime for their correct management.
Scientific literature has lately been enriched with studies focused on the hydraulic and thermal modeling of aquifers and on assessment of their thermal conditions in order to understand which conditions are the main drivers for the sustainability of SGS, both for local or regional applications.
In particular Fujii et al. (2007) developed regional and local numerical models to provide the basis for estimating heat exchange performances of BHE; Händel et al. (2013) built a regional model to evaluate the sustainability of SGS exploitation in Leibnitz (Austria). Bridger and Allen (2014) reconstructed a numerical model to assess the influence of aquifer heterogeneity on the operation of an ATES system (Aquifer Thermal Energy Storage) while Sheldon et al. (2015) investigated the suitability and sustainability of groundwater extraction to cool a supercomputer by using numerical models. Boon et al. (2019) investigated the urban aquifer of Cardiff (UK) using a coupled hydro/thermal model to ensure the suitability for shallow geothermal exploitation. The work of Epting et al. (2013) highlighted the need to create tools to perform sustainable management of an urban aquifer and to assess the potential presence of a subsurface urban heat island (SUHI) effect (Zhu et al. 2015) compared to a "natural state", which is useful to understand the anthropogenic impact on groundwater quality, significantly produced by SGS. The common feature of all these scientific works is the detailed hydraulic and thermal characterization of regional aquifers and the study of their boundary conditions to properly simulate the influence of current or future SGS.
Cantone Ticino, a region in southern Switzerland, has witnessed in the past 20 years a constant growth of subsurface heat exploitation through the use of SGS, constituted by closed and open-loops (Perego et al. 2019; Repubblica e Cantone Ticino 2016). On the Cantonal territory, currently 677 open-loop doublet SGS are present, 389 of which are installed in the six main aquifers of Cantone Ticino, with an average density of 4 open-loop doublets/km 2 (eight perforations considering a representative scheme made up of a production and an injection well). More than 5,000 BHEs are also mapped, with their density in the six main aquifers reaching values of 10 BHE/km 2 . Short distances between adjoining SGS, both currently installed and planned, could influence ground temperatures and system performances in the long term. This phenomenon can be proficiently studied in the delta of river Maggia, a torrential regime river that flows into Lake Maggiore and hosts the cities of Locarno in the east and Ascona in the west. In this area there is a great density of both closed and open SGS (40 SGS/km 2 ), installed at relatively small distances (Fig. 1); however, no harmonic monitoring is currently adopted, with lack of detailed subsurface and SGS data. The absence of adequate monitoring data does not allow a better management of the shallow geothermal resource, which should achieve adequate systems efficiency with thermal sustainability and a harmonic decision process (Vienken et al. 2015).
To properly study the sustainability of SGS exploitation in the long term at urban scale, it is necessary to firstly characterize the aquifer thermal regime, understanding the influence of hydrogeological boundary conditions, which is the main objective of the present report. The proposed procedure involves hydrogeological characterization, groundwater monitoring and modeling, in order to shed light on the processes governing regional flow and heat transport in an aquifer affected by the presence of several SGS, through the characterization of the "potential natural state" (Epting et al. 2017a). Literature data, lowcost and noninvasive geophysical investigations were useful to describe the shape and thickness of the aquifer, while groundwater monitoring allowed characterization of its hydraulic and thermal regimes. All the collected information is then conveyed into numerical models useful to represent a "potential natural state" of the aquifer and to unravel the regional and local processes governing flow and heat transport, without considering anthropogenic factors. The steady-state and transient simulations of both piezometric levels and groundwater temperatures are useful tools to understand if the measured groundwater temperatures are consistent (or not) with a natural scenario. Moreover, they provide the starting point, represented by initial conditions, for future local-scale analyses through numerical models, which will be focused on assessing the impact of both closed and open-loop SGS towards ground temperature and to evaluate their mutual interferences. This last aspect is not included in the present work but will be the focus of further studies. A summary of abbreviations and symbols is given in the Appendix.

Study area
The region was subject to the presence of the Ticino and Maggia glaciers that during the Quaternary remodeled the deep valleys dug by pre-Quaternary fluvial erosion, lately filled by water of Lake Maggiore. The lake reached its maximum level approximately 20,000 years ago (Scapozza 2016), flooding the study area and reaching the convergence between stream Melezza and river Maggia (Fig. 1). The progressive recession of the lake was due to the delta progradation, with the continuous deposition of coarse sediments by both river Maggia and stream Melezza, characterized by a marked torrential regime. This torrential regime reflects on the heterogeneity and texture of the deposited sediments, mainly gravel, sand and a small portion of silt, where thick sandy lenses are often observed. Within this geological framework, gravel «channels» can be detected in the subsurface, representing the paleo-riverbeds of river Maggia before the river-training actions took place at the end of the nineteenth century, constraining the river in an artificial rectilinear channel. The thick sediments pack of fluvio-glacial origin lies over metamorphic bedrock. The area is located in proximity of the Insubric Line, which is an important tectonic division between the Central and Southern Alps.
To characterize the shallow portion of the subsurface of the study area, geological profiles were executed relying on borehole stratigraphy from GESPOS (SUPSI-IST) geodatabase (SUPSI -IST GESPOS 2020), which also contains hydrogeological information for the management of boreholes, wells and springs. From the analysis of the stratigraphic logs ( Fig. 1, lower part), it can be observed that the aquifer is composed of a pack of heterogeneous gravels and sands with no evidence of confinement. There is a continuous phreatic aquifer, constituted of gravels, gravelly sand, sandy gravel and sands, with local presence of silt in the northern portion of the delta. The major vertical discontinuity in this geological framework is the transition between the homogeneous fluvio-glacial cover and the underlying deep crystalline bedrock, composed mainly of metamorphic rocks as gneiss. This bedrock is observed at shallow depths only in the proximity of Mt. Verità (West) and near the reliefs, while it stands at greater depths all over the investigated area. The geophysical surveys performed by Fondazione Lerici of the Polytechnic of Milan in the 1970s (Fondazione Lerici 1970), consisting of electric and seismic surveys, were digitized and used to characterize the bedrock shape in the stream Melezza valley. These geophysical surveys were performed for different aquifers of Cantone Ticino in order to assess the thickness of groundwater bodies for drinking-water exploitation purposes. In this location the bedrock is found at great depths (maximum 250 m), while it rises to a depth of 50-100 m proceeding southward. Moreover, from Fondazione Lerici data, the bedrock cannot be detected anymore proceeding from the city of Locarno towards Lake Maggiore.

Passive microseismic
Fondazione Lerici surveys (Fondazione Lerici 1970) provided preliminary information on the thickness of the sediment pack in the northern portion of the study area. The continuous representation of the bedrock depth throughout the study area was therefore performed through the noninvasive technique called horizontal to vertical spectral ratio (HVSR), which allows for the collection of signals related to environmental micro-tremors or seismic "noise" (Nogoshi and Igarashi 1970). Environmental micro-tremors (seismic noise) are micro-movements of the ground (order of magnitude in the range of 10 −2 and 10 −6 mm) observed everywhere and consist mainly of Rayleigh (Rayleigh 1885) and Love (Love 1906;Bonnefoy-Claudet et al. 2008) waves produced by the interference between S and P waves at the surface. These microvibrations are produced by natural sources such as wind, tides, sea waves or climatic phenomena of regional interest, and anthropic noise-so-called "cultural noise" (Longuet-Higgins 1950). The core of the method is based on the study of the spectral ratio between the horizontal and vertical components of the seismic movement. Assuming horizontally polarized shear waves, the site resonance frequency represents the frequency at which the portion of investigated ground is subject to the maximum horizontal movement and to the minimum vertical movement. By creating a H/V spectral ratio, which is the ratio between the averaged horizontal component and the vertical component-HVSR technique (Nogoshi and Igarashi 1970;Nakamura 1989; Nakamura 2019)-the resonance frequency corresponds to the observed spectral peak.
In the presence of impedance contrast in the subsoil, e.g. sediment layer lying on stiff bedrock, the site resonance frequency f 0 can be calculated by the relation: where V s is the velocity of shear waves and H is the thickness of the soft top layer. This technique can identify the major resonance frequency of the shallower horizons, especially when considering a two-layer conceptual ground model, composed of an unconsolidated deposits cover above thick bedrock.
Having detailed information of the subsurface at a certain location (calibration point) allows for correlation of the resonance frequency with the depth at which the discontinuity is observed. Through this passage, it is possible to estimate the V s velocity of the sediments pack and therefore derive the discontinuity depth in an unknown location acquiring another signal. The theory, specific guidelines for the acquisition, elaboration and analysis of passive seismic signals are contained in the outcomes of the SESAME project (SESAME 2004).
Firstly, the GESPOS database was used in order to select stratigraphies from boreholes and geophysical surveys detecting bedrock. This allowed finding five calibration points to understand the range of V s velocity of the sediments pack, while other 17 locations were chosen as measurement points. A first measurement campaign was performed during 9-10 July 2018, a second during 18-19 September, 2018, and a third on 24 October 2018 and 15 November 2018. The locations mainly involved rural, leisure or agricultural areas in the northern portion of the investigated area, while they mainly involved urbanized areas in the southern portion.
Passive microseismic signal acquisition was performed by using simultaneously two tomographs, GEMINI and TROMINO devices, to compare results of both instruments and to more accurately characterize the investigated location. For these devices, a sampling frequency of 200 Hz and acquisition duration of 20 min were chosen. A total of 70 recordings in 22 locations were collected. For an 11-km 2 wide area such as the investigated one, this means an average of two locations per km 2 and six measurements per km 2 . Passive seismic measurements were performed even where the depth of the bedrock was known to be deep enough to be detectable, for example in Lido Locarno and Bosco Isolino; this however was an intrinsic result, since it offered the ability to understand where the discontinuity between sediments and bedrock could not be detected anymore.
The acquired signals from both GEMINI and TROMINO devices were analysed with both WINMASW (Eliosoft 2012) licensed software and with the freeware GEOPSY (Wathelet et al. 2020): an example of the acquired horizontal and vertical components of the passive seismic signal and of the H/V spectra for both software are reported in Fig. 2. It is possible to see in green, blue and red the acquired signals respectively in the north-south, east-west and up-down directions. The "eye-shaped" spectrum between approximately 1 and 4 Hz testifies a clear and ideal result of a good signal acquisition and processing with the presence of a clear discontinuity between the pack of sediments and the underlying bedrock.
The signal processing was performed by using the Konno-Ohmachi (Konno and Ohmachi 1998) smoothing technique with a constant of 40 and the squared average of horizontal components. For some points, the values are averaged over different measurement campaigns, since they were particularly affected by anthropic activities. The results of the signals analysis are reported in Table 1. V s velocities derived from the acquired signals in calibration points ranged between 170 and 790 m/s, which are respectively representative values for stiff and very dense sediments (FEMA P-1050 2015).
Maps of both peak frequency (f 0 ) and peak amplitude (A 0 ) were created and are respectively reported in Figs. 3 and 4. Figure 3 shows a range of f 0 between 0.6 and 13 Hz, with frequencies progressively increasing when proceeding southward, going from the deep valley dug by the glacier to the metamorphic rocks, mainly heterogeneous gneiss, constituting the nearby rocky reliefs. The increase in resonance frequencies testifies to the progressive thinning of the Quaternary sediments cover and the presence of shallower bedrock, which is in agreement with the current knowledge of the subsurface from Fondazione Lerici reports.
The A 0 map reported in Fig. 4 instead provides information of the weathering degree of the bedrock since the amplitude value increases when there is a strong impedance contrast between the sediments cover and the underlying bedrock (Galgaro et al. 2014). It can therefore be assumed that the bedrock beneath the glacial valley is more weathered than the one present near the rocky reliefs. In particular, a high A 0 value is observed in a field measurement right adjacent to Monte Verità, a hill composed of gneiss and partially covered with glacial deposits. After f 0 and A 0 isomaps reconstruction, a comparison between the reconstructed bedrock top surface with different V s and the reconstruction performed by Fondazione Lerici with active seismic was performed. A V s velocity of 430 m/s could adequately describe both the pit and the peaks reported in Fondazione Lerici reports. The hypothesized velocity of 430 m/s for the unconsolidated top cover, coherent with the investigated ground type-dense ground according to (FEMA P-1050 2015)-and the interpolated f 0 map was then used to reconstruct on a GIS platform the shape of the bedrock surface. In the southern and final part of the delta, the bedrock was reconstructed by hypothesizing constant dipping. The hypothesis was that the bedrock depth of Piano di Magadino (Fondazione Lerici 1970; Fig. 5), an adjacent glacial valley located eastern from the Maggia delta, could be comparable to the bedrock depth under the delta, since no discontinuity between sediments and bedrock was detected by the HVSR method; therefore, a similar steepness index was hypothesized. The reconstruction of the bedrock top was doubly useful-for the regional numerical model it represents the bottom edge and provides information on the thickness of sediments forming the Maggia delta aquifer.

Groundwater monitoring
To characterize the hydrogeological conditions across the selected study area, the first step was to create a groundwater monitoring network from scratch. The created network included a total of 16 piezometers/wells investigated with a wellcamera and evaluated as appropriate due to their good conditions-three for Melezza plain, six for the study of river Maggia/groundwater interactions and seven distributed along the delta. The eastern side of the delta (Locarno) is however the most investigated portion because it includes the future local-scale case study.
From the first piezometric level measurements, initially carried out between February and July 2019, it can be stated that the groundwater flow in the regional aquifer proceeds Horizontal to vertical spectral ratio (HVSR) in WINMASW; c HVSR in GEOPSY from the west to the south-east, with a steeper gradient in the Melezza aquifer and in the northern part of the Maggia delta which gently flattens when approaching Lake Maggiore. The water table has an elevation of approximately 235 m asl in the north-western portion of the study area and reaches elevations of approximately 193-194 m asl, when approaching Lake Maggiore in the south-east. The piezometric level on the delta is strongly affected by Lake Maggiore level and by its daily variations. During May-July 2019, spot groundwater temperature profiles were performed by using three different temperature sensors to assess average thermal conditions of the area, with a procedure similar to previous works (Farr et al. 2017). Temperature profiles (Fig. 6) testify to a great variability in groundwater temperature and a nearly complete absence of vertical temperature variations. Intuitively the shallower piezometers are affected by the thermal wave of air temperature, which propagates downward until a depth of approximately 10 m from ground surface, below which groundwater temperature reaches a constant value comparable to the mean annual air temperature (MAAT) of the location. Piezometers installed in the Melezza aquifer (409.2, 425.52 and 425.80) showed lower average groundwater temperature than the southern ones. Moreover deep piezometers located in the city of Locarno showed higher (423.268) and much higher (423.471) groundwater temperatures, probably related to anthropogenic alteration.
During the period October 2019-March 2020 a second measurement campaign was designed and executed to characterize the aquifer dynamic behaviour and to obtain observations to calibrate the transient model. The piezometric level and temperature measurements were executed monthly by using a PASI WMS-02 phreatimeter equipped with a PT100 temperature sensor, having ±0.1% accuracy.
The second set of temperature profiles confirmed the previous hypotheses, revealing that piezometer 425.52, located near stream Melezza, is strongly affected by river exfiltration and groundwater temperature is altered along the entire investigated depth. This piezometer also showed a significant vertical variability of groundwater temperature, probably highlighting a higher exfiltration velocity  (172, 173, 174 and 407), which seem to partially detect river exfiltration and for piezometers 423.471, demonstrating the presence of a highly-altered urban groundwater temperature, with values constantly between 18 and 20°C during the year. Shallow piezometers 423.189, 423.119 and 423.275 are, as expected, affected by the downward thermal propagation of air temperature delayed with depth, which seasonally contributes to increase groundwater temperature in autumn and decrease it in spring, as shown in the shift of the vertical profiles. The zone of seasonal fluctuation is greatly variable, depending on the analysed location and on the thermal diffusivity of the medium, but on average it has a depth between 9 and 12 m.
Thermal profiles were conveyed, by averaging the temperature value over the entire water column, into the reconstruction of the baseline preliminary groundwater temperature map reported in Fig. 7 for the May-July 2019 period. The created groundwater temperature map shows that on average the groundwater temperature ranges around 13.5-13.6°C, which corresponds to the MAAT in the investigated area.

Aquifer properties
The following step consisted of hydrogeologically characterizing the study area by subdividing it into macro-zones of hypothesized homogeneous properties, according to the hydrogeological map of Locarno (Carta idrogeologica del Cantone Ticino, Foglio Locarno 1989). The use of hypothesized homogeneous zones in numerical modeling has the advantage of a simpler calibration and shorter computational times; considering the aim of the proposed regional model, the hypothesized zones were considered as a reasonable approximation given also the almost undifferentiated lithologies observed in the subsurface down to 50 m depth and due to macro vertical homogeneity. To assign numerical values of hydraulic conductivities to the delimited zones, pumping test information was used. There are places where the information on hydraulic conductivity obtained by pumping tests was accurate and well distributed, while in other places (e.g. Solduno, Ascona, alluvial fans, southernmost part of the delta), there was a complete lack of any type of information. The hydraulic conductivity values of pumping tests were transformed in log 10 k, and averaged for each identified zone. Three different types of hypothesized homogeneous zones were identified, with four distinct values of hydraulic conductivity ( Table 2): The estimated hydraulic conductivity values, reported in Table 2, were used as initial values in the numerical model and adjusted during the calibration procedure.
The delta, as previously said, is composed of fluvio-glacial sediments, mainly gravelly sand and sandy gravel. The subsurface of the Maggia delta can be roughly considered, from a regional modeling standpoint, as thermally similar, assuming that gravels and sands could form an aquifer with thermally homogeneous characteristics.
In order to further characterize the subsurface of the investigated area, it was therefore decided to collect a shallow ground sample (4-5 m depth) near Lake Maggiore to identify the thermal conductivity of the sediments, which, as a rough approximation, could be generalized to the entire aquifer volume. Two buckets of sediments were collected and analysed: the sample was composed of gravelly sand, with a small component of silt, as can be seen in Fig. 8.
Thermal conductivity measurements were performed at FISTEC laboratory (Environmental Technical Physics Laboratory) of the University of Venice (Italy), by using an experimental method described in detail in Dalla Santa et al. (2017): this method allows taking into account the mineralogical composition variability and is based on a steady-state assumption, since it works by imposing and maintaining a constant temperature difference between three parallel plates. Two samples of the same material are placed in between: the volume of the sample-box allows involving a representative volume of unconsolidated material, including all the different lithologies and minerals. The temperature difference imposed between the hot and the cold plate was lowered from the standard value to 5°C, in order to minimize the convective contribution on the measured thermal values. Water was nebulized over the samples, in order to measure a dry and wet thermal conductivity.
Results of the experimental measurements are reported in Table 3. The experiment was replicated two times, in order to assess the variability of the results and to have a second set of values to compare against the first set. Results show that the measured thermal conductivity of the ground sample ranges between 0.19 W/mK in dry conditions to 0.69 W/mK in wet conditions. The resulting values are quite low, even for wet conditions, since, as reported in Dalla Santa et al. (2020), the recommended thermal conductivity value for wet heterometric gravel with sand should be around 1.08 W/ mK-in the regional numerical model, both steady-state and transient, a precautionary value of 0.7 W/mK for the thermal conductivity of the solid medium was therefore set. Total porosity of the samples was also extrapolated from thermal conductivity lab measurements. Knowing that the sample was analysed in a box 0.52 m wide, 0.52 m long and 0.06 m high, with a total volume of 0.016224 m 3 and knowing that the sample reached wet conditions when 3.4 L of water were poured into the box using a spray bottle, resulted in a total porosity of 0.2. The effective porosity, given the gravel and sand matrix, can be assumed as equal or slightly lower than 0.2, according to the literature (Heath 1983).

FE numerical model
All the previous information was implemented in the numerical model, which was firstly built in FEFLOW (Diersch 2014) with steady-state plurennial average conditions and then in transient conditions. The first solution was adopted to create a representation of the average hydrogeological and thermal conditions of the aquifer, while the transient model was created to understand time-dependencies of the analysed phenomena (Anderson et al. 2015). FEFLOW is a commercial numerical code designed for the simulation of flow, mass and contaminant or heat transport in groundwater, considering porous or fractured media in both saturated and unsaturated conditions. The finite element mesh had an extension of 11.2 km 2 ( Fig. 9) and was refined in proximity of the drinking water wells of regional importance and in proximity of the stream Melezza and river Maggia. In particular, around the drinking water wells, the mesh was optimized by creating a series of auxiliary points at a distance proportional on the water well radius, according to Diersch et al. (2010).
The vertical discretization was designed and applied in order to have a more discretized shallower portion of the model, since the totality of both the hydraulic head and temperature observations resided between 0 and 30 m depth and the deepest drinking water wells are drilled down to a maximum of 35 m. A digital elevation model (DEM) with a spatial resolution of 25 m was assigned to the top slice of the numerical model, while for the bottom slice elevation bedrock shape reconstruction described in section 'Passive microseismic' was assigned. The numerical model was then cut at a depth

Material properties
The groundwater recharge was calculated by considering the land use from the CORINE project (Büttner et al. 2004), and a division of the landscape into urban and nonurban was performed. The effective recharge rates were therefore set according to the land use-for non-urbanized zones the effective recharge was set to 1/3 of the mean annual precipitation rate, while for urbanized areas this factor was set to 1/30 according to a previous work (Epting et al. 2013). The average cumulative precipitation in the study area is 1,530 mm/year (Federal Office of Meteorology and Climatology 2020), which is 4.2 mm/day. The effective recharge for nonurbanized area was therefore estimated to be 510 mm/year (0.0014 m/day), while the effective recharge for urbanized areas was estimated to be 51 mm/ year (0.00014 m/day). These values were not further calibrated. The drainable/fillable porosity of the porous medium was set to 0.2 or 20%, according to thermal conductivity measurements. The initial value for the river Maggia and stream Melezza conductance (vertical hydraulic conductivity/ riverbed thickness) was set to 1.43 × 10 −5 1/s: this value was extrapolated from a numerical model created for another aquifer in Canton Ticino. The flow transfer rates of both streams were however calibrated in FePEST, with resulting different values for both Maggia and Melezza.
The imposed thermal conductivity of the solid matrix was 0.7 W/mK, while the thermal conductance of the entire model was set to 5 W/m 2 K. The use of an arbitrary value, equal for both rivers and for the air transfer condition, arises from the fact that both thermal conductivity and thickness of the clogging layers (the riverbeds) were not known. In the steady-state model this value has a relative importance (as effective porosity for example), while for the transient model, the heat transfer coefficients for Maggia, Melezza and air temperature were progressively calibrated to adequately fit field measurements. Heat porosity required by FEFLOW, which is the effective porosity of the medium, was set to 0.2.

Steady-state model
To build the steady-state regional model, time-averaged average hydraulic and thermal conditions were considered. The Fig. 8 a Heap of sediments located near Locarno lido, composed of gravelly sand/sandy gravel with a small component of silt, constituting the main lithology of the river Maggia delta; b samples of gravelly sands and sandy gravels of Maggia delta within buckets and the container for the measurement test; c TAURUS device, a hot guarded plate tool to measure the thermal conductivity of insulation materials. In this case, the device was modified on purpose to be capable of measuring the thermal conductivity of gravel A hydraulic head of 237.5 m asl was applied upstream at the western edge of the numerical model, with an average gradient of 0.0093. The eastern flow boundary condition corresponded to Lake Maggiore level; therefore, a constant head BC was set to the eastern edge of the model, with a 2016-2018 average hydrometric level of 193.5 m asl. No flow rate or hydrometric level measurements were available for stream Melezza, therefore a DEM of 2 m was used in order to assign the elevation of the hydrometric level, while its riverbed elevation was assumed to be 1 m below the water level elevation, since also no riverbed data were available. The continuous shape of the hydrometric level of river Maggia was created through the software HEC-RAS (USACE 2010), by implementing the riverbed shape through 38 sections and setting a flow rate of 20 m 3 /s, which is the plurennial average for 2016-2018. The continuous river hydrometric level was then obtained by inverse distance interpolation of the hydrometric level calculated for each section, from downstream (lake) to upstream (confluence of Maggia/Melezza).
Drinking water wells of regional importance were also implemented in the regional numerical model-they provide drinking water to Terre di Pedemonte, Ascona and Locarno municipalities and represent important factors in the groundwater budget. The characteristics of these water wells, including the screened depth, were taken from GESPOS database, while the declared realistic pumping rates applied to the numerical model were retrieved from previous works and reports (Studio Ammann 1991;Studio Baumer 2005).
The plurennial mean air temperature registered at the MeteoSwiss (Swiss Federal Office of Meteorology and Climatology) station of Locarno-Monti between January 2016 and December 2018 of 13.6°C was applied to the uppermost slice of the model as the third kind of boundary condition, while the plurennial average temperature of the shallow portion of Lake Maggiore (0-20 m depth) between January 2016 and December 2018 was equal to 13.7°C. A constant groundwater temperature of 13.6°C was applied as the first kind of boundary condition to the western edge of the model. The plurennial registered average temperatures of river Maggia and Melezza were respectively 11.7 and 11.5°C, assigned as the third kind of heat transfer boundary conditions (BCs).
The amount of geogenic flow under the Maggia delta region was taken from the heat flux map of Switzerland (Federal Office of Toporaphy Swisstopo 2011) and recorded as 0.055 W/m 2 ; it was applied upwards from the bottom slice. All the main boundary conditions implemented in the numerical model are reported in Fig. 10.
The numerical model was firstly run in steady-state conditions. Hydraulic head targets were all used within the calibration procedure performed with FePEST, the FEFLOW integrated parameter estimation code PEST (Doherty 2010), while 425.52 and 423.471 temperature measurements were Fig. 9 a Three-dimensional FE mesh, cut at a depth of 150 m from the ground surface; b mesh of the FE model, discretized around the major elements affecting the thermo-hydrogeological status of the aquifer excluded from the calibration since they were identified as outliers and do not represent an average thermal state of the groundwater; instead, they represent specific local phenomena/thermal anomalies. The groundwater temperatures used in the calibration process refer to an average temperature over the entire water column in order to have a representative value of the monitored site, since (as shown in Fig. 6) no significant vertical temperature variations were discovered. Shallower piezometers 423.189, 119 and 275 were also removed from thermal calibration, since they were influenced by seasonal air temperature fluctuations. The model was automatically run by FePEST 37 times in order to decrease the value of the objective function and to find the solution that best reproduced the measured head and temperature targets by minimizing the total error. The calibrated hydraulic conductivity values for the four zones were 504.5 m/day (zone 1), 192.7 m/day (zone 2), 117.4 m/day (zone 3) and 212.6 m/day (zone 4). The calibrated flow transfer rates for Maggia and Melezza were respectively 2.39 × 10 −6 and 4.7 × 10 −6 1/s. The heat transfer coefficients for streams Maggia and Melezza were respectively set to 0.5 and 5 W/m 2 K but they are relatively significant in a steady-state model. They were further calibrated in the transient model.

Transient model
The steady-state model was a valid preliminary tool to understand the average groundwater flow and to determine which were the main thermal elements influencing groundwater temperature. The steady-state model was also used to calculate the initial hydraulic heads and temperatures for the transient model, which can in turn help in identifying the major processes governing the interactions of surface/groundwater, providing further insights into the transient behaviour of the aquifer and understanding thermal influences of the rivers or Lake Maggiore towards groundwater temperature.
The three-dimensional (3D) structure of the regional transient model was maintained unchanged from the steady-state one, while transient BCs were applied to the regional model. In particular, the following BCs of the period 30 June 2015/31 March 2020 were set-effective recharge rate for both urban and non-urban area (monthly), daily Lake Maggiore hydrometric level, daily air temperature, monthly shallow Lake Maggiore temperature (0-20 m; courtesy of CNR-IRSA; National Research Council of Italy 2019), Maggia daily river temperature, Maggia monthly transient hydrometric levels A sensitivity analysis was performed by varying the heat conductance of the 3rd kind of boundary condition and thermal dispersivity of the medium. Heat conductances were varied between 0.1/1/5 W/m 2 K, and longitudinal and transversal thermal dispersivities were varied between 1/0.1, 5/0.5 and 10/1 m.
This was useful to understand that air heat conductance greatly influences the temperature of the shallow piezometers 119, 189 and 275, while thermal dispersivity greatly influences the simulated temperatures in piezometers located in the proximity of the river Maggia, on the eastern side of the delta. Reasonably low overall MAE and RMSE (respectively 1.2 and 1.3°C) were obtained by setting a longitudinal thermal dispersivity of 10 m and a transversal thermal dispersivity of 1 m: the scenario with thermal dispersivity values of 1/0.1 m led to MAE and RMSE increase of 14%, while setting the thermal dispersivity to 5/0.5 m did not produce a better model fit. Therefore, after manual calibration, the Melezza heat transfer boundary condition was deepened, the heat transfer rate for river Maggia was set to 0.5 W/m 2 K, heat transfer rate for air was set to 0.5 W/m 2 K and thermal dispersivity (longitudinal/transversal) was set to 10/ 1 m; this was revealed to be an adequate value for the simulated scale or advective transport distances, in this case equal to hundreds of meters (Zech et al. 2015).

Steady-state modeling
From the scatter plots of Fig. 11 it can be seen that the fit between measured and computed heads is good, with a MAE (Mean Absolute Error) of 0.50 m, a RMSE (root mean squared error) of 0.9 m and a R 2 of 0.99. The fit between the measured and computed temperatures without considering the two identified outliers is equally good, with a RMSE of 0.57°C.
A validation of the steady-state model created with February-July 2019 measurements was performed against groundwater levels and temperatures collected between October 2019 and March 2020. This was useful to evaluate the representativeness of the model not only of summer conditions, but during the year: the comparison is reported in Table 4.
The validation reports an average MAE and RMSE of respectively 0.94 and 1.13 m for the groundwater levels and of 0.62 and 0.83°C for groundwater temperatures without the two outliers, which can be seen as a good result considering the complexity of the hydrogeological setting, with the water table shaped by the presence of several surface water bodies, as reported in Fig. 12. Lake Maggiore drives the regional groundwater level, while the hydrodynamic flow is locally altered by the presence of stream Melezza and river Maggia.
In particular, river Maggia shows a draining behaviour in the northern part of its path and becomes mainly losing in its southern part, when reaching the lake. The local hydrodynamic flow near Locarno regional drinking water wells is affected by the losing behaviour of Maggia, which recharges the aquifer with its water.
The simulated groundwater temperature over the entire domain is reported in Fig. 13. The mean annual air temperature is the parameter that mostly affects the undisturbed groundwater temperature, since it oscillates between 12 and 14°C during the year, not considering places where an anthropogenic alteration is observed such as in the city of Locarno. The local losing behaviour of river Maggia is also collaborated by Fig. 11 Scatter plots of a heads and b temperatures after automatic calibration with FePEST. b Refers to the comparison of measured vs computed temperatures without points 423.471 and 425.52, which represent specific local phenomena and do not describe an "average" thermal status temperature measurements in piezometer 423.173 and clearly observed in Fig. 13; this behaviour is also probably facilitated by the pumping of Locarno drinking water wells, which seem to partially accelerate the natural exfiltration. This particular aspect is not caught by groundwater level and temperature measurements in the western side of the delta (Ascona municipality). On this side the influence of the river on groundwater temperature is minimal to negligible, since the groundwater flow is affected by the presence of the Mt. Verità rock relief which constrains the local groundwater flow direction. Therefore this asymmetric groundwater flow does not allow propagation towards the west of the exfiltrated river water, leaving the monitored piezometers unaffected by river temperature.  Figures 14 and 15 show the results of the transient regional model. Piezometers located within the Melezza Valley were not reported, since the calculated groundwater level from the model is less reliable in the northern portion of the study area due to the inaccurate estimates of the river hydrometric level and riverbed elevation. Simulated groundwater levels for the southern portion of the delta are in agreement with the measured ones (data available only between February 2019 and March 2020), with a global MAE and RMSE of respectively 0.43 and 0.47 m. Considering the percentage ratio between the RMSE (0.47 m) and the range between the maximum and minimum target elevation (200.62-193.2 = 7.38 m), a result of 6% is obtained, which is more than acceptable according to Anderson et al. (2015). Groundwater levels in piezometers 42 and 46 seem to be poorly affected by Maggia hydrometric level oscillations or by its temperature, confirming that on this side the river influence is less significant than in the left side, due to the asymmetric groundwater flow direction (south-east-east) constrained by Mount Verità relief. The groundwater levels in piezometers located near the lake (288,119,189,275) are strongly affected by the Lake Maggiore hydrometric level. Urban deeper piezometers 268 and 471 are also strongly affected by Lake Maggiore hydrometric level, but in this case the model tends to constantly underestimate the groundwater level, showing that the hydrodynamic field could be locally altered; other factors should therefore be considered, e.g. groundwater withdrawal and reinjection by large openloop SGSs operating for decades, which could alter both the groundwater level and temperature compared to a natural state. Figure 15 shows the simulated groundwater temperatures in the considered period. Temperature monitoring showed that piezometers located near river Maggia show different behaviours depending on their hydraulic side. Piezometers on the left side of the delta (considering Maggia flow direction) are affected by river Simulated urban groundwater temperatures are often lower (119, 189, 275) and much lower (471) than the measured ones, testifying that the baseline groundwater temperature under the city of Locarno is altered and is far from being only the result of natural processes. In undisturbed conditions, the groundwater temperature should be comparable to the MAAT: a significant anthropogenic contribution, unambiguous for piezometer 471 (18-20°C), can therefore be the explanation for these high groundwater temperatures. Figure 16 shows a cross-section, located as reported in Fig. 13, including stream Melezza and piezometer 425.52, which is the representation of yearly river-groundwater interactions. In particular the transient model adequately explains the measurements in piezometer 425.52, where a sinusoidal oscillation of great amplitude is observed during the year, with temperature values that range from 6 to 18°C. The section represents the exfiltration of river water, producing a thermal plume reaching piezometer 425.52. The comparison between measured and simulated temperatures in the piezometer qualitatively testifies to a very strong influence of the river towards groundwater, proved by the observed seasonal groundwater temperature oscillations even at great depth as shown in the thermal profile of Fig. 6. The river, standing at higher elevation in that location and showing a strongly losing behaviour, exfiltrates into the groundwater and locally alters its It was observed that, with a good approximation and considering that groundwater temperature is not a conservative parameter, the simulated temperatures and the measured ones are in agreement as observed in the graph of Fig. 16 (with a slight overestimate of lower spring groundwater temperatures) and they could be used as a tracer (Anderson 2005;Constantz 2008). The lag time between the maximum river temperature and the peak of groundwater temperature could be assessed through a crosscorrelation function (Lee et al. 2013) between the simulated river and groundwater temperature signals. To perform this correlation, simulated daily stream Melezza temperature and groundwater temperature in piezometer 425.52 were detrended and normalized, creating time series with amplitudes comprised between 0 and 1. The cross-correlation function returned statistically significant amplitude, with a maximum value of 0.94 after 110 days. Considering that the linear distance between the hypothesized exfiltration area of the river and the piezometer is estimated to be 300 m and considering that from crosscorrelation the lag between river temperature and groundwater temperature signals amounts to 110 days, the empirical velocity of the temperature signal should be approximately 2.7 m/day, not considering the effect of porosity and thermal properties of the sediments that retard the heat pulse.

Discussion
The steady-state and transient models, compared against real mea sur emen ts, g av e se ver al ind i ca tions on the hydrogeological setting of the investigated study area. Some discrepancies between simulated and measured temperatures can be observed, especially in piezometers close to the river Maggia, which can be explained by the fact that transient heat transport simulation requires detailed information on hydrogeological and thermal properties, which could be very heterogeneous, even locally. These parameters contribute to retarding the heat transport but are hard to determine in the field (Vogt et al. 2010); therefore, some values were taken from literature and may differ from real values.
However real measurements and simulations showed that Lake Maggiore level fluctuations have a major influence on the overall piezometric surface of the delta, since the measured groundwater levels are in complete agreement with lake hydrometric levels pattern, especially proceeding southward. The oscillations of Lake Maggiore hydrometric level therefore affect the hydraulic gradient of groundwater, resulting in modified groundwater velocities during the year which in turn could modify the shape of urban SGS thermal plumes, similar to previous works (García-Gil et al. 2014).
Stream Melezza locally influences the amount of heat in the aquifer, causing important oscillations in the observed proximal groundwater temperatures of piezometer 425.52, but not in the two other monitored piezometers (409.2, 425.80), suggesting a different river/groundwater connection. The great temperature oscillation is caused by the consistent exfiltration of stream Melezza into the groundwater, recordable even at 30 m depth. A preferential pathway for river exfiltration is observed, due to South from piezometer 425.52, there is an industrial area, with factories that in the future could potentially ask for the concession of new wells devoted to the thermal use of the subsurface. In this case, the designer should be very aware of the particular hydrogeological conditions and must dimension the open-loop system by taking into account the seasonality of groundwater temperature, far from being constant during the year and even characterized by a 110-day lag from river temperature. In this specific case, an open-loop system would find optimal conditions for an increased efficiency, since groundwater temperature would show peaks of 17°C in fall-winter, when heating needs are greater, and values of 7°C in spring-summer when cooling needs are dominant. Compared to an undisturbed groundwater temperature of 13.6°C (MAAT), this would mean an improved coefficient of performance (COP) in winter and an improved energy efficiency ratio (EER) in summer, which represents a synergic effect that could be profitably exploited (García-Gil et al. 2020b). For large energy demands, this would represent significant energy and economic savings; it must be however considered that the installation of large GWHP systems with significant pumping rates could contribute to altering the natural exfiltration velocity and therefore the arrival time of the hot/cold thermal plume.
The thermal influence of the river Maggia was revealed to be asymmetric towards the eastern part of the delta and this is currently explained by the particular shape of the aquifer, affected by the presence of Monte Verità on the right side which acts as a hydraulic barrier and conveys the groundwater flow towards the river Maggia, despite its otherwise losing behaviour observed on the eastern side of the delta.
In the city of Locarno, some of the measured and calculated groundwater temperatures could not be explained by simply considering natural effects such as air temperature, river or lake influence. There are regions where groundwater is at least constantly 1°C higher than in an undisturbed scenario (groundwater temperature = MAAT) and this testifies to the presence of a potential subsurface urban heat island effect (SUHI; Zhu et al. 2010;Huang et al. 2009) with underground structures releasing heat in the aquifer (Epting et al. 2017b). The thermal regime of the urban aquifer is strongly altered in specific locations, testifying to a significant presence of anthropogenic activities affecting groundwater temperature.
After detailed investigations, it was discovered that the anomalous registered high temperatures of 18-20°C observed in piezometer 423.471 throughout the analysed period could be partially related to the infiltration through wells of the water drained from a nearby underground tunnel. In this tunnel, the groundwater is pumped, and used for fire-extinguishing and cleaning purposes, or to cool electro-mechanical devices. The infiltration rate of these wells, placed in the evacuation duct of the tunnel, was investigated through a saline injection test using a SalinoMADD device, discovering a flow rate of 1,500 L/ min. The water is infiltrated constantly during the year at a temperature of 16-18.5°C and seems to reach the monitored piezometer 423.471, partially justifying the anomalously high temperature. However, the higher temperature observed in 423.471 is still not completely explained: other large SGS are present between the tunnel and the piezometer which could further alter groundwater temperature. Here the measured groundwater temperature, as shown in the thermal profiles, reaches values between 18 and 20°C, approximately 4.5-6.5°C higher than the MAAT (13.6°C); therefore, testifying to a significant thermal use of the subsurface, particularly by SGS. This results in an increased amount of available heat in winter, but also results in excessive heat in summer. This aspect is extremely important because it affects the production temperatures and therefore efficiencies of the installed SGS: installing a new open-loop SGS for heating purposes in a location where groundwater temperature is constantly at 18-20°C instead of 13.6°C could increase the system efficiency, while installing it for cooling purposes would significantly decrease its efficiency, forcing the designer to choose another heat sink rather than groundwater. For large SGS with great energy demand, this would imply significant economic benefits for heating and severe issues for systems designed to provide cool.

Conclusions
This work presents an approach to hydrogeologically and thermally characterizing the potential natural state of an aquifer affected by the presence of several SGS (not simulated in the present work), in order to understand the natural processes governing the flow and heat transport field, and to assessing the influence of hydraulic and thermal boundary conditions on the aquifer's hydrogeological natural thermal regime. The method was applied to a fluvio-glacial aquifer located in southern Switzerland: through the use of both literature (in terms of active seismic surveys) and collected field data, with the support of 70 experimental passive microseismic surveys, the shape and the litho-textural setting of the aquifer were characterized. A groundwater network, composed of 16 piezometers/ wells, was designed to monitor the hydrogeological and thermal status of the aquifer. Moreover, two buckets of sediments were collected and analysed in the laboratory to estimate their thermal conductivity by using an experimental device, designed for coarse textural terms, which considers a representative sample. The collected data were conveyed into steady-state and transient numerical models useful to understand both the static and dynamic hydrogeological/thermal regime.
Both real measurements and numerical modeling revealed very different behaviours of the aquifer depending on the analysed location. In particular, the northern portion of the study area is affected by the strong exfiltration of stream Melezza, which produces a delayed hot/cold thermal plume in the groundwater, depending on the season. This thermal alteration of groundwater in the range of ±6°C could positively affect the design and installation of SGS downgradient depending on their distance, given the 110-day delay of groundwater temperature from river (and air) temperature: this would strongly improve the efficiency of future systems since the amount of heat stored in the groundwater is higher during the buildings' heating season and lower during the buildings' cooling season.
Altered groundwater temperatures in the range of 18-20°C were also detected in city of Locarno subsurface: further investigations revealed that this anomalous temperature could be partially related to the presence of an underground anthropic structure that infiltrates groundwater at higher temperatures during the year, affecting locally, along with SGS, the aquifer thermal regime. The greater groundwater temperature would result in better SGS heating performances in winter, but it would result in an efficiency decrease for a SGS devoted to cooling purposes. This of course would completely change the design solutions and the efficiency of future SGS, since different heat sinks rather than groundwater are needed.
The proposed approach should be identified as a good practice, consisting of the preliminary assessment of the hydrothermal status of the aquifer, which greatly influences local groundwater conditions. As a matter of fact, a proper assessment of SGS operation on both subsurface temperature and their mutual interferences at local scale cannot be performed without knowing the broader hydrogeological setting in which they are installed and the interconnections between parameters, namely a "holistic" approach. This would contribute to optimizing the design and sizing of these systems, helping the authorities responsible for SGS authorization to define sustainable management plans from an environmental and energy standpoint, given that potential positive or negative thermal interference between SGS could also be increased or mitigated by the regional hydrogeological setting. This proposed approach would be useful not only in countries with a current high density of SGS, but also in countries with low SGS penetration to preemptively optimize the allocation of the shallow geothermal resource, since the current "EoE" (electrification of everything) trend is leading to the recurring use of heat pumps for building conditioning, especially based on SGE. statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.