A New Simulation Model for Vertical Spiral Ground Heat Exchangers Combining Cylindrical Source Model and Capacity Resistance Model

Considering the heat capacity inside vertical spiral ground heat exchanger (VSGHEX) in the simulation is one of the most noteworthy challenge to design the ground source heat pump (GSHP) system with VSGHEXs. In this paper, a new simulation model for VSGHEXs is developed by combining the ICS model with the CaRM. The developed simulation model can consider the heat capacity inside VSGHEX and provide dynamic calculation with high speed and appropriate precision. In order to apply the CaRM, the equivalent length was introduced. Then, the equivalent length was approximated by comparing the results of the CaRM and the numerical calculation. In addition, the calculation model of the VSGHEX was integrated into the design and simulation tool for the GSHP system. The accuracy of the tool was verified by comparing with the measurements. The error between supply temperatures of the measurements and calculation is approximately 2 °C at the maximum. Finally, assuming GSHP systems with VSGHEXs, whose spiral diameter was 500 mm and depth was 4 m, were installed in residential houses in Japan, the required numbers of VSGHEXs were estimated. The results showed a strong correlation between the total heating or cooling load and the required number. Therefore, the required number can be estimated by using the simplified approximate equation.


Introduction
In recent years, global warming due to increase of CO 2 emission has become a worldwide environmental issue. This problem has created interest in applying ground source heat pump (GSHP) systems and borehole thermal energy storage (BTES) systems with vertical ground heat exchangers (GHEXs). However, these systems are not widely installed in Japan because of high installation cost. Especially, the boring cost is more expensive than in other countries. GSHP systems using vertical spiral ground heat exchangers (VSGHEXs) can reduce the installation costs, because the borehole for the VSGHEX can be drilled by using the vehicle that is used for installing the electric poles. However, the diameter of VSGHEXs is much larger than the diameter of conventional GHEXs with U-tubes. Therefore, it is important to consider the heat capacity inside the VSGHEX in the performance prediction of the GSHP system with VSGHEXs. In order to evaluate the VSGHEX, several simulation Figure 1 shows the concept of applying the CaRM to VSGHEX and Figure 2 is the plane view. The borehole of VSGHEX is regarded as a hollow cylinder in the infinite solid. The ground temperature surrounding VSGHEX can be calculated by applying the ICS model [12,13] and the temperatures inside VSGHEX are calculated by CaRM [10].

Calculation Method for Inside Spiral Ground Heat Exchangers
As shown in Figure 2, the grout between Pipe1 and Pipe2 is called as the Core and the grout between Pipe2 and the borehole surface is called as the Shell [10]. Then, the nodes are built at the Core and Shell. In this paper, the inside of the VSGHEX is regarded as a multilayer cylinder as shown in Figure 3 when the CaRM is applied. The heat balance of the heat carrier fluid in Pipe1, Pipe2 (The point , in Figure 2), the surface of Pipe1, Pipe2 (The point , in Figure 2), Core, and Shell are expressed by Equations (1)-(6), respectively.
Heat carrier fluid in Pipe1 Heat carrier fluid in Pipe2

Calculation Method for Inside Spiral Ground Heat Exchangers
As shown in Figure 2, the grout between Pipe1 and Pipe2 is called as the Core and the grout between Pipe2 and the borehole surface is called as the Shell [10]. Then, the nodes are built at the Core and Shell. In this paper, the inside of the VSGHEX is regarded as a multilayer cylinder as shown in Figure 3 when the CaRM is applied. The heat balance of the heat carrier fluid in Pipe1, Pipe2 (The point , in Figure 2), the surface of Pipe1, Pipe2 (The point , in Figure 2), Core, and Shell are expressed by Equations (1)-(6), respectively.
Heat carrier fluid in Pipe1 Heat carrier fluid in Pipe2

