Next Article in Journal
Redox Dependent Arsenic Occurrence and Partitioning in an Industrial Coastal Aquifer: Evidence from High Spatial Resolution Characterization of Groundwater and Sediments
Next Article in Special Issue
Numerical Simulation of Shallow Geothermal Field in Operating of a Ground Source Heat Pump System—A Case Study in Nan Cha Village, Ping Gu District, Beijing
Previous Article in Journal
Simulation of Water-Use Efficiency of Crops under Different Irrigation Strategies
Previous Article in Special Issue
Distortion of the Estimated Hydraulic Conductivity from a Hydraulic Test in Fractured Rock Due to Excessive Injection or Extraction
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Evaluation of Ground Temperature Changes by the Operation of the Geothermal Heat Pump System and Climate Change in Korea

Groundwater Research Center, Korea Institute of Geoscience and Mineral Resources, Daejeon 34132, Korea
*
Author to whom correspondence should be addressed.
Water 2020, 12(10), 2931; https://doi.org/10.3390/w12102931
Submission received: 30 August 2020 / Revised: 15 October 2020 / Accepted: 16 October 2020 / Published: 20 October 2020

Abstract

:
To evaluate long-term temperature changes caused by the operation of a geothermal heat pump (GHP) system, temperatures near borehole heat exchangers (BHEs) of the GHP system in Korea were measured. The temperature measurements showed increasing rates of 0.135 °C/year at a depth of 10 m and 0.118 °C/year at a depth of 50 m for approximately 10 years. Simulations for the analysis of climate change effects on measured temperature fluctuations showed that a rate of temperature increase was 0.010 °C/year at a depth of 50 m owing to changes in surface air temperatures (SATs). From two-dimensional heat transfer simulations, the discharged heat measuring 16.7 W/m in the cooling season and extracted heat measuring 12.4 W/m in the heating season could cause an annual mean temperature increase of 0.109 °C over approximately 10 years. Additionally, results of simulations for future prediction of ground temperatures assuming that the GHP system retains its level of operation showed that in 2050, temperature at a depth of 50 m will increase by approximately 3.00 °C from that in 2005. Thus, balancing the heat discharged into and extracted from the ground by considering climate change to minimize long-term changes in the ground temperature is necessary.

1. Introduction

The use of renewable energy sources is rapidly increasing worldwide with increasing awareness regarding global climate change. The geothermal heat pump (GHP) system is one of the promising renewable energy technologies; it has drawn increased attention recently because it is more environmentally friendly and cost-effective compared to other systems that use electricity or fossil fuels [1,2,3,4,5,6,7].
For the sustainable use of the GHP system, long-term (>~30 years) efficiency should be evaluated in consideration of the lifespan of the system and the thermal response time of the ground. The long-term efficiency of the GHP system depends on several factors, such as thermal interference between boreholes, thermal conductivity of the ground, groundwater flow, ground temperature change caused by climate change, and the operation of the GHP system [8]. Therefore, to elucidate the overall performance of the GHP system, evaluation of the long-term efficiency of the GHP system should consider these factors.
In previous studies, the long-term efficiency of the GHP system was mainly evaluated via numerical modeling including factors such as thermal interference, thermal conductivity, and groundwater flow [1,9,10,11,12]. However, changes in the ground temperature due to climate change have rarely been considered in the evaluation of the GHP system.
When the surface air temperature (SAT) changes because of climate change, the ground surface temperature (GST) changes accordingly, and these changes are transferred to the ground, eventually changing the ground temperature [13,14,15]. Therefore, the identification of changes in ground temperature caused by climate change is essential to the elucidation and evaluation of the long-term efficiency of the GHP system.
In addition, long-term monitoring of changes in ground temperature resulting from the operation of the GHP system and climate change are uncommon. In this study, ground temperature changes caused by both operation of the GHP and climate change were measured via long-term temperature monitoring and evaluated through numerical simulations of the GHP system. The objective of this study is to provide valuable information on the design of the GHP system for long-term sustainable use and predict future environmental changes in the ground because of the operation of the GHP system.

2. Materials and Methods

2.1. Study Area

The study site is located at the Korea Institute of Geoscience and Mineral Resources (KIGAM), Daejeon, Korea (Figure 1a). Of the 31 buildings at KIGAM, 3 have GHP systems installed for cooling and heating purposes. The GHP system installed in the Earthquake Research Center (ERC) building of KIGAM, is the oldest among the three, and has been in operation since November 2005 (Figure 1b). The ERC building features three stories above ground and one basement level, and 28 borehole heat exchangers (BHEs) are buried under the front yard located on the east side of the building. The dimensions of the BHE field are 35 m from east to west and 42 m from north to south. A total of 79 heat pumps were installed individually for labs and offices in the building. Four circulating pumps with closed circuits connected to BHEs and heat pumps are located in the mechanical room on the basement floor. Heat pumps installed on each floor are connected to each BHE group. Three groundwater monitoring wells were drilled at the center and in the southern part of the BHE field, and the temperature of the groundwater was measured at depths of 10 and 50 m in a 300 m deep monitoring well located in the south (monitoring well A). Figure 1c shows the layout of the GHP systems and its monitoring system for the ERC building. A detailed description of this site was presented in a previous study [8]. The geothermal gradient of 20 °C/km was measured at 0.5 m intervals using fiber optic cable sensors in monitoring well A. The basal heat flow measuring 59.7 mW/m2 was estimated from the mean thermal conductivity of 2.98 W/mK obtained by taking 61 core samples from monitoring well A and the geothermal gradient [16].

2.2. Monitoring Data

