Long-Term Temperature Evaluation of a Ground-Coupled Heat Pump System Subject to Groundwater Flow

: The performance of ground-coupled heat pump systems (GCHPs) operating under signiﬁcant groundwater ﬂow can be di ﬃ cult to predict due to advective heat transfer in the subsurface. This is the case of the Carignan-Sali è res elementary school located on the south shore of the St. Lawrence River near Montr é al, Canada. The building is heated and cooled with a GCHP system including 31 boreholes subject to varying groundwater ﬂow conditions due to the proximity of an active quarry being irregularly dewatered. A study with the objective of predicting the borehole temperatures in order to anticipate potential operational problems was conducted, which provided an opportunity to evaluate the impact of groundwater ﬂow. For this purpose, a numerical model was calibrated using a full-scale heat injection test and then run under di ﬀ erent scenarios for a period of twenty years. The heat exchange capacity of the GCHP system is clearly enhanced by advection when the Darcy ﬂux changes from 6 × 10 − 8 m s − 1 (no dewatering) to 8 × 10 − 7 m s − 1 (high dewatering). This study further suggests that even the lowest groundwater ﬂow condition can be beneﬁcial to avoid a progressive cooling of the subsurface due to the unbalanced building loads, which can have important impacts for design of new systems.


Introduction
Groundwater flow can have a significant impact on the long-term performance of vertical ground heat exchangers (GHEs), especially when the Darcy flux is greater than 1 × 10 −7 m s −1 [1,2]. Flow in the subsurface can actually improve long-term GHE performance by dissipating heat injected into or extracted from the ground [3]. However, the design of ground-coupled heat pump (GCHP) systems is commonly based on the assumption of conductive heat transfer in the subsurface [4]. Although numerical tools are available to optimize the operation of a GCHP system under the influence of groundwater flow [5,6], an estimate of the specific groundwater flux affecting a GHE field can be difficult to define accurately. Alternatively, different authors have performed thermal response tests (TRTs) on a single GHE subject to groundwater flow to evaluate the equivalent subsurface thermal conductivity impacted by advection [7,8]. Practitioners commonly using such a heat conduction approach to simulate the long-term operating temperature of GHEs, based on an equivalent subsurface thermal conductivity assumption, have tried to cope with systems influenced by groundwater flow. This simplified approach can be useful, but neglects the fact that flow conditions can change over time.
The Carignan-Salières elementary school is an example of such a building, which is heated and cooled with a GCHP system operating under varying groundwater flow conditions because of its location within a kilometer of two quarries, one of which is being actively dewatered. In addition, because of the significant groundwater flow rates, upon installation the GHEs were backfilled with sand instead of using a geothermal grout. The grouting procedure was attempted but was not successful because the high flow rate most likely dispersed the fine grout particles into the fractures of the geological formations. This change of borehole filling material from the initial design plans, as well as the variable groundwater flow conditions, has likely affected the long-term operating temperature of the GHEs and their performance. The objective of this study is, therefore, to better understand the heat transfer mechanisms affecting the long-term operating temperature of the GHEs, including variable groundwater flow conditions, to anticipate potential operational problems.
Previous studies have been conducted to evaluate the impact of groundwater flow on the operating temperature of GHEs. For example, different authors simulated the operation of thermal energy storage systems with GHEs influenced by groundwater flow [9][10][11], while similar modelling was completed for GCHP systems by Dehkordi et al. [1]. However, most available studies are based on theoretical modelling exercises and have not been validated with operational data from systems that are influenced by significant groundwater flow conditions. This study provides a unique field case, where the GCHP system is located near an active quarry, at which groundwater flow conditions and ground thermal properties could be assessed. The GHE operating temperature could thus be accurately predicted based on a large-scale heat injection test carried out for the whole bore field, rather than a single GHE used for a TRT. Long-term numerical simulations of the system temperature under different field-based scenarios were carried out for a period of twenty years. A comparison between the simulated scenarios allowed a quantitative evaluation of the ground physical parameters affecting the GHE temperature, including the site hydraulic gradient. The knowledge gained can be used to improve the design and construction of new GCHP systems under the influence of groundwater flow.

Site Description
The Carignan-Salières elementary school is located on the south shore of the St. Lawrence River near Montreal ( Figure 1) and was built above the Nicolet Formation. This sedimentary rock formation belongs to the Loraine Group, known to have a thermal conductivity on the order of 3 W m −1 K −1 [12][13][14], and is part of the St. Lawrence Lowlands sedimentary basin [15]. From a structural perspective, the St. Lawrence Lowlands are characterized by a series of normal faults extending from the southwest to the northeast, and which dip toward the southeast [16,17].

Geological and Hydrogeological Setting
The Nicolet Formation, hosting the school GHE field, consists of sequences of silty gray shale, with interbedded sandstone, siltstone and limestone [18]. Gabbro dykes, which are observed in the school area, are oriented EW and cut the stratigraphic sequence [19][20][21]. Most of the dykes are sub-vertical and have a varying thickness of 0.5 to more than 20 m [22]. One of the two quarries near the school is currently active and water is pumped irregularly (Figure 2). Regional aquifers in this area of the south shore of the St. Lawrence River are hosted in fractured rocks. The direction of groundwater flow at the site is oriented toward the active quarry where water is being pumped. The average hydraulic conductivity of the host rock reported in the area is 10 −5 m s −1 , varying from 3.2 × 10 −7 to 1.6 × 10 −4 m s −1 , and the annual net recharge of the aquifer is approximately 100 mm y −1 [23]. The town of Carignan has a humid continental climate with an annual average air temperature of 5.9 • C and an amplitude of 30 • C [24]. The heating period is mostly from October to June while the cooling period is from July to September.