Calculation Method for Inside Spiral Ground Heat Exchangers
As shown in Figure 2, the grout between Pipe1 and Pipe2 is called as the Core and the grout between Pipe2 and the borehole surface is called as the Shell [10]. Then, the nodes are built at the Core and Shell. In this paper, the inside of the VSGHEX is regarded as a multilayer cylinder as shown in Figure 3 when the CaRM is applied. The heat balance of the heat carrier fluid in Pipe1, Pipe2 (The point T f 1 , T f 2 in Figure 2), the surface of Pipe1, Pipe2 (The point T sp1 , T sp2 in Figure 2), Core, and Shell are expressed by Equations (1)-(6), respectively. ) + (6) Here, = , = , = . As explained above, the inside of the VSGHEX is regarded as a multilayer cylinder in this paper. However, the calculation error occurred because the actual configuration of the heat exchanger is spiral as shown at the left in Figure 3. Therefore, the new model is proposed by introducing the equivalent length as shown in Equations (4) and (5) to improve the calculation accuracy. In this paper, the equivalent length is expressed as the following equation.
l pitch Heat carrier fluid in Pipe1 Heat carrier fluid in Pipe2 Surface of Pipe1 Surface of Pipe2 Core c grout ρ grout V core dT core dt = T sp1 − T core 1 2πl b λ grout ln r core r sp1 Shell Here, T f 1in = T pin , T f 2in = T f 1out , T pout = T f 2out . As explained above, the inside of the VSGHEX is regarded as a multilayer cylinder in this paper.
However, the calculation error occurred because the actual configuration of the heat exchanger is spiral as shown at the left in Figure 3. Therefore, the new model is proposed by introducing the Energies 2020, 13, 1339 5 of 17 equivalent length l p2 as shown in Equations (4) and (5) to improve the calculation accuracy. In this paper, the equivalent length l p2 is expressed as the following equation.
Here, coefficient c can be expressed as the function of parameters d p−out , d b , l pitch and the parameters are indicated in Figure 3. The coefficient c is determined by comparing the detailed numerical calculation, in which the finite volume method is applied. In order to compare the two calculation methods, the simplified problem shown in Figure 4 was provided. In the simplified problem, the fluid inside the Pipe1 and Pipe2 are not considered. Also, only a part of grout and Pipe2 shown in Figure 4 is considered and the part of Pipe2 is regarded as a ring. The initial and boundary condition are indicated in Figure 4. The parameters d p−out , d b , l pitch (in the numerical calculation), or l p2 (in the developed method) are changed as shown in Table 1 and the heat transmissions are repeatedly calculated by using the numerical calculation and the developed method. The commercially available numerical software stream Ver. 13 was used for numerical calculation. Then, when the heat transmissions at steady state agree with each other as shown in Figure 5, the value of l p2 is determined. The values of l p2 and c are indicated in Table 1. Furthermore, the coefficient c is approximated as the function of parameters d p−out , d b , l pitch by using the multiple linear regression analysis. As the result, the approximate equation of c can be represented in the following equation.
Energies 2020, 13, x FOR PEER REVIEW 5 of 17 = Here, coefficient c can be expressed as the function of parameters , , and the parameters are indicated in Figure 3. The coefficient c is determined by comparing the detailed numerical calculation, in which the finite volume method is applied. In order to compare the two calculation methods, the simplified problem shown in Figure 4 was provided. In the simplified problem, the fluid inside the Pipe1 and Pipe2 are not considered. Also, only a part of grout and Pipe2 shown in Figure 4 is considered and the part of Pipe2 is regarded as a ring. The initial and boundary condition are indicated in Figure 4. The parameters , , (in the numerical calculation), or (in the developed method) are changed as shown in Table 1 and the heat transmissions are repeatedly calculated by using the numerical calculation and the developed method. The commercially available numerical software stream Ver. 13 was used for numerical calculation. Then, when the heat transmissions at steady state agree with each other as shown in Figure 5, the value of is determined. The values of and c are indicated in Table 1. Furthermore, the coefficient c is approximated as the function of parameters , , by using the multiple linear regression analysis. As the result, the approximate equation of c can be represented in the following equation.     Table 2 shows the calculation results of heat transmission at steady state obtained by the numerical calculation and the developed method (Applying the CaRM and ). At the maximum, the relative error is 2.1 %.