Approximately two and a half years after the operation of the GHP system, groundwater temperature monitoring began at depths of 10 and 50 m in monitoring well A at 1 h intervals (Figure 2). Groundwater temperatures were measured and recorded using two TD-Divers manufactured by Van Essen Instruments (Delft, Netherlands), which has specifications of temperature range of −20–+80 °C, accuracy of ±0.1 °C, and resolution of 0.01 °C. The monitoring period was divided into two parts: from April 30, 2008 to August 8, 2010 (first monitoring period) and from December 8, 2015 to May 28, 2018 (second monitoring period). Data from October 6 to December 7, 2019 were not recorded. Further monitoring was not possible after May 2018 because the monitoring well was damaged.
Although groundwater temperatures for each monitoring period did not clearly show a steadily increasing trend, comparisons between the changes in the mean annual temperature of groundwater over the two monitoring periods revealed an increasing trend in groundwater temperatures at the study site. In addition, temperature fluctuations in the form of a sinusoid function with a period of one year were observed at depths of 10 and 50 m. The temperature fluctuations at a depth of 10 m can be amplified or diminished owing to the seasonal variation of SATs. As the groundwater temperature at a depth of 50 m is hardly affected by the seasonal variation of SATs, temperature fluctuations occurring every year at this depth can be considered to be caused mainly by the operation of the GHP system. The annual mean rates of temperature increase calculated using linear regression analysis of all measured data are 0.135 °C/year and 0.118 °C/year at depths of 10 and 50 m, respectively. As no data were recorded in October and November 2009, and only data from the cooling season were recorded from May to August 2010, the annual mean rate of temperature increase may rise. There are three periods of continuous and uninterrupted temperature measurement during the year (May 2008 to April 2009, May 2016 to April 2017, and May 2017 to April 2018). The annual mean rates of temperature increase calculated using the temperatures of these periods were 0.132 °C/year (coefficient of determination, r2 = 0.996) and 0.119 °C/year (r2 = 0.991) at depths of 10 and 50 m, respectively.

2.3. Simulation Model and Physical Background

In this study, the following three steps of simulation were performed:
  • Simulation 1: One-dimensional vertical heat transfer simulations to analyze the effect of SAT change on temperature fluctuations in monitoring well A.
  • Simulation 2: Two-dimensional horizontal heat transfer simulations to evaluate the contribution of the GHP (i.e., amount of heat discharged into or extracted from the ground through the GHP system) to the increase in the ground temperature of monitoring well A without considering the SAT change effect.
  • Simulation 3: Three-dimensional heat transfer simulations taking into account future climate change of the ground temperature assuming the GHP system continues to operate as it does now.
Transient simulations of heat fluxes were performed using TOUGH3 (Lawrence Berkeley National Laboratory, Berkeley, CA, USA), which has been developed for multi-dimensional fluid and heat flows of multiphase, multicomponent fluid mixtures in porous and fractured media [17]. The fluid flow was not considered in the simulations because using the concept of effective thermal conductivity would yield a significantly different result only if the groundwater flow rate is very high. The general form of the basic energy balance equations in a porous medium is expressed as follows:
d d t V n M d V n = Γ n F n d Γ n + V n q d V n ,
where Vn is an arbitrary subdomain bounded by the closed surface Γn, and n is a normal vector on the surface element dΓn pointing inward into Vn. M denotes the energy per unit volume. F represents the heat flux, n is a normal vector on the surface dΓn, and q represents heat sources or sinks (J/s) [18]. The volume, Vn, should be large enough to be a “representative elementary volume” including many pores and mineral grains, to ensure the validity of the continuum approximation for the porous medium.
The heat accumulation term (M) is expressed as follows:
M = ρ R c R T ,
where ρR is the rock density, cR is the specific heat of the rock, and T is the temperature. Within each subdomain Vn, the fluid and rock are assumed to have the same temperature.
The conductive heat flux (F) is estimated as follows:
F = λ T ,
where λ is the thermal conductivity.
Figure 3 shows two- and three-dimensional model domains. A developed mesh generator applying an adaptive gridding technique, known as centroidal Voronoi tessellation (CVT), was used for the discretization of the model domains [19,20]. The connections between two adjacent grid blocks (a straight line connecting the center of each grid block) in a TOUGH3 grid should be orthogonal to their connection interface. The mesh produced using the CVT always satisfies the orthogonal condition of the TOUGH3 grid. The model domain was designed to be dense around boreholes where the heat source or sink existed and sparsely toward the side boundary. Mesh density can affect discretization errors and simulation results. In general, the denser the mesh, the less discretization error and the higher the computational amount. If the mesh size is not small enough in the transient heat transfer simulation, the temperature rises more slowly than it actually is when heat is injected. The mesh size used in this study was determined by comparing the increased temperature when heat was injected into the center of the well for 30 days. As the mesh size was gradually reduced (especially around the well), the temperature at the center of the well increased gradually when 30 days was reached. The mesh size was reduced until the difference in temperature value disappeared at the second decimal place.

3. Results

3.1. Simulation 1—The Effect of Climate Change on Temperature Fluctuations in Monitoring Well A

