Simulating thermal pollution caused by a hypothetical groundwater heat pump system under different climate, operation and hydrogeological conditions

Abstract

Thermal–hydraulic analysis was conducted to examine the propagation of temperature anomalies under different climatic, subsurface and operation conditions. A series of scenarios simulated by numerical modeling indicate that the injection rate, injection depth and anisotropy of the aquifer are important parameters that need to be properly estimated to avoid critical thermal interference in shallow subsurface areas. The potential impact due to climate change is predominant only in shallow subsurface areas, but the situation may become severe if thermal anomalies that propagate upward due to injection of warm water mix with climate change anomalies that penetrate downward. Groundwater recharge and horizontal groundwater flows retard upward migration of temperature anomalies. Therefore, a proper analysis of groundwater flow in the area would facilitate harvesting a large amount of geothermal energy while limiting its impact to acceptable levels. Furthermore, a set of type curves was developed to select operational conditions with a reasonable accuracy while maintaining thermal pollution in an acceptable range.

Introduction

Subsurface temperatures, from the ground surface to thousands of meters depth, have been measured in many parts of the world for commercial purposes such as exploitation of geothermal energy and for non-commercial purposes such as ecological research. These measurements indicate that aquifer temperatures are 1–5 °C higher in urban areas than in rural areas (Ferguson and Woodbury 2007; Wang et al. 2009; Gunawardhana et al. 2011). This difference can be attributed primarily to urbanization effects that increase aquifer temperatures by land use changes, warm industrial water injection and sewage leaks. Gunawardhana and Kazama (2015) found that land use change alone may account for approximately three-fourths of the total warming of the ground surface in five urban areas in Japan. More specifically, Menberg et al. (2013) found that heat losses from basement of buildings are one of the dominant factors to influence subsurface thermal regime. Moreover, climate change is expected to increase surface air temperature and change the precipitation regime leading to frequent and intense extreme events. The Fifth Assessment Report (AR5) of the United Nations Intergovernmental Panel on Climate Change (IPCC) concluded that the rise in the mean global surface air temperature by the end of the twenty-first century relative to the pre-industrial period is likely to be between 1.5 and 4.5 °C (Symon 2013). These projected changes with continuing urban growth will significantly alter the subsurface thermal regime (Kurylyk et al. 2014).

Although an increase in subsurface temperature above the natural range is considered pollution, enhancing geothermal potential is indisputably important for energy independence and coping with global warming (Haehnlein et al. 2010). For example, Zhu et al. (2010) showed that the minimum extractable geothermal yield is at least 2.5 times greater than the total annual residential heating demand in Cologne, Germany.

Due to the low thermal diffusivity of aquifers relative to air, subsurface temperatures experience less fluctuation than surface air temperatures. Therefore, geothermal energy systems use aquifers as heat sources during the winter months and as heat sinks during the summer months. In principle, there are two types of geothermal energy systems: the ground source heat pump (GSHP) and the groundwater heat pump (GWHP). The GSHP is a closed system with a heat-carrying fluid circulating through a tube system installed underground. Heat is exchanged between the aquifer and the GSHP system, but the groundwater level remains unchanged. The GWHP is an open system, within which groundwater is extracted and sent through the heat pump installed at the ground surface. Heat is released or extracted depending on the season, and groundwater is re-injected within a short time period back into the aquifer (Haehnlein et al. 2010). Even though the cones of depression from groundwater withdrawal wells are likely to affect local-scale groundwater flow, the water budget does not change on a regional scale. However, the issue of changing the subsurface temperature around an injection well is crucial in some respects. Increasing the number of geothermal energy systems within a small area can affect the efficiency of the individual systems (Ferguson and Woodbury 2006). Moreover, effects of temperature changes on the groundwater ecosystems such as estuaries, wetlands and ponds can be critical (Brielmann et al. 2009; Bonte et al. 2011).

