Development of Simulation Tool for Ground Source Heat Pump Systems Influenced by Ground Surface

The authors developed a ground heat exchanger (GHE) calculation model influenced by the ground surface by applying the superposition theorem. Furthermore, a simulation tool for ground source heat pump (GSHP) systems affected by ground surface was developed by combining the GHE calculation model with the simulation tool for GSHP systems that the authors previously developed. In this paper, the outlines of GHE calculation model is explained. Next, in order to validate the calculation precision of the tool, a thermal response test (TRT) was carried out using a borehole GHE with a length of 30 m and the outlet temperature of the GHE calculated using the tool was compared to the measured one. The relative error between the temperatures of the heat carrier fluid in the GHE obtained by measurement and calculation was 3.3% and this result indicated that the tool can reproduce the measurement with acceptable precision. In addition, the authors assumed that the GSHP system was installed in residential houses and predicted the performances of GSHP systems using the GHEs with different lengths and numbers, but the same total length. The result showed that the average surface temperature of GHE with a length of 10 m becomes approximately 2 ◦C higher than the average surface temperature of a GHE with a length of 100 m in August.


Introduction
In the recent years, global warming due to increase of CO 2 emission has become a worldwide environmental issue. Therefore, ground source heat pump (GSHP) systems, which use the underground thermal energy as a heat sink/source and supply cooling, heating, and/or hot water with high efficiency, have attracted widespread attention. Several million units have been installed in the world [1]. However, the number of GSHP system installed in Japan is still low although it is gradually increasing. The main reason is that the installation cost of ground heat exchangers (GHEs) is expensive due to the complex underground formations. The utilization of shallow underground thermal energy such as short GHEs, energy piles and horizontal GHEs is effective to reduce the installation cost especially for the small building.
Recently, several research reports related to the use of short GHEs have been provided. A short vertical ground heat exchanger with a borehole length of 3 m and a helical shaped pipe was evaluated by means of a network of interconnected thermal resistance and capacitance [2]. Zarrella et al. analyzed short helical and double U-tube borehole heat exchangers and concluded that the helical heat exchanger had better performance than the double U-tube [3]. The operation modes of GSHP systems with short helical GHEs were also investigated [4]. Furthermore, there are several research reports on the use of piles. Hamada et al. evaluated the performance of energy piles installed in a building [5]. The authors have reported heating tests of GSHP systems using steel pipe foundation piles as GHEs in cold regions [6] and the effectiveness of using foundation piles in residential houses in moderate climate regions using simulations [6]. Furthermore, a simulation of a GSHP system using energy piles with short lengths and a large diameter was introduced [7]. With regard to research on horizontal GHEs, Fujii et al. developed a numerical model for slinky-coil horizontal GHEs and evaluated the performance [8,9]. Li et al. developed a simulation model for horizontal slinky GHEs by applying the analytical solution of ring source [10,11]. Xiong et al. developed a response function (g-function) for horizontal slinky GHEs [12].
Regarding the utilization of shallow ground thermal energy, the performance of GHEs may decrease due to the influence of the ambient air temperature and solar radiation. In a heating test using steel pipe foundation piles as a GHE actually performed in the cold district mentioned above [6], no building was constructed above the GHE, and the ground surface was in contact directly with the ambient air. As a result, the decrease of ground temperature become larger than the decrease calculated using the simulation. From this, it is necessary to predict the performance of GHEs affected by the ambient air temperature and the solar radiation on the ground surface. The ground temperature considering the influence of the ambient air temperature and solar radiation can be obtained by applying a numerical analysis [13] or an analytical solution of semi-infinite solid, surface temperature periodic with time [14]. For the latter (applying the analytical solution of semi-infinite solid), it is possible to superpose it with the ground temperature variation caused by the heat injection/extraction via GHEs, which is calculated by applying the analytical solution of line source or cylindrical heat source. Bandos et al. calculated the ground temperature variation in a thermal response test considering the influence of ground surface temperature by superposing the analytical solutions of finite line source and semi-infinite solid, surface temperature periodic with time [15].
In this study, the authors applied the superposition theorem and considered the influence of the ground surface on the GHE calculation model, which can be applied for a short GHE and a horizontal GHE. Then, a simulation tool for a GSHP system influenced by the ground surface was developed by combining the GHE calculation model with the simulation tool that the authors developed [6,7,16,17]. The outlines of the calculation model for a short GHE is firstly explained (the calculation model for a horizontal GHE will be introduced in the future). Next, a thermal response test (TRT) was carried out using a borehole GHE with a length of 30 m, and the outlet temperature of the GHE was calculated using the developed tool and compared with the measured one. In addition, the authors assumed that a GSHP systems using GHEs of different numbers and different lengths, but the same total length, was installed in residential houses and evaluated the performance using a simulation tool.