The ground temperature profile with respect to the depth can be determined using a basal heat flux, thermal conductivity of the ground, GST, and the inherent heat source. Assuming that the heat source within a depth of 50 m from the surface is negligible, and that the basal heat flux and thermal conductivity of the ground remain almost constant, changes in the ground temperature depend only on the GST. The GST around the study area (Figure 1b) has been measured by the Daejeon Regional Office of Meteorology of the Korea Meteorological Administration since 1969. The observation station is located 1.3 km away from the study area. The ground temperature data measured at depths of 1, 1.5, 3, and 5 m were downloaded from the website of the Korea Meteorological Administration (data.kma.go.kr). As the EOS3 module of the TOUGH3 exhibits a numerical problem in the calculation of physical parameters for sub-zero temperatures, the ground temperature at a depth of 1 m was used instead of the GST. The GST data beginning from 1969 are not extensive enough for the calculation of the present-day temperature change at a depth of 50 m. As the GST data before 1969 were not available in Daejeon, the SAT and GST data from other regions in Korea (Gangneung, Seoul, Daegu, Jeonju, Busan, and Mokpo) were used to estimate them (Figure 4). The linear regression slopes of the annual mean SAT and GST for each city are summarized in Table 1. The annual mean GST of Daejeon before 1969 was estimated using the following procedure:
  • The annual mean SAT values of other regions in Korea before 1969 were calculated.
  • The difference between the annual mean SAT of other regions in Korea and the annual mean SAT of Daejeon after 1969 was calculated.
  • The difference between the annual mean SAT and the annual mean GST of Daejeon after 1969 was calculated.
  • The results yielded by the three steps above were summed up. Then, GSTs were estimated from the SATs of Daejeon.
  • The annual mean GST of other regions in Korea before 1969 was calculated.
  • The difference between the annual mean GST of other regions in Korea and the annual mean GST of Daejeon after 1969 was calculated.
  • The results yielded by the two preceding steps were summed up. The GSTs of Daejeon were estimated from the GSTs of other regions before 1969.
Figure 5 shows the annual mean GSTs of Daejeon estimated from the SATs (procedures 1 to 4) and GSTs (procedures 5 to 7). The overall means were 13.88 and 14.17 °C, respectively.
Change in ground temperature was simulated using the estimated GSTs before 1969 and the measured GSTs in Daejeon after 1969. The GSTs were used as the time-dependent upper boundary condition, while the basal heat flux (59.7 mW/m2) measured at monitoring well A was used as the lower boundary condition. Assuming that there was no significant change in temperature before 1916, the depth-dependent initial conditions were determined according to the setting shown in Table 2. In the one-dimensional simulation, the meshmaker software of TOUGH3 was used to develop the model domain. The 200 m deep model domain consisting of 55 elements and 53 connections was divided into four zones according to thermal conductivities. The first zone (1 m to 6 m depth) was located above the groundwater table. The second zone (6 m to 14 m depth) was the alluvium where the casing of the well was installed. Third (14 m to 40 m depth) and fourth zones (40 m to 200 m depth) represent the fractured zone and bedrock, respectively. Although the thermal property data of each zone was measured using the core samples in several sections of each zone, the upper three zones are highly heterogeneous; thus, the thermal properties in those zones could not be used as representative values. Accordingly, simulations were performed to determine the thermal conductivities of three zones, which was consistent with temperature changes in the ground due to GST changes. The measured values of the specific heat and density of three zones were used in the simulation since their variations in the realm of nature were relatively less compared to the thermal conductivities (Table 2). 1134 simulations were performed for cases of combinations of GSTs and thermal conductivities as described in Table 2.
Figure 6 shows the simulation results of four cases for changes in temperature at a depth of 10 and 50 m from 2008 to 2018 due to changes in the GST. There was no significant difference in the trend of temperature changes for each case. Table 3 shows the mean thermal conductivity of each zone for each case and the rates of temperature increase obtained from linear regression slopes. The rates of temperature increase in four cases were 0.029 °C/year and 0.010 °C/year, respectively, at depths of 10 and 50 m. Therefore, the annual mean rates of temperature increase at depths of 10 and 50 m in monitoring well A only due to the influence of the GHP system are 0.103 °C/year and 0.109 °C/year for about 10 years, respectively.
Among all the tested simulation cases, only four cases of the simulation were in good agreement with the temperature measurements at a depth of 50 m in monitoring well A in November 2005 and ground temperatures measured at depths of 1.5, 3, and 5 m since 1969.

3.2. Simulation 2–The Effect of the GHP System on Temperature Fluctuations in Monitoring Well A

The BHE discharges heat into the ground during the cooling season and absorbs heat from the ground during the heating season. If there is an annual imbalance between the heat discharged into and absorbed from the ground over many years, the temperature of the ground around BHEs may continuously increase or decrease depending on the energy imbalance. In this study site, the increasing temperature of the BHE field could be attributed to the fact that the cooling loads of the ERC building were greater than the heating loads.
Two-dimensional heat transfer simulations were performed to determine the amount of heat discharged into or extracted from each BHE during the cooling or heating season. A constant amount of heat was discharged into or extracted from 28 geothermal wells in the two-dimensional model domain during the cooling or heating season over 12 years, assuming 120 days of cooling and 120 days of heating per year. The two-dimensional model domain of 200 m-by-200 m discretized by unstructured meshes consisted of 3621 elements and 10,004 connections (interfaces) between them. Changes in the GST, basal heat flux, and geothermal gradient were not taken into account in the two-dimensional model because the effect of temperature increase in the monitoring well due to the GHP system was the purpose of the simulation. The BHE is simulated as a cylindrical heat source or sink without taking into account the circulation of the fluid in the closed circuit. No flux boundary conditions were specified on the upper and lower boundaries, and fixed temperature (15 °C) boundary conditions were specified on the side boundaries. The initial temperature of the domain was also 15 °C. The thermal conductivity, specific heat capacity, and density of the ground were set to be 3.0 W/mK, 0.82 kJ/kgK, and 2670 kg/m3, respectively.
Through several simulations performed with different values of heat sources or sinks, the annual mean temperature was found to increase by 0.109 °C in the case of 28 heat sources of 16.7 W/m released into the ground through each BHE in the cooling season and 28 heat sinks of 12.4 W/m absorbed from the ground through each BHE in the heating season. Figure 7 shows temperature fluctuations at monitoring well A. The annual mean rate of temperature increase decreased gradually over time. Therefore, if the GHP system continues to operate under the same conditions and there is enough ground space for the heat to transfer, the ground temperature around the BHE field would not increase infinitely.