Haehnlein et al. (2010) reviewed the legal status of the use of geothermal energy in 46 countries and found that in 35 of those countries, there were no regulations or recommendations for subsurface energy use. Some other countries impose legal regulations or recommendations based on the minimum distance to the next property line (for example, 2.5 m in Austria and 10 m in Finland), geothermal system (for example, 5–10 m in Germany) or building (for example, 5 m in Greece) as well as the temperature threshold for environmental sensitivity. For example, the allowable groundwater temperature change at a distance of 50 m from a geothermal system is less than 2 °C in the city of Stuttgart in southern Germany (Paly et al. 2012), and a minimum distance of 300 m is required between a geothermal system and a drinking groundwater well in Denmark. Haehnlein et al. (2010) also noted that the legal status of the use of geothermal energy can vary considerably even among the different states of a single country and that the defined threshold values are generally not based on reliable scientific analysis but rather appear to be arbitrary. Several recent studies have attempted to address these issues by conducting numerical modeling of geothermal systems under specific site conditions. Ferguson and Woodbury (2006) studied the use of groundwater for cooling purposes in four buildings in Canada and found that the increases in temperature observed were due to the spacing between injection wells being smaller than the optimal spacing for such systems. Kim et al. (2010) evaluated the performance of geothermal energy systems for different site conditions and found that their performance can be significantly affected by varying the distance between two wells and the pumping/injection rates. However, to date, most studies of geothermal systems have been conducted under specific site conditions and have addressed the influence of limited parameters. The objective of this study was to understand the distribution of thermal anomalies under various subsurface (e.g., horizontal and vertical groundwater flow), operation and climatic conditions. One of the unique aspects addressed in this study is the effect of ground surface warming that has been attributed to urbanization and climate change in geothermal studies. The results of this study will provide a set of type curves to select operational conditions of geothermal energy systems under different site conditions in a manner that minimizes their potential impact on the environment.

Numerical model setup

The VS2DH numerical code (Healy and Ronan 1996), which can simulate heat transport for unsteady, non-uniform water movement with variable boundary conditions, was used to simulate a site with two injection wells (IW1 and IW2, as shown in Fig. 1a) located in two neighboring private properties. The governing equation for subsurface temperature distribution can be written as follows:

$${\partial \mathord{\left/ {\vphantom {\partial {\partial t\left[ {\uptheta C_{W} + \left( {1 - \varphi } \right)C_{S} } \right]}}} \right. \kern-0pt} {\partial t\left[ {\uptheta C_{W} + \left( {1 - \varphi } \right)C_{S} } \right]}}T = \nabla { \cdot }K_{T} \left( \uptheta \right)\nabla T + \nabla{\cdot}\uptheta C_{W} D_{H} \nabla T - \nabla \uptheta C_{W} vT + qC_{W} T^{ * }$$

(1)

where $t$ is the time in sec, $\uptheta$ is the volumetric moisture content, $C_{W}$ is the heat capacity of water in J/m³ °C, $\varphi$ is the porosity, $C_{S}$ is the heat capacity of the dry soil in J/m³ °C, $T$ is the temperature in  °C, $K_{T}$ is the thermal conductivity of the water–solid matrix in W/m °C, $D_{H}$ is the hydrodynamic dispersion tensor in m²/sec, $v$ is the water velocity in m/sec, $q$ is the rate of the fluid source in sec⁻¹, and $T^{ * }$ is the temperature of the fluid source in  °C.

Fig. 1

a Hypothetical site conditions, b Initial and boundary conditions of the numerical model, c Surface air temperature derived from 6 GCM models and 2 scenarios for time-varying ground surface temperatures at the top boundary (Source: Gunawardhana and Kazama 2012)

Gunawardhana and Kazama (2012) obtained good agreement between temperature–depth (T–D) profiles simulated using the VS2DH numerical model and T–D profiles observed in the Sendai plain in Japan. In this study, the same initial, boundary and geological conditions were used in a two-dimensional VS2DH numerical model to examine various geothermal input scenarios and site conditions (Fig. 1b). Depending on the season, cold or warm water is extracted through the production well (PW), and thermally altered water is injected back into the aquifer through the IW. For the sake of simplicity, we simulate the effect of injecting thermally altered water only. The injection rate and the temperature of the injection water are represented by the parameters q and T*, respectively, in Eq. (1).

Boundary conditions

Injection rate and timing