Calculation Method for Underground Temperature
The ground temperature surrounding the VSGHEX is calculated by applying the ICS model. The GHEX is considered as a hollow cylinder in an infinite medium, and the temperature variation is calculated as a transient heat transfer problem of a two-dimensional cylindrical coordinate system. The temperature change ∆ ( , ) can be calculated by the following equation, which applies the principle of superposition based on the analytical solution of ICS [13] and Duhamel's theorem.  Table 2 shows the calculation results of heat transmission at steady state obtained by the numerical calculation and the developed method (Applying the CaRM and l p2 ). At the maximum, the relative error is 2.1%.

Calculation Method for Underground Temperature
The ground temperature surrounding the VSGHEX is calculated by applying the ICS model. The GHEX is considered as a hollow cylinder in an infinite medium, and the temperature variation is calculated as a transient heat transfer problem of a two-dimensional cylindrical coordinate system. The temperature change ∆T s (r b , t) can be calculated by the following equation, which applies the principle of superposition based on the analytical solution of ICS [13] and Duhamel's theorem. Now, Also, if the heat injection/extraction from the surface of the underground heat exchanger is considered to occur in a step-wise manner, Equation (9) can be simplified as in Equation (10) [18].
Furthermore, Equation (10) is translated to the following equation if the dimensionless quantities introduced [17].
This property has been used to simplify the computation of ∆I(r b , t − t l ) and further speed up the calculation. Having determined q l for each instant from the formula q = Also, the superposition of the temperature field in space was applied when calculating the underground temperature due to the heat injection/extraction into/from multiple GHEXs. The detail of calculation method for underground temperature due to the heat injection/extraction into/from multiple GHEXs is described in earlier reports [15][16][17]. Furthermore, when the VSGHEXs are applied, there is quite a lot of case where the diameter of borehole is much larger and the borehole depth is much smaller. In this case, in order to consider the influence of heat transfer from the ground surface and the edge of GHEX, the method that calculates the average temperature on the surface of GHEX affected by the ground surface and the edge is applied [15,17].

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 GHEX [12]. 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 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 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.
Energies 2020, 13, 1339 8 of 17 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 GHEX T pin , the outlet temperature of GHEX T pout can be calculated by using Equations (1)- (6). If these calculations are carried out repeatedly, the hourly temperature variation can be obtained.

Outlines of Measurement
Considering the case of a GSHP system with VSGHEXs installed in a residential house in Miyagi Prefecture, Japan, the temperature change of heat carrier fluid at the outlet of VSGHEXs was calculated, and its reproducibility was verified by comparing with the measured values. The outlines of the residential house are shown in Table 3. The floor area of the residential house is 57 m 2 , and the Q value (Heat loss coefficient per floor area and temperature difference) is 2.92 W/m 2 ·K. The GSHP system with VSGHEX was installed in the house.  Figure 6 shows the arrangement of VSGHEXs. The number of VSGHEXs was six. Three of the six VSGHEXs had a diameter of 400 mm, and the other three had a diameter of 600 mm. The depth of each VSGHEX was 4 m. The interval of the adjacent VSGHEXs was 3-5 m to prevent a large influence from the neighboring VSGHEX's injection or extraction. Figure 7 shows the description of GSHP system. A commercially available GSHP unit was installed, and a fan-coil unit was used for heating and cooling. There is an inhabitant in the residential house. In case the resident was absent from home, the temperature setting was fixed, as shown in Table 4. When the inhabitant stayed in the house, the temperature setting depends on the inhabitant.    Figure 8 shows the integrated value of the measured electric energy, heat extraction, and heating output. The system COP (SCOP = heating or cooling output/electric energy of heat pump and circulation pump) is also indicated in Figure 8. The system COPs for the cooling season and heating season were approximately 4.0 and 3.3, respectively. The hourly heating and cooling load is shown in the upper portion of Figure 9, and the temperature variations of heat carrier fluid in the primary side are indicated at the bottom of Figure 9. The temperatures of the heat carrier fluid were changed     Figure 8 shows the integrated value of the measured electric energy, heat extraction, and heating output. The system COP (SCOP = heating or cooling output/electric energy of heat pump and circulation pump) is also indicated in Figure 8. The system COPs for the cooling season and heating season were approximately 4.0 and 3.3, respectively. The hourly heating and cooling load is shown in the upper portion of Figure 9, and the temperature variations of heat carrier fluid in the primary side are indicated at the bottom of Figure 9. The temperatures of the heat carrier fluid were changed   Figure 8 shows the integrated value of the measured electric energy, heat extraction, and heating output. The system COP (SCOP = heating or cooling output/electric energy of heat pump and circulation pump) is also indicated in Figure 8. The system COPs for the cooling season and heating season were approximately 4.0 and 3.3, respectively. The hourly heating and cooling load is shown in the upper portion of Figure 9, and the temperature variations of heat carrier fluid in the primary side are indicated at the bottom of Figure 9. The temperatures of the heat carrier fluid were changed in response to the heating and cooling load. The supply temperature decreased to approximately 2 • C during the heating period and increased to approximately 30 • C during the cooling period.