3.3. Simulation 3—Future Changes in Ground Temperature

If the operation of the GHP system is maintained as is, the change in the temperature distribution around the BHE field continues to be a concern. Based on the previous simulation results, three-dimensional heat transfer simulations were performed to predict changes in the temperature of the surrounding ground via the operation of the GHP system until 2050. Unlike in the two-dimensional model, variations of the GST, the geothermal gradient due to the basal heat flux, and change in the ground thermal conductivity with respect to depth were considered in the three-dimensional model. Model calibration was performed by comparing simulated temperature variations with those produced in actual temperature measurements. The same boundary conditions and thermal properties acquired from Simulation 1 were used in the three-dimensional model except for zone 4, and thermal properties were further discretized depending on the depth. The boundary conditions, initial conditions, and thermal properties used in the simulations for the calibration are summarized in Table 4. Three-dimensional simulations of the GHP system were performed using the same amount of heat discharge and extraction calculated from the two-dimensional simulation.
Figure 8 shows temperature fluctuations at the depth of 50 m in the groundwater monitoring well. The increasing trends of the temperature yielded by the simulation results and measured data were similar, although it was difficult to fit measured data for each cooling or heating season well with the simulation results because the actual operating pattern of the GHP system varied every year. Similar to the two-dimensional simulation, the increasing rates of the mean annual temperature tend to decrease over time. To identify the case that best matches the measured temperature data, the linear regression analysis of temperature fluctuations at a depth of 50 m in monitoring well A of each case was performed. The slope of the linear regression of all cases including measured data was almost identical, and the difference between the y-intercept of the linear regression of measured data and each case was 0.02, 0.09, 0.18, and 0.24 °C, respectively. Therefore, the linear regression of case 1 coincided best with the linear regression of measured data (Figure 8). Therefore, simulations of long-term operation of the GHP system from 2018 to 2050 were performed using the result of case 1 as the initial condition. Boundary conditions, except the upper boundary, and thermal properties used in simulations for future temperature prediction were the same as those used in case 1. Possible SATs obtained from climate change scenarios [21] (RCP (Representative Concentration Pathway) 2.6, RCP 4.5, RCP 6.0, and RCP 8.5) were converted to GSTs and used as time-dependent upper boundary conditions (Figure 9a). The conversion from SATs to GSTs was accomplished using the same method as that used in the Section 3.1.
Figure 10 shows temperature distributions around the BHE field at the depth of 50 m under the RCP 8.5 (the worst-case climate change scenario). The ground temperature around the BHE field increases over time. Figure 11 shows temperature changes at a depth of 50 m in groundwater monitoring well A for each RCP scenario. The temperature in RCP 8.5 recorded the greatest increase in temperature, while in RCP 6.0, the temperature exhibited the least increase at a depth of 50 m in groundwater monitoring well A. The least increase in temperature under RCP 6.0 could be because the predicted annual mean temperature in Daejeon from 2020 to 2050 was lowest (15.84 °C) in RCP 6.0 and highest (16.93 °C) under RCP 8.5. In May 2050, the difference between temperatures at the depth of 50 m in the groundwater monitoring well A in under RCP 8.5 and 6.0 will be approximately 0.17 °C. The highest temperatures at the depth of 50 m in the groundwater monitoring well A calculated for each RCP scenario will be 19.10 °C for RCP 2.6, 19.06 °C for RCP 4.5, 18.96 °C for RCP 6.0, and 19.12 °C for RCP 8.5 in November 2049. These are the mean increase of 1.22 °C and 3.00 °C from the maximum measured temperature of 17.84 °C at the depth of 50 m in the groundwater monitoring well A on 9 November 2017 and the estimated temperature of 16.06 °C at the depth of 50 m in the groundwater monitoring well A immediately before the operation of the GHP system.

4. Discussion