Overview of Ground Temperature Calculation
Soil is treated as an infinite isotonic constant solid, and the ground surface is considered as a flat solid surface without water evaporation and condensation. Then, the heat transfer in the ground surrounding the GHE can be treated as an axisymmetric three-dimensional unsteady heat transfer problem in a cylindrical coordinate system, and it is expressed as Equation (1).
Here, the initial condition is set as T s (r, z, 0) = T s0 , and the boundary condition is set as T s (r, 0, t) = ψ(t c ). By applying the superposition theorem, the temperature variation in the ground Energies 2020, 13, 4491 3 of 18 surrounding the GHE caused by the influence of the ground surface and the heat injection/extraction via the GHEs can be calculated as in Equation (2) [15].
Here, the gradient of ground temperature is not considered because the authors targeting GHEs up to about 100 m in depth. However, if the deep boreholes with several hundred meters are targeted, it is required to consider the gradient of ground temperature [18]. This challenge will be addressed in the future. Where, the first term on the right side is the initial temperature. The second term is the temperature variation caused by the heat injection/extraction via the GHEs, and the initial condition is set as ∆T s1 (r, z, 0) = 0 and the boundary condition is set as ∆T s1 (0, t) = 0. If the heat injection rate via the GHE is constant, ∆T s1 (r, z, t) can be calculated by using the analytical solution of finite line source. The equation can be expressed as the following [19]. Here, The third term is the temperature variation affected by the influence of the ground surface, and the initial condition is set as ∆T s2 (z, 0) = 0 and the boundary condition is set as Here, ψ(t c ) is the equivalent ambient air temperature and is calculated using Equation (4) [20].
Here, the first term on the right side is the ambient temperature, the second term and the third term are temperature variations affected by the solar radiation and the effective radiation, respectively. The effective radiation is obtained as the difference between the downward atmospheric radiation to the ground surface and the upward ground surface radiation to the atmosphere. In order to calculate the underground temperature distribution, the equivalent ambient air temperature is approximated to a periodic function by harmonic analysis of monthly temperature and solar radiation at an arbitrary point [20]. The monthly temperature and solar radiation can be obtained from the extended AMeDAS weather data [21]. Then, the underground temperature considering the influence of ground surface T s2 (r, z, t c ) can be calculated using Equation (5) [20].
where T T is the ground temperature affected by the ambient air temperature variation, T I is the temperature variation affected by the solar radiation and T J is the temperature affected by the effective radiation. For example, according to the extended AMeDAS weather data of Yahata (Kitakyushu), the outdoor temperature and solar radiation values for a standard year can be obtained, therefore, T T and T I can be approximated to Equations (6) and (7).
The calculated results are compared with the measured data. The 20 m borehole used in this experiment is installed on the campus of The University of Kitakyushu (Kitakyushu city, Japan). The schematic diagram and comparison results are shown in Figure 1. Here, T J is given as −3.5 • C [19]. Figure 1 compares the temperature distributions in different months (February, May and August) and at different depths. The calculation result shows the acceptable precision. Therefore, it can be determined that T J = −3.5 • C is used for simulation. The calculated results are compared with the measured data. The 20 m borehole used in this experiment is installed on the campus of The University of Kitakyushu (Kitakyushu city, Japan). The schematic diagram and comparison results are shown in Figure 1. Here, is given as −3.5 °C [19]. Figure 1 compares the temperature distributions in different months (February, May and August) and at different depths. The calculation result shows the acceptable precision. Therefore, it can be determined that = −3.5 °C is used for simulation.