Geological and Hydrogeological Setting
The Nicolet Formation, hosting the school GHE field, consists of sequences of silty gray shale, with interbedded sandstone, siltstone and limestone [18]. Gabbro dykes, which are observed in the school area, are oriented EW and cut the stratigraphic sequence [19][20][21]. Most of the dykes are subvertical and have a varying thickness of 0.5 to more than 20 m [22]. One of the two quarries near the school is currently active and water is pumped irregularly (Figure 2). Regional aquifers in this area of the south shore of the St. Lawrence River are hosted in fractured rocks. The direction of groundwater flow at the site is oriented toward the active quarry where water is being pumped. The average hydraulic conductivity of the host rock reported in the area is 10 −5 m s −1 , varying from 3.2 × 10 −7 to 1.6 × 10 −4 m s −1 , and the annual net recharge of the aquifer is approximately 100 mm y -1 [23]. The town of Carignan has a humid continental climate with an annual average air temperature of 5.9 °C and an amplitude of 30 °C [24]. The heating period is mostly from October to June while the cooling period is from July to September.

Ground Heat Exchanger Characteristics
The GCHP system consists of thirty-one GHEs connected to fifty heat pumps distributed in the school building. Each heat pump has a net heating capacity of 3.62 to 44.2 kW. The system has sufficient capacity to supply the entire heating and cooling loads of the building. The GHEs are spaced by 6 m and each borehole is ~152 m deep. A high-density polyethylene single U-pipe with omega-shaped spacers constitutes the GHEs. During the installation, the boreholes could not be sealed with grout because groundwater flow along the intersecting fractures flushed the fine particles from the grout mixture. The boreholes were consequently filled with olivine sand instead of a common thermally-enhanced grout made of bentonite and sand. The olivine sand has a thermal conductivity of 1.75 W m −1 K −1 , a value that was measured in the laboratory during a previous study [25]. The heat carrier fluid is a mixture of water and propylene glycol at 25 vol. % and circulates in the closed loops with an average flow rate of 1017 m 3 d −1 within the entire GHE field. The GHEs are connected in parallel so the flow rate in each borehole is similar and around 3.8 × 10 −4 m s −1 . Two TRTs were performed in two different boreholes, by injecting heat in order to evaluate the in-situ thermal conductivity of the subsurface. The first TRT was carried out before the GHE installation and revealed a bulk subsurface thermal conductivity of 2.58 W m −1 K −1 [26]. The second TRT was conducted after the GHE installation and indicated a lower value of 2.27 W m −1 K −1 [27]. The difference between the two tests is believed to be due to water pumping in the neighbouring and active quarry, inducing changes in the groundwater flow regime near the school.

Methodology
Taking advantage of the two quarries, field activities were undertaken to estimate the thermal conductivity of the subsurface and the groundwater flow conditions near the school site using representative values characteristics of a heterogeneous system. Twelve rock samples (6 shales, 1 calcarenite and 5 gabbros) were collected in the quarries to be analyzed in the laboratory. Water levels

Ground Heat Exchanger Characteristics
The GCHP system consists of thirty-one GHEs connected to fifty heat pumps distributed in the school building. Each heat pump has a net heating capacity of 3.62 to 44.2 kW. The system has sufficient capacity to supply the entire heating and cooling loads of the building. The GHEs are spaced by 6 m and each borehole is~152 m deep. A high-density polyethylene single U-pipe with omega-shaped spacers constitutes the GHEs. During the installation, the boreholes could not be sealed with grout because groundwater flow along the intersecting fractures flushed the fine particles from the grout mixture. The boreholes were consequently filled with olivine sand instead of a common thermally-enhanced grout made of bentonite and sand. The olivine sand has a thermal conductivity of 1.75 W m −1 K −1 , a value that was measured in the laboratory during a previous study [25]. The heat carrier fluid is a mixture of water and propylene glycol at 25 vol. % and circulates in the closed loops with an average flow rate of 1017 m 3 d −1 within the entire GHE field. The GHEs are connected in parallel so the flow rate in each borehole is similar and around 3.8 × 10 −4 m s −1 . Two TRTs were performed in two different boreholes, by injecting heat in order to evaluate the in-situ thermal conductivity of the subsurface. The first TRT was carried out before the GHE installation and revealed a bulk subsurface thermal conductivity of 2.58 W m −1 K −1 [26]. The second TRT was conducted after the GHE installation and indicated a lower value of 2.27 W m −1 K −1 [27]. The difference between the two tests is believed to be due to water pumping in the neighbouring and active quarry, inducing changes in the groundwater flow regime near the school.

Methodology
Taking advantage of the two quarries, field activities were undertaken to estimate the thermal conductivity of the subsurface and the groundwater flow conditions near the school site using representative values characteristics of a heterogeneous system. Twelve rock samples (6 shales, 1 calcarenite and 5 gabbros) were collected in the quarries to be analyzed in the laboratory. Water levels in the two quarry lakes were measured with a global positioning system to calculate the local hydraulic gradient. A heat injection test was further carried out for the entire GHE field with the system operating in cooling mode during the summer. Field-based observations allowed developing a representative numerical model of the GHE field, taking into account hourly building heating and cooling loads transferred to the ground.

