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

: 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 e ﬀ ects on measured temperature ﬂuctuations 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.


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.

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/m 2 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].

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, r 2 = 0.996) and 0.119 • C/year (r 2 = 0.991) at depths of 10 and 50 m, respectively.

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: where V n 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 V n . 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, V n , 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: where ρ R is the rock density, c R is the specific heat of the rock, and T is the temperature. Within each subdomain V n , the fluid and rock are assumed to have the same temperature. The conductive heat flux (F) is estimated as follows: 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.

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: 1.
The annual mean SAT values of other regions in Korea before 1969 were calculated. 2.
The difference between the annual mean SAT of other regions in Korea and the annual mean SAT of Daejeon after 1969 was calculated. 3.
The difference between the annual mean SAT and the annual mean GST of Daejeon after 1969 was calculated.

4.
The results yielded by the three steps above were summed up. Then, GSTs were estimated from the SATs of Daejeon. 5.
The annual mean GST of other regions in Korea before 1969 was calculated. 6.
The difference between the annual mean GST of other regions in Korea and the annual mean GST of Daejeon after 1969 was calculated. 7.
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.  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 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/m 2 ) 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.

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/m 3 , 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.

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.

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.

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.