Due to the nature of the heat pump that converts the power consumed by the compressor into thermal energy, even if the cooling and heating loads of the building are the same, the amount actually discharged into or extracted from the ground varies [22]. For example, if a GHP system with a coefficient of performance, known as COP, value of 3.0 for both cooling and heating is used to emit heat amounting to 100 W from the room in the cooling season and bring heat of 100 W in the heating season, a heat of 133 W is discharged into the ground during the cooling season and heat measuring 67 W is extracted from the ground during the heating season. Therefore, heat may continuously accumulate in the ground even in a building with the same cooling load as the heating load. Moreover, since the global warming will increase both the annual mean GST and cooling load of the building in the future, the ground temperatures around the GHP system are likely to increase continuously. In relatively cold regions, there could be continuous decreases in the ground temperatures around the GHP system. However, the effect of the global warming on the GST and cooling and heating loads of the building will be able to change decreases in the ground temperature to increases in the ground temperature in the future. In relatively warm regions where the GHP systems are used in an environment where heating is less important than cooling, heat is highly likely to accumulate in the ground. This may cause a problem that may easily be overlooked because this phenomenon does not have to be considered in conventional air-source heat pump systems. If the temperature in the ground increases in this way, it will change not only the performance of the heat pump, but also the physical, chemical [23], and microbiological environments in the ground [24,25].
Heat transfer characteristics and geothermal heat exchanger configurations can also influence the behavior of the ground temperature change due to the long-term operation of the GHP system. If the same heat as the study area is discharged or extracted through BHE in areas with higher effective thermal conductivity than in the study area, the increase or decrease in the ground temperature occurs less than in the study area. Therefore, in the case of the GHP system installed in such an area, not only the performance of each season is better, but also the temperature changes in the ground due to thermal imbalance are less during the long-term operation. On the contrary to this, the opposite occurs in areas with lower effective thermal conductivity than in the study area. In the case of GHP systems with horizontal heat exchangers, heat accumulation by the GHP system may be less than that of a vertical system because heat can easily transfer from the ground surface to the atmosphere. However, their cooling efficiency may decrease with increasing the GST during long-term operation because they are much more affected by global warming than vertical systems.
Temperature data measured at a depth of 50 m were mainly used for these analyses because those obtained from a depth of 10 m were affected by thermal conductivities of alluvial layers and weathered rocks with high heterogeneity and thermal disturbance by rainfall, which are not considered in the model used in this study. However, it would be better to consider the groundwater flow to further analyze the thermal behavior at depths that were shallower than 10 m [26]. This is because the hydraulic conductivity of alluvial aquifers at that depth is likely to be relatively higher than that of bedrock aquifers. In addition, at that depth, changes in the ground temperature can result from advection due to the groundwater flow in the vertical direction, rapid thermal property changes due to the movement of the groundwater table, and so on [27].
The operation pattern of the GHP system and thermal process occurring inside the BHE were simplified as constant heat discharge or extraction and the line source, respectively, which was an appropriate assumption in the determination of the cause of the temperature increase in the monitoring well. However, to analyze the effect of the increase in ground temperature on the performance of the GHP system, an operation pattern that is similar to the actual operation pattern, the thermal process occurring inside the BHE, and the groundwater flow should be considered.

5. Conclusions

In this study, the causes of temperature increase were analyzed and evaluated by conducting simulations using measured temperatures in the groundwater monitoring well around the GHP system operated for nearly 15 years. In addition, the prediction of the future temperature distribution of the BHE field was quantitatively estimated through the simulation model. The annual mean rate of temperature increase calculated using continuously measured temperatures of three one-year periods without interruption (May 2008 to April 2009, May 2016 to April 2017, and May 2017 to April 2018) were 0.132 °C/year and 0.119 °C/year at depths of 10 and 50 m, respectively. In total, 1134 simulations of the one-dimensional vertical heat transfer were performed to analyze the effect of climate change on measured temperature fluctuations. Only four cases of the simulation agreed well with temperature measurements at a depth of 50 m in monitoring well A in November 2005 and ground temperatures measured at depths of 1.5, 3, and 5 m since 1969. The simulation results showed that temperature increments due to changes in SATs were 0.029 °C/year and 0.010 °C/year at 10 m and 50 m depths, respectively. The remaining temperature increments were probably due to the operation of the GHP system.
Through several two-dimensional horizontal heat transfer simulations to evaluate the amount of heat discharged into or extracted from the ground through the GHP system, the discharged heat measuring 16.7 W/m in the cooling season and extracted heat measuring 12.4 W/m in the heating season were found to cause an increase of 0.109 °C in the annual mean temperature at a depth of 50 m in monitoring well A, where the thermal conductivity of the ground was 3.0 W/mK. Based on the previous simulation results, three-dimensional heat transfer simulations were performed to calculate changes in the temperature of the surrounding ground through the operation of the GHP system over next 30 years. If the operating conditions of the GHP system are retained in the future, results showed that in 2050, the temperature at a depth of 50 m in monitoring well A would reach a maximum of 19.10, 19.06, 18.96, and 19.12 °C based on RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5 scenarios, respectively.
As temperatures are predicted to increase because of climate change, the demand for cooling is very likely to increase gradually. The use of GHP systems, which is considered eco-friendly, has the potential to adversely affect the environment because of the acceleration of the temperature increase in the ground. Therefore, during the design of the GHP system, the heat discharged into and extracted from the ground should be balanced to minimize long-term changes in the ground temperature.

Author Contributions

Conceptualization, S.-K.K.; formal analysis, S.-K.K. and Y.L.; investigation, S.-K.K. and Y.L.; writing—original draft preparation, S.-K.K. and Y.L.; writing—review and editing, S.-K.K. and Y.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Basic Research Project (20-3411) of the Korea Institute of Geoscience and Mineral Resources (KIGAM) funded by the Ministry of Science, ICT, and Future Planning of Korea.

Acknowledgments