The range of injection rates assigned to the model boundaries (q _i1 and q _i2, as shown in Fig. 1b) was 50–150 m³/day (Kim et al. 2010). The injection timing was varied in the different scenarios studied. Cyclic operation of the geothermal system was studied by simulating the injection of warm water at 35 °C in the summer (4 months) and cold water at 10 °C in the winter (4 months) (Haehnlein et al. 2010). The system operation during the rest of the year was not considered. In some cases, geothermal systems are operated during the summer months only (Ferguson and Woodbury 2006). To simulate such a scenario, it was assumed that the geothermal system is in operation for 4 months only. Another scenario, in which injection of warm water throughout the year was simulated, was considered to represent geothermal energy use for cooling purposes only (e.g., commercial buildings).

Horizontal and vertical water flow

Heat carried by groundwater recharge and regional groundwater flow is important for heat distribution in subsurface layers (Lu and Ge 1996; Bense and Beltrami 2007). In this study, the effect of groundwater recharge (8 × 10⁻⁹ m/sec) was simulated by specifying constant pressure head conditions at the top and bottom boundaries of the model to induce a hydraulic gradient. Similarly, the effect of regional groundwater flow was simulated by specifying constant head boundaries of 50 m height between the left and right sides of the model (Fig. 1b). Two horizontal water flow (HWF) scenarios, 3 × 10⁻⁵ and 3 × 10⁻⁷ m/sec, were generated by varying the pressure head difference between two boundaries and the hydraulic conductivity of the porous medium. To investigate the individual effects of recharge and horizontal water flow on the subsurface temperature distribution, the vertical anisotropy of the hydraulic conductivity was assigned values of 10² and 10⁻², respectively.

Ground surface temperature change

The ground surface temperature is mostly affected by land use changes (due to changes in albedo, evapotranspiration and thermal insulation) and surface air temperature changes (due to non-local effects such as global climate change) (Ferguson and Beltrami 2006; Bense and Beltrami 2007; Gunawardhana et al. 2011). Therefore, variations in the ground surface temperature generally reflect changes in the surface air temperature and vice versa, and in long-term comparisons, variations in the ground surface temperature follow trends that are similar to those of variations in surface air temperature (Huang et al. 2000). However, exceptions may occur due to the changes in timing and duration of snow cover (Kurylyk et al. 2013). Due to heat capacity differences, surface air temperatures exhibit more variation than ground surface temperatures. In this study, surface air temperature was used as a proxy for ground surface temperature estimation and was filtered using a 5-year moving average (Ferguson and Woodbury 2005; Majorowicz et al. 2006). This filtering reduces the ground surface temperature variation with respect to the surface air temperature variation. According to meteorological records from the Sendai plain, surface air temperatures exhibited no significant trend until the middle of the twentieth century but started increasing rapidly beginning in 1947. Therefore, 1947 surface air temperatures were taken to compute initial conditions in the numerical simulation performed in this study. Until 2000, records of surface air temperature observations were used, and from 2000 to 2046, surface air temperature results from 6 general circulation models (GCM) and 2 emission scenarios (12 scenarios in total) were used to approximate variations in ground surface temperatures (Fig. 1c). More details of the numerical model assumptions, parameter values and methodology for downscaling GCM output to the local scale can be found in Gunawardhana and Kazama (2012).

Model structure and initial conditions

Two IWs were placed sufficiently far apart (1000 m) to detect their individual effects on temperature distribution (Fig. 1b). A general injection depth (ID) of 150 m was specified for thermally altered water injection, and this depth was later changed to investigate the effect of IW depth on shallow subsurface temperature change. The horizontal and vertical grid spacing was 1 m near the injection point and increased away from the injection point up to maximum grid spacing of 10 m. The time step was started at 1 day and was allowed to increase to 30 days. The total simulation time was 100 years. The initial temperatures at the top and bottom of the model were assigned by extrapolating the undisturbed linear portion of the observed T–D profile. The temperature at the top boundary was allowed to change with time (as the ground surface temperature changed), while the temperature at the bottom of the model was held constant throughout the simulation time.

Results and discussion

In this study, the subsurface temperature change due to injection of thermally altered water from geothermal energy use was assessed by considering different scenarios of operation conditions, subsurface characteristics and ground surface climate change. The degree of thermal interference in each case was compared to the movement of a 2 °C temperature anomaly (Paly et al. 2012) in the horizontal and vertical directions after 10 and 100 years. Due to the symmetric nature and large lateral extent of the model boundaries (Fig. 1b), propagation of temperature anomalies in only half of the model is shown for the following scenarios.

Ground surface warming due to global climate change

According to the results for 12 GCM scenarios, the surface air temperature may increase by 1.3–4.7 °C during the 2060–2099 time period, compared to the observed averages between 1967 and 2006. Figure 2a shows the propagation of subsequent temperature anomalies downward into the subsurface as predicted for the highest temperature ground surface warming scenario. Approximately 1.5–2.0 °C warming can be observed within a depth of 4 m from the ground surface 100 years after the simulation starting year of 1947. However, the magnitude of the temperature change attenuation increases with depth, and the temperature change becomes insignificant (>0.1 °C) beyond a depth of 50 m. According to some lower impact warming scenarios extend over 100 years, signatures of ground surface warming dissipate at depths less than 35 m. Stallman (1965) explained the theoretical relationship between the attenuation of a harmonic surface temperature signal and the depth and applied this concept in developing an analytical solution for subsurface temperature distribution with transient ground surface temperature change. Since then, researchers in many parts of the world have confirmed similar phenomena by monitoring aquifer temperatures and have found greater temperature anomalies closer to the ground surface than at greater depths, primarily due to urbanization and climate change effects (Taniguchi et al. 1999; Huang et al. 2000; Ferguson and Woodbury 2005; Majorowicz et al. 2006).

Fig. 2

Development of temperature anomalies: a after 100 years due to ground surface temperature change as predicted by the highest impact scenario, b after 10 years due to 35 °C water injection at a 50-m³/day rate and a 150-m depth, c after 100 years due to 35 °C water injection at a 50-m³/day rate and a 150-m depth, d after 10 years due to 35 °C water injection at a 150-m³/day rate and a 150-m depth, e after 100 years due to 35 °C water injection at a 150-m³/day rate at a 150-m depth, and f after 100 years due to 35 °C water injection at a 150-m³/day rate and a 100-m depth

Injection rate

Two injection rates (50 and 150 m³/day) were selected to represent domestic and small-scale industrial geothermal use (Haehnlein et al. 2010). Figure 2b–f shows temperature distributions in the subsurface after 10 and 100 years of injection of warm water. Injection water temperatures are initially distributed in the shape of a symmetrical bubble. A hotspot generated with overheated convective flow appears to be stronger in the central part of the bubble and expand with time. Based on field data from over 30 years of observation of an area that has been subjected to geothermal developments, Ferguson and Woodbury (2006) showed that aquifer temperatures gradually increase over time and that temperature anomalies propagate in a manner similar to that presented here. In this study, when the injection rate was 50 m³/day, the thermal front propagated over 50 m from the injection well after 10 years and over 110 m after 100 years. Beyond these distances, the temperature change was less than 2 °C compared to the initial temperature. When the injection rate was increased to 150 m³/day, for the same conditions of hydraulic conductivity and injection depth, the thermal front propagated over 75 and 160 m in 10 and 100 years, respectively. Notably, the thermal bubble started to move upward (e.g., approximately 10 and 40 m, as observed in Fig. 2d, e, respectively). This phenomenon is expected because energy flows from a high-temperature region to a low-temperature region, which in this case is from a greater depth toward the ground surface. Once temperature anomalies generated by injection of warm water start to mix with anomalies resulting from ground surface temperature changes in shallow subsurface areas (depths <10 m), the thermal bubble begins to lose its initial symmetry, resulting in anomalies concentrating upward while expanding at the bottom (Fig. 2e). This phenomenon implies that the thermal front propagates more with increasing injection rate and time, and may cause additional stress in shallow subsurface areas if mixed with the effects of anthropogenic changes on ground surface temperature. However, the impact on shallow subsurface areas remains very localized. The subsurface temperature does not change by more than 2 °C at distances greater than 60 m from the injection point.

Injection depth