Overview of Ground Temperature Calculation in Simulation Tool
The average surface temperature of the GHE, , , , can be obtained using Equation (8).
The second term on the right side is a temperature variation affected by the heat injection/extraction via GHEs. Here, the initial condition is ∆ , , 0 = 0 and the boundary condition is ∆ ( , 0, ) = 0. The average surface temperature variation , , can be obtained using Equation (9), which can consider multiple GHEs [7].

Overview of Ground Temperature Calculation in Simulation Tool
The average surface temperature of the GHE, 1 L L 0 T s r p−out , z, t dZ, can be obtained using Equation (8).
The second term on the right side is a temperature variation affected by the heat injection/extraction via GHEs. Here, the initial condition is ∆T s1 r p , z, 0 = 0 and the boundary condition is ∆T s1 (r, 0, t) = 0. The average surface temperature variation 1 L L 0 ∆T s1 r p−out , z, t dZ can be obtained using Equation (9), which can consider multiple GHEs [7].
where the subscript i indicates the GHE under consideration, j indicates a neighboring GHE, and j i. Thus, the first term on the right side means the average temperature variation due to the heat injection/extraction via the GHE under consideration and the second term is the average temperature variation due to the heat injection/extraction via the neighboring GHEs. Applying the superposition in the time of the ICS solution is the most suitable to calculate the surface temperature variation due to the heat injection/extraction via the considered GHE [22,23]. Then, the average surface temperature variation due to the heat injection/extraction via the considered GHE ∆T s,i,i r p−out,i , t is expressed as the following equation [7].
In addition, the average surface temperature variation due to the heat injection/extraction via the neighboring GHE, ∆T s,i,j r d,i, j , t , can be obtained as the following equation [7].
The average temperature, T * L−ave , can be obtained by modifying the temperature response, T * L , which can be obtained using the solution of infinite line source. The equation is expressed by the following [7]. Here, The average temperature, T * C−ave , can also be evaluated by substituting T * L = T * C in Equation (12). The reason is that the temperature response T * C , which can be calculated by the solution of infinite cylindrical source, is almost equal to T * L for large values of t * [7]. In addition, the modification coefficient, C, on the constant temperature condition on the ground surface can be expressed by the following equations [7].
The third term on the right side of Equation (8) is the average temperature variation caused by the ground surface distribution described in Section 2.1. This temperature variation is calculated as the Energies 2020, 13, 4491 6 of 18 difference between the initial ground temperature and the underground temperature at each depth affected by the ground surface, as shown in Equation (14).

Calculation Method for Ground Source Heat Pump System
The operation of the GSHP system was simulated by using the method for calculating the heat carrier fluid and underground temperature described in the previous section. The GSHP system mainly comprises three parts, namely, the indoor unit, the GSHP unit and the GHE [16,24]. The calculation formulas used for each element are shown below. It is assumed that there is no heat loss in the piping connecting the various parts.
(1) Indoor unit In the indoor unit, it is supposed that hot water with the temperature T 2out is sent to the generated heat load (heating load) Q 2 to process the load. It is further assumed that only those air conditioners capable of processing the load by the supply of hot water at the temperature T 2out are being used.
(2) GSHP unit Assuming that the coefficient of performance (COP) of the GSHP unit is determined by the primary inlet temperature, T 1in , and the secondary outlet temperature, T 2out , it can be expressed as follows.
Furthermore, the power consumption, E, of the heat pump can be obtained from the following equation.
Next, the heat extraction quantity (heat exchange quantity), Q 1 , in the primary side evaporator of the GSHP unit can be calculated by the following equation, using Q 2 and E.
Then, the outlet temperature, T 1out , in the primary side of the GSHP unit can be calculated using the following equation.
If T 1out is given as the inlet temperature of the GHE T pin , the fluid temperature in the GHE, T f , can be calculated using Equation (19).
Then, the outlet temperature of the GHE, T pout , is evaluated as T pout = T f . If these calculations are carried out repeatedly, the hourly temperature variation can be obtained. Figure 2 shows a calculation flow of the simulation tool. This tool calculates the approximate expression of the ground temperature distribution from the weather data before simulating the performance of the GSHP system and calculates the average temperature variation caused by the underground temperature distribution from the approximate expression. Consequently, this simulation tool can consider the influence of ground surface.