Result of Measurement
Energies 2020, 13, x FOR PEER REVIEW 10 of 17 in response to the heating and cooling load. The supply temperature decreased to approximately 2 °C during the heating period and increased to approximately 30 °C during the cooling period.

Validation of Calculation
The calculation model of the VSGHEX was integrated into the design and simulation tool for the GSHP system [12,[15][16][17]. The accuracy of the tool was verified by comparing with the measurements. The hourly heating and cooling load was given as the condition. Then, the temperature of the heat carrier fluid was simulated by using the tool. Figure 10 shows the supply temperatures of the measurement and calculation, and Figure 11 shows a comparison of the temperatures between

Validation of Calculation
The calculation model of the VSGHEX was integrated into the design and simulation tool for the GSHP system [12,[15][16][17]. The accuracy of the tool was verified by comparing with the measurements. The hourly heating and cooling load was given as the condition. Then, the temperature of the heat carrier fluid was simulated by using the tool. Figure 10 shows the supply temperatures of the measurement and calculation, and Figure 11 shows a comparison of the temperatures between

Validation of Calculation
The calculation model of the VSGHEX was integrated into the design and simulation tool for the GSHP system [12,[15][16][17]. The accuracy of the tool was verified by comparing with the measurements. The hourly heating and cooling load was given as the condition. Then, the temperature of the Energies 2020, 13, 1339 11 of 17 heat carrier fluid was simulated by using the tool. Figure 10 shows the supply temperatures of the measurement and calculation, and Figure 11 shows a comparison of the temperatures between measurement and calculation. The error between supply temperatures of the measurements and calculation is approximately 2 • C at the maximum. It can be said that the calculated temperature reproduces the measured one with appropriate accuracy.
Energies 2020, 13, x FOR PEER REVIEW 11 of 17 measurement and calculation. The error between supply temperatures of the measurements and calculation is approximately 2 °C at the maximum. It can be said that the calculated temperature reproduces the measured one with appropriate accuracy.

Calculation Conditions
By repeating the calculation of heat carrier fluid with changing the conditions of the area, the insulation performance of the residential house, and the thermal conductivity of soil, the required numbers of VSGHEXs were estimated. The floor plan is the same as the standard residential house determined by the Architectural Institute of Japan. The floor area of the house is 125 m 2 . All the prefectural capital cities in Japan except for Okinawa (46 cities) were selected. It was assumed that the houses have two types of insulation performance, expressed as Q values (heat loss coefficient per floor area and temperature difference) or average U values. Some values were the same as the next-

Calculation Conditions
By repeating the calculation of heat carrier fluid with changing the conditions of the area, the insulation performance of the residential house, and the thermal conductivity of soil, the required numbers of VSGHEXs were estimated. The floor plan is the same as the standard residential house determined by the Architectural Institute of Japan. The floor area of the house is 125 m 2 . All the prefectural capital cities in Japan except for Okinawa (46 cities) were selected. It was assumed that the houses have two types of insulation performance, expressed as Q values (heat loss coefficient per floor area and temperature difference) or average U values. Some values were the same as the next-

Calculation Conditions
By repeating the calculation of heat carrier fluid with changing the conditions of the area, the insulation performance of the residential house, and the thermal conductivity of soil, the required numbers of VSGHEXs were estimated. The floor plan is the same as the standard residential house determined by the Architectural Institute of Japan. The floor area of the house is 125 m 2 . All the prefectural capital cities in Japan except for Okinawa (46 cities) were selected. It was assumed that the houses have two types of insulation performance, expressed as Q values (heat loss coefficient per floor area and temperature difference) or average U values. Some values were the same as the next-generation energy standard values in Japan, and the others were the lower values, as shown in Table 5. The set temperature is 22 • C for heating and 26 • C for cooling. The hourly heat load was calculated by using a commercially available heat load calculation tool, AE-CAD/AE-Simheat [19]. As example of thermal load calculation, Figure 12 shows the total heating load, the maximum heating load, the total cooling load, and the maximum cooling load when the Q values of next-generation energy standard were given. generation energy standard values in Japan, and the others were the lower values, as shown in Table  5. The set temperature is 22 °C for heating and 26 °C for cooling. The hourly heat load was calculated by using a commercially available heat load calculation tool, AE-CAD/AE-Simheat [19]. As example of thermal load calculation, Figure 12 shows the total heating load, the maximum heating load, the total cooling load, and the maximum cooling load when the Q values of next-generation energy standard were given.  The required number of VSGHEXs was defined as the minimum number of VSGHEXs that can satisfy the temperature of heat carrier fluid at the inlet of GSHP's primary side in the range from −5 °C (Heating) to 35 °C (Cooling). Then, the required number of VSGHEXs was estimated by using the design and simulation tool for the GSHP system integrating the VSGHEX model. The GSHP system description and soil data, which were given as the calculated conditions are shown in Figure 13. It was assumed that the same heat pump as the measurement one was used. The spiral diameter of the VSGHEXs was 500 mm. The depth of the VSGHEXs was 4 m and The VSGHEXs were installed at 1 to 5 m depth. Grout and soil had the same thermal conductivities, which were set at 1.2 and 1.8 W/m⋅K (these are the thermal conductivities of common unsaturated and saturated sands in Japan). The undisturbed temperatures (=average ambient temperature at each city +1.5 °C) indicated in Figure 14 were used. The required number of VSGHEXs was defined as the minimum number of VSGHEXs that can satisfy the temperature of heat carrier fluid at the inlet of GSHP's primary side in the range from −5 • C (Heating) to 35 • C (Cooling). Then, the required number of VSGHEXs was estimated by using the design and simulation tool for the GSHP system integrating the VSGHEX model. The GSHP system description and soil data, which were given as the calculated conditions are shown in Figure 13. It was assumed that the same heat pump as the measurement one was used. The spiral diameter of the VSGHEXs was 500 mm. The depth of the VSGHEXs was 4 m and The VSGHEXs were installed at 1 to 5 m depth. Grout and soil had the same thermal conductivities, which were set at 1.2 and 1.8 W/m·K (these are the thermal conductivities of common unsaturated and saturated sands in Japan). The undisturbed temperatures (=average ambient temperature at each city +1.5 • C) indicated in Figure 14 were used.   Figures 15 and 16 show the required number of VSGHEXs. Here, in Japan, there is quite a lot of case where the cooling is not operated according to the schedule. In this case, the operation time of cooling is smaller. Therefore, the required numbers of VSGHEXs, that are determined by only the constrains of heating period (T1in > −5 °C), are also shown in Figures 15 and 16. When the thermal conductivity changes from 1.2 W/m⋅K to 1.8 W/m⋅K, the required numbers of VSGHEXs can be reduced. In addition, when the Q values were changed from the value of the next-generation energy standard to smaller ones, the numbers of VSGHEXs can be reduced, especially in the case of heating. Figure 17 shows the relationship between the number of VSGHEXs and the total heating load ∑ and Figure 18 shows the relationship between the number of VSGHEXs and the cooling load ∑ . This figure indicates that the number of VSGHEXs correlates with the total load regardless of the insulation performance of residential house. Also, it seems that the effective thermal conductivity and the undisturbed temperature influence the number of VSGHEXs. Therefore, the required numbers of VSGHEXs are regarded as the function of ∑ or ∑ (The units are MWh), λs and undisturbed (Initial) ground temperature Ts0. The approximate equation of the required number of VSGHEXs is determined by using multiple regression analysis. As a result, Equations (16) and (17) are obtained.

Result and Discussion
Heating: Cooling: = 20.86 + 3.28 ∑ − 6.99λ − 0.27 By using Equations (12) and (13), the required number of VSGHEXs can be estimated more easily.   Figures 15 and 16 show the required number of VSGHEXs. Here, in Japan, there is quite a lot of case where the cooling is not operated according to the schedule. In this case, the operation time of cooling is smaller. Therefore, the required numbers of VSGHEXs, that are determined by only the constrains of heating period (T1in > −5 °C), are also shown in Figures 15 and 16. When the thermal conductivity changes from 1.2 W/m⋅K to 1.8 W/m⋅K, the required numbers of VSGHEXs can be reduced. In addition, when the Q values were changed from the value of the next-generation energy standard to smaller ones, the numbers of VSGHEXs can be reduced, especially in the case of heating. Figure 17 shows the relationship between the number of VSGHEXs and the total heating load ∑ and Figure 18 shows the relationship between the number of VSGHEXs and the cooling load ∑ . This figure indicates that the number of VSGHEXs correlates with the total load regardless of the insulation performance of residential house. Also, it seems that the effective thermal conductivity and the undisturbed temperature influence the number of VSGHEXs. Therefore, the required numbers of VSGHEXs are regarded as the function of ∑ or ∑ (The units are MWh), λs and undisturbed (Initial) ground temperature Ts0. The approximate equation of the required number of VSGHEXs is determined by using multiple regression analysis. As a result, Equations (16) and (17)

Result and Discussion
Cooling: = 20.86 + 3.28 ∑ − 6.99λ − 0.27 By using Equations (12) and (13), the required number of VSGHEXs can be estimated more easily.  Figures 15 and 16 show the required number of VSGHEXs. Here, in Japan, there is quite a lot of case where the cooling is not operated according to the schedule. In this case, the operation time of cooling is smaller. Therefore, the required numbers of VSGHEXs, that are determined by only the constrains of heating period (T 1in > −5 • C), are also shown in Figures 15 and 16. When the thermal conductivity changes from 1.2 W/m·K to 1.8 W/m·K, the required numbers of VSGHEXs can be reduced. In addition, when the Q values were changed from the value of the next-generation energy standard to smaller ones, the numbers of VSGHEXs can be reduced, especially in the case of heating.     Table 5).   Table 5). Figure 17 shows the relationship between the number of VSGHEXs and the total heating load Q h and Figure 18 shows the relationship between the number of VSGHEXs and the cooling load Q c . This figure indicates that the number of VSGHEXs correlates with the total load regardless of the insulation performance of residential house. Also, it seems that the effective thermal conductivity and the undisturbed temperature influence the number of VSGHEXs. Therefore, the required numbers of VSGHEXs are regarded as the function of Q h or Q c (The units are MWh), λ s and undisturbed (Initial) ground temperature T s0 . The approximate equation of the required number of VSGHEXs is determined by using multiple regression analysis. As a result, Equations (16) and (17) Table 5).

