Case Studies of Geothermal System Response to Perturbations in Groundwater Flow and Thermal Regimes

Global demands for energy‐efficient heating and cooling systems coupled with rising commitments toward net zero emissions is resulting in wide deployment of shallow geothermal systems, typically installed to a depth of 100 to 200 m, and in the continued growth of the global ground source heat pump (GSHP) market. Ground coupled heat pump (GCHP) systems take up to 85% of the global GSHP market. With increasing deployment of GCHP systems in urban areas coping with limited regulations, there is growing potential and risk for these systems to impact the subsurface thermal regime and to interact with each other or with nearby heat‐sensitive subsurface infrastructure. In this paper, we present three numerical modeling case studies, from the UK and Canada, which examine GCHP systems' response to perturbation of the wider hydrogeological and thermal regimes. The studies demonstrate how GCHP systems can be impacted by external influences and perturbations arising from subsurface activities that change the thermal and hydraulic regimes in the area surrounding these systems. Additional subsurface heat loads near existing schemes are found to have varied impacts on system efficiency with reduction ranging from <1% to 8%, while changes in groundwater flow rates (due to a nearby groundwater abstraction) reduced the effective thermal conductivity at the study site by 13%. The findings support the argument in favor of regulation of GCHP systems or, to a minimum, their registration with records of locations and approximate heat pump capacity—even though these systems do not abstract/inject groundwater.