Thermal Response Test
To validate the accuracy of the simulation tool, a thermal response test (TRT) was performed as shown in Figure 3. The test was conducted at a private company's field in Karatsu City, Saga Prefecture, for two days from March 11 to 13, 2013. Using an electric boiler and a circulating pump, water was circulated at a constant flow rate and heating rate (average 2 kW) to inject heat into the ground. The temperatures of the U-tube surface inside the borehole GHE were measured using thermocouples at 9 depths: 2 m, 4 m, 6 m, 8 m, 10 m, 15 m, 20 m, 25 m and 30 m. In addition, the inlet and outlet temperature of the GHE were measured using Pt-100 temperature sensors and the flow rate was measured using an electromagnetic flowmeter.

Thermal Response Test
To validate the accuracy of the simulation tool, a thermal response test (TRT) was performed as shown in Figure 3. The test was conducted at a private company's field in Karatsu City, Saga Prefecture, for two days from March 11 to 13, 2013. Using an electric boiler and a circulating pump, water was circulated at a constant flow rate and heating rate (average 2 kW) to inject heat into the ground. The temperatures of the U-tube surface inside the borehole GHE were measured using thermocouples at 9 depths: 2 m, 4 m, 6 m, 8 m, 10 m, 15 m, 20 m, 25 m and 30 m. In addition, the inlet and outlet temperature of the GHE were measured using Pt-100 temperature sensors and the flow rate was measured using an electromagnetic flowmeter. Table 1 shows the calculation conditions. The conditions of the GHE and soil were set so as to be the same as the measurement. The underground geology is sandy silt up to a depth of 4 m, sand from 4 m to 18 m and weathered granite from 18 m to 30 m. The values of soil density and specific heat were estimated for these geological conditions as shown in Figure 4, and the weighted average was calculated based on the thickness. The value of effective thermal conductivity from a previously performed TRT was given. Silica sand was used as grouting material. As a condition for the simulation, the grout (silica sand) thermal conductivity is considered to be about 1.8 W/m.K. However, there is a possibility that convection may occur, thus in consideration of the effect, the value was doubled to 3.6 W/m.K. In the simulation, the measured values of the heating rate and the flow rate were given hourly, and the temperature of the heat carrier fluid in the GHE and the temperature of the GHE at a depth of 2 m, 6 m and 10 m were calculated.  Table 1 shows the calculation conditions. The conditions of the GHE and soil were set so as to be the same as the measurement. The underground geology is sandy silt up to a depth of 4 m, sand from 4 m to 18 m and weathered granite from 18 m to 30 m. The values of soil density and specific heat were estimated for these geological conditions as shown in Figure 4, and the weighted average was calculated based on the thickness. The value of effective thermal conductivity from a previously performed TRT was given. Silica sand was used as grouting material. As a condition for the simulation, the grout (silica sand) thermal conductivity is considered to be about 1.8 W/m.K. However, there is a possibility that convection may occur, thus in consideration of the effect, the value was doubled to 3.6 W/m.K. In the simulation, the measured values of the heating rate and the flow rate were given hourly, and the temperature of the heat carrier fluid in the GHE and the temperature of the GHE at a depth of 2 m, 6 m and 10 m were calculated.   In addition, the weather data at Esashiki (Saga Prefecture), which is the closest to the test site, was used to calculate the underground temperature distribution. When comparing the borehole inside temperatures, the borehole surface temperature was calculated in the simulation, while the U-tube surface temperature were measured by placing the temperature sensor in the middle of the inlet and the outlet pipe in the TRT, as shown in Figure 4. Therefore, the surface temperature was calculated using Equation (20) to compare and validate.

Calculation Conditions
Here, Q p is the heat injection rate (negative value of the heat extraction rate), K p−out is the overall heat transfer coefficient from the borehole surface to the U-tube surface, and A p−out is the surface area of the borehole. K p−out was obtained using the boundary element method shown in the previous research [16].