Conclusions
(1) A new simulation model for VSGHEXs was developed by combining the ICS model with the CaRM. In this model, the ground temperature surrounding VSGHEX can be calculated by applying the ICS model and the temperatures inside VSGHEX are calculated by CaRM. The developed simulation model can consider the heat capacity inside VSGHEX and provide dynamic calculation with high speed and appropriate precision. (2) In order to apply the CaRM, the equivalent length was introduced. Then, the equivalent length was approximated by comparing the results of CaRM and the numerical calculation. As the result, the relative error of heat transmission (thermal resistance) at steady sate between the numerical calculation and the developed method (Applying the CaRM and equivalent length) was 2.1% at the maximum. (3) The calculation model of the VSGHEX was integrated into the design and simulation tool for the GSHP system. The accuracy of the tool was verified by comparing with the measurements. The error between supply temperatures of the measurements and calculation is approximately 2 °C at the maximum annually. It can be said that the calculated temperature reproduces the measured one with appropriate accuracy. (4) Assuming the GSHP system with VSGHEXs, whose spiral diameter was 500 mm and depth was 4 m, was installed in the residential house in Japan, the required numbers of VSGHEXs were estimated. The results showed a strong correlation between the total heating or cooling load and the required number. Therefore, the required number can be regarded as the function of the total heating or cooling load, the effective thermal conductivity of soil, and the undisturbed ground temperature. Then, the required number can be estimated by using the simplified approximate equation.