Introduction
Ground source heat pump (GSHP) systems provide an efficient and clean technology for heating and cooling of buildings in the worldwide energy market using a renewable energy resource. With increasing deployment of these systems in the subsurface of urban areas, there is growing risk for these systems to impact the subsurface thermal regime and interact with other heat-sensitive subsurface infrastructure, such as tunnels, building foundations, or other shallow energy systems including underground thermal energy storage (Bidarmaghz et al. 2019(Bidarmaghz et al. , 2020. This impact can be positive or negative, depending on the system operating mode and the type of thermal interference. When used in heating mode, GSHPs extract heat from the ground, that is, they use the heat resource that accumulates in the urban subsurface, partly due to anthropogenic activities. Such heat accumulation (also referred to as the Subsurface Urban Heat Island-SUHI) is a widely observed phenomenon, which increases the urban technical potential of geothermal use by up to 40% when compared to rural conditions (Rivera et al. 2017). The SUHI phenomenon is attributed to land use changes associated with urbanization, specifically heat losses from building basements (Ferguson and Woodbury 2005), pavements (Taylor and Stefan 2009), and buried infrastructure, such as tunnels or sewers (Menberg et al. 2013a), which have resulted in elevated groundwater temperatures beneath many cities (Taniguchi et al. 2005;Banks et al. 2009;Headon et al. 2009;Farr et al. 2017;Rivera et al. 2017;Zhu et al. 2017) causing temperature perturbations to depths of 100 m or more. When used in cooling mode, ground source heat pump systems reject heat to the subsurface, hence have considerable potential to contribute to and enhance the SUHI effect.
In either case, these systems benefit from the thermal inertia and storage capacity of the subsurface, which permits its use for both heating and cooling. However, there is a fragile equilibrium between the heat pump system's thermal loads and the rate of thermal renewal in the subsurface. This equilibrium needs to be maintained over the life of the system to ensure sufficient energy savings. However, many factors can affect this thermal equilibrium, such as unbalanced ground loads, groundwater flow and interference with other energy systems or subsurface infrastructure.
There are three main types of GSHP systems (Banks 2012;Self et al. 2013;Dowling et al. 2016): (1) ground-coupled heat pump (GCHP) or closed-loop systems which use a borehole heat exchanger (BHE) installed in the subsurface through which a heat exchanger fluid is circulated, (2) groundwater heat pump (GWHP) or open-loop systems which use pumped groundwater for the heat exchange, for example, via an intermediate plate heat exchanger installed in building at the surface, and (3) surface water sourced open and closed loop systems. This paper focuses on GCHP systems becoming one of the most widely used, currently occupying 85% of the ground source heat pump market share worldwide (Gupta and Singh Bais 2018), due in part to their reduced potential environmental impact. Regulations for these systems vary between countries, ranging from (1) no regulations (systems do not require any permits or registration), to (2) notification schemes (no permits but need to be registered/reported to authorities) and (3) permitting schemes (systems require relevant permits). A review of guidelines can be found in Dehkordi and Schincariol (2014) and Haehnlein et al. (2010). In Canada, regulations exist in provinces to protect groundwater resources from potential threats that can occur during geothermal system installation and operation. There is, however, no regulation to help maintain the thermal equilibrium of the subsurface. In the UK, GCHP systems are unregulated, requiring neither a permit nor registration.
The design and operation of GCHP systems is commonly based on the assumption of conductive heat transfer in the subsurface (Bernier 2001), even though thermal advection (i.e., via groundwater flow) can have a significant impact on the subsurface thermal equilibrium and the long-term performance of the systems, especially when the Darcy flux is greater than 1 × 10 −7 m/s (Dehkordi and Schincariol 2014;Ferguson 2015). While improving long-term performance of systems by dissipating heat/ cold injected into the ground (Raymond et al. 2011;Zanchini et al. 2012), it also enlarges the system's footprint as well as its sensitivity to far-field boundary conditions. Numerical tools are available to optimize the operation of GCHPs under the influence of groundwater flow (Fujii et al. 2005), but estimation of site-specific groundwater fluxes can be difficult to define accurately. Thermal response tests (TRTs) can be performed on a single BHE to evaluate an effective subsurface thermal conductivity affected by the groundwater flow (Signorelli et al. 2007;Bozdag et al. 2008). Peclet number analysis made with downhole temperature measurements during a TRT can further help distinguish between conductive and advective heat transfer (Koubikana Pambou et al. 2019). This approach can be useful, but neglects the fact that flow conditions, and hence effective thermal properties, can change over time (Abesser et al. 2020).
GCHPs have been studied in urban environments with respect to their effects on the subsurface thermal regime (Rivera et al. 2015), and the impacts of a building and groundcover on a single BHE system (Rivera et al. 2016) have been assessed. Various modeling studies (e.g., Hecht-Méndez et al. 2013;Casasso and Sethi 2014;Hein et al. 2016) have investigated system sensitivity to key hydrogeological and operational parameters to identify the controls on GCHP functional efficiency. The focus of these studies has been on isolated systems, where flow conditions and background subsurface temperature are assumed to be constant, impacted by the modeled system only during its operation. Extensive monitoring and modeling studies have been undertaken on large BHE systems at the EPIC systems site and the Ball State University system to investigate internal interferences between BHEs within the same design field (Florea et al. 2017). However, less attention has been paid to the effects on GCHP functional efficiency from external influences, such as perturbations in the wider hydrogeological and thermal regime, for example, due to urbanization, groundwater abstraction, multiple BHEs within tight residential clusters or competing subsurface uses.
This paper details three modeling-based cases studies that investigate the changes in the performance of typical GCHP installations (different designs and operational pattern) in response to perturbations in the hydrogeological and/or thermal regimes. The specific modeling objectives vary for the different case studies, but the overall aims of this paper are to compare (1) GCHP systems' response to changing state or process variables within different hydrogeological and thermal systems, and (2) the impact of interferences with other subsurface uses on GCHP functional efficiency. In doing so, we will identify general factors that need to be considered in the planning and design of different, potentially competing, subsurface uses.

Research Methodology
Modeling within all three case studies is performed using FEFLOW ® , a three-dimensional finite-element, fully coupled variable-density groundwater flow and transport code. FEFLOW offers different approaches for simulating heat transport around the GCHPs (Diersch et al. 2010;Diersch et al. 2011) through implementation of the BHE: (1) via a Heat Nodal Sink/Source Boundary Condition within a fully discretized two-dimensional (FD2DM) or three-dimensional model (FD3DM) (this approach simulates BHE thermal exchange with the surrounding soil/rock, while thermal transfers within the BHE configuration are not explicitly considered); (2) by discretizing all borehole elements and assigning flow and thermal material properties on a nodal/element basis in a FD3DM; or (3) via built-in modules, based on numerical (Al-Khoury et al. 2005;Al-Khoury and Bonnier 2006) or analytical (Eskilson and Claesson 1988) methods, where the BHE is represented by a simplified one-dimensional (1D) element, inserted at the center node of the BHE and coupled with the rest of the model domain. FEFLOW solves the governing flow and heat transport equations for the area surrounding the BHE; a BHE solution is coupled with the rest of the model domain through the temperatures at borehole nodes. Modeling Studies 1 and 3 applied the built-in, discrete-element BHE solution (Approach 3), while Modeling Study 2 simulated heat exchange via a nodal boundary condition (Approach 1).

Modeling Study I: University of Western Ontario (UWO) Campus, Canada Objectives
This case study has three main objectives: (1) to assess how a functioning GCHP system, serving a small portion of a building on the University of Western Ontario (UWO) campus, could be expanded within the space available between buildings, (2) to investigate the effects of a future upgrading of a BHE field installation on the efficiency of the existing system, and (3) to assess the importance of fully accounting for near surface thermal disturbances in the modeling process.

Study Site
The study site is a 450 m by 250 m area aligned with regional groundwater flow toward the south (Figure 1).
The site contains two active vertical BHEs (90 m) and two horizontal ground heat exchangers (Figure 1). Three monitoring boreholes with thermistors at 30 m, 45 m, 60 m, 75 m, and 90 m depth are adjacent to the vertical BHE. The vertical BHEs and monitoring boreholes extend to a depth of 90 m, through 34 m of glacial till and into Paleozoic limestone and dolostone formations (Armstrong and Carter 2010). The upper portion of the till is clayey silt, stone poor, has a relatively low hydraulic conductivity, and often acts as a confining layer or aquitard (Matrix Solutions Inc. 2014). While the underlying silty-sandy till can locally act as an aquifer (Schwartz 1974), regionally it is considered an aquitard. The upper few meters of weathered bedrock surface, where fractured and/or karstic, is considered an aquifer (Schlumberger 2011;Matrix Solutions Inc. 2014). Overall, the limestone and dolostone members vary from fossiliferous to crystalline, and massive to bedded; generally, they can be considered aquifers (Matrix Solutions Inc. 2014).
The pipe dimensions and operational requirements for the vertical geothermal system can be found in Table S1. The functioning BHE system operates in conjunction with a shallow horizontal geothermal system, following an alternating 7-day cycle. Most BHE systems do not operate on an intermittent cycle. The BHE system is nearly balanced with six cooling months (May to October), two transitionary months (April and November) where the system may alternate between heating and cooling, and four heating months.

Methodology
A 3D model was developed of the study site. The model hydraulic head boundary conditions ( Figure 1) were determined from overburden and bedrock aquifer potentiometric maps (Matrix Solutions Inc. 2014). For steady-state and transient simulations, the lateral and basal boundary conditions remained constant. The temperature boundary at the model base (9.4 • C) was derived from a geothermal study performed by Judge (1972). Temperatures at depths of 200 m are relatively stable, and unaffected by climate shifts within the last 200 years (Pollack and Huang 2000; Kukkonen et al. 2011). The ground surface boundary condition was derived from measured air temperatures using relations developed by Taylor and Stefan (2009). For the steady-state simulations, the surface boundary was set at a constant 10.1 • C for grass and 13.2 • C for asphalt and concrete. As discussed later, the transient model fully accounted for monthly fluctuations in these surface temperatures. At the lateral boundaries, a zero-gradient (adiabatic) temperature condition was used. Buildings and sewer systems were not included in the steady-state model but were added to a transient model as part of the spin-up process.
Physical model properties (Table 1) were estimated from regional studies (Matrix Solutions Inc. 2014) in correlation with the site BHE borehole logs, and borehole logs (Judge and Beck 1967). Thermal conductivities for the bedrock units were measured by Judge (1972), overburden values were estimated from Banks (2008). Volumetric heat capacity values were estimated from Banks (2008). Porosity was measured by Judge (1972).
First, an expansion to a BHE field with a spacing of 10 m (18 BHEs) and a field with a spacing of 5 m (69 BHEs) were assessed. Second, to investigate potential upgradient influences, a similarly sized 18 BHE field was added adjacent to a nearby building (Figure 1). Energy for BHE systems is injected for cooling and extracted for heating. In balanced systems, the energy difference between injection and extraction is close to 0. This study defines the total energy exchanged as the absolute value of the sum of injected and extracted energy. It was used to assess the energy that the systems were able to produce over a lifespan of 20 years.

Model Input and Spin-up Process
For the steady-state model, initial modeling attempts using present-day infrastructure resulted in greater heat accumulation in the upper 80 m than shown from the monitoring borehole data prior to BHE activation. It was determined that the model spin-up (i.e., a set of repeated runs to determine the initial model conditions that best represent the system's thermal-dynamic balance) needed to be completed in phases. There were three main phases; (1) the initial steady-state spin-up, (2) a multistep transient spin-up, and (3) the final working model which would serve as initial conditions for predictive models. The transient model spin-up was started in 1942, when little infrastructure was present at the site, and moved through 12 phases, each bringing in buildings, sewage systems, roads, and parking lots as they appeared in the aerial photo and building records until 2011. Air temperature climate records were used to adjust the average annual temperature for grass cover at the start (9.1 • C, 1942) to the end of the multi-step transient spinup (August 30, 2011, 10.1 • C). Building basements were represented by a 20 • C temperature boundary condition (Ferguson and Woodbury 2004;Menberg et al. 2013b) set at a depth of 2.5 m. The heat from the basement walls are not represented as Thomas and Rees (1998) and Emery et al. (2007) showed that heat loss through basement walls was mostly connected to the atmosphere. Asphalt, concrete and grass cover temperatures were represented by their respective annual average values except for the final 10 years of spin up when the average monthly temperatures were used. A temperature of 18.5 • C was used for sewage pipe temperatures which correlated with nearby sewage treatment facility data and Menberg et al. (2013b). Mesh spacing was optimized around the BHE at 0.46 m, following Diersch et al. (2010), and increased laterally. Vertical discretization was 0.1 m for the first 1.0 m, and then followed at 0.25 m, 0.5 m, 1 m, and 5 m for depths up to 3.5 m, 6 m, 10 m, and 200 m. Mesh sensitivity analyses were completed to ensure that the thermal transport solution was mesh-independent, that is, not influenced by further discretization. The functioning BHE system was incorporated in the model as discrete linear elements (Diersch et al. 2010) representing the vertical U-tube. BHE inlet temperatures and flow rates varied depending on operational cycles and were recorded and implemented in the model calibration and steady-state initial conditions phases. An average inlet temperature was applied for each month (Table S2). The shallow horizontal geothermal system was represented as a specified temperature boundary condition matching operating cycles and temperatures.

Results
The modeled thermal profile of the subsurface prior to geothermal system activation compares well to the monitoring data except for a minor divergence (maximum 0.3 • C) centered around 60 m depth ( Figure S1). The cause of the temperature difference may be due to localized groundwater flow in fractures in the limestone, or a shifting of thermistors during installation. The geoexchange systems were then added to the model. Field data from the monitoring boreholes were used to further calibrate the BHE model.
The average annual energy exchanged, over a 20-year simulation period, for the BHE field scenarios is shown in Table 2. By comparison, the active two BHE systems exchanged an average of 46 MWh per year. Expanding the system from 2 to 18 (10 m spacing) BHEs resulted in a loss of efficiency of 3% (energy exchanged per BHE), while increasing the density to 69 (5 m spacing) BHEs increased this loss to 6.9%. The addition of a similarly sized (18 BHE) upgradient installation (Figure 1) only had a minor effect (0.3% loss in average annual energy exchange). The depths to which the thermal perturbations from the building, asphalt, and grass cover extend are clearly seen in the thermal difference plot ( Figure 2). The effect of the upgradient system and groundwater flow is seen in the thermal difference plots after 5 and 20 years of operation ( Figure 3). Model simulations where the effects of infrastructure were removed by conducting the model spin-up with only grass conditions on surface, showed a small increase (2.2%) in annual energy exchange ( Table 2). The small change in energy exchange was increased to 3.1% when proper accounting for the unsaturated zone through the application of the Richards Equation in FEFLOW was removed (i.e., phreatic option). Finally, when groundwater flow was set to a zero-gradient (i.e., no flow), a much larger decrease in energy exchange (15.8%) was noticed.

Discussion
The current BHE system is expandable within the tight inter-building space with little loss in efficiency per borehole. This in itself is a routine investigative outcome. An equivalent BHE system operating 100 m upgradient was also found to have minimal impact on the downgradient field. This also is an expected outcome as the depth-averaged specific discharge across the BHE field is approximately 7 × 10 −8 m/s, which by using the screening tool developed by Ferguson (2015), puts the system into the boundary area where advective effects become more important over conduction. It also correlates with the findings of Dehkordi and Schincariol (2014) who found that groundwater influence on ground loop temperatures becomes significant at fluxes of approximately 10 −7 m/s and higher.
A noteworthy finding of this study is the importance of applying the correct initial conditions by assessing the appropriate level of model spin-up in relation to BHE functional efficiency. Here, to adequately match nearsurface temperatures during initial model calibration, a multi-step spin-up process over a 70-year period of infrastructure development was required. However, removing infrastructure effects and using a simple unimpacted subsurface temperature distribution only affected BHE energy exchange by approximately 2%. More significant was accounting for the effect of the unsaturated zone on thermal transport (3%). Finally, for the site conditions, removing groundwater flow effects had the most significant impact on the BHE energy exchange, reducing it by 16%.
Thermal impacts from infrastructure are known to extend over 100 m deep as shown by Ferguson and Woodbury (2004) and this study ( Figure 2). However, the effects of infrastructure appear to impact minimally BHE energy exchange over 90 m borehole depth which is typical for these systems. Additional simulations reducing BHE depth by 50% to 45 m (not typical), showed, as expected, an increased effect of infrastructure with a 2.7% increase in energy exchange. However, this is still considered minimal in light of the uncertainty in other model parameters such as hydraulic and thermal conductivities. Overall, it can be concluded that properly accounting for surface infrastructure in BHE modeling is an onerous process, but had no significant impact on the outcome of this study; and this is expected to be the case for most investigations of a similar kind.

Modeling Study II: London Road, Reading, UK Objectives
A modeling case study was conducted to (1) assess interactions between systems in high-density deployment of GCHP systems in an urban setting typical for the south of the UK, where a large increase in use of these systems is predicted (Committee on Climate Change 2017), (2) investigate the impact of changing hydrogeological conditions and heating loads on the subsurface temperature field (thermal footprint) and system's performance.

Study Site
The study site, a residential area in the city of Reading (UK), about 60 km west of London (Figure 4), comprises two blocks of semi-detached houses built in the 1930s with frontage width varying between 5 m and 18 m and an approximate distance of 65 m between the blocks. The modeling exercise assumed that each of the 58 properties in the two housing blocks is fitted with a separate, verticalborehole BHE system used to provide seasonal heating only (i.e., unbalanced system). The houses are located about 100 m south of the River Thames (Figure 4a). The bedrock geology is Cretaceous Chalk, which in places, is overlain by Paleogene (clay with fine-grained sand) or superficial deposits (sand and gravels) and by river valley alluvium along the River Thames (Figure 4b). The Chalk is an important, dual-permeability aquifer of considerable thickness (∼400 m) that is generally productive due to the elevated secondary porosity/permeability provided by fractures. The heterogeneity of these natural fracture systems is a significant control on the distribution of groundwater flow rates and flow paths within the aquifer (Bloomfield 1996). The general groundwater flow direction at the study site is from the higher grounds in the SSE toward the river in the NNE. Water levels at the study site are at around 2-3 m below the ground surface.

Methodology
A 2D-model was set up of the study site in FEFLOW (Figure 4c), simulating a fully saturated aquifer with an initial thickness of 100 m and a groundwater gradient of 0.005 m/m (Darcy flux = 2.9 × 10 −7 m/s) representing regional groundwater flow. Hydraulic conductivities (K) within the Chalk are controlled by the distribution and properties of the inherent secondary fracture systems which vary considerably with depths as well as between different boreholes (Williams et al. 2006). In the absence of site-specific K data, a simpler 2D modeling approach was deemed sufficient for this study, integrating hydraulic variables over the vertical thickness of the aquifer, for example, by assigning an aquifer transmissivity of 500 m 2 /d (equivalent to K = 5.8 × 10 −5 m/s for an aquifer  thickness of 100 m), estimated from statistical analyses of pumping test data (Allen et al. 1997), rather than assigning speculative vertical K distributions. The approach is consistent with the model objectives to understand the risk of interactions between adjacent systems, which can be assessed from lateral temperature distributions provided by the 2D model. A temperature of 12 • C was assigned to the entire model area as the initial condition, consistent with measured groundwater temperature in the region (Shand et al. 2003), and also to inflowing groundwater via a heat transport boundary condition (BC) along the southern boundary (Figure 4c). Heat extraction at 58 nodes (corresponding to GCHPs in individual dwellings) was defined via a nodal sink/source heat transport BC (Figure 4c). Heat loads were calculated for each node (i.e., dwelling) by estimating the heat demand (HD) for a single dwelling. Estimations were based on published degree day data, available for the period August 20, 2007, to February 16, 2015(Environmental Change Institute 2015, and building parameter values in Table S3 to derive minimum, maximum, and median heat demand envelopes for each day of the year ( Figure 5). Monthly average air temperatures were assigned to the top boundary. Heat losses from buildings to the subsurface were ignored, as considerable losses are assumed to only occur from basements in direct contact with the underlying aquifer (Menberg et al. 2013a)-which is not the situation here as the properties do not have basements. Furthermore, high permeability settings within an extensive saturated zone (as assumed in this study) were found to promote (horizontal) heat dissipation away from the basements, thereby reducing the impact on vertical temperature disturbances beneath the buildings (Epting et al. 2017b;Bidarmaghz et al. 2019). Other model parameters are given in Table S3. The model was run for a period of 25 years for the three thermal load scenarios representing years with above average, average, and below average air temperatures, which correspond to total heat abstractions of 3.3 MWh, 6.2 MWh, and 10.1 MWh per dwelling per year. Model calibration and validation were not undertaken within this study as it relates to hypothetical installations for which there are no actual data. Instead, to assess model performance, parameter sensitivity was tested for thermal heating loads (L H ), transmissivity (T ), thermal conductivity (K th ), groundwater gradient (dl /dh), thermal dispersivity (α x ,α y ) and subsurface temperature (t ss ), and corresponding normalized sensitivity coefficients (SC) were calculated as the ratio of relative changes in model output over relative changes in parameter input.

Results
Periodic ground temperature variations in Figure 5 (dotted lines) within the BHE field are typical for seasonal BHE schemes with ground temperatures with decreasing temperature during the winter period (heat abstraction) and increasing during the summer period (recovery). In the absence of groundwater flow, annual mean ground temperatures decrease to 6.5 • C, 1 • C, and -5 • C for the minimum, median and maximum heat demand scenario, respectively. The system does not reach a steady-state condition for heat transfer during the 25 modeled heating seasons ( Figure 5) even for average or low heating loads. The spatial footprint of the thermally affected zone under these no-flow conditions is limited to a few (<10 m) meters around the installed BHE systems as heat transport is dominantly conductive, constrained by the subsurface thermal conductivity. In the presence of groundwater flow ( Figure 5, plots 1-3, Figure 6a), ground temperatures at individual BHEs during the heating season drop by up to 2.5 • C, 3.8 • C, and 7.0 • C for the minimum, median, and maximum heat demand scenario, respectively. Ground temperatures recover during the summer (no heating and higher surface temperatures) period, but remain below the background temperature of 12 • C by about 1 • C, 1.8 • C, and 3 • C, for minimum, median, and maximum heat demand scenarios, respectively. Mean annual ground temperatures (solid lines) stabilize after about 10 heating seasons at 11 • C (minimum HD), 10 • C (median HD), and 8.5 • C (maximum HD), suggesting that the system has reached a seasonal equilibrium or dynamic balance, even for high heat demands. The thermally affected zone around the BHE field is markedly dispersed in the direction of groundwater flow, extending to the northern model boundary, which represents the River Thames located approximately 100 m north/north east and down-gradient of the site. Temperature reductions along the river of up to 0.8 • C, 1.6 • C, and 2.6 • C (for the minimum, median, and maximum HD scenarios) highlight the potential impact that the modeled BHE schemes could have on nearby energy installations or heat-sensitive (eco)systems.
System efficiency is assessed via the seasonal performance factors (SPF), which, in this study, is calculated from the coefficient of performance (COP) of the heat pump averaged over the heating season (Singh et al. 2019). The SPF declines by about 0.1 for every 1 • C reduction in ground temperatures, hence higher reductions in efficiency in Figure 6b are associated with higher overall heating loads-as would be expected.
Efficiency reductions result in rising energy consumption (Figure 6c), and these were used to compare the impact of different operational and interference scenarios (Table 3). In the presence of groundwater flow, for the median HD scenario, the SPF stabilizes at an average value of 2.75 and an energy consumption of 396 MWh/year for the entire BHE field (or an average of 6.8 MWh/year per system). Corresponding CO 2 emissions are 138 t CO 2 /year per BHE field (and 2.4 t CO 2 /year per system), assuming a conversion factor of 0.35 kg CO 2 per kWh electricity (BEIS 2017a). Within each HD scenario, consumption of individual systems varies depending on their position within the borehole field, with differences of 3%, 5%, and 9% in daily consumption between the least and the most efficient systems at low, median, and high HD, respectively. In the absence of groundwater flow, as ground temperatures continue to decline, there is a dramatic decline in efficiency and an associated rise in energy consumption. After 25 years, the annual energy consumption of the BHE field is 533 MWh/year with corresponding CO 2 emissions of 187 t CO 2 /year (3.2 t CO 2 /year per system). For comparison, generating the equivalent amount of heating energy using gas-boilers would produce a total of 273 t CO 2 /year (or 4.7 t CO 2 /year per system), assuming a boiler efficiency of 80% (BEIS 2017b) and a conversion factor of 0.2 kg CO 2 per kWh for natural gas (BEIS 2017a). These CO 2 emissions are considerably (46-98%) higher than those produced by the BHE field even under suboptimal conditions, that is, in the absence of groundwater flow. The potential impact of interference with other nearby installations on system efficiency was assessed for the median HD by simulating a decrease in groundwater temperatures of 1-2 • C, which can be expected from the operation of a similar scheme 100-200 m upstream of the site. The results indicate reductions in system efficiency of 4-8% and an overall increase in CO 2 emissions of 5-10 t CO 2 /year (0.09-1.8 t CO 2 /year per system).
Model sensitivity ( Figure S2) is mostly associated with the hydraulic head gradient dl /dh, thermal loads L H to the ground (i.e., heat demand) and aquifer transmissivity T . The importance of correctly estimating heat loads to the subsurface are obvious. High model sensitivities at low transmissivities (T < 500 m 2 /d) and low hydraulic gradients (dl /dh < 0.005 m/m) can also be expected due to the decreasing effect of advection at lower groundwater flow velocities resulting in reduced dissipation of heat and hence larger temperature increases in response to heating loads. Thermal Peclet numbers Pe t were calculated after Bear (1972) for the Darcy fluxes listed in Table 3 to assess the influence of advection on heat transport for the different transmissivity/hydraulic head gradient settings. Peclet numbers Pe t were >1 in all cases, except where dl /dh = 0 m/m (no groundwater flow), suggesting that advection of heat by flowing groundwater is a significant process contributing to heat transfer in the ground (Chiasson et al. 2000).
Other parameters that impact ground temperatures at the BHE are longitudinal (α x ) and transverse (α y ) dispersivity, as demonstrated in more detail by other modeling studies (Molina-Giraldo et al. 2011a;Piga et al. 2017;Pophillat et al. 2020). Model outputs are relatively insensitive to thermal conductivity values, which controls conductive heat transport (Liuzzo-Scorpo et al., 2015), again confirming the dominance of advective (rather than conductive) heat transport within the modeled systems.

Discussion
An initial heat balance estimate for a single GCHP over 1 year for the median heat load scenarios suggests that, in the absence of groundwater flow, heat abstraction at the proposed rate is not sustainable in the long term, not even for a single system. The modeling confirms this and highlights that sustainability and system efficiency over the anticipated operational lifespan is largely controlled by the presence of groundwater flow. Ground temperature Although efficiency reductions lead to an increase in CO 2 emissions, and are therefore undesirable and should be minimized, it is interesting to observe that overall CO 2 emissions of the simulated systems after 25 years of operation remain below those that would have been produced if gas boilers had been used to provide the heating. This applies even for scenarios where groundwater flows were absent, and the systems are considered unsustainable (on the basis that ground temperature continues to drop due to an imbalance in thermal extraction and recharge). Finally, model sensitivity has highlighted key controls on model performance, confirming the importance of processes linked to groundwater flow, that is, thermal advection and dispersion.

Modeling Study III: Carignan-Salières Elementary School, South of Montréal, Canada Objectives
The objective of this third study was to (1) predict the long-term performance of an entire BHE field installed for a school building and affected by variable groundwater flow in order to (2) anticipate potential operational interference with dewatering of a nearby quarry and (3) evaluate the effect of groundwater flow on the thermal plume around the BHEs. To this purpose, a numerical model was calibrated with a large-scale heat injection test using the whole borefield and then simulations were run under different scenarios for a period of 20 years. The distribution of the thermal plume around the BHE field is newly addressed in this article as a complement to results given in a previous study (Jaziri et al. 2020).

Description of the Case Study
The Carignan-Salières elementary school is located on the south shore of the St. Lawrence River near Montréal, Canada, about 1 km away from an active quarry which is irregularly dewatered to facilitate excavations (Figure 7). The building lies on the Nicolet Formation, a sedimentary rock unit belonging to the Loraine Group, which is part of the St. Lawrence Lowlands sedimentary basin (Brisebois and Brun 1994). The formation consists of sequences of silty gray shale, with interbedded sandstone, siltstone, and limestone (Globensky 1987). Gabbro dykes, observed in the school area, are oriented EW and cut the stratigraphic sequence (Foster and Symons 1979;Foland et al. 1986;Feninger and Goodacre 1995). The direction of the groundwater flow is locally oriented toward the active quarry (SW) due to dewatering (Figure 7a and 7b). The school building constructed in 2013 is heated and cooled with a GCHP system experiencing varying groundwater flow conditions.
The GCHP system consists of 31 BHEs connected to 50 heat pumps, with net heating capacity from 3.62 to 44.2 kW, depending on the size of rooms to heat and cool. The BHEs are 152 m deep, spaced by 6 m and made with high-density polyethylene single U-pipes (outer diameter 32 mm, thermal conductivity 0.39 W/m/K) with omega-shaped spacers. During the installation, the boreholes could not be sealed with thermally enhanced grout made of sand and bentonite, which is commonly used in Canada to fill boreholes, because groundwater along the intersecting fractures flushed the fine particles from the grout mixture. As a consequence, boreholes were filled with olivine sand having thermal conductivity of 1.75 W/m/K (Côté et al. 2012). The heat carrier fluid is a mixture of water and propylene glycol at 25% vol. and circulates in the BHE loops at a total average flow rate of 1017 m 3 /d (0.38 L/s in each BHE). Heating and cooling annual energy consumption of the school building was determined with an eQuest simulation using the DOE2.2 algorithm (Hirsch 2004) and is 290 MWh/year, with peak heating and cooling loads of 494 kW and 253 kW occurring in January and July, respectively. This induces significant unbalanced ground conditions that can affect the long-term thermal response of the system. Initial site characterization included two TRTs, carried out before and after the BHE field installation, which revealed a bulk subsurface thermal conductivity of 2.58 W/m/K and 2.27 W/m/K, respectively (Jaziri et al. 2020 and references therein). Thermal conductivity values are in agreement with literature data (Bédard et al. 2017;Raymond et al. 2017;Raymond et al. 2019); and the difference between the two tests is assumed to be due to changes in the groundwater flow regime near the school. Rock samples of the gabbro dykes, shales, and calcarenite, collected from the quarries, showed an average thermal conductivity of 1.85, 2.64, and 3.5 W/m/K, when respectively measured in the laboratory (Jaziri et al. 2016). Hydraulic conductivity and recharge were assessed by reproducing the hydraulic head measured in the abandoned quarry (h) considered as an observation point and using an analytical solution for steady-state flow in an unconfined aquifer (Fetter 2001). The conditions that best match the observation point (h = 20.7 m a.s.l.) are a hydraulic conductivity of 1.26 × 10 −5 m/s and a net recharge of 100 mm/year (Jaziri et al. 2016), both in agreement with the available regional groundwater flow assessment (Carrier et al. 2013).

Methodology
A large-scale heat injection test enclosing the 31 BHEs was carried out during hot summer days in July 2015. The test was carried out by using the cooling system at its full capacity for 16.9 days (305 kW total; 9.8 kW per BHE). This was achieved by opening the school windows during summer vacation to allow the outdoor heat to enter the building while it was not used. The cooling system was then stopped and the heat carrier fluid was kept circulating in the loop to monitor the thermal re-equilibration during an additional 13.3 days similar to a TRT with monitoring of the recovery period. The flow rate and the inlet/outlet temperature of the GHE field were monitored during the whole test by means of flowmeters (accuracy ±1.5%) and temperature sensors (accuracy 0.1 • C) at a 30-s sampling interval (Figure 8a).
Numerical simulations were run to calibrate the FEFLOW model with data from the large-scale heat injection test. The size of the 3D model was 500 × 500 × 300 m and spatially discretized in six layers of 50 m each, resulting in 195,720 triangular prismatic elements and 114,450 nodes ( Figure 7c). Before the transient simulations, the initial temperature was achieved with a stationary simulation based on a local geothermal gradient of 23.1 • C/km and heat flow of 35 mW/m (Bédard et al. 2017;Raymond et al. 2017;Nasr et al. 2018). Steady-state groundwater flow in a simplified unconfined aquifer system with surface recharge was considered. Calibration parameters (hydraulic conductivity, porosity, thermal conductivity of BHE's grout, and host rock) were adjusted manually, one at a time until the model reproduced the observed BHE outlet temperature with a maximum error of 2% (Figure 8b). For simplicity, the geology surrounding the BHEs was assumed to be uniform with dominating silty gray shales and the same material properties were assigned to all six layers: a hydraulic conductivity of 10 −4 m/s in the x and y direction and in the z direction 10 −6 m/s; a matrix and fluid thermal conductivity of 2.4 and 0.6 W/m/K, respectively, and a porosity of 0.03. After the calibration, the long-term performance of the BHE field was evaluated by means of 20-year simulations in order to predict the long-term impact of groundwater flow on the GCHP operation. Two different scenarios were evaluated, simulating conditions of low groundwater flow (Scenario 1 with a hydraulic gradient of 0.0006 m/m), and high groundwater flow associated with dewatering activities in the quarry (Scenario 2 with a hydraulic gradient of 0.008 m/m). Constant hydraulic heads, with different values according to chosen simulation scenarios, were then imposed on the eastern and western boundaries of the model. The bottom surface was set as an impermeable (no flow) boundary, and an annual net recharge of 100 mm/year was imposed at the top surface (see Figure 7c). Lateral heat transfer boundaries were set adiabatic. Average and constant coefficients of performance (COP) of 4.7 and 4.1 in heating and cooling mode, respectively, were assumed for all the school heat pumps when calculating the ground loads to be used as inputs to the BHE model in FEFLOW.

Results
The 20-year simulation results conducted using the calibrated model show significant differences in BHE fluid temperature of the two scenarios (Figure 8c and 8d). Despite the differences, both scenarios show an adequate thermal response of the subsurface although ground loads are unbalanced. Lower groundwater flow (Scenario 1), which represents a case where pumping in the active quarry is stopped or reduced, has a clear negative impact on the whole system temperature. In heating mode, the BHE inlet temperature drops to −5 • C and 3 • C in Scenarios 1 and 2, while the outlet temperature reaches 3 • C and 7 • C, respectively. As expected, Scenario 2 provides a better operating temperature, and therefore better GCHP performance. In Scenario 1, the minimum BHE outlet temperature is sufficiently higher than the minimal operating temperature of the heat pump system recommended by the manufacturer (−9.62 • C). However, the minimum inlet temperature is within 5 • C of the freezing point of the heat carrier fluid, here −10 • C. After 1 year, the thermally affected zone caused by cooling the building is little affected in Scenario 1, while it is markedly dispersed and follows the groundwater flow direction in Scenario 2 ( Figure 9). Groundwater flow appears to have an important impact on the dispersion of the thermal plume around the BHE field that is at least 25 m wider in Scenario 2. Dispersion of the hot and cold front around the BHEs due to heat transfer enhanced by advection is believed to be the mechanism responsible for better operation temperatures obtained with Scenario 2. Therefore, under low groundwater flow conditions, in the event that the quarry stops dewatering activity, care should be taken to follow the system minimum operating temperature during winter periods to avoid potential freezing problems at the GHE inlet.

Discussion
This case study illustrates the 20-year performance of a GCHP system with temperature simulations affected by dewatering activities in a nearby active quarry (less than 1 km from the BHE). The GCHP system of the Carignan-Salières School provides a unique field case with BHEs interfering with the groundwater drawdown around the quarry and where the local thermal and hydraulic conditions of the GCHP system have uncommonly been assessed at a large scale. The subsurface heat exchange capacity of the GCHP system is clearly enhanced by groundwater advection when the specific Darcy velocity increases from 6 × 10 −8 m/s (no dewatering) to 8 × 10 −7 m/s (high dewatering). This study provides further evidence that even the lowest groundwater flow conditions expected at the site can be beneficial to avoid a progressive cooling of the underground over the expected life of the system due to the unbalanced heating and cooling loads. In a previous study, (Jaziri et al. 2016) simulated the operation of the GCHP system with a heat conduction approach considering an equivalent subsurface thermal conductivity up to 3 W/m/K and assumed affected by groundwater flow. The BHE operating temperatures at the beginning of the simulations were similar to those obtained with FEFLOW and presented in this paper for Scenario 1 (low groundwater flow), but decreased by 4 to 6 • C over the 20 years of system operation. BHE simulations considering advection did not show a significant decrease of the minimum outlet temperature over the life of the system, even with low groundwater flow (Figure 8), which is believed to be due to dispersion of the warm and cold front around the BHEs (Figure 9). The fact that low groundwater flow can help dissipate heat in the ground to help cope with unbalanced ground loads has important implications for GCHP system design, especially for systems subject to interference.

Summary and Conclusions
The three case studies highlight that GCHP systems can be impacted by perturbations arising from subsurface activities that change the thermal and hydraulic regimes in the surrounding areas.
Changes in the thermal regime arising from additional subsurface heat loads near existing schemes were found to have varied impacts on system efficiency with reductions ranging from <1% to 8%. A clear difference was observed between impacts of additional loads on balanced (Case Study 1) compared to unbalanced (Case Studies 2 and 3) systems, with overall efficiency reduction being much smaller for balanced schemes (<1%) compared to unbalanced schemes (3-8%) despite similar (or higher) subsurface temperature changes.
For unbalanced systems, thermal interference is unavoidable where individual systems are installed in close proximity, on the order of tens of meters. However, interaction within the field between the individual BHE had less impact on the efficiency of individual systems than interaction with large heat loads from neighboring systems, for example when an additional borehole field with similar heat loads is installed upstream of the existing systems. Such thermal interferences between GSHP systems have long been predicted (Fascì et al., 2019;Ferguson and Woodbury, 2007), but evidence of system interference in published case studies remains rare. By analyzing temperatures of the pumped groundwater for an open-loop system in London, Herbert et al. (2013) identified thermal interference that was attributed to operations of a nearby GCHP system. However, the source of the interference could not be confirmed as there are no requirements for licensing or monitoring of GCHPs in the UK, not even for recording their location.
Thermal losses from near-surface infrastructure were found to result in significant temperature changes (up to 10 • C) in the zone of 0-20 m below ground surface, with observable impacts ( T > 1 • C) up to a depth of 75 m. While such temperature increases can be expected to benefit the performance of a heating-only system, it had only a minor impact (+2.2% increase) on the efficiency of the balanced system.
Changes in the hydrogeological regime were confirmed as the main control on GCHP performance in all three studies. Even low groundwater flow rates were found to improve the performance of the system overall, and vice versa: small reductions in groundwater flow reduced system efficiencies by a considerable margin. This is especially true in the case of unbalance ground loads where groundwater flow can decrease the temperature effect of the unbalance loads as shown in Case Study 2 and 3.
The effects of groundwater flow on the design and performance of GCHP systems has been demonstrated by various analytical and numerical modeling studies (Chiasson et al. 2000;Diao et al. 2004;Molina-Giraldo et al. 2011b). However, it is difficult to have welldocumented field cases with calibrated models to show the groundwater impact of BHE fields. The current study offers such field cases confirming that advection is important to consider in the system design, especially where Darcy fluxes of greater than 1 × 10 −7 m/s are expected (Dehkordi and Schincariol 2014;Ferguson 2015).
The modeling highlights the importance of considering subsurface activities that can change subsurface groundwater flows in the design and operation of BHEs as they have potential to impact the efficiency of nearby GCHP systems. This is particularly evident in Case Study 3 where quarry dewatering activities (at ∼1 km distance from the study site) showed a clear interaction with GCHP's operating temperatures and system efficiency, through influencing groundwater gradients and hence flow rates. A ∼13% difference in thermal conductivity was obtained from two thermal response tests (TRT) at the site which is attributed to changes in groundwater flow related to dewatering activities. In situ measurements of thermal conductivity using TRTs (Raymond et al. 2011) are now widely recommended as part of the design process (e.g., GSHPA 2017), but this study demonstrates that these need to be considered within the context of groundwater flow conditions at the time of the test to ensure accurate sizing of the BHE installation. The school system of Case Study 3 was designed for subsurface conditions with active dewatering and can obviously experience decreasing performances if the quarry is shut down and the dewatering is stopped. Understanding the thermal regime of a given site under varying flow conditions should be a priority for GCHP design, specifically where groundwater flows are expected to vary. Modeling studies can help to evaluate the impact of changing thermal and groundwater flow conditions on the BHE operating temperature and need to be systematically considered for GCHP design. While groundwater resources are regulated in most countries, thermal abstractions/ discharges to the subsurface are largely unregulated, although approaches for the regulations of GSHPs vary greatly between different countries (Tsagarakis et al. 2020;Haehnlein et al. 2010;Dehkordi and Schincariol 2014;García-Gil and Moreno 2020). The modeling studies presented here support the argument in favor of regulation to, as a minimum, register GCHP systems with records of locations and approximate heat pump capacity-even though these systems do not abstract/inject groundwater. Additional regulation can be put in place to ensure that subsurface thermal equilibrium is maintained around the properties with GCHP systems using a threshold temperature yet to be defined. This is currently not the case in the UK or in Canada, even for large systems with high heating/cooling loads. As others have pointed out (Herbert et al. 2013), this poses an increasing risk for interference problems as numbers of installations increase in densely populated areas.
More comprehensive data on the actual system location as well as their cumulative heating and cooling loads is also required if the underground thermal resource is to be managed sustainability. In some countries, this may require the designation of heat as a natural resource in order to legislate its use (Abesser et al. 2018). The management of the subsurface thermal resource requires some assessment of where systems should be deployed and how. In the city of Zürich, for example, active regeneration of the underground thermal resource is mandatory for GCHP systems in areas of high-density GCHP deployment (Knüsel 2015;Stadt Zürich 2014). Various tools and approaches have been developed for assessing and managing subsurface thermal resources (García- Gil et al. 2020aGil et al. , 2020b(Epting et al. 2013), but operational subsurface temperature data are rarely available for calibration and validation of such models-although exceptions exist (e.g., Zaragoza, Basel- Epting et al. 2017a). The temperature monitoring during the heat injection test for Case Study 3 was done in the scope of a research project to anticipate GCHP operational problems but is certainly not a requirement in Canada.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Supporting Information is generally not peer reviewed.
Table S1. BHE model specifications Table S2. Average monthly inlet temperatures Table S3. Parameter values used for heat demand calculation, heat pump design, and heat transport modeling Figure S1. Temperature versus depth profile of the UWO monitoring wells compared to simulated temperatures. The prefix "F" denotes field data, and the prefix "M" denotes modeled data. The minor divergence (maximum 0.3 • C) centered around a 60 m depth may be due to localized fractures in the limestone, or a shifting of thermistors during installation. Figure S2. Cumulative time curves of the normalized sensitivity coefficient for selected model input parameters