Result and Discussion
The temperatures of heat carrier fluid at inlet and outlet of the GHE in the TRT were measured and calculated, and they were compared as shown in Figure 4. In addition, the variations of U-tube surface temperature at different depths obtained using the measurement and calculation are shown in Figure 5. According to the results, the average relative error of the temperatures of the heat carrier fluid between the measurement and calculation is 3.3%. This value shows acceptable precision. When comparing the U-tube surface temperature during the TRT, the difference between the temperatures was hardly observed at a depth of 2 m. The differences between the temperatures at the depth of 6 m and 10 m are approximately 1~2 • C. Therefore, this simulation tool can reproduce the measurement within an allowable range.
In addition, it is considered that the calculated value of the temperature variation after the TRT well reproduces the practically measured value.
( , ) = , , − / Here, is the heat injection rate (negative value of the heat extraction rate), is the overall heat transfer coefficient from the borehole surface to the U-tube surface, and is the surface area of the borehole.
was obtained using the boundary element method shown in the previous research [16].

Result and Discussion
The temperatures of heat carrier fluid at inlet and outlet of the GHE in the TRT were measured and calculated, and they were compared as shown in Figure 4. In addition, the variations of U-tube surface temperature at different depths obtained using the measurement and calculation are shown in Figure 5.
According to the results, the average relative error of the temperatures of the heat carrier fluid between the measurement and calculation is 3.3%. This value shows acceptable precision. When comparing the U-tube surface temperature during the TRT, the difference between the temperatures was hardly observed at a depth of 2 m. The differences between the temperatures at the depth of 6 m and 10 m are approximately 1~2 °C. Therefore, this simulation tool can reproduce the measurement within an allowable range.
In addition, it is considered that the calculated value of the temperature variation after the TRT well reproduces the practically measured value.

Calculation Conditions
The influence of the ground surface on the temperature variation during the TRT, as shown in Figure 6, was evaluated using the developed tool. GHEs with borehole lengths of 20 m and 100 m were used and three meteorological conditions for February, May, and August were set. The temperature variations of the heat carrier fluid in the GHE were calculated and the results were compared. Table 2 shows the calculation conditions. As other calculation conditions, a constant flow rate of 0.2 L/min/m (4 L/min for GHE with 20 m long and 20 L/min for GHE with 100 m long) and a heating rate of 50 W/m (1 kW for GHE with 20 m long and 5 kW for GHE with 100 m long) were given. The values of the GHE and soil were assumed to be general values in Japan. The calculation period was 100 h. Yahata (Kitakyushu) was selected as the area to give the weather data.

Calculation Conditions
The influence of the ground surface on the temperature variation during the TRT, as shown in Figure 6, was evaluated using the developed tool. GHEs with borehole lengths of 20 m and 100 m were used and three meteorological conditions for February, May, and August were set. The temperature variations of the heat carrier fluid in the GHE were calculated and the results were compared. Table 2 shows the calculation conditions. As other calculation conditions, a constant flow rate of 0.2 L/min/m (4 L/min for GHE with 20 m long and 20 L/min for GHE with 100 m long) and a heating rate of 50 W/m (1 kW for GHE with 20 m long and 5 kW for GHE with 100 m long) were given. The values of the GHE and soil were assumed to be general values in Japan. The calculation period was 100 h. Yahata (Kitakyushu) was selected as the area to give the weather data.
The influence of the ground surface on the temperature variation during the TRT, as shown in Figure 6, was evaluated using the developed tool. GHEs with borehole lengths of 20 m and 100 m were used and three meteorological conditions for February, May, and August were set. The temperature variations of the heat carrier fluid in the GHE were calculated and the results were compared. Table 2 shows the calculation conditions. As other calculation conditions, a constant flow rate of 0.2 L/min/m (4 L/min for GHE with 20 m long and 20 L/min for GHE with 100 m long) and a heating rate of 50 W/m (1 kW for GHE with 20 m long and 5 kW for GHE with 100 m long) were given. The values of the GHE and soil were assumed to be general values in Japan. The calculation period was 100 h. Yahata (Kitakyushu) was selected as the area to give the weather data.    Figures 7 and 8 show the variations of the heat carrier fluid in the GHE with the different borehole lengths. When the borehole length of the GHE was set to 20 m, there was a temperature difference of approximately 2 • C between winter and summer, and a difference in the temperature gradient for the first several hours was observed due to the influence of the ground surface. However, when the gradients of logarithmic approximation obtained from the temperature variations were compared, it was found that there was no significant difference. In addition, when the borehole length of the GHE was set to 100 m, a difference in the temperature variation was hardly observed. As a result, it is confirmed that the influence of the ground surface on the TRT results such as the effective thermal conductivity is small even when a GHE with a borehole length of 20 m is used. borehole lengths. When the borehole length of the GHE was set to 20 m, there was a temperature difference of approximately 2 °C between winter and summer, and a difference in the temperature gradient for the first several hours was observed due to the influence of the ground surface. However, when the gradients of logarithmic approximation obtained from the temperature variations were compared, it was found that there was no significant difference. In addition, when the borehole length of the GHE was set to 100 m, a difference in the temperature variation was hardly observed. As a result, it is confirmed that the influence of the ground surface on the TRT results such as the effective thermal conductivity is small even when a GHE with a borehole length of 20 m is used.