We thank the reviewers of Water for their insightful comments, which have helped us to improve the manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Banks, D. An Introduction to Thermogeology: Ground Source Heating and Cooling, 2nd ed. ; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2012; ISBN 978-0-470-67034-7. [Google Scholar]
  2. Self, S.J.; Reddy, B.V.; Rosen, M.A. Geothermal heat pump systems: Status review and comparison with other heating options. Appl. Energy 2013, 101, 341–348. [Google Scholar] [CrossRef]
  3. Rybach, L. Geothermal energy: Sustainability and the environment. Geothermics 2003, 32, 463–470. [Google Scholar] [CrossRef]
  4. Lee, W.-N.; Kim, H.-J.; Park, J.-B.; Cho, K.-S.; Roh, J.H.; Son, S.-Y. Economic analysis of heating and cooling systems from the various perspectives: Application to EHP and GHP in Korea. Renew. Sust. Energ. Rev. 2012, 16, 4116–4125. [Google Scholar] [CrossRef]
  5. Chang, Y.; Gu, Y.; Zhang, L.; Wu, C.; Liang, L. Energy and environmental implications of using geothermal heat pumps in buildings: An example from north China. J. Clean. Prod. 2017, 167, 484–492. [Google Scholar] [CrossRef]
  6. Marrasso, E.; Roselli, C.; Sasso, M.; Tariello, F. Global and local environmental and energy advantages of a geothermal heat pump interacting with a low temperature thermal micro grid. Energ. Convers. Manag. 2018, 172, 540–553. [Google Scholar] [CrossRef]
  7. Deng, J.; Wei, Q.; Liang, M.; He, S.; Zhang, H. Field test on energy performance of medium-depth geothermal heat pump systems (MD-GHPs). Energy Build. 2019, 184, 289–299. [Google Scholar] [CrossRef]
  8. Kim, S.-K.; Bae, G.-O.; Lee, K.-K.; Song, Y. Field-scale evaluation of the design of borehole heat exchangers for the use of shallow geothermal energy. Energy 2010, 35, 491–500. [Google Scholar] [CrossRef]
  9. Fan, R.; Jiang, Y.; Yao, Y.; Shiming, D.; Ma, Z. A study on the performance of a geothermal heat exchanger under coupled heat conduction and groundwater advection. Energy 2007, 32, 2199–2209. [Google Scholar] [CrossRef]
  10. Dehkordi, S.E.; Schincariol, R.A. Effect of thermal-hydrogeological and borehole heat exchanger properties on performance and impact of vertical closed-loop geothermal heat pump systems. Hydrogeol. J. 2014, 22, 189–203. [Google Scholar] [CrossRef]
  11. Kurevija, T.; Vulin, D.; Krapec, V. Effect of borehole array geometry and thermal interferences on geothermal heat pump system. Energ. Convers. Manag. 2012, 60, 134–142. [Google Scholar] [CrossRef]
  12. Han, C.; Yu, X.B. Sensitivity analysis of a vertical geothermal heat pump system. Appl. Energy 2016, 170, 148–160. [Google Scholar] [CrossRef] [Green Version]
  13. Lachenbruch, A.H.; Marshall, B.V. Changing Climate: Geothermal Evidence from Permafrost in the Alaskan Arctic. Science 1986, 234, 689–696. [Google Scholar] [CrossRef] [PubMed]
  14. Deming, D. Climatic warning in North America: Analysis of borehole temperatures. Science 1995, 268, 1576–1577. [Google Scholar] [CrossRef]
  15. Pollack, H.N. Climate Change Record in Subsurface Temperatures: A Global Perspective. Science 1998, 282, 279–281. [Google Scholar] [CrossRef] [Green Version]
  16. Shim, B.-O.; Lee, Y.; Kim, H.-C.; Song, Y. Investigation of Thermal and Hydraulic Characteristics for the Performance Analysis of a Borehole Heat Exchanger. J. Korean Soc. Miner. Energy Resour. 2006, 43, 97–105. [Google Scholar]
  17. Jung, Y.; Pau, G.S.H.; Finsterle, S.; Doughty, C.A. TOUGH3 User’s Guide, Version 1.0.; Lawrence Berkeley National Lab. (LBNL): Berkeley, CA, USA, 2018. [Google Scholar]
  18. Pruess, K.; Oldenburg, C.M.; Moridis, G. TOUGH2 User’s Guide Version 2; Lawrence Berkeley National Lab. (LBNL): Berkeley, CA, USA, 1999; pp. 1–197. [Google Scholar]
  19. Kim, S.-K.; Bae, G.-O.; Lee, K.-K. Improving accuracy and flexibility of numerical simulation of geothermal heat pump systems using Voronoi grid refinement approach. Geosci. J. 2015, 19, 527–535. [Google Scholar] [CrossRef]
  20. Kim, S.-K. TOUGH2 Mesh Generation Based on Constrained Centroidal Voronoi Tessellations for Simulation of Geothermal Heat Pump Systems. J. Korean Soc. Miner. Energy Resour. 2019, 56, 639–644. [Google Scholar] [CrossRef]
  21. Rozenberg, J.; Fay, M.; Riahi, K.; van Vuuren, D.; Kriegler, E.; Edmonds, J.; O’Neill, B.; Fujimori, S.; Bauer, N.; van Vuuren, D. The Representative Concentration Pathways: An Overview. Clim. Chang. 2011, 109, 5. [Google Scholar]
  22. Zhao, Z.; Shen, R.; Feng, W.; Zhang, Y.; Zhang, Y. Soil Thermal Balance Analysis for a Ground Source Heat Pump System in a Hot-Summer and Cold-Winter Region. Energies 2018, 11, 1206. [Google Scholar] [CrossRef] [Green Version]
  23. Drake, J.J.; Wigley, T.M.L. The effect of climate on the chemistry of carbonate groundwater. Water Resour. Res. 1975, 11, 958–962. [Google Scholar] [CrossRef]
  24. Yates, M.V.; Gerba, C.P.; Kelley, L.M. Virus persistence in groundwater. Appl. Environ. Microbiol. 1985, 49, 778–781. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. John, D.E.; Rose, J.B. Review of Factors Affecting Microbial Survival in Groundwater. Environ. Sci. Technol. 2005, 39, 7345–7356. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Perego, R.; Guandalini, R.; Fumagalli, L.; Aghib, F.S.; De Biase, L.; Bonomi, T. Sustainability evaluation of a medium scale GSHP system in a layered alluvial setting using 3D modeling suite. Geothermics 2016, 59, 14–26. [Google Scholar] [CrossRef]
  27. Anderson, M.P. Heat as a ground water tracer. Ground Water 2005, 43, 951–968. [Google Scholar] [CrossRef]