Figure 2f shows the distribution of temperature anomalies when the injection depth is 100 m, which is 50 m shallower than the scenario presented in Fig. 2e, with all other parameters held constant. Comparison of these two figures reveals that the shallow subsurface temperature change is stronger when the injection is conducted at a shallow depth. In the 100-m injection depth scenario (Fig. 2f), a temperature anomaly as high as 14 °C reaches a depth of 20 m after 100 years, but over the same time period, the 150-m injection depth scenario results in a temperature change of only 8 °C at a depth of 20 m. However, when water is injected at 50 m³/day rate at a 150-m depth, the temperature change at a depth of 20 m is less than 2 °C (Fig. 2c). These results indicate that the injection depth is an important parameter that needs to be selected by considering the injection rate, the aquifer properties and the lifetime of the geothermal system.

Groundwater recharge

Figures 2e and 3b show the temperature anomalies associated with the effect of groundwater recharge, with other input and boundary conditions held constant. The scenario in Fig. 3b was produced with groundwater recharge at a rate of 8 × 10⁻⁹ m/sec, which is typical of recharge rates observed in the study region (Gunawardhana et al. 2011), while the scenario of Fig. 2e was produced without a groundwater recharge effect. Figure 2e shows that a temperature anomaly of 4 °C reaches upward to a depth of 5 m after 100 years. However, with the presence of significant groundwater recharge, a temperature anomaly of 4 °C reaches a 25-m depth, and a 2 °C anomaly reaches a depth of only 22 m after 100 years. This happens because the downward heat advection due to groundwater recharge is opposing the upward conductive heat flux from the injection point. On the other hand, if the warm water is injected into the very shallow subsurface, groundwater recharge could warm up aquifers below the injection point more quickly.

Fig. 3

Distribution of temperature anomalies in the presence of groundwater recharge (GWR) and horizontal water flow (HWF). Injection rate qi = 150 m³/day and injection depth ID = 150 m for all scenarios

Horizontal groundwater flow

Vertical hydraulic gradients and temperature–depth profiles abruptly change in the presence of regional groundwater flow, which typically occurs at greater aquifer depths as horizontal flow. The distribution of temperature anomalies was evaluated when injection water was allowed to mix with natural horizontal flow. Figure 3c–f shows that the thermal front propagates farther downstream with stronger groundwater flow. A groundwater flow of 3 × 10⁻⁵ m/sec carries a temperature anomaly of 2 °C up to 300 m after 100 years, but when the groundwater flow rate is 3 × 10⁻⁷ m/sec, the same anomaly extends only 250 m after 100 years. This is because strong groundwater flow enhances lateral advective heat transport, which takes excess heat away from the source. Therefore, Fig. 3f shows less excess heat presence near the injection point than does Fig. 3d. The lateral extents of temperature anomalies due to horizontal water flow determined in this study were generally in agreement with those reported by Bense and Beltrami (2007), who found that temperature anomalies caused by increased surface temperature due to deforestation coupled with groundwater flow at a rate of 1 × 10⁻⁷–1 × 10⁻⁸ m/sec can extend several hundred meters away from the deforested area after 100 years. Moreover, unlike the asymmetric temperature bubbles produced in the conduction-dominated heat transport scenarios (Fig. 2b–f), temperature bubbles due to advection-dominated heat transport scenarios seem to stretch extensively in the lateral direction (Fig. 3a–d) while posing absolutely no risk of reaching shallow subsurface areas.

Operational methods

All of the scenarios considered thus far assumed that warm water was injected continuously throughout year. However, domestic geothermal systems generally operate in a cyclic mode. Therefore, simulation was performed by assuming water at two different temperatures (35 °C in the summer and 10 °C in the winter) was injected for 4 months in each case (Fig. 4a). In this scenario, the 2 °C temperature anomaly propagated laterally 30 m less than in the scenario that considered system operation throughout the year (Fig. 2e). However, the symmetric shape was deformed due to the cyclic nature of the water injection at the two temperatures. During the warm water injection phase, heat travels preferentially along the direction of the geothermal gradient toward the ground surface. In contrast, during the cold water injection phase, cool water mixes with naturally warm water deep in the aquifer. Similarly, the scenario that assumes warm water injection for 4 months only extends its 2 °C temperature anomaly approximately 120 m in the lateral direction (Fig. 4b), which is 40 m less than in the scenario that considered system operation throughout the year (Fig. 2e) and 10 m less than in the scenario with water injected at two temperatures (Fig. 4a).

Fig. 4

Distribution of temperature anomalies when subjected to different operational methods: a every year, water at 35 °C was injected for four summer months and water at 10 °C was injected for four winter months. System operations for the rest of the months of the year were not considered. b Warm water at 35 °C was injected for four summer months every year. System was assumed to be shut down for rest of the months. c Warm water at 35 °C was injected at a 50 m³/day rate throughout the year for 50 years in the presence of comparatively small HWF. d Conditions are similar to scenario c, but injection was stopped after 50 years and the temperature was allowed to be distributed with HWF. The injection depth ID = 150 m for all scenarios

The horizontal component of the groundwater flow is comparatively small in extended flat areas and, therefore, its effect on attenuating the impact of injection water in shallow subsurface areas is weak. Operation of a geothermal system under such conditions was simulated by considering a reasonably small horizontal water flow (1 × 10⁻⁸ m/sec, as in Bense and Beltrami 2007). Figure 4c shows the resulting temperature distribution after 50 years of continuous injection of warm water at a temperature of 35 °C and a rate of 50 m³/day. The 2 °C temperature anomaly propagates up to 140 m in the lateral direction. After 50 years of simulated operation under these conditions, warm water injection was stopped and simulation was continued for another 50 years to understand how temperature anomalies may attenuate with the specified groundwater flow. As observed in Fig. 4d, a hotspot thermal bubble generated by 50 years of operation is diluted and its recorded maximum temperature change decreases from 13 to 7 °C. Such temperature anomalies appear to pose no significant impact even if mixed with a significant groundwater discharge flow and spread in shallow subsurface areas.

Sensitivity for anisotropy of porous medium

The anisotropy of the porous medium (the ratio of vertical to horizontal hydraulic conductivity) influences vertical and horizontal groundwater flow rates. When the anisotropy is significant, advective heat travels preferentially along the direction of higher hydraulic conductivity. A practical example is a carbonate aquifer, in which hydraulic conductivity is greater when fractures are aligned with the predominant flow directions (Knochenmus and Robinson 1996). This type of situation was modeled using values of 0.01 and 100 for the vertical anisotropy (Fig. 5). The hydraulic conductivity in the lateral direction was fixed for both scenarios at the same value (1 × 10⁻⁸ m/sec) used in the scenarios presented earlier (Fig. 1a–f). According to Fig. 5a, temperature anomalies propagate more in the vertical direction when the vertical hydraulic conductivity is larger, which limits the lateral extent of a temperature anomaly of 2 °C to only 120 m from the injection point. Although the impact in the lateral direction is comparatively small, the temperature changes in shallow subsurface areas become significant: for example, at a depth of 10 m, the temperature change is approximately 3 °C higher than in the scenario that assumes the same hydraulic conductivity in both the vertical and horizontal directions (Fig. 1e). Kim et al. (2010) presented a similar finding: the propagation of a thermal front increases with decreasing hydraulic resistance. Similarly, when the hydraulic conductivity in the lateral direction is greater than that in the vertical direction, temperature anomalies propagate more in the lateral direction (Fig. 5b). Consequently, a temperature anomaly of 2 °C extends up to 260 m away from the injection point, which is approximately 100 m farther than in the scenario that assumes the same hydraulic conductivity in both directions (Fig. 1e). Therefore, the anisotropy of the porous medium needs to be determined with reasonable accuracy to avoid affecting the performance of neighboring geothermal systems (Ferguson and Woodbury 2006).

Fig. 5

Distribution of temperature anomalies in the presence of significant anisotropy. a Vertical hydraulic conductivity is 100 times larger than that in horizontal direction. b Horizontal hydraulic conductivity is 100 times larger than that in vertical direction

Application of results

In the absence of generally accepted guidelines and recommendations for subsurface energy use, we developed a set of type curves to select appropriate operational conditions while maintaining thermal pollution in an acceptable range. The allowable degree of thermal interference at a depth of 20 m after 100 years of operation time is considered to be 2 °C (Paly et al. 2012). The numerical model was then run with different combinations of injection temperatures, injection rates and injection depths under four common geothermal gradients to obtain the predetermined 2 °C warming at the 20 m depth. Developed curves are shown in Fig. 6. When the injection depth is as shallow as 50 m, injection water at a lower temperature of about 17 °C with a smaller rate of 50 m³/sec would cause 2 °C warming at the 20 m depth. Potential injection temperature increases exponentially as the injection depth increases. However, comparatively smaller injection temperature differences were estimated when the model was run under different geothermal gradients. When all combinations were compared, maximum of 1.5 °C injection temperature difference was calculated between 1 °C/100 m and 4 °C/100 geothermal gradients. One reason for this comparatively small temperature differences could be the shallow depth we considered, where subsurface temperature difference at the 20 m depth between two geothermal gradients is only 0.6 °C.

Fig. 6

Different installation and operational conditions of a geothermal system when allowable temperature change at a depth of 20 m is 2 °C

When these type curves were developed, certain hydrogeological and climatic conditions were assumed. Firstly, thermal and hydraulic properties of the aquifer are assumed to be homogeneous. Secondly, groundwater recharge was neglected. Finally, climate is assumed to be static. According to our previous discussion, groundwater recharge could slow down the propagation of temperature anomalies. Similarly, anisotropy of the porous medium proved to have a significant effect on propagation of temperature anomalies. Moreover, climate change impacts could be significant at shallow subsurface layers. These type curves are, therefore, intended to use in places where computational resources are not available to simulate thermal pollution under specific site conditions. Furthermore, they can be used as initial conditions if the model is developed with different site conditions, which eventually, could reduce the computational effort.

Conclusions

The distribution of temperature anomalies within the subsurface during geothermal energy production depends on the rate of injection, operational methods and subsurface conditions, which if not properly assessed, may not only influence the performance of neighboring geothermal systems but also pose a serious threat to surrounding ecosystems. In this study, we simulated the propagation of temperature anomalies due to warm water injection under different scenarios of site, operational and climatic conditions.

Over 100 years of simulation time, climate change impacts below the ground surface will be limited to shallow subsurface areas only. Domestic geothermal systems could easily impact neighboring systems if operated at high production rates for long periods of time. In an urban area, this range of influence may cross the boundaries of several neighboring systems. Geothermal systems that are operated at high production rates (e.g., commercial and industrial users) cause temperature anomalies to extend farther laterally and propagate toward shallow subsurface areas. Switching a geothermal system from continuous operation to cyclic operation or operation only in the summer months reduce the lateral extent of temperature anomalies but still poses a risk of propagation to shallow subsurface areas. Furthermore, greater shallow subsurface temperature change occurs when the injection is conducted at shallower depths.

Groundwater recharge and horizontal groundwater flows play positive roles in retarding upward migration of thermal anomalies. Similarly, strong horizontal groundwater flows tend to dissipate slow-moving excess heat bubbles in such a way that there is no risk of them reaching shallow subsurface areas. Therefore, our results suggest that installation of a geothermal system after a proper analysis of groundwater recharge and horizontal water flow in the area would facilitate harvesting a large amount of geothermal energy while limiting the impact of the geothermal system to acceptable levels.

Anisotropy of the porous medium also plays an important role when injection heat flux is distributed within the subsurface. A critical situation may occur in shallow subsurface areas if hydraulic resistance in the vertical direction is lower than that in the horizontal direction. A severe situation may arise in shallow subsurface areas if the thermal anomalies that propagate upward mix with climate change anomalies that penetrate downward.

The set of type curves developed in this study can be used to select operational conditions with a reasonable accuracy while maintaining thermal pollution in an acceptable range. However, the injection rate and depth of a geothermal system should be designed based on site conditions, and the system should be regulated when installed in ecologically sensitive areas. The installation of a periodic monitoring system is recommended, especially in groundwater discharge areas where ecologically critical temperature anomalies could easily transport from deeper areas to surface ecosystems by the advective effect of groundwater discharge.