Calculation Conditions
Assuming that the GSHP system is installed in a residential house in a moderate climate area, the effect of the short GHE on the system performance is investigated. Table 3 shows the calculation conditions. The building is a standard house set by the Japan Architectural Association. The airconditioning section is set to the temperature and humidity shown in Table 3. The heating and cooling load under this condition is calculated using the building energy simulation tool BEST [25].

Calculation Conditions
Assuming that the GSHP system is installed in a residential house in a moderate climate area, the effect of the short GHE on the system performance is investigated. Table 3 shows the calculation conditions. The building is a standard house set by the Japan Architectural Association. The air-conditioning section is set to the temperature and humidity shown in Table 3. The heating and cooling load under this condition is calculated using the building energy simulation tool BEST [25]. Figure 9 shows the time-dependent fluctuation of the heating and cooling load of the building. The heat source equipment used in this house is assumed to be a small-capacity water-cooled package air-conditioning heat pump. The total borehole length of the GHE was set to be 100 m, and the calculation was performed under four conditions of 100 m × 1 borehole, 50 m × 2 boreholes, 20 m × 5 boreholes and 10 m × 10 boreholes. In the case of multiple GHEs, the GHEs are arranged in a horizontal row at intervals of 5 m and connected in series to each other. Single U-tube was used, and general values were given for soil conditions. Two types of soil effective thermal conductivity, 1.0 W/m.K for considering unsaturated soil and 1.5 W/m.K for considering saturated soil, are given. Yahata (Kitakyushu) was selected as the area to give the weather data. The calculation period is two years. After changing the length condition as described above, the monthly ground temperature variation, power consumption, etc. in the second year were calculated and compared.   Figure 11. Figures 10 and 11 also show the monthly variation of the average surface temperature of the GHE when the boundary condition of the ground surface is set as the constant temperature condition (17.8 • C). When the soil effective thermal conductivity is lower (1.0 W/m.K), the average surface temperature of the GHE during the cooling season is approximately from 2 to 4 • C higher compared to a temperature for the condition of soil effective thermal conductivity of 1.5 W/m.K. In addition, the average surface temperature of the GHE during the heating season is approximately

Result and Discussion
Under the condition that the borehole length and number of the GHE are 100 m and 1 borehole (with or without considering the temperature variation of the ground surface), the results show that the temperatures in the two cases have a difference of from 0.3 to 0.5 • C during cooling and from 0.2 to 0.3 • C during heating, respectively.
The comparison of temperature variations with the different borehole length conditions shows that the GHE with the shorter borehole length and the larger number yields larger temperature variation. If the influence of ground surface (temperature variation on the ground surface) is considered, the average surface temperature of the GHE with the borehole length of 10 m becomes approximately 2 • C higher than the average surface temperature of the GHE with the borehole length of 100 m in August, regardless of the effective thermal conductivity. In January, the average surface temperature of the GHE with the borehole length of 10 m decreases 2 • C more, while the influence of ground surface on the condition of borehole length of 20 m is less than the case of 10 m. air-conditioning heat pump. The total borehole length of the GHE was set to be 100 m, and the calculation was performed under four conditions of 100 m × 1 borehole, 50 m × 2 boreholes, 20 m × 5 boreholes and 10 m × 10 boreholes. In the case of multiple GHEs, the GHEs are arranged in a horizontal row at intervals of 5 m and connected in series to each other. Single U-tube was used, and general values were given for soil conditions. Two types of soil effective thermal conductivity, 1.0 W/m.K for considering unsaturated soil and 1.5 W/m.K for considering saturated soil, are given. Yahata (Kitakyushu) was selected as the area to give the weather data. The calculation period is two years. After changing the length condition as described above, the monthly ground temperature variation, power consumption, etc. in the second year were calculated and compared. Figure 9. Total heating and cooling load of a residential house.  Figure 11. Figures 10 and 11 also show the monthly variation of the average surface temperature of the GHE when the boundary condition of the ground surface is set as the constant temperature condition (17.8 °C). When the soil effective thermal conductivity is lower (1.0 W/m.K), the average surface temperature of the GHE during the cooling season is approximately from 2 to 4 °C higher compared to a temperature for the condition of soil effective thermal conductivity of 1.5 W/m.K. In addition, the average surface temperature of the GHE during the heating season is approximately 1 °C lower. Under the condition that the borehole length and number of the GHE are 100 m and 1 borehole (with or without considering the temperature variation of the ground surface), the results show that the temperatures in the two cases have a difference of from 0.3 to 0.5 °C during cooling and from 0.2 to 0.3 °C during heating, respectively.

Result and Discussion
The comparison of temperature variations with the different borehole length conditions shows that the GHE with the shorter borehole length and the larger number yields larger temperature variation. If the influence of ground surface (temperature variation on the ground surface) is considered, the average surface temperature of the GHE with the borehole length of 10 m becomes approximately 2 °C higher than the average surface temperature of the GHE with the borehole length of 100 m in August, regardless of the effective thermal conductivity. In January, the average surface temperature of the GHE with the borehole length of 10 m decreases 2 °C more, while the influence of ground surface on the condition of borehole length of 20 m is less than the case of 10 m.   Figure 13 shows the results for a soil Under the condition that the borehole length and number of the GHE are 100 m and 1 borehole (with or without considering the temperature variation of the ground surface), the results show that the temperatures in the two cases have a difference of from 0.3 to 0.5 °C during cooling and from 0.2 to 0.3 °C during heating, respectively.
The comparison of temperature variations with the different borehole length conditions shows that the GHE with the shorter borehole length and the larger number yields larger temperature variation. If the influence of ground surface (temperature variation on the ground surface) is considered, the average surface temperature of the GHE with the borehole length of 10 m becomes approximately 2 °C higher than the average surface temperature of the GHE with the borehole length of 100 m in August, regardless of the effective thermal conductivity. In January, the average surface temperature of the GHE with the borehole length of 10 m decreases 2 °C more, while the influence of ground surface on the condition of borehole length of 20 m is less than the case of 10 m.   Figure 13 shows the results for a soil Finally, Figures 14 and 15 show the annual energy consumption for the different borehole lengths and the ratios when 100 m × 1 borehole is set to 1. Annual energy consumption increased by 5 to 6% for 10 m GHEs compared to a 100 m GHE, however, it increased only approximately 2% for 20 m GHEs. From this, it was confirmed that the performance of the GSHP system using GHEs with the borehole length of 20 m was not significantly decreased compared with the GSHP system using a GHE with a borehole length of 100 m. thermal conductivity of 1.5 W/m.K, and Figure 14 shows the results for 1.0 W/m.K. When the electric power consumption for the conditions of borehole length of 10 m (×10 boreholes) and 100 m (×1 borehole) are compared, the electric power consumption for the condition of borehole length of 10 m increases 8.1% in the cooling season. Similarly, the electric power consumption for the condition of borehole length of 10 m increases 2.4% in the heating season. When the COPs are compared, it was found that the COP in the case of a 10 m borehole length decreases by 0.5 during the cooling season compared to the COP in the case of a 100 m borehole length. Finally, Figures 14 and 15 show the annual energy consumption for the different borehole lengths and the ratios when 100 m × 1 borehole is set to 1. Annual energy consumption increased by 5 to 6% for 10 m GHEs compared to a 100 m GHE, however, it increased only approximately 2% for 20 m GHEs. From this, it was confirmed that the performance of the GSHP system using GHEs with the borehole length of 20 m was not significantly decreased compared with the GSHP system using a GHE with a borehole length of 100 m.    thermal conductivity of 1.5 W/m.K, and Figure 14 shows the results for 1.0 W/m.K. When the electric power consumption for the conditions of borehole length of 10 m (×10 boreholes) and 100 m (×1 borehole) are compared, the electric power consumption for the condition of borehole length of 10 m increases 8.1% in the cooling season. Similarly, the electric power consumption for the condition of borehole length of 10 m increases 2.4% in the heating season. When the COPs are compared, it was found that the COP in the case of a 10 m borehole length decreases by 0.5 during the cooling season compared to the COP in the case of a 100 m borehole length. Finally, Figures 14 and 15 show the annual energy consumption for the different borehole lengths and the ratios when 100 m × 1 borehole is set to 1. Annual energy consumption increased by 5 to 6% for 10 m GHEs compared to a 100 m GHE, however, it increased only approximately 2% for 20 m GHEs. From this, it was confirmed that the performance of the GSHP system using GHEs with the borehole length of 20 m was not significantly decreased compared with the GSHP system using a GHE with a borehole length of 100 m.    thermal conductivity of 1.5 W/m.K, and Figure 14 shows the results for 1.0 W/m.K. When the electric power consumption for the conditions of borehole length of 10 m (×10 boreholes) and 100 m (×1 borehole) are compared, the electric power consumption for the condition of borehole length of 10 m increases 8.1% in the cooling season. Similarly, the electric power consumption for the condition of borehole length of 10 m increases 2.4% in the heating season. When the COPs are compared, it was found that the COP in the case of a 10 m borehole length decreases by 0.5 during the cooling season compared to the COP in the case of a 100 m borehole length. Finally, Figures 14 and 15 show the annual energy consumption for the different borehole lengths and the ratios when 100 m × 1 borehole is set to 1. Annual energy consumption increased by 5 to 6% for 10 m GHEs compared to a 100 m GHE, however, it increased only approximately 2% for 20 m GHEs. From this, it was confirmed that the performance of the GSHP system using GHEs with the borehole length of 20 m was not significantly decreased compared with the GSHP system using a GHE with a borehole length of 100 m.

Conclusions
(1) By applying the superposition theorem, the authors developed a GHE calculation model considering the influence of the ground surface. In addition, the GHE calculation model was combined with a simulation tool for a GSHP system. The outlines of the calculation model for a short GHE were explained. (2) A TRT was carried out using a GHE with a borehole length of 30 m, and the heat carrier fluid temperature in the GHE and the underground temperature were obtained using measurements. The results were then compared with the results from the simulation. The relative error of the temperature of the heat carrier fluid was 3.3%. This result indicated that the tool can reproduce the measurements with acceptable precision. (3) The influence of the ground surface on the TRT result was investigated using the developed tool.
The result showed that the difference of the effective thermal conductivities obtained from TRTs for different seasons is small even when a GHE with a borehole length of 20 m is used. (4) It was assumed that the GSHP system is installed in a residential house to discuss the influence of the ground surface temperature on the temperature variation of heat carrier fluid in GHEs and the performance of the GSHP system. Comparing the calculation results of GHEs with borehole lengths of 10 m and 100 m under the same total length of GHEs and the same effective thermal conductivity, it was found that the average surface temperature of a GHE with a borehole length of 10 m becomes approximately 2 °C higher than the average surface temperature of a GHE with a borehole length of 100 m in August. Comparing the annual power consumption, a GSHP system using 10 GHEs with the borehole length of 10 m increases the consumption by 5 to 6% compared to a GSHP system using one GHE with the borehole length of 100 m. In the case of the GHE with a borehole length of 20 m, it was found that the increase of annual power consumption was only about 2%.  Heat pump primary side 1f Heat carrier fluid in the primary side 1in Inlet in the primary side 1out Outlet in the primary side 2 Heat pump secondary side a Ambient air C Temperature response of cylindrical heat sources d Distance f Heat carrier fluid (water or antifreeze liquid) I Solar radiation J Effective radiation L Temperature response of line heat source p Ground heat exchanger pin Inlet of ground heat exchanger pout Outlet of ground heat exchanger p-in Inside of ground heat exchanger p-out Outside of ground heat exchanger s Soil s0 Soil initial s1 Soil affected by the heat injection/extraction via GHEs s2 Soil affected by the underground temperature distribution t Affected by ambient air temperature U-out U-tube outside -ave Average