Laboratory Measurement of Thermal Conductivity
The thermal conductivity assessments were made with two different methods: (1) a needle probe for hard rocks [12,[28][29][30] and (2) the modified transient plane source (MTPS) method for friable rocks [30][31][32]. The KD2 Pro unit [33], part of a standardized method under the ASTM D5334 norm, was used for the needle probe analysis with gabbro samples, while the C-therm heating plate following the ASTM D7984 norm was used for the MTPS analysis with shales and calcarenites [34]. Before making the measurements, the samples were immersed for 24 h in distilled water to saturate the samples and obtain a value representative of subsurface conditions.
The samples analyzed with the needle probe were cut and a hole was drilled in each sample to insert a 6 cm long and 4 mm thick needle injecting heat at a rate of 6 W m −1 ; the probe was covered with thermal grease to reduce contact resistance. After completing a measurement with a reference polyethylene sample, each rock sample was analyzed in a controlled-temperature room with at least ten measurements. An interval of 1 h was adopted for thermal equilibrium to be restored and a correction factor was taken into account to consider the reference sample measurements. With the MTPS method, a polished surface of the soft rock samples was placed on the heating plate to evaluate its thermal conductivity according to the transient increase of temperature. The heating plate has a disk shape and the electric signal used for the heat source serves as a proxy for temperature. The power level applied to the heat source is 90 V. Thermal conductivity of each sample was measured five times and an average was then calculated to obtain the final value.

In-Situ Heat Injection Test
The heat injection test to evaluate the thermal response of the entire GHE field was carried out during a hot summer period in July 2015. The test was completed by using the cooling system at its full capacity over 16.9 days, while the school windows had been opened to allow the outdoor heat to enter the building. The cooling system was then stopped, and the heat carrier fluid was kept circulating in the loop to monitor the thermal recovery during an additional 13.3 days. The flow rate of the heat carrier fluid and the temperature at the inlet and outlet of the entire GHE field were measured during the test by using the temperature sensors and flowmeters installed in the mechanical room of the GCHP system near the heat pump inlets. The temperature sensors and flowmeter have an accuracy of ±0.1 • C and ±1.5%, respectively. The instruments allowed taking measurements every 30 s. The heat injection rate was calculated from the temperature and flow rate measurements.

Building and Ground Load Evaluation
The school building loads were simulated using the program eQuest, a graphical interface for the DOE-2 program [35]. The simulations were used to establish a thermal energy budget of the building, which depends essentially on the building dimensions, the construction and insulation materials, the size and number of windows and doors, the operation schedule, as well as the internal and external temperature. Heat losses and gains of the indoor spaces were evaluated hourly to determine heating and cooling loads every hour over a full year. Simulated building loads, used as inputs in the numerical GHE model, were converted to ground loads according to: and where P ground and P building [W] are the loads for the ground and the building, respectively, and COP heating and COPcooling [-] are the heat pump coefficients of performance in heating and cooling modes, respectively. An average and constant COP for all heat pumps of the building throughout the simulation time was assumed for simplification.

GCHP System Simulation
Numerical simulations to calibrate the model using the short-term heat injection test and to subsequently evaluate the long-term operating temperature of the GHEs over 20 years were carried out with the finite element program FEFLOW [36]. This program was used because it allows simulating transient heat transfer and groundwater flow in 3D porous media hosting GHEs embedded as 1D elements. A general iterative finite element strategy is used to solve the overall flow and heat transfer problem coupling the ground with the GHEs [37]. The numerical approach of Eskilson et al. [38], which represents GHEs with equivalent resistances, was selected instead of the more general numerical approach of Al-Khoury et al. [39,40]; both being available in FEFLOW. This choice was made because the Eskilson approach requires less computational effort and has been demonstrated to be accurate for long-term predictions [37].

Governing Equations
The global fluid flow and heat transfer problem is solved in the form of fluid, mass and thermal energy balances for the subsurface s and the groundwater gw. The flow equation is described by [36]: is a source/sink term for flow and EOB refers to the Extended Oberbeck-Boussineq approximation. In Equation (3), q [m s −1 ] is the Darcy flux in the porous medium and is expressed with Darcy's law: where K [m s −1 ] is the hydraulic conductivity tensor. The heat transport equation with conductive and advective terms is described by [41]: is the volumetric heat capacity and H [W m −3 ] represents in this case heat sources and sinks that can be, for example, tied to the GHEs. The GHE equations used in the present model apply for single U-pipe heat exchangers, that, according to the thermal resistance and capacity model of Bauer et al. [42], consist of four components, namely an inlet-pipe denoted with subscript il, an outlet-pipe with subscript ol and two filling material zones (grout) with exponent or subscript g, and can be written as [36]: where ρ [kg m −3 ] and c [J kg −1 K −1 ] are the density and the specific heat of the heat carrier fluid r and grout g, respectively, is the vector of the heat carrier fluid velocity, 2 r [W m −2 K −1 ] is the thermal hydrodynamic dispersion tensor for the heat carrier fluid, λ g [W m −1 K −1 ] is the thermal conductivity of the grout and ε g [−] is the volume fraction of the grout or filling material. The relations are used to express the heat exchange between the borehole with 1 U-pipe and the subsurface. The effect of the 1 U-pipe components is lumped into effective heat transfer coefficients, which represent the sum of thermal resistances between the different components of the GHE elements [37].

Model Geometry and Properties
The surface area of the numerical GHE model is 500 m × 500 m and extends from ground surface to a depth of 300 m, divided into 6 layers of 50 m each. The same lithological characteristics of the subsurface, such as hydraulic conductivity, thermal conductivity of both fluid and solids, as well as porosity are assigned to the six layers that have a different initial temperature based on the geothermal gradient. The model mesh consists of 3D triangular prismatic elements for a total of 195,720 elements and 114,450 nodes. A 2D horizontal-plane mesh was built and then extended in 3D to cover the entire domain ( Figure 3). Each GHE was surrounded by 6 nodes. Thirty-one GHEs were added to the domain, following the same layout as at the Carignan-Salières elementary school. In the plan view, the triangular mesh dimensions are up to approximately 4 m at the boundaries of the model (Figure 3), decreasing to 0.06 m around each GHE.
is the thermal conductivity of the grout and εg [−] is the volume fraction of the grout or filling material. The relations are used to express the heat exchange between the borehole with 1 U-pipe and the subsurface. The effect of the 1 U-pipe components is lumped into effective heat transfer coefficients, which represent the sum of thermal resistances between the different components of the GHE elements [37].

Model Geometry and Properties
The surface area of the numerical GHE model is 500 m × 500 m and extends from ground surface to a depth of 300 m, divided into 6 layers of 50 m each. The same lithological characteristics of the subsurface, such as hydraulic conductivity, thermal conductivity of both fluid and solids, as well as porosity are assigned to the six layers that have a different initial temperature based on the geothermal gradient. The model mesh consists of 3D triangular prismatic elements for a total of 195,720 elements and 114,450 nodes. A 2D horizontal-plane mesh was built and then extended in 3D to cover the entire domain ( Figure 3). Each GHE was surrounded by 6 nodes. Thirty-one GHEs were added to the domain, following the same layout as at the Carignan-Salières elementary school. In the plan view, the triangular mesh dimensions are up to approximately 4 m at the boundaries of the model (Figure 3), decreasing to 0.06 m around each GHE.

Initial and Boundary Conditions
Type-1 constant hydraulic heads with different values according to chosen simulation scenarios were imposed on the eastern and western boundaries of the model. This head gradient was applied to represent local groundwater flow conditions due to pumping in the active quarry, which is 1 km west (down-gradient) of the school location. The bottom surface was set impermeable and an annual net recharge of 100 mm y −1 was imposed at the top surface [23]. The initial ground temperature was assigned to each layer according to the geothermal gradient of the area (Figure 4), which is equal to 23.1 °C km −1 [12]. A fixed temperature boundary condition was used at the top of the model, while a

Initial and Boundary Conditions
Type-1 constant hydraulic heads with different values according to chosen simulation scenarios were imposed on the eastern and western boundaries of the model. This head gradient was applied to represent local groundwater flow conditions due to pumping in the active quarry, which is 1 km west (down-gradient) of the school location. The bottom surface was set impermeable and an annual net recharge of 100 mm y −1 was imposed at the top surface [23]. The initial ground temperature was assigned to each layer according to the geothermal gradient of the area (Figure 4), which is equal to 23.1 • C km −1 [12]. A fixed temperature boundary condition was used at the top of the model, while a constant heat flux was imposed at the bottom of the model in order to represent the Earth's heat flux in the St. Lawrence lowlands, which is assumed equal to 35 mW m −1 [43]. The lateral side boundaries were considered as Type-2 for heat transfer, with a zero temperature gradient assuming no conductive heat transfer. The GHEs in FEFLOW are Type-4 boundary conditions for which two approaches were used for the simulations. First, the inlet temperature was specified with an average flow rate of 1076 m 3 d −1 for the short-term calibration simulations to reproduce the outlet temperature during the heat injection test. Second, a total flow rate of 1017 m 3 d −1 , anticipated for the long-term operation of the system, and ground loads defined from the building simulation and the assumed heat pump COPs (Equations (1) and (2)), were assigned to predict the operating temperature of the GHE field for twenty years. For all numerical simulations, including model calibration, hourly time steps were used to ensure accurate results. during the heat injection test. Second, a total flow rate of 1017 m 3 d −1 , anticipated for the long-term operation of the system, and ground loads defined from the building simulation and the assumed heat pump COPs (Equations (1) and (2)), were assigned to predict the operating temperature of the GHE field for twenty years. For all numerical simulations, including model calibration, hourly time steps were used to ensure accurate results.

Model Calibration
The model was calibrated to reproduce the outlet temperature recorded during the heat injection test. Parameters with the highest uncertainty, which are the thermal conductivity of the subsurface solids and the borehole backfill material, as well as the subsurface hydraulic conductivity, were manually adjusted to provide the best match between measured and simulated temperature. The objective was to identify possible ranges for uncertain parameters, which have an influence on the operating GHE temperature. All other parameters were kept constant during the calibration (Table  1), which are basically linked to the subsurface and the GHE entities. The constant parameters are the GHE pipe spacing, inlet and outlet pipe diameters, pipe thermal conductivity, pipe thickness, heat carrier fluid thermal conductivity and density, volumetric heat capacity of the fluid and subsurface solids, thermal conductivity of the fluid and specific storage of the subsurface.

Model Calibration
The model was calibrated to reproduce the outlet temperature recorded during the heat injection test. Parameters with the highest uncertainty, which are the thermal conductivity of the subsurface solids and the borehole backfill material, as well as the subsurface hydraulic conductivity, were manually adjusted to provide the best match between measured and simulated temperature. The objective was to identify possible ranges for uncertain parameters, which have an influence on the operating GHE temperature. All other parameters were kept constant during the calibration (Table 1), which are basically linked to the subsurface and the GHE entities. The constant parameters are the GHE pipe spacing, inlet and outlet pipe diameters, pipe thermal conductivity, pipe thickness, heat carrier fluid thermal conductivity and density, volumetric heat capacity of the fluid and subsurface solids, thermal conductivity of the fluid and specific storage of the subsurface.

Groundwater Flow Conditions
Field observations combined with GHE drilling reports allowed to define a geological cross-section used as a conceptual model ( Figure 5). The 2800 m cross-section is bounded by a stream and the active quarry that are considered constant-head hydraulic boundaries defining the water table at these locations. The thickness of the unconsolidated sediment cover was determined from well logs of the bore field [26]. The abandoned quarry closest to the school is 189 m down-gradient from the GHE field, while the active quarry is approximately 1 km down-gradient. The water level depth in the abandoned quarry was 21 m above sea level (asl) during the field work. Water is pumped irregularly in the active quarry to maintain the groundwater level below the excavation, which affects the local groundwater flow regime. The water level in the active quarry was measured at −3 m asl, but can vary by a few meters.

Groundwater Flow Conditions
Field observations combined with GHE drilling reports allowed to define a geological crosssection used as a conceptual model ( Figure 5). The 2800 m cross-section is bounded by a stream and the active quarry that are considered constant-head hydraulic boundaries defining the water table at these locations. The thickness of the unconsolidated sediment cover was determined from well logs of the bore field [26]. The abandoned quarry closest to the school is 189 m down-gradient from the GHE field, while the active quarry is approximately 1 km down-gradient. The water level depth in the abandoned quarry was 21 m above sea level (asl) during the field work. Water is pumped irregularly in the active quarry to maintain the groundwater level below the excavation, which affects the local groundwater flow regime. The water level in the active quarry was measured at −3 m asl, but can vary by a few meters. Steady-state groundwater flow in a simplified unconfined aquifer system with surface recharge was considered to calculate the level of the water table and the hydraulic gradient near the school site. A total distance of 2800 m separates the active quarry from the stream beyond the school, in the direction of the groundwater flow. Under an assumed constant hydraulic head maintained in the quarry and the stream, the hydraulic head near the GHE field was calculated using the Dupuit  Steady-state groundwater flow in a simplified unconfined aquifer system with surface recharge was considered to calculate the level of the water table and the hydraulic gradient near the school site. A total distance of 2800 m separates the active quarry from the stream beyond the school, in the direction of the groundwater flow. Under an assumed constant hydraulic head maintained in the quarry and the stream, the hydraulic head near the GHE field was calculated using the Dupuit formulation according to [44]:

Subsurface Thermal Conductivity
Gabbro samples taken from a dyke in the quarry and analyzed with both the needle probe and the MTPS have a thermal conductivity below 2.0 W m −1 K −1 , which is lower than the shales and calcarenites, due to the mineralogy of the intrusive rock containing abundant feldspars ( Table 2). Thermal conductivity of gabbro was measured in another case study and ranged between 1.65 and 2.29 W m −1 K −1 [45], which is consistent with our results. Shale samples have an average thermal conductivity ranging between 2.4 and 2.9 W m −1 K −1 ( Table 2), which is in agreement with a value of 2.8 W m −1 K −1 measured for other samples of similar lithology in the St. Lawrence lowlands [10,12,44,46]. The calcarenite sample has the highest thermal conductivity among all the samples collected in the quarry, with a value of 3.5 W m −1 K −1 . Results obtained from the laboratory analyses were considered to constrain the range of the solids thermal conductivity between 2 and 3 W m −1 K −1 in the subsequent numerical model since the thermal response tests reported in Section 2.2 were believed to be affected by groundwater flow.

In-Situ Heat Injection Test
The heat injection test was carried out at an average heat injection rate of 305 kW for the entire bore field or 9.8 kW per borehole ( Figure 6). The total average flow rate was 39.9 L min −1 . The injection period lasted 406 h with GHE water temperature fluctuating between 15 • C to over 35 • C. Monitoring of the temperature recovery, where the fluid is kept circulating in the GHE field without heat injection, lasted 318 h and the GHE water temperature eventually recovered to about 13 • C.

In-Situ Heat Injection Test
The heat injection test was carried out at an average heat injection rate of 305 kW for the entire bore field or 9.8 kW per borehole ( Figure 6). The total average flow rate was 39.9 L min −1 . The injection period lasted 406 h with GHE water temperature fluctuating between 15 °C to over 35 °C. Monitoring of the temperature recovery, where the fluid is kept circulating in the GHE field without heat injection, lasted 318 h and the GHE water temperature eventually recovered to about 13 °C.

Ground Loads
The building simulation using eQuest revealed that the total annual heating energy consumption of the school building is 171 MWh and the total cooling energy consumption is 119 MWh for a regular year of operation. The peak heating load of 494 kW occurs in January while the peak cooling load of 253 kW occurs during July. The building loads were converted to ground loads using Equations (1) and (2) assuming a COP of 4.7 and 4.1 in heating and cooling modes, respectively (Figure 7). The COPs were chosen according to the heat pump specifications provided by the manufacturer and represent average conditions. Heating loads are greater than cooling loads, which makes the ground loads unbalanced and can affect the long-term thermal response of the system.

Ground Loads
The building simulation using eQuest revealed that the total annual heating energy consumption of the school building is 171 MWh and the total cooling energy consumption is 119 MWh for a regular year of operation. The peak heating load of 494 kW occurs in January while the peak cooling load of 253 kW occurs during July. The building loads were converted to ground loads using Equations (1) and (2) assuming a COP of 4.7 and 4.1 in heating and cooling modes, respectively (Figure 7). The COPs were chosen according to the heat pump specifications provided by the manufacturer and represent average conditions. Heating loads are greater than cooling loads, which makes the ground loads unbalanced and can affect the long-term thermal response of the system.

Model Calibration
The total duration for the calibration simulation was 724 h, which corresponds to the duration of the heat injection test, including the heat injection and the thermal recovery periods. Calibration parameters were adjusted one at a time (Table 3) until the calibrated model best reproduced the GHE outlet temperature with a maximum error of 2% (Figure 8). Calculation time using a workstation with a 3.3 GHz-core processor and 32Gb RAM was about 30 min for the calibration simulations.

Model Calibration
The total duration for the calibration simulation was 724 h, which corresponds to the duration of the heat injection test, including the heat injection and the thermal recovery periods. Calibration parameters were adjusted one at a time (Table 3) until the calibrated model best reproduced the GHE outlet temperature with a maximum error of 2% (Figure 8). Calculation time using a workstation with a 3.3 GHz-core processor and 32Gb RAM was about 30 min for the calibration simulations.

System Operation-Twenty-Year Simulations
The eight simulation scenarios used to evaluate the long-term operating GHE temperature at the Carignan-Salières elementary school were based on field observations and calibration results, varying key parameters one at a time to determine the impact of associated heat transfer mechanisms ( Table 4). The start time of the simulation was the month of September, when the GHE system was put in operation. Each simulation took 5 days using the same workstation as reported above for the calibration simulations. Scenario A is a base case scenario, considering measured data for the thermal conductivity of the solids and the borehole filling material as well as hydraulic conductivity of the host rock reported in the literature [23] and a moderate hydraulic gradient, when compared to the hydraulic gradient measured during the field work. Scenario B was defined according to the best-fit parameters identified with the calibration and has a small change of solids thermal conductivity when compared to Scenario A. Only the thermal conductivity of the borehole filling material was changed in Scenarios C and D with respect to the base case. Scenarios E and F focused on thermal conductivity of the host rock solids, while Scenarios G and H were run to verify the influence of groundwater flow with a change of hydraulic head at the eastern and western boundaries of the model (Figure 4). Simulation results for Scenarios A and B do not show significant differences in the simulated GHE fluid temperature at both the outlet and inlet (Figure 9 and Table 5). The small decrease of the subsurface thermal conductivity in Scenario B did not have a significant impact on the results. Scenario D shows a GHE fluid temperature that is 0.4 °C higher than the maximum and lower than

System Operation-Twenty-Year Simulations
The eight simulation scenarios used to evaluate the long-term operating GHE temperature at the Carignan-Salières elementary school were based on field observations and calibration results, varying key parameters one at a time to determine the impact of associated heat transfer mechanisms ( Table 4). The start time of the simulation was the month of September, when the GHE system was put in operation. Each simulation took 5 days using the same workstation as reported above for the calibration simulations. Scenario A is a base case scenario, considering measured data for the thermal conductivity of the solids and the borehole filling material as well as hydraulic conductivity of the host rock reported in the literature [23] and a moderate hydraulic gradient, when compared to the hydraulic gradient measured during the field work. Scenario B was defined according to the best-fit parameters identified with the calibration and has a small change of solids thermal conductivity when compared to Scenario A. Only the thermal conductivity of the borehole filling material was changed in Scenarios C and D with respect to the base case. Scenarios E and F focused on thermal conductivity of the host rock solids, while Scenarios G and H were run to verify the influence of groundwater flow with a change of hydraulic head at the eastern and western boundaries of the model (Figure 4). Simulation results for Scenarios A and B do not show significant differences in the simulated GHE fluid temperature at both the outlet and inlet (Figure 9 and Table 5). The small decrease of the subsurface thermal conductivity in Scenario B did not have a significant impact on the results. Scenario D shows a GHE fluid temperature that is 0.4 • C higher than the maximum and lower than the minimum temperature when compared to Scenario C, which is due to the decrease of the thermal conductivity of the backfilling material. A low thermal conductivity of the subsurface (Scenario E) affects the maximum and the minimum GHE fluid temperature by up to 0.9 • C, when compared to the case with a high thermal conductivity (Scenario F). the minimum temperature when compared to Scenario C, which is due to the decrease of the thermal conductivity of the backfilling material. A low thermal conductivity of the subsurface (Scenario E) affects the maximum and the minimum GHE fluid temperature by up to 0.9 °C, when compared to the case with a high thermal conductivity (Scenario F).   Simulations with differences in hydraulic gradient imposed across the model domain show the most significant differences in GHE fluid temperature. In Scenario G, the GHE inlet temperature reaches 40 • C and drops to −3 • C and the GHE outlet temperature reaches more than 30 • C and drops to less than 2 • C. However, in Scenario H, the GHE inlet temperature reaches 32 • C and falls to only 3 • C while the GHE outlet temperature is between 23 • C and 8 • C. The maximum GHE fluid temperature in Scenario H is 8 • C lower than the maximum of Scenario G, while the minimum GHE fluid temperature is 6 • C higher than minimum of Scenario G. This trend can be explained by the higher imposed hydraulic gradient of 0.008, which was considered for Scenario H that represents conditions with significant pumping in the active quarry. This last scenario provides better heat exchange with the subsurface and therefore better GHE temperature or performance. The lower hydraulic gradient can represent a case where pumping in the active quarry is stopped or reduced, which has a negative impact on the GHE temperature and consequently a potential negative impact on the system performance. Overall, the thermal conductivity of the grouting material has a small impact on twenty-year predictions of the GHE inlet and outlet temperature. The thermal conductivity of the subsurface has a moderate impact on the maximum and minimum GHE temperature. The hydraulic gradient, which is thought to be affected by pumping water in the active quarry, has the greatest impact on heat exchange with the subsurface. Therefore, long-term GHE performance mostly depends on the local groundwater flow conditions. As the minimum GHE outlet temperature drops, the performance of the heat pump can decrease. Despite the significant differences between the simulated temperature of the heat carrier fluid in Scenarios G and H, the twenty-year simulations show an appropriate thermal response of the subsurface with constant temperature changes from year to year, although ground loads are unbalanced. A low groundwater flow rate appears enough to reduce the effect of unbalanced loads that can potentially cool the subsurface since this is not noticed in the long-term simulation. The main impact of groundwater flow is on yearly temperature changes. The minimum GHE outlet temperature dropped to 2 • C in Scenario G, which is still far from the minimal operating temperature of the heat pump system recommended by the manufacturer (−9.62 • C). However, the minimum GHE inlet temperature dropped to −3 • C. The freeze protection provided by the 25 vol. % propylene glycol solution circulating in the GHE is −10 • C and geothermal system designers recommend a minimum fluid temperature 5 to 7 • C higher than the freezing point of the solution. Therefore, under low groundwater flow conditions, care should be taken to follow the minimum operating system temperature during winter periods to avoid potential freezing problems at the GHE inlet.

Discussion
Previous studies that involved the numerical evaluation of GHE performance under the influence of groundwater flow identified the conditions where advection becomes the dominant heat transfer mechanism affecting GHE operation, which has been compared in terms of Darcy flux in the following discussion. For example, Chiasson et al. [47] concluded that groundwater flow can be expected to affect GHE performance, with finite element simulations of groundwater flow and heat transfer that indicated a higher minimum GHE water temperature when increasing Darcy flux from 1.9 × 10 −6 to 1.9 × 10 −5 m s −1 . In addition, Fujimoto et al. [48] studied a GCHP system where significant groundwater flow was shown to have more effect on the thermal response of GHEs when compared to the ground thermal conductivity, even under different building load scenarios. In less permeable subsurface layers, which had a groundwater Darcy flux of 2 × 10 −10 m s −1 , Fujimoto et al. [48] showed that the operating temperature of a GHE exhibits more changes than that of the permeable layer having a Darcy flux of 6 × 10 −7 m s −1 . Another study that used numerical simulations has shown that increasing the Darcy flux in a theoretical aquifer system from 10 −9 to 10 −7 m s −1 resulted in a longer and narrower thermal plume propagating from the GHE [1]. Further increasing the Darcy flux to 10−6 m s −1 made the thermal plume dramatically smaller. Dehkordi et al. [1] concluded that groundwater flow can substantially improve heat transfer by enhancing advection. This was confirmed by Ferguson [2], who provided a screening tool indicating that the boundary between GCHP systems with heat conduction versus advection as a dominant heat transfer mechanism is when the Darcy flux is on the order of 1 × 10 −7 m s −1 .
The GCHP system at the Carignan-Salières elementary school has been operating near this boundary. Constant hydraulic heads imposed at the lateral boundaries of the model presented in this study were varied to evaluate the operating temperature of the GHE under specific Darcy fluxes varying from 8 × 10 −7 to 6 × 10 −8 m s −1 , which represented conditions with high to low groundwater flow rates in situations with and without pumping, respectively, in the neighbouring active quarry. The subsurface heat exchange capacity of the GCHP system appears to be enhanced by groundwater flow, as a consequence of pumping in the active quarry. The simulation carried out under the most likely operating conditions, assuming a thermal conductivity of the grouting material and the subsurface solids of 1.75 and 2.4 W m −1 K −1 , respectively, and with groundwater flow conditions that induce a Darcy flux of 2 × 10 −7 m s −1 at the GHE field, indicates a reasonable operating GHE temperature for the next twenty years. However, if pumping activity ceases in the active quarry and the specific groundwater flux drops below 6 × 10 −8 m s −1 , the minimum GHE outlet temperature can become critical and may drop closer to the minimum heat pump working temperature.
Underground thermal energy storage systems can also be affected by groundwater flow when the Darcy flux is on the same order of magnitude as that reported for GCHP systems (1 × 10 −7 m s −1 ) [11]. However, the impact of groundwater flow on underground thermal energy storage systems is negative because it increases heat loss around GHEs used as a heat storage medium. Giordano and Raymond [10] demonstrated that~10% of the energy lost by the bore field of an underground energy storage system can be due to advection when the Darcy flux is 8 × 10 −7 m s −1 , while this heat loss can be reduced bỹ 60 % if the connection and layout of the GHEs is designed considering the magnitude and direction of groundwater flow.
The present study further highlights the potential of groundwater flow to help operating GHEs with unbalanced ground loads, which were affected in this study by~30 % more heat extraction due to greater heating than cooling of the building. The minimum GHE water temperature is expected to decrease year after year when simulating heating-dominated GCHP systems using a heat conduction approach to represent heat transfer in the subsurface. Jaziri et al. [49] previously simulated the operation of the GCHP system at the Carignan-Salières elementary school with such a heat conduction approach based on the HyGCHP program [50], where GHEs were represented with the duct storage model with no groundwater flow contribution [51]. The GHE minimum operating temperature at the beginning of the simulations was similar to simulations made with FEFLOW for the case with a low groundwater flow rate (Scenario G) but decreased by 4 to 6 • C over the twenty-years of system operation. On the other hand, GHE simulations made with FEFLOW and considering advection did not show a significant decrease of the minimum GHE water temperature over the expected life-span of the system, even with a low groundwater flow rate (Figure 9). This finding has important implications for GCHP system design. When calculating the required GHE length according to building energy needs, and considering the thermal state and properties of the subsurface [52], practitioners will ensure that the desired minimum GHE temperature is not reached before 10 to 20 years, assuming heat conduction only in the subsurface. This approach can overestimate the required GHE length if advection becomes significant. The present study suggests that this can be the case, even at a low groundwater flow rate with a Darcy flux on the order of 10−8 m s −1 . Possible approaches to design GCHP systems under the influence of groundwater flow are to assume heat conduction only, but calculate the GHE length for a one-year heat pulse [53] or consider advection with a moving line-source solution to complete the sizing calculation [54,55].

Conclusions
Field characterization and modelling work conducted at the Carignan-Salières elementary school has shown that the installed GCHP system can be influenced by the varying groundwater flow conditions due to dewatering of a neighbouring quarry located approximately 1 km down-gradient of the borehole field. The specific groundwater Darcy flux in the host rock, in which the GHE system is installed, has been estimated near 10 −7 m s −1 at the time when the water levels in the nearby quarries were measured. The school building is constructed above fractured shale and limestone of the Nicolet Formation in the Lorraine Group of the St. Lawrence Lowlands and the groundwater flow is oriented westward toward the active quarry. Measurement of thermal conductivity of the rock samples collected in the active quarry revealed that the most abundant rock types have a thermal conductivity ranging from 1.8 to 3.6 W m −1 K −1 . A TRT conducted at the site, which reported an equivalent subsurface thermal conductivity of 2.58 W m −1 K −1 [26], was likely affected by groundwater flow. A heat injection test was therefore carried out with the purpose of evaluating the thermal response of the entire GHE field subject to groundwater flow. The data collected were used to calibrate a groundwater flow and heat transfer model built using FEFLOW to simulate the operating temperature of the GHEs. The numerical model was developed according to the characteristics of the subsurface and the GHE field to reproduce the heat injection test and then predict the long-term response according to different scenarios.
The parameters that impact the GHE operating temperature in order of increasing importance are: the thermal conductivity of the backfilling material, the thermal conductivity of the subsurface and the groundwater flow rate. The impact of the thermal conductivity of the borehole filling material and the subsurface (Scenarios B to F) on the inlet and outlet temperature of the GHEs was actually shown to be minor when compared to the possible change of groundwater flow conditions at the site (Scenarios G and H). The Carignan GCHP system is a unique field case with GHEs 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 flow inducing significant advective heat transfer when the specific Darcy flux increases from 6 × 10 −8 m s −1 (no dewatering) to 8 × 10 −7 m s −1 (high dewatering). This study further suggests that even the lowest groundwater flow rates expected at the site can be beneficial to avoid a progressive cooling of the system over its expected lifetime due to the unbalanced heating and cooling loads.
Further research is necessary to define the conditions where groundwater flow can help to cope with unbalance ground loads. While it is clear that advection plays an important role on heat transfer when the Darcy flux is 10 −7 m s −1 or higher, a lower groundwater flow rate may be sufficient to ensure a constant minimum GHE temperature year after year, depending on the magnitude of the unbalanced loads.
Author Contributions: N.J. did the field, laboratory and numerical modelling work and wrote the manuscript. J.R. and J.M. supervised the work and improved the original manuscript prepared by N.J., N.G. helped with numerical modelling and contributed to the revision of the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by a BMP-Innovation scholarship given to Nehed Jaziri from the Fonds de recherche du Québec-Nature et technologies, the Natural Sciences and Engineering Research Council (NSERC) and the company GBi as well as an NSERC Discovery grant given to Jasmin Raymond.

Acknowledgments:
The support of GBi, particularly Maxime Boisclair and Jérôme Plante, and the Commission scolaire des Patriotes, through Jean-François Rondeau, who facilitated this study, are acknowledged.