Figure 1. Location of the study area (a), a satellite view of the Earthquake Research Center (ERC) building and borehole heat exchanger (BHE) field (b) and the layout of the geothermal heat pump (GHP) system in the ERC building (c). The vertical closed-loop GHP system comprises 79 heat pumps, 4 fluid pumps, and 28 BHEs. The dimensions of the BHE field are 35 m from east to west and 42 m from north to south. Red, blue, and green boxes in the ERC building represent heat pumps for each floor.
Figure 1. Location of the study area (a), a satellite view of the Earthquake Research Center (ERC) building and borehole heat exchanger (BHE) field (b) and the layout of the geothermal heat pump (GHP) system in the ERC building (c). The vertical closed-loop GHP system comprises 79 heat pumps, 4 fluid pumps, and 28 BHEs. The dimensions of the BHE field are 35 m from east to west and 42 m from north to south. Red, blue, and green boxes in the ERC building represent heat pumps for each floor.
Water 12 02931 g001aWater 12 02931 g001b
Figure 2. Groundwater temperatures measured at 10 and 50 m depths in the 300 m deep monitoring well for approximately 10 years. The blue and red dots represent the temperature at depths of 10 and 50 m, respectively. The green triangular and yellow circular marks are the annual mean temperature at depths of 10 and 50 m, respectively. The green and yellow dotted lines are the results of the linear regression analysis of the annual mean temperature at depths of 10 and 50 m, respectively.
Figure 2. Groundwater temperatures measured at 10 and 50 m depths in the 300 m deep monitoring well for approximately 10 years. The blue and red dots represent the temperature at depths of 10 and 50 m, respectively. The green triangular and yellow circular marks are the annual mean temperature at depths of 10 and 50 m, respectively. The green and yellow dotted lines are the results of the linear regression analysis of the annual mean temperature at depths of 10 and 50 m, respectively.
Water 12 02931 g002
Figure 3. Two- and three-dimensional model domains. The two-dimensional model domain (a) is the same as a plan view of the three-dimensional model domain (b). The domain sizes of two- and three-dimensional models are 200 × 200 m (length × width) and 200 × 200 × 300 m (length × width × depth), respectively. There are 3621 and 29,150 elements and 10,004 and 108,894 connections in the two- and three-dimensional model domains, respectively.
Figure 3. Two- and three-dimensional model domains. The two-dimensional model domain (a) is the same as a plan view of the three-dimensional model domain (b). The domain sizes of two- and three-dimensional models are 200 × 200 m (length × width) and 200 × 200 × 300 m (length × width × depth), respectively. There are 3621 and 29,150 elements and 10,004 and 108,894 connections in the two- and three-dimensional model domains, respectively.
Water 12 02931 g003aWater 12 02931 g003b
Figure 4. The annual mean (a) surface air temperatures (SAT) and (b) GST in Gangneung, Seoul, Daegu, Jeonju, Busan, Mokpo, and Daejeon observed by the Korea Meteorological Administration. The SATs in Gangneung, Seoul, Daegu, Jeonju, Busan, Mokpo, and Daejeon have been measured since 1912, 1907, 1909, 1919, 1904, 1906, and 1969, respectively. The GSTs in Gangneung, Seoul, Daegu, Jeonju, Busan, Mokpo, and Daejeon have been measured since 1917, 1916, 1918, 1921, 1917, 1917, and 1969, respectively.
Figure 4. The annual mean (a) surface air temperatures (SAT) and (b) GST in Gangneung, Seoul, Daegu, Jeonju, Busan, Mokpo, and Daejeon observed by the Korea Meteorological Administration. The SATs in Gangneung, Seoul, Daegu, Jeonju, Busan, Mokpo, and Daejeon have been measured since 1912, 1907, 1909, 1919, 1904, 1906, and 1969, respectively. The GSTs in Gangneung, Seoul, Daegu, Jeonju, Busan, Mokpo, and Daejeon have been measured since 1917, 1916, 1918, 1921, 1917, 1917, and 1969, respectively.
Water 12 02931 g004
Figure 5. The annual mean GSTs of Daejeon used as the time-dependent upper boundary condition of the simulation model. Before 1969, the GSTs of Daejeon were estimated from the SATs of Daejeon (blue line) and the GSTs of other regions (orange line). They have been measured since 1969 (grey line).
Figure 5. The annual mean GSTs of Daejeon used as the time-dependent upper boundary condition of the simulation model. Before 1969, the GSTs of Daejeon were estimated from the SATs of Daejeon (blue line) and the GSTs of other regions (orange line). They have been measured since 1969 (grey line).
Water 12 02931 g005
Figure 6. The simulation results of four cases considering temperature changes at depths of 10 and 50 m due to changes in the GST from 2008 to 2018.
Figure 6. The simulation results of four cases considering temperature changes at depths of 10 and 50 m due to changes in the GST from 2008 to 2018.
Water 12 02931 g006
Figure 7. Temperature fluctuations at monitoring well A obtained from the result of the two-dimensional heat transfer simulation. The temperature of the ground around the BHE field increases gradually because the heat released to the ground through the BHE in the cooling season (16.7 W/m) is greater than the heat absorbed from the ground through the BHE in the heating season (12.4 W/m). The black dotted line is the result of linear regression, and the slope is 0.109 °C/year.
Figure 7. Temperature fluctuations at monitoring well A obtained from the result of the two-dimensional heat transfer simulation. The temperature of the ground around the BHE field increases gradually because the heat released to the ground through the BHE in the cooling season (16.7 W/m) is greater than the heat absorbed from the ground through the BHE in the heating season (12.4 W/m). The black dotted line is the result of linear regression, and the slope is 0.109 °C/year.
Water 12 02931 g007
Figure 8. Comparisons of simulated and measured temperatures at the depth of 50 m in the monitoring well for 12 years of operation.
Figure 8. Comparisons of simulated and measured temperatures at the depth of 50 m in the monitoring well for 12 years of operation.
Water 12 02931 g008
Figure 9. Boundary conditions of the simulations for future changes in ground temperature: (a) the annual mean GSTs of Daejeon used as the time-dependent upper boundary condition of the simulations for future changes in ground temperature. Possible SATs obtained from each RCP scenario were converted to GSTs. (b) the initial temperature distribution for lower and sides boundary conditions.
Figure 9. Boundary conditions of the simulations for future changes in ground temperature: (a) the annual mean GSTs of Daejeon used as the time-dependent upper boundary condition of the simulations for future changes in ground temperature. Possible SATs obtained from each RCP scenario were converted to GSTs. (b) the initial temperature distribution for lower and sides boundary conditions.
Water 12 02931 g009aWater 12 02931 g009b
Figure 10. Temperature distribution around the BHE field at the depth of 50 m in 2018, 2030, 2040, and 2050 in order from the top: (a) after the cooling season and (b) after the heating season.
Figure 10. Temperature distribution around the BHE field at the depth of 50 m in 2018, 2030, 2040, and 2050 in order from the top: (a) after the cooling season and (b) after the heating season.
Water 12 02931 g010aWater 12 02931 g010b
Figure 11. Temperature changes at the depth of 50 m in the groundwater monitoring well due to the operation of the GHP system under possible GSTs obtained from RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5 as the time-dependent upper boundary conditions.
Figure 11. Temperature changes at the depth of 50 m in the groundwater monitoring well due to the operation of the GHP system under possible GSTs obtained from RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5 as the time-dependent upper boundary conditions.
Water 12 02931 g011
Table 1. Linear regression slopes of the annual mean SAT and GST (°C/year).
Table 1. Linear regression slopes of the annual mean SAT and GST (°C/year).
Gang-NeungSeoulDaeguJeonjuBusanMokpoDaejeon
Before 1969 (SAT)0.01100.01830.01270.02550.01180.0137-
1969–2019 (SAT)0.03150.03470.03930.03020.02980.01520.0393
Before 1969 (GST)0.02530.01290.03480.00250.01450.0089-
1969–2019 (GST)0.02090.03770.03170.03620.02650.02990.0456
Table 2. Summary of GST and thermal property data used for simulation.
Table 2. Summary of GST and thermal property data used for simulation.
GST(°C)from SATs13.88 (before 1916)
estimated GST (1916–1969)
measured GST (1969–2018)
from GSTs14.17 (before 1916)
estimated GST (1916–1969)
measured GST (1969-–2018)
Thermal Conductivity (W/mK)1–6 m1.00, 1.25, 1.50, 1.75, 2.00, 2.25, 2.50
6–14 m1.00, 1.25, 1.50, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00
14–40 m1.00, 1.25, 1.50, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00
40–200 m3.00
Specific Heatand Density1–6 m0.90 kJ/kgK, 2000 kg/m3
6–14 m0.82 kJ/kgK, 2670 kg/m3
14–40 m0.82 kJ/kgK, 2670 kg/m3
40–200 m0.85 kJ/kgK, 2670 kg/m3
Table 3. Summary of thermal conductivity data and rates of temperature increase (°C/year) for each case for approximately 10 years.
Table 3. Summary of thermal conductivity data and rates of temperature increase (°C/year) for each case for approximately 10 years.
Upper Boundary ConditionThermal Conductivity (W/mK)Rates of Temperature Increase (°C/year)
GSTZone 1Zone 2Zone 310 m50 m
Case 1From GSTs2.251.001.500.02760.0102
Case 2From GSTs2.251.251.500.02930.0110
Case 3From SATs2.251.001.250.02800.0099
Case 4From SATs2.251.251.250.02960.0106
Mean----0.0290.010
Table 4. Summary of boundary and initial conditions, and thermal properties used for simulation.
Table 4. Summary of boundary and initial conditions, and thermal properties used for simulation.
Boundary ConditionsUpper
(Time-Dependent)
before 1916
1916–1969
1969–2018
13.88 °C (SATs), 14.17 °C (GSTs)
estimated GST
measured GST
Lower and Sidesconstant temperature (Figure 9b)
(temperature profile of initial conditions)
Initial Conditionstemperature profile calculated from the basal heat flux (59.7 mW/m2), the initial upper boundary temperature, and the thermal conductivity of each case
Sources and Sinks at the BHEsCooling season: 28 heat sources of 16.7 W/m released through each BHE
Heating Season: 28 heat sinks of 12.4 W/m absorbed through each BHE
Thermal PropertiesDepthThermal Conductivity (W/mK)Specific Heat (kj/kgK)Density (kg/m3)
1–6 m2.250.902000
6–14 m1.00 (case 1, 2)
1.25(case 3, 4)
0.822670
14–40 m1.25 (case 3, 4)
1.50 (case 3, 4)
0.822670
40–60 m3.000.842660
60–100 m3.040.772680
100–200 m2.980.882670
200–300 m3.100.842680
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Kim, S.-K.; Lee, Y. Evaluation of Ground Temperature Changes by the Operation of the Geothermal Heat Pump System and Climate Change in Korea. Water 2020, 12, 2931. https://doi.org/10.3390/w12102931

AMA Style

Kim S-K, Lee Y. Evaluation of Ground Temperature Changes by the Operation of the Geothermal Heat Pump System and Climate Change in Korea. Water. 2020; 12(10):2931. https://doi.org/10.3390/w12102931

Chicago/Turabian Style

Kim, Seong-Kyun, and Youngmin Lee. 2020. "Evaluation of Ground Temperature Changes by the Operation of the Geothermal Heat Pump System and Climate Change in Korea" Water 12, no. 10: 2931. https://doi.org/10.3390/w